# TL;DR

This code block is not evaluated. Need a breakdown? Look at the following sections.

    suppressWarnings(suppressMessages(require(netDx)))
    suppressMessages(library(curatedTCGAData))

    # fetch RNA, methylation and proteomic data for TCGA BRCA set
    brca <- suppressMessages(
       curatedTCGAData("BRCA",
                   c("mRNAArray","RPPA*","Methylation_methyl27*"),
        dry.run=FALSE))

    # process input variables
    # prepare clinical variable - stage
    staget <- sub("[abcd]","",sub("t","",colData(brca)$pathology_T_stage))
    staget <- suppressWarnings(as.integer(staget))
    colData(brca)$STAGE <- staget
    # exclude normal, HER2 (small num samples)
    pam50 <- colData(brca)$PAM50.mRNA
    idx <- union(which(pam50 %in% c("Normal-like","HER2-enriched")),
        which(is.na(staget)))
    idx <- union(idx, which(is.na(pam50)))
    pID <- colData(brca)$patientID
    tokeep <- setdiff(pID, pID[idx])
    brca <- brca[,tokeep,]
    pam50 <- colData(brca)$PAM50.mRNA
    colData(brca)$pam_mod <- pam50

    # remove duplicate names
    smp <- sampleMap(brca)
    for (nm in names(brca)) {
        samps <- smp[which(smp$assay==nm),]
        notdup <- samps[which(!duplicated(samps$primary)),"colname"]
        brca[[nm]] <- suppressMessages(brca[[nm]][,notdup])
    }

    # colData must have ID and STATUS columns
    pID <- colData(brca)$patientID
    colData(brca)$ID <- pID
    colData(brca)$STATUS <- gsub(" ","_",colData(brca)$pam_mod)

    # create grouping rules
    groupList <- list()
    # genes in mRNA data are grouped by pathways
    pathList <- readPathways(fetchPathwayDefinitions("January",2018))
    groupList[["BRCA_mRNAArray-20160128"]] <- pathList[1:3]
    # clinical data is not grouped; each variable is its own feature
    groupList[["clinical"]] <- list(
          age="patient.age_at_initial_pathologic_diagnosis",
           stage="STAGE"
    )
    # for methylation generate one feature containing all probes
    # same for proteomics data
    tmp <- list(rownames(experiments(brca)[[2]]));
    names(tmp) <- names(brca)[2]
    groupList[[names(brca)[2]]] <- tmp

    tmp <- list(rownames(experiments(brca)[[3]]));
    names(tmp) <- names(brca)[3]
    groupList[[names(brca)[3]]] <- tmp

    # create function to tell netDx how to build features
    # (PSN) from your data
    makeNets <- function(dataList, groupList, netDir,...) {
        netList <- c() # initialize before is.null() check
        # correlation-based similarity for mRNA, RPPA and methylation data
        # (Pearson correlation)
        for (nm in setdiff(names(groupList),"clinical")) {
           # NOTE: the check for is.null() is important!
            if (!is.null(groupList[[nm]])) {
            netList <- makePSN_NamedMatrix(dataList[[nm]],
                         rownames(dataList[[nm]]),
                       groupList[[nm]],netDir,verbose=FALSE,
                         writeProfiles=TRUE,...)
            }
        }

        # make clinical nets (normalized difference)
        netList2 <- c()
        if (!is.null(groupList[["clinical"]])) {
        netList2 <- makePSN_NamedMatrix(dataList$clinical,
            rownames(dataList$clinical),
            groupList[["clinical"]],netDir,
            simMetric="custom",customFunc=normDiff, # custom function
            writeProfiles=FALSE,
            sparsify=TRUE,verbose=TRUE,...)
        }
        netList <- c(unlist(netList),unlist(netList2))
        return(netList)
    }

    # run predictor
    set.seed(42) # make results reproducible
    outDir <- paste(tempdir(),"pred_output",sep=getFileSep())
    # To see all messages, remove suppressMessages()
    # and set logging="default".
    # To keep all intermediate data, set keepAllData=TRUE
    numSplits <- 2L
    out <- suppressMessages(
       buildPredictor(dataList=brca,groupList=groupList,
          makeNetFunc=makeNets,outDir=outDir,
          numSplits=numSplits, featScoreMax=2L, featSelCutoff=1L,
           numCores=1L)
    )

    # collect results for accuracy
    st <- unique(colData(brca)$STATUS)
    acc <- matrix(NA,ncol=length(st),nrow=numSplits)  # accuracy by class
    colnames(acc) <- st
    for (k in 1:numSplits) {
        pred <- out[[sprintf("Split%i",k)]][["predictions"]];
        tmp <- pred[,c("ID","STATUS","TT_STATUS","PRED_CLASS",
                         sprintf("%s_SCORE",st))]
        for (m in 1:length(st)) {
           tmp2 <- subset(tmp, STATUS==st[m])
           acc[k,m] <- sum(tmp2$PRED==tmp2$STATUS)/nrow(tmp2)
        }
    }

    # accuracy by class
    print(round(acc*100,2))

    # confusion matrix
    res <- out$Split1$predictions
    print(table(res[,c("STATUS","PRED_CLASS")]))

    sessionInfo()

# Introduction

In this example, we will use clinical data and three types of 'omic data - gene expression, DNA methylation and proteomic data - to classify breast tumours as being one of three types: Luminal A, Luminal B, or Basal. This example is nearly identical to the one used to build a binary classifier.

We also use several strategies and definitions of similarity to create features:

- Clinical variables: Each variable is its own feature (e.g. age); similarity is defined as normalized difference.
- Gene expression: Features are defined at the level of pathways; i.e. a feature groups genes corresponding to the pathway. Similarity is defined as pairwise Pearson correlation
- Proteomic and methylation data: Features are defined at the level of the entire data layer; a single feature is created for all of proteomic data, and the same for methylation. Similarity is defined by pairwise Pearson correlation

# Setup

Load the `netDx` package.

    suppressWarnings(suppressMessages(require(netDx)))

# Data

For this example we pull data from the The Cancer Genome Atlas through the BioConductor `curatedTCGAData` package. The fetch command automatically brings in a `MultiAssayExperiment` object.

    suppressMessages(library(curatedTCGAData))

We use the `curatedTCGAData()` command to look at available assays in the breast cancer dataset.

    curatedTCGAData(diseaseCode="BRCA", assays="*",dry.run=TRUE)

    ##                                         Title DispatchClass
    ## 31                       BRCA_CNASeq-20160128           Rda
    ## 32                       BRCA_CNASNP-20160128           Rda
    ## 33                       BRCA_CNVSNP-20160128           Rda
    ## 35             BRCA_GISTIC_AllByGene-20160128           Rda
    ## 36                 BRCA_GISTIC_Peaks-20160128           Rda
    ## 37     BRCA_GISTIC_ThresholdedByGene-20160128           Rda
    ## 39  BRCA_Methylation_methyl27-20160128_assays        H5File
    ## 40      BRCA_Methylation_methyl27-20160128_se           Rds
    ## 41 BRCA_Methylation_methyl450-20160128_assays        H5File
    ## 42     BRCA_Methylation_methyl450-20160128_se           Rds
    ## 43                 BRCA_miRNASeqGene-20160128           Rda
    ## 44                    BRCA_mRNAArray-20160128           Rda
    ## 45                     BRCA_Mutation-20160128           Rda
    ## 46              BRCA_RNASeq2GeneNorm-20160128           Rda
    ## 47                   BRCA_RNASeqGene-20160128           Rda
    ## 48                    BRCA_RPPAArray-20160128           Rda

In this call we fetch only the gene expression, proteomic and methylation data; setting `dry.run=FALSE` initiates the fetching of the data.

    brca <- suppressMessages(
       curatedTCGAData("BRCA",
                   c("mRNAArray","RPPA*","Methylation_methyl27*"),
        dry.run=FALSE))

This next code block prepares the TCGA data. In practice you would do this once, and save the data before running netDx, but we run it here to see an end-to-end example.

    # prepare clinical variable - stage
    staget <- sub("[abcd]","",sub("t","",colData(brca)$pathology_T_stage))
    staget <- suppressWarnings(as.integer(staget))
    colData(brca)$STAGE <- staget

    # exclude normal, HER2 (small num samples)
    pam50 <- colData(brca)$PAM50.mRNA
    idx <- union(which(pam50 %in% c("Normal-like","HER2-enriched")),
        which(is.na(staget)))
    idx <- union(idx, which(is.na(pam50)))
    pID <- colData(brca)$patientID
    tokeep <- setdiff(pID, pID[idx])
    brca <- brca[,tokeep,]

    pam50 <- colData(brca)$PAM50.mRNA
    colData(brca)$pam_mod <- pam50

    # remove duplicate names
    smp <- sampleMap(brca)
    for (nm in names(brca)) {
        samps <- smp[which(smp$assay==nm),]
        notdup <- samps[which(!duplicated(samps$primary)),"colname"]
        brca[[nm]] <- suppressMessages(brca[[nm]][,notdup])
    }

    ## harmonizing input:
    ##   removing 59 sampleMap rows with 'colname' not in colnames of experiments

    ## harmonizing input:
    ##   removing 19 sampleMap rows with 'colname' not in colnames of experiments

    ## harmonizing input:
    ##   removing 26 sampleMap rows with 'colname' not in colnames of experiments

The important thing is to create `ID` and `STATUS` columns in the sample metadata slot. netDx uses these to get the patient identifiers and labels, respectively.

    pID <- colData(brca)$patientID
    colData(brca)$ID <- pID
    colData(brca)$STATUS <- gsub(" ","_",colData(brca)$pam_mod)

# Rules to create features (patient similarity networks)

Our plan is to group gene expression data by pathways and clinical data by single variables. We will treat methylation and proteomic data each as a single feature, so each of those groups will contain the entire input table for those corresponding data types.

In the code below, we fetch pathway definitions for January 2018 from (http://download.baderlab.org/EM_Genesets) and group gene expression data by pathways. To keep the example short, we limit to only three pathways, but in practice you would use all pathways meeting a size criterion; e.g. those containing between 10 and 500 genes.

Grouping rules are accordingly created for the clinical, methylation and proteomic data.

    groupList <- list()

    # genes in mRNA data are grouped by pathways
    pathList <- readPathways(fetchPathwayDefinitions("January",2018))

    ## ---------------------------------------

    ## Fetching http://download.baderlab.org/EM_Genesets/January_01_2018/Human/symbol/Human_AllPathways_January_01_2018_symbol.gmt

    ## File: df02e2b1043_Human_AllPathways_January_01_2018_symbol.gmt

    ## Read 3028 pathways in total, internal list has 3009 entries

    ##  FILTER: sets with num genes in [10, 200]

    ##    => 971 pathways excluded
    ##    => 2038 left

    groupList[["BRCA_mRNAArray-20160128"]] <- pathList[1:3]
    # clinical data is not grouped; each variable is its own feature
    groupList[["clinical"]] <- list(
          age="patient.age_at_initial_pathologic_diagnosis",
           stage="STAGE"
    )
    # for methylation generate one feature containing all probes
    # same for proteomics data
    tmp <- list(rownames(experiments(brca)[[2]]));
    names(tmp) <- names(brca)[2]
    groupList[[names(brca)[2]]] <- tmp

    tmp <- list(rownames(experiments(brca)[[3]]));
    names(tmp) <- names(brca)[3]
    groupList[[names(brca)[3]]] <- tmp

## Define patient similarity for each network

We provide `netDx` with a custom function to generate similarity networks (i.e. features). The first block tells netDx to generate correlation-based networks using everything but the clinical data. This is achieved by the call:

    makePSN_NamedMatrix(..., writeProfiles=TRUE,...)`

The second block makes a different call to `makePSN_NamedMatrix()` but this time, requesting the use of the normalized difference similarity metric. This is achieved by calling:

       makePSN_NamedMatrix(,...,
                           simMetric="custom", customFunc=normDiff,
                           writeProfiles=FALSE)

`normDiff` is a function provided in the `netDx` package, but the user may define custom similarity functions in this block of code and pass those to `makePSN_NamedMatrix()`, using the `customFunc` parameter.

    makeNets <- function(dataList, groupList, netDir,...) {
        netList <- c() # initialize before is.null() check
        # correlation-based similarity for mRNA, RPPA and methylation data
        # (Pearson correlation)
        for (nm in setdiff(names(groupList),"clinical")) {
           # NOTE: the check for is.null() is important!
            if (!is.null(groupList[[nm]])) {
            netList <- makePSN_NamedMatrix(dataList[[nm]],
                         rownames(dataList[[nm]]),
                       groupList[[nm]],netDir,verbose=FALSE,
                         writeProfiles=TRUE,...)
            }
        }

        # make clinical nets (normalized difference)
        netList2 <- c()
        if (!is.null(groupList[["clinical"]])) {
        netList2 <- makePSN_NamedMatrix(dataList$clinical,
            rownames(dataList$clinical),
            groupList[["clinical"]],netDir,
            simMetric="custom",customFunc=normDiff, # custom function
            writeProfiles=FALSE,
            sparsify=TRUE,verbose=TRUE,...)
        }
        netList <- c(unlist(netList),unlist(netList2))
        return(netList)
    }

# Build predictor

Finally we make the call to build the predictor.

    set.seed(42) # make results reproducible

    # location for intermediate work
    # set keepAllData to TRUE to not delete at the end of the
    # predictor run.
    # This can be useful for debugging.
    outDir <- paste(tempdir(),"pred_output",sep=getFileSep())

    numSplits <- 2L
    out <- suppressMessages(
       buildPredictor(dataList=brca,groupList=groupList,
          makeNetFunc=makeNets,outDir=outDir,
          numSplits=numSplits, featScoreMax=2L, featSelCutoff=1L,
           numCores=1L)
    )

    ## function(dataList, groupList, netDir,...) {
    ##  netList <- c() # initialize before is.null() check
    ##  # correlation-based similarity for mRNA, RPPA and methylation data
    ##  # (Pearson correlation)
    ##  for (nm in setdiff(names(groupList),"clinical")) {
    ##     # NOTE: the check for is.null() is important!
    ##      if (!is.null(groupList[[nm]])) {
    ##      netList <- makePSN_NamedMatrix(dataList[[nm]],
    ##                   rownames(dataList[[nm]]),
    ##                    groupList[[nm]],netDir,verbose=FALSE,
    ##                   writeProfiles=TRUE,...)
    ##      }
    ##  }
    ##
    ##  # make clinical nets (normalized difference)
    ##  netList2 <- c()
    ##  if (!is.null(groupList[["clinical"]])) {
    ##  netList2 <- makePSN_NamedMatrix(dataList$clinical,
    ##      rownames(dataList$clinical),
    ##      groupList[["clinical"]],netDir,
    ##      simMetric="custom",customFunc=normDiff, # custom function
    ##      writeProfiles=FALSE,
    ##      sparsify=TRUE,verbose=TRUE,...)
    ##  }
    ##  netList <- c(unlist(netList),unlist(netList2))
    ##  return(netList)
    ## }
    ##             IS_TRAIN
    ## STATUS       TRAIN TEST
    ##   Basal-like    77   20
    ##   Luminal_A    184   46
    ##   Luminal_B    101   26
    ##
    ## Luminal_A   nonpred      <NA>
    ##       184       178         0
    ##
    ## Basal-like    nonpred       <NA>
    ##         77        285          0
    ##
    ## Luminal_B   nonpred      <NA>
    ##       101       261         0
    ##             IS_TRAIN
    ## STATUS       TRAIN TEST
    ##   Basal-like    77   20
    ##   Luminal_A    184   46
    ##   Luminal_B    101   26
    ##
    ## Luminal_A   nonpred      <NA>
    ##       184       178         0
    ##
    ## Basal-like    nonpred       <NA>
    ##         77        285          0
    ##
    ## Luminal_B   nonpred      <NA>
    ##       101       261         0

# Examine output

The results are stored in the list object returned by the `buildPredictor()` call. This list contains:

- `inputNets`: all input networks that the model started with.
- `Split<i>`: a list with results for each train-test split
  - `featureScores`: feature scores for each label (list with `g` entries, where `g` is number of patient labels). Each entry contains the feature selection scores for the corresponding label.
  - `featureSelected`: vector of features that pass feature selection. List of length `g`, with one entry per label.
  - `predictions`: real and predicted labels for test patients
  - `accuracy`: percent accuracy of predictions

    summary(out)

    ##           Length Class  Mode
    ## inputNets 14     -none- character
    ## Split1     4     -none- list
    ## Split2     4     -none- list

    summary(out$Split1)

    ##                 Length Class      Mode
    ## featureScores      3   -none-     list
    ## featureSelected    3   -none-     list
    ## predictions     2693   data.frame list
    ## accuracy           1   -none-     numeric

Compute accuracy for three-way classificationL

    # Average accuracy
    st <- unique(colData(brca)$STATUS)
    acc <- matrix(NA,ncol=length(st),nrow=numSplits)
    colnames(acc) <- st
    for (k in 1:numSplits) {
        pred <- out[[sprintf("Split%i",k)]][["predictions"]];
        tmp <- pred[,c("ID","STATUS","TT_STATUS","PRED_CLASS",
                         sprintf("%s_SCORE",st))]
        for (m in 1:length(st)) {
           tmp2 <- subset(tmp, STATUS==st[m])
           acc[k,m] <- sum(tmp2$PRED==tmp2$STATUS)/nrow(tmp2)
        }
    }
    print(round(acc*100,2))

    ##      Luminal_A Basal-like Luminal_B
    ## [1,]     46.43        100     35.71
    ## [2,]     44.83        100     55.00

Also, examine the confusion matrix. We can see that the model perfectly classifies basal tumours, but performs poorly in distinguishing between the two types of luminal tumours.

    res <- out$Split1$predictions
    print(table(res[,c("STATUS","PRED_CLASS")]))

    ##             PRED_CLASS
    ## STATUS       Basal-like Luminal_A Luminal_B
    ##   Basal-like         14         0         0
    ##   Luminal_A           5        13        10
    ##   Luminal_B           4         5         5

# sessionInfo

    sessionInfo()

    ## R version 4.0.2 (2020-06-22)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server 2012 R2 x64 (build 9600)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] C
    ## system code page: 1252
    ##
    ## attached base packages:
    ## [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    ## [8] methods   base
    ##
    ## other attached packages:
    ##  [1] rhdf5_2.32.2                BiocFileCache_1.12.0
    ##  [3] dbplyr_1.4.4                curatedTCGAData_1.10.0
    ##  [5] MultiAssayExperiment_1.14.0 SummarizedExperiment_1.18.2
    ##  [7] DelayedArray_0.14.1         matrixStats_0.56.0
    ##  [9] Biobase_2.48.0              GenomicRanges_1.40.0
    ## [11] GenomeInfoDb_1.24.2         IRanges_2.22.2
    ## [13] S4Vectors_0.26.1            BiocGenerics_0.34.0
    ## [15] netDx_1.0.3
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] bitops_1.0-6                  bit64_0.9-7.1
    ##  [3] doParallel_1.0.15             RColorBrewer_1.1-2
    ##  [5] httr_1.4.2                    tools_4.0.2
    ##  [7] R6_2.4.1                      HDF5Array_1.16.1
    ##  [9] DBI_1.1.0                     colorspace_1.4-1
    ## [11] tidyselect_1.1.0              bit_1.1-15.2
    ## [13] curl_4.3                      compiler_4.0.2
    ## [15] glmnet_4.0-2                  graph_1.66.0
    ## [17] labeling_0.3                  scales_1.1.1
    ## [19] rappdirs_0.3.1                stringr_1.4.0
    ## [21] digest_0.6.25                 R.utils_2.9.2
    ## [23] XVector_0.28.0                pkgconfig_2.0.3
    ## [25] htmltools_0.5.0               highr_0.8
    ## [27] fastmap_1.0.1                 rlang_0.4.7
    ## [29] RSQLite_2.2.0                 shiny_1.5.0
    ## [31] farver_2.0.3                  shape_1.4.4
    ## [33] generics_0.0.2                combinat_0.0-8
    ## [35] dplyr_1.0.0                   R.oo_1.23.0
    ## [37] RCurl_1.98-1.2                magrittr_1.5
    ## [39] GenomeInfoDbData_1.2.3        Matrix_1.2-18
    ## [41] Rhdf5lib_1.10.1               Rcpp_1.0.5
    ## [43] munsell_0.5.0                 lifecycle_0.2.0
    ## [45] R.methodsS3_1.8.0             stringi_1.4.6
    ## [47] yaml_2.2.1                    RJSONIO_1.3-1.4
    ## [49] zlibbioc_1.34.0               plyr_1.8.6
    ## [51] AnnotationHub_2.20.0          grid_4.0.2
    ## [53] blob_1.2.1                    promises_1.1.1
    ## [55] ExperimentHub_1.14.0          bigmemory.sri_0.1.3
    ## [57] crayon_1.3.4                  lattice_0.20-41
    ## [59] splines_4.0.2                 knitr_1.29
    ## [61] pillar_1.4.6                  igraph_1.2.5
    ## [63] markdown_1.1                  reshape2_1.4.4
    ## [65] codetools_0.2-16              XML_3.99-0.5
    ## [67] glue_1.4.1                    BiocVersion_3.11.1
    ## [69] evaluate_0.14                 BiocManager_1.30.10
    ## [71] vctrs_0.3.2                   httpuv_1.5.4
    ## [73] foreach_1.5.0                 gtable_0.3.0
    ## [75] purrr_0.3.4                   assertthat_0.2.1
    ## [77] ggplot2_3.3.2                 xfun_0.16
    ## [79] mime_0.9                      xtable_1.8-4
    ## [81] pracma_2.2.9                  later_1.1.0.1
    ## [83] RCy3_2.8.0                    survival_3.2-3
    ## [85] tibble_3.0.3                  iterators_1.0.12
    ## [87] AnnotationDbi_1.50.3          memoise_1.1.0
    ## [89] bigmemory_4.5.36              ellipsis_0.3.1
    ## [91] interactiveDisplayBase_1.26.3 ROCR_1.0-11