diff --git a/README.md b/README.md
index 0c622bcb..7f665587 100644
--- a/README.md
+++ b/README.md
@@ -11,21 +11,28 @@ The `examples/` folder contains R code that should just run once both `netDx/` a
For more information and FAQ, visit **http://netdx.org**
* [Install netDx](#install-netdx)
+ * [Prerequisites](#prerequisites)
* [Test functionality](#test-functionality)
- * [Known issues with compiling pdfs](#known-issues-with-compiling-pdfs)
-* [Run breastcancer LumA example](#run-breastcancer-luma-example)
+ * [Run Medulloblastoma example](#mblastoma)
+ * [Run breastcancer LumA example](#brca)
+ * [See full list of examples](#other-examples)
+ * [Known issues with compiling pdfs](#pdfissue)
+**Other useful information:**
+* [Read the netDx preprint at bioRXiv](https://doi.org/10.1101/084418): Pai et al. (2016). netDx: Patient classification using integrated patient similarity networks. https://doi.org/10.1101/084418
+* Once you have run the included netDx examples, [read the user manual](http://netdx-manual.readthedocs.io/en/latest/) to learn how to design features or predictors.
## Install netDx
### Prerequisites
**netDx has been tested on Mac OS/X and on Linux systems. For now we recommend you run netDx on these operating systems.** Future versions of netDx will have Windows support.
-You must have Java, Python and R installed. Within R, you must have BioConductor installed. This section helps you figure out which
-of these you need to install. If you already have all these, skip to the next section.
+You must have Java, Python and R installed. Within R, you must have BioConductor installed. To plot the results of the predictor run, including network visualizations such as the EnrichmentMap and integrated patient similarity network, you will need Cytoscape with the latest EnrichmentMap and AutoAnnotate apps installed. This section helps you figure out which of these you need to install.
+
+If you already have all these, skip to the next section.
#### Java (1.8+ recommended, but will probably work on 1.6+)
-netDx uses the GeneMANIA algorithm to integrate patient networks and recommend patients by similarity (Mostafavi and Morris (2008). *Genome Biol* 9:Suppl 1). GeneMANIA is currently implemented in Java, making this interpreter a requirement for netDx.
+The engine netDx uses to integrate patient networks and recommend patients by similarity is implemented in Java, making this interpreter a requirement for netDx.
At command line, run `java --version`. You should see output like this:
```
@@ -36,7 +43,7 @@ Java HotSpot(TM) 64-Bit Server VM (build 25.31-b07, mixed mode)
If you don't see this kind of output, you may need to first [install java](https://java.com/en/).
#### Python (2.7 recommended)
-netDx uses legacy Python scripts in creating the GeneMANIA database, so for now a Python interpreter is required to run netDx. Future versions of netDx will not have this requirement.
+A legacy script to create the database of input patient networks was implemented in Python. For now, netDx requires Python to run but future versions of netDx will not have this requirement.
At command line, run `python --version`. You should see output like this:
@@ -70,22 +77,58 @@ $ R
```
Say `yes` to all dependencies that need to be installed.
+#### Cytoscape (3.5.1+ recommended)
+In order to generate and visualize patient similarity networks and enrichment maps from the data generated by netDx we will need Cytoscape. If Cytoscape is not located in your Applications directory, install the latest version from http://www.cytoscape.org.
+
+With Cytoscape open, install the following Apps:
+* AutoAnnotate v1.2 (http://apps.cytoscape.org/apps/autoannotate) *Earlier versions may not work.*
+* Enrichment Map v3.0.0 (http://apps.cytoscape.org/apps/enrichmentmap) *Earlier versions may not work.*
+
## Install `netDx` and `netDx.examples`
This section assumes you have Java, Python, R and Bioconductor installed. From command-line, download the git repo for these packages and install them. In the code below, output from intermediate steps is omitted for clarity.
*Note: For now, R package dependencies must be separately installed using the install.packages() call as shown below. netDx will be submitted to CRAN following publication; thereafter, dependencies can be automatically installed with the call to install netDx.*
+On Unix systems you may need one or more of these packages as these are dependencies for R packages:
+
+```
+sudo apt-get install zlib1g-dev libssl-dev libssh2-1-dev libcurl-devel
+```
+
+If you are on an RPM system you may be need to run this:
+```
+sudo apt-get install libcurl4-openssl-dev
+```
+
+Now we install the necessary R packages:
```
$ cd netDx-master/
$ R
> install.packages(c('devtools','curl'))
> install.packages(c("bigmemory","foreach","combinat","doParallel","ROCR","pracma","RColorBrewer","reshape2"))
> devtools::install_github("igraph/rigraph") # install from CRAN has a bug and can fail (31 Aug 2017).
+> devtools::install_github('cytoscape/cytoscape-automation/for-scripters/R/r2cytoscape')
+> devtools::install_github('BaderLab/Easycyrest/EasycyRest@0.1')
> install.packages("netDx",type="source",repos=NULL)
> install.packages("netDx.examples",type="source",repos=NULL)
> install.packages("knitr") # needed to run examples
```
+On Unix systems you may need to install the libraries below at command-line:
+```
+$ sudo apt-get install libssl-dev # for openssl & httr
+$ sudo apt-get install xml2-config # for XML & r2cytoscape
+```
+
+```
+>install.packages(c("openssl","httr","RJSONIO"))
+>devtools::install_github('cytoscape/cytoscape-automation/for-scripters/R/r2cytoscape')
+>devtools::install_github('BaderLab/Easycyrest/EasycyRest@0.1')
+```
+Note: On Unix systems, installing `httr` requires a prior install of the `openssl` package. If the `openssl` install fails with a message like: `ERROR: configuration failed for package 'openssl'`, you will need to install openssl for your system. e.g. for Debian and Ubuntu, install `libssl-dev`
+
+
+
## Test functionality
Run the medulloblastoma vignette to make sure the netDx pipeline works from end to end.
Each vignette is in Sweave format (`.Rnw`) . To run these, you need to have both `netDx` and `netDx.examples` installed. You will also need to install the R package `knitr` to compile the Sweave file. If you have [Rstudio](https://www.rstudio.com/home/) installed (highly recommended), you should be able to open the `Rnw` file and click `Compile PDF`. Alternately, you may run the vignette through an interactive R session:
@@ -98,6 +141,20 @@ $ R
```
This should generate `Medulloblastoma.pdf` in the `examples/` directory.
+## Run BreastCancer LumA example
+This vignette is presented in the netDx manuscript. Here we start with 348 primary tumours from the Cancer Genome Atlas, and build a predictor for Luminal A subtype classification (The Cancer Genome Atlas (2012). *Nature.* **490**:61-70). This example illustrates feature selection using a simple design in which networks are scored out of 10 based on a single round of 10-fold cross validation. On a MacBook Air laptop (late 2014), this vignette takes ~1.5 hours to run to completion. You may speed it up by running it on a machine with more processors and changing the `numCores` variable in the vignette.
+
+**We do not recommend running it on a Mac with less than 8Gb RAM. A Unix machine manages memory differently and may require as much as 32Gb RAM. If such a machine is not available, set `numCores=2L` in the `.Rnw` file before running.**
+
+```
+$ cd netDx/examples/
+$ R
+> require(knitr)
+> knit2pdf("BreastCancer.Rnw")
+```
+**NOTE:** The vignette will generate a pdf file. All intermediate files will be stored in the `TCGA_BRCA/` subdirectory of the examples directory.
+
+
## Known issues with compiling pdfs
#### (Linux)
@@ -119,13 +176,11 @@ ln -s /Library/TeX/Distributions/.DefaultTeX/Contents/Programs/texbin /usr/texbi
Step 3. Now, in the terminal check the value of `echo $PATH`. Make sure that `/usr/texbin` is present. If it isn't present, then you need to add `/usr/texbin` to your PATH variable. This can be done by updating the `PATH` variable in `~/.bashrc`.
However, if you find yourself having to mess with the PATH variable, try reinstalling the [MacTex](http://tug.org/mactex/) package.
-## Run BreastCancer LumA example
-This vignette is presented in the netDx manuscript. Here we start with 348 primary tumours from the Cancer Genome Atlas, and build a predictor for Luminal A subtype classification (The Cancer Genome Atlas (2012). *Nature.* **490**:61-70). This example illustrates feature selection using a simple design in which networks are scored out of 10 based on a single round of 10-fold cross validation. On a MacBook Air laptop (late 2014), this vignette takes ~1.5 hours to run to completion. You may speed it up by running it on a machine with more processors and changing the `numCores` variable in the vignette. We do not recommend running it on a machine with less than 8Gb RAM.
+
-```
-$ cd netDx/examples/
-$ R
-> require(knitr)
-> knit2pdf("BreastCancer.Rnw")
-```
-**NOTE:** The vignette will generate a pdf file. All intermediate files will be stored in the `TCGA_BRCA/` subdirectory of the examples directory.
+## Other examples
+The `examples/` directory contains R notebooks (`.Rmd`) that teach basic functionality useful for any predictor building. These include plotting predictor results and running nested cross-validation, the design we recommend to use after feature design.
+
+The R notebooks must be run from within Rstudio. Install [Rstudio](https://www.rstudio.com/products/rstudio/download/) if necessary.
+
+We have also posted the results of running all the examples [here](http://netdx.org/index.php/examples/).
diff --git a/README.txt b/README.txt
deleted file mode 100755
index 5f9823c4..00000000
--- a/README.txt
+++ /dev/null
@@ -1,131 +0,0 @@
-# netDx
-netDx is an algorithm for building patient classifiers by using patient similarity networks as features.
-
-This repo contains two R packages, each of which need to be separately installed:
-
-1. `netDx/`: Software implementing the netDx method
-2. `netDx.examples/`: Data for use cases
-
-The `examples/` folder contains R code that should just run once both `netDx/` and `netDx.examples/` are installed.
-
-For more information and FAQ, visit **http://netdx.org**
-
-* [Install netDx](#install-netdx)
-* [Test functionality](#test-functionality)
- * [Known issues with compiling pdfs](#known-issues-with-compiling-pdfs)
-* [Run breastcancer LumA example](#run-breastcancer-luma-example)
-
-
-## Install netDx
-
-### Prerequisites
-**netDx has been tested on Mac OS/X and on Linux systems. For now we recommend you run netDx on these operating systems.** Future versions of netDx will have Windows support.
-
-You must have Java, Python and R installed. Within R, you must have BioConductor installed. This section helps you figure out which
-of these you need to install. If you already have all these, skip to the next section.
-
-#### Java (1.8+ recommended, but will probably work on 1.6+)
-netDx uses the GeneMANIA algorithm to integrate patient networks and recommend patients by similarity (Mostafavi and Morris (2008). *Genome Biol* 9:Suppl 1). GeneMANIA is currently implemented in Java, making this interpreter a requirement for netDx.
-
-At command line, run `java --version`. You should see output like this:
-```
-java version "1.8.0_31"
-Java(TM) SE Runtime Environment (build 1.8.0_31-b13)
-Java HotSpot(TM) 64-Bit Server VM (build 25.31-b07, mixed mode)
-```
-If you don't see this kind of output, you may need to first [install java](https://java.com/en/).
-
-#### Python (2.7 recommended)
-netDx uses legacy Python scripts in creating the GeneMANIA database, so for now a Python interpreter is required to run netDx. Future versions of netDx will not have this requirement.
-
-At command line, run `python --version`. You should see output like this:
-
-```
-Python 2.7.9 :: Anaconda 2.1.0 (x86_64)
-```
-
-#### R (3.2+ recommended)
-At command line, run `R --version`. You should see output like this:
-```
-R version 3.2.4 (2016-03-10) -- "Very Secure Dishes"
-Copyright (C) 2016 The R Foundation for Statistical Computing
-Platform: x86_64-apple-darwin13.4.0 (64-bit)
-
-R is free software and comes with ABSOLUTELY NO WARRANTY.
-You are welcome to redistribute it under the terms of the
-GNU General Public License versions 2 or 3.
-For more information about these matters see
-http://www.gnu.org/licenses/.
-```
-
-If not, [install R](https://www.r-project.org/).
-
-#### BioConductor (2.3+ recommended)
-[BioConductor](http://bioconductor.org/) is a system of R objects and software specifically for biological applications.
-Once you have R installed, install the `biobase` and `GenomicRanges` packages via the BioConductor installer:
-```
-$ R
-> source("http://bioconductor.org/biocLite.R")
-> biocLite(c("Biobase","GenomicRanges"))
-```
-Say `yes` to all dependencies that need to be installed.
-
-## Install `netDx` and `netDx.examples`
-This section assumes you have Java, Python, R and Bioconductor installed. From command-line, download the git repo for these packages and install them. In the code below, output from intermediate steps is omitted for clarity.
-
-*Note: For now, R package dependencies must be separately installed using the install.packages() call as shown below. netDx will be submitted to CRAN following publication; thereafter, dependencies can be automatically installed with the call to install netDx.*
-
-```
-$ cd netDx-master/
-$ R
-> install.packages(c('devtools','curl'))
-> devtools::install_github("igraph/rigraph")
-> install.packages(c("bigmemory","foreach","combinat","doParallel","ROCR","pracma","RColorBrewer","reshape2"))
-> install.packages("netDx",type="source",repos=NULL)
-> install.packages("netDx.examples",type="source",repos=NULL)
-> install.packages("knitr") # needed to run examples
-```
-
-## Test functionality
-Run the medulloblastoma vignette to make sure the netDx pipeline works from end to end.
-Each vignette is in Sweave format (`.Rnw`) . To run these, you need to have both `netDx` and `netDx.examples` installed. You will also need to install the R package `knitr` to compile the Sweave file. If you have [Rstudio](https://www.rstudio.com/home/) installed (highly recommended), you should be able to open the `Rnw` file and click `Compile PDF`. Alternately, you may run the vignette through an interactive R session:
-
-```
-$ cd netDx-master/examples/
-$ R
-> require(knitr)
-> knit2pdf("Medulloblastoma.Rnw")
-```
-This should generate `Medulloblastoma.pdf` in the `examples/` directory.
-
-## Known issues with compiling pdfs
-
-#### (Linux)
-If you are getting an error saying that R cannot find `/usr/bin/texti2dvi`, install the `texinfo` and 'texlive' packages in linux. You can use the commands:
-apt-get install texinfo
-apt-get install texlive
-
-#### (OS/X): "`pdfLaTex` not found" error
-When compiling the pdf, you may get a message saying that `pdfLaTex` is not installed. We have had one such report on OS/X and it is known to occur after an upgrade to OS X Mavericks. The following steps resolved the issue (paraphrased from [this post](http://stackoverflow.com/questions/22081991/rmarkdown-pandoc-pdflatex-not-found)):
-Step 1. Check that `/usr/texbin` exists. In terminal, type:
-```
-cd /usr/texbin
-```
-Step 2. If you get a "No such file or directory" message, create a symbolic link to your installation's `texbin` file. It may be in `/Library/TeX/Distributions/.DefaultTeX/Contents/Programs/texbin`.
-In terminal type:
-```
-ln -s /Library/TeX/Distributions/.DefaultTeX/Contents/Programs/texbin /usr/texbin
-```
-Step 3. Now, in the terminal check the value of `echo $PATH`. Make sure that `/usr/texbin` is present. If it isn't present, then you need to add `/usr/texbin` to your PATH variable. This can be done by updating the `PATH` variable in `~/.bashrc`.
-However, if you find yourself having to mess with the PATH variable, try reinstalling the [MacTex](http://tug.org/mactex/) package.
-
-## Run BreastCancer LumA example
-This vignette is presented in the netDx manuscript. Here we start with 348 primary tumours from the Cancer Genome Atlas, and build a predictor for Luminal A subtype classification (The Cancer Genome Atlas (2012). *Nature.* **490**:61-70). This example illustrates feature selection using a simple design in which networks are scored out of 10 based on a single round of 10-fold cross validation. On a MacBook Air laptop (late 2014), this vignette takes ~1.5 hours to run to completion. You may speed it up by running it on a machine with more processors and changing the `numCores` variable in the vignette. We do not recommend running it on a machine with less than 8Gb RAM.
-
-```
-$ cd netDx/examples/
-$ R
-> require(knitr)
-> knit2pdf("BreastCancer.Rnw")
-```
-**NOTE:** The vignette will generate a pdf file. All intermediate files will be stored in the `TCGA_BRCA/` subdirectory of the examples directory.
diff --git a/examples/BreastCancer.Rnw b/examples/BreastCancer.Rnw
index cf7e4b09..d73e72a9 100644
--- a/examples/BreastCancer.Rnw
+++ b/examples/BreastCancer.Rnw
@@ -2,6 +2,7 @@
\usepackage{caption}
\begin{document}
+\SweaveOpts{concordance=TRUE}
\title{netDx use case \linebreak Integrate gene expression and CNV for \linebreak binary classification of breast tumour}
\author{Shraddha Pai}
@@ -29,7 +30,7 @@ The workflow is shown in Figure 1. The algorithm proceeds in two steps:
\begin{figure}[ht]
\begin{center}
-\includegraphics[width=\textwidth]{tcga_brca.png}
+\includegraphics[width=\textwidth]{images/tcga_brca.png}
\caption{netDx workflow for a binary tumour classifier from gene expression and CNV data. \newline A. Two sets of patient similarity networks are generated: the first based on correlation of gene expression in cellular pathways (magenta), and the second based on shared overlap of CNVs in cellular pathways (teal). Each datatype generates ~1,000-2,000 networks, and these are integrated into a single database by GeneMANIA. \newline B. Feature selection is separately carried out for the `LumA' class for the `other' class. A GeneMANIA query is run on the integrated database is queried 10 times; each time a different 9/10th of the training ``+'' samples is used as query. A network's score is the frequency with which GeneMANIA marks it as being informative. Networks scoring 9 or 10 out of 10 are feature selected. Before patient classification, two enriched databases are constructed (orange and grey cylinders); each contain feature selected nets and train as well as test samples. \newline C. Patient similarity to a class is ranked by running a query against the class database; this is done once per class. Test patients are assigned to the class for which they have the highest-ranking similarity.}
\end{center}
\end{figure}
diff --git a/examples/Medulloblastoma.Rnw b/examples/Medulloblastoma.Rnw
index 7642cd0a..95998f0f 100644
--- a/examples/Medulloblastoma.Rnw
+++ b/examples/Medulloblastoma.Rnw
@@ -1,6 +1,7 @@
\documentclass{article}
\begin{document}
+\SweaveOpts{concordance=TRUE}
\title{netDx use case: 4-way classification: Medulloblastoma tumour subtype}
\author{Shraddha Pai}
@@ -15,7 +16,7 @@ The netDx workflow is shown in Figure 1. We use gene signatures identified by pr
\begin{figure}[ht]
\begin{center}
-\includegraphics[width=\textwidth]{mblastoma.png}
+\includegraphics[width=\textwidth]{images/mblastoma.png}
\caption{netDx workflow for a 4-way classification of medulloblastoma tumour from known gene signatures. \newline A. A GeneMANIA database is built for each subtype. The networks are at the gene-level, with a custom similarity metric (see Section 3); each subtype-specific database contains networks for the genes in the corresponding subtype signature. The networks and database contain all training and test patients. \newline B. A sample is classified by running four GeneMANIA queries, one per subtype. In each case, the query comprises of training samples for that subtype. The query results in subtype-specific similarity ranking. After such ranks have been obtained for all four subtypes, the sample is assigned to the class for which it has the highest rank.}
\end{center}
\end{figure}
diff --git a/examples/NestedCV_MultiData.Rmd b/examples/NestedCV_MultiData.Rmd
new file mode 100644
index 00000000..28cd16e8
--- /dev/null
+++ b/examples/NestedCV_MultiData.Rmd
@@ -0,0 +1,227 @@
+---
+title: "Nested Cross-Validation with Multiple Datatypes"
+output:
+ html_document:
+ toc: true
+ toc_depth: 2
+ html_notebook:
+ toc: true
+ toc_depth: 2
+---
+
+# Introduction
+
+In this example, we will use netDx to run **nested cross-validation design** using multiple datatypes and custom similarity metrics.
+
+In the nested design, cross-validation is performed repeatedly over multiple random splits of the data into train and blind test partitions. Feature selected networks are those that consistently score highly across the multiple splits.
+
+Conceptually, this is what the higher-level logic looks like for a nested cross-validation design with 10-fold CV in the inner loop, and 100 splits in the outer loop:
+
+*(Note: these aren't real function calls; this block just serves to illustrate the concept of the nested CV design for our purposes)*
+```{r,eval=FALSE}
+outerLoop <- 100 # num times to split data into train/blind test samples
+innerLoop <- 10 # num folds for cross-validation, also max score for a network
+netScores <- list() # collect set of netScores
+perf <- list() # collect set of test evaluations
+
+for k in 1:outerLoop
+ [train, test] <- splitData(80:20) # split data using RNG seed
+ netScores[[k]] <- runCV(train)
+ perf[[k]] <- collectPerformance(netScores[[k]], test)
+end
+```
+
+# Setup
+
+```{r}
+suppressWarnings(suppressMessages(require(netDx)))
+suppressWarnings(suppressMessages(require(netDx.examples)))
+```
+
+# Data
+
+In this example, we use data from The Cancer Genome Atlas (http://cancergenome.nih.gov/), downloaded from the PanCancer Survival project (https://www.synapse.org/#!Synapse:syn1710282). We integrate gene expression profiles and clinical variables from renal clear cell carcinoma tumours to predict poor and good survival (Refs 1-2). The data are from 150 tumours and the classes are `SURVIVEYES` and `SURVIVENO`.
+
+```{r,eval=T}
+load(sprintf("%s/extdata/nestedCV_input.rda",
+ path.package("netDx.examples")))
+head(pheno)
+```
+
+# Design custom patient similarity networks (features)
+
+netDx allows the user to define a custom function that takes patient data and variable groupings as input, and returns a set of patient similarity networks (PSN) as output. The user can customize what datatypes are used, how they are grouped, and what defines patient similarity for a given datatype. When running the predictor (next section), the user simply passes this custom function as an input variable (the `makeNetFunc` parameter).
+
+***Note:*** While netDx provides a high degree of flexibility in achieving your design of choice, it is up to the user to ensure that the design, i.e. the similarity metric and variable groupings, is appropriate for your application. Domain knowledge is almost likely required for good design.
+
+netDx requires that this function take some generic parameters as input. These include:
+
+* `dataList`: the patient data
+* `groupList`: sets of input data that would correspond to individual networks (e.g. genes grouped into pathways)
+* `netDir`: the directory where the resulting PSN would be stored.
+
+This section provides more details on the `dataList` and `groupList` variables.
+
+### dataList
+This contains the input patient data for the predictor.
+Each key is a datatype, while each value is the corresponding data matrix. Note that columns are patients and rows are unit names (e.g. genes for rna, or variable names for clinical data).
+
+**Important:** The software expects the patient order in the columns to match the row order in the `pheno` table.
+
+The names are datatypes:
+
+```{r,eval=TRUE}
+names(dats)
+```
+
+The `rna` container has the patient gene expression matrix:
+
+```{r,eval=TRUE}
+head(dats[["rna"]][,1:6])
+```
+
+The `clinical` entry has the clinical data matrix:
+
+```{r,eval=TRUE}
+head(dats[["clinical"]][,1:6])
+```
+
+
+Look at the number of units per datatype:
+```{r,eval=TRUE}
+lapply(dats, nrow)
+```
+
+
+### groupList
+This object tells the predictor how to group units when constructing a network. For examples, genes may be grouped into a network representing a pathway. This object is a list; the names match those of `dataList` while each value is itself a list and reflects a potential network.
+
+```{r,eval=TRUE}
+names(groupList)
+```
+
+For example, here are the networks to be created with RNA data. Genes corresponding to pathways are to be grouped into individual network. Such a groupList would create pathway-level networks:
+
+```{r,eval=TRUE}
+groupList[["rna"]][1:3]
+```
+
+
+For clinical data, we want to keep each variable as its own network:
+
+```{r,eval=TRUE}
+head(groupList[["clinical"]])
+```
+
+To speed up this example, reduce the number of input networks created by limiting to three pathways *(prediction may not be great!)*:
+```{r,eval=TRUE}
+groupList[["rna"]] <- groupList[["rna"]][1:3]
+```
+
+## Define patient similarity for each network
+This function is defined by the user and tells the predictor how to create networks from the provided input data.
+
+**This function must take `dataList`,`groupList`, and `netDir` as input variables.** The residual `...` parameter is to pass additional variables to `makePSN_NamedMatrix()`, notably `numCores`.
+
+In this particular example, the custom similarity function does the following:
+
+1. Creates *pathway-level networks from RNA* data using the default Pearson correlation measure `makePSN_NamedMatrix(writeProfiles=TRUE,...)`
+2. Creates *variable-level networks from clinical* data using a custom similarity function of normalized difference: `makePSN_NamedMatrix(writeProfiles=FALSE,simMetric="custom",customFunc=normDiff)`.
+
+
+```{r}
+KIRC_makeNets <- function(dataList, groupList, netDir,...) {
+ netList <- c()
+ # make RNA nets: group by pathway
+ if (!is.null(groupList[["rna"]])) {
+ netList <- makePSN_NamedMatrix(dataList$rna,
+ rownames(dataList$rna),
+ groupList[["rna"]],netDir,verbose=FALSE,
+ writeProfiles=TRUE,...)
+ netList <- unlist(netList)
+ cat(sprintf("Made %i RNA pathway nets\n", length(netList)))
+ }
+
+ # make clinical nets,one net for each variable
+ netList2 <- c()
+ if (!is.null(groupList[["clinical"]])) {
+ netList2 <- makePSN_NamedMatrix(dataList$clinical,
+ rownames(dataList$clinical),
+ groupList[["clinical"]],netDir,
+ simMetric="custom",customFunc=netDx::normDiff, # custom function
+ sparsify=TRUE,verbose=TRUE,append=TRUE,...)
+ }
+ netList2 <- unlist(netList2)
+ cat(sprintf("Made %i clinical nets\n", length(netList2)))
+ netList <- c(netList,netList2)
+ cat(sprintf("Total of %i nets\n", length(netList)))
+ return(netList)
+}
+```
+
+
+**Note:** `dataList` and `groupList` are generic containers that can contain whatever object the user requires to create PSN. **The custom function gives the user complete flexibility in feature design**. For instance, a dataset with mutation data may pass `GenomicRanges` objects in `dataList`, to represent patient mutations, and `groupList` would contain `GenomicRanges` corresponding to loci to be grouped together in the resulting network.
+
+# Run nested cross-validation
+
+Finally we call the function that runs the netDx predictor. We provide:
+
+* number of train/test splits (outer loop of nested CV: `numSplits` parameter),
+* number of folds for cross-validation (inner loop: `nFoldCV`),
+* cutoff for networks within an inner loop (`CVcutoff`),
+* and the information to create the PSN, including patient data (`dataList`), how variables are to be grouped into networks (`groupList`) and the custom function to generate features (`makeNetFunc`).
+
+Running the below takes a lot of time so we have commented it out. Feel free to uncomment and run. Change `numCores` to match the number of cores available on your machine for parallel processing.
+
+The call below runs nested two-fold cross-validation over 2 splits. Within each split, it:
+
+* splits data into train/blind test using the default split of 80:20
+* runs 2-fold cross-validation using the training samples
+* uses networks that score >=1 out of 2 (`CVcutoff`) to classify blind test samples.
+
+**!!! Note:**Setting `numCores`>1 will not work in Rstudio. To increase `numCores`, use `knitr::purl("NestedCV_MultiData.Rmd")` to extract the R code from this notebook into its own file, change `numCores` and run in an interactive R client (i.e. by calling "R" on command-line). **!!!**
+
+
+```{r,eval=TRUE}
+runPredictor_nestedCV(pheno,
+ dataList=dats,groupList=groupList,
+ makeNetFunc=KIRC_makeNets, ### custom network creation function
+ outDir=sprintf("%s/nestedCV_output",getwd()), ## absolute path
+ numCores=1L,nFoldCV=2L, CVcutoff=1L,numSplits=2L)
+```
+
+# Examine results directory
+The results directory contains the log file (`log.txt`) and results for each split are stored in its own directory (here, `rng1/` and `rng2/`).
+
+```{r,eval=TRUE}
+outDir <- sprintf("%s/nestedCV_output",getwd())
+dir(outDir)
+```
+
+Look into one of the split directories. `predictionResults.txt` contains blind test prediction labels for that split. Each class has its own directory with corresponding cross-validation results.
+
+```{r}
+dir(sprintf("%s/rng1",outDir))
+pred <- read.delim(sprintf("%s/rng1/predictionResults.txt",outDir),h=T,as.is=T)
+head(pred)
+```
+
+The cross-validation results are in the `/GM_results/` directory. The set of network scores for the split are in `/GM_results/_pathway_CV_score.txt`.
+
+```{r}
+dir(sprintf("%s/rng1/SURVIVEYES",outDir))
+sc <- read.delim(sprintf("%s/rng1/SURVIVEYES/GM_results/SURVIVEYES_pathway_CV_score.txt",outDir),
+ sep="\t",h=T,as.is=T)
+head(sc)
+```
+
+
+# sessionInfo
+
+```{r}
+sessionInfo()
+```
+
+# References
+1. Yuan, Y. *et al.* (2014) Assessing the clinical utility of cancer genomic and proteomic data across tumor types. *Nat Biotechnol* **32**, 644-52.
+2. The Cancer Genome Atlas Research Network (2013). Comprehensive molecular characterization of clear cell renal cell carcinoma. *Nature* **499**, 43-9.
\ No newline at end of file
diff --git a/examples/mblastoma.png b/examples/images/mblastoma.png
similarity index 100%
rename from examples/mblastoma.png
rename to examples/images/mblastoma.png
diff --git a/examples/tcga_brca.png b/examples/images/tcga_brca.png
similarity index 100%
rename from examples/tcga_brca.png
rename to examples/images/tcga_brca.png
diff --git a/examples/KIRCsurvivalPred_nestedCV_multiData.R b/examples/other/KIRCsurvivalPred_nestedCV_multiData.R
similarity index 100%
rename from examples/KIRCsurvivalPred_nestedCV_multiData.R
rename to examples/other/KIRCsurvivalPred_nestedCV_multiData.R
diff --git a/examples/TCGA_BRCA.R b/examples/other/TCGA_BRCA.R
similarity index 100%
rename from examples/TCGA_BRCA.R
rename to examples/other/TCGA_BRCA.R
diff --git a/examples/other/example_multiDataset.R b/examples/other/example_multiDataset.R
new file mode 100644
index 00000000..99092c65
--- /dev/null
+++ b/examples/other/example_multiDataset.R
@@ -0,0 +1,396 @@
+#' This is an example of advanced predictor design with netDx. It illustrates
+#' 1) nested cross-validation predictor design
+#' 2) creating nets with multiple similarity measures
+#' 3) a function to compute mean normalized difference when a datasource has
+#' 2-5 variables. See function normDiff2()
+#'
+#' Context
+#' ----------------------
+#' Predict binarized survival for kidney renal clear cell carcinoma (KIRC)
+#' using 5 types of data (clinical, RNA, DNA methylation, proteomic,
+#' miRNA).
+#' Feature design:
+#' 1) For each datatype create a single similarity network
+#' 2) In addition, using just RNA data, create pathway-level networks.
+#'
+#' Predictor design:
+#' Uses nested cross validation. The outer loop splits data 100 times
+#' (numSplits variable) into
+#' train and blind test. The inner loop first performs 10-fold cross validation
+#' on the specific train, and classifies the test.
+#' Therefore the output is:
+#' 1) 100 predictionResults.txt files, each with classification for the 100
+#' test samples
+#' 2) 100 sets of network scores, each out of 10.
+#'
+#' The rule can then be set to call networks that score consistently well
+#' as being "feature selected". e.g. Networks that score 10/10 100% of the time
+#'
+#' Output: See Example 3 of outputs.md in the main folder of the github repo
+
+# ----------------------------------------------------------------------
+# Setup
+# ----------------------------------------------------------------------
+
+#' 10-fold CV predictor design
+rm(list=ls())
+require(netDx)
+require(netDx.examples)
+
+numCores <- 8L # parallel processing cores
+GMmemory <- 4L # memory in Gigabytes for each fold of cross-validation to use
+trainProp <- 0.8 # % split for train/test
+cutoff <- 9 # score for inner cross-validation fold that decides which
+ # nets get used within a fold for classifying a test sample
+numSplits <- 100
+
+### CHANGE THIS to your ROOT FOLDER
+rootDir <- "/mnt/data2/BaderLab"
+inDir <- sprintf("%s/PanCancer_KIRC/input",rootDir)
+outRoot <- sprintf("%s/PanCancer_KIRC/output",rootDir)
+
+dt <- format(Sys.Date(),"%y%m%d")
+megaDir <- sprintf("%s/AllPlusPathways_%s",outRoot,dt)
+
+# ----------------------------------------------------------------
+# helper functions
+# takes average of normdiff of each row in x
+normDiff2 <- function(x) {
+ # normalized difference
+ # x is vector of values, one per patient (e.g. ages)
+ normDiff <- function(x) {
+ #if (nrow(x)>=1) x <- x[1,]
+ nm <- colnames(x)
+ x <- as.numeric(x)
+ n <- length(x)
+ rngX <- max(x,na.rm=T)-min(x,na.rm=T)
+
+ out <- matrix(NA,nrow=n,ncol=n);
+ # weight between i and j is
+ # wt(i,j) = 1 - (abs(x[i]-x[j])/(max(x)-min(x)))
+ for (j in 1:n) out[,j] <- 1-(abs((x-x[j])/rngX))
+ rownames(out) <- nm; colnames(out)<- nm
+ out
+ }
+
+ sim <- matrix(0,nrow=ncol(x),ncol=ncol(x))
+ for (k in 1:nrow(x)) {
+ tmp <- normDiff(x[k,,drop=FALSE])
+ sim <- sim + tmp
+ rownames(sim) <- rownames(tmp)
+ colnames(sim) <- colnames(tmp)
+ }
+ sim <- sim/nrow(x)
+ sim
+}
+
+# -----------------------------------------------------------
+# Process input
+# -----------------------------------------------------------
+inFiles <- list(
+ clinical=sprintf("%s/KIRC_clinical_core.txt",inDir),
+ survival=sprintf("%s/KIRC_binary_survival.txt",inDir)
+ )
+datFiles <- list(
+ rna=sprintf("%s/KIRC_mRNA_core.txt",inDir),
+ prot=sprintf("%s/KIRC_RPPA_core.txt",inDir),
+ mir=sprintf("%s/KIRC_miRNA_core.txt",inDir),
+ dnam=sprintf("%s/KIRC_methylation_core.txt",inDir),
+ cnv=sprintf("%s/KIRC_CNV_core.txt",inDir)
+)
+
+pheno <- read.delim(inFiles$clinical,sep="\t",h=T,as.is=T)
+colnames(pheno)[1] <- "ID"
+
+#======transform clinical data=========
+pheno$grade <- as.vector(pheno$grade)
+pheno$grade[pheno$grade=="G1"] <- "G2"
+pheno$grade[pheno$grade=="GX"] <- "G2"
+pheno$grade <- as.factor(pheno$grade)
+pheno <- pheno[, -which(colnames(pheno)=="gender")]
+#======================================
+
+surv <- read.delim(inFiles$survival,sep="\t",h=T,as.is=T)
+colnames(surv)[1:2] <- c("ID","STATUS_INT")
+survStr <- rep(NA,nrow(surv))
+survStr[surv$STATUS_INT<1] <- "SURVIVENO"
+survStr[surv$STATUS_INT>0] <- "SURVIVEYES"
+surv$STATUS <- survStr
+pheno <- merge(x=pheno,y=surv,by="ID")
+pheno$X <- NULL
+# pheno$gender <- ifelse(pheno$gender=="FEMALE",1, 0)
+pheno_nosurv <- pheno[1:4]
+
+cat("Collecting patient data:\n")
+dats <- list() #input data in different slots
+cat("\t* Clinical\n")
+clinical <- pheno_nosurv
+rownames(clinical) <- clinical[,1];
+clinical$grade <- as.numeric(factor(clinical$grade))
+clinical$stage <- as.numeric(factor(clinical$stage))
+clinical$ID <- NULL
+clinical <- t(clinical)
+dats$clinical <- clinical; rm(clinical)
+
+# create master input net
+for (nm in names(datFiles)) {
+ cat(sprintf("\t* %s\n",nm))
+ tmp <- read.delim(datFiles[[nm]],sep="\t",h=T,as.is=T)
+ if (colnames(tmp)[ncol(tmp)]=="X") tmp <- tmp[,-ncol(tmp)]
+ rownames(tmp) <- tmp[,1]
+ tmp <- t(tmp[,-1])
+ class(tmp) <- "numeric"
+ dats[[nm]] <- tmp
+}
+
+cat("\t Ordering column names\n")
+# include only data for patients in classifier
+dats <- lapply(dats, function(x) { x[,which(colnames(x)%in%pheno$ID)]})
+dats <- lapply(dats, function(x) {
+ midx <- match(pheno$ID,colnames(x))
+ x <- x[,midx]
+ x
+})
+
+# confirm patient order the same for all input nets
+pname <- colnames(dats[[1]])
+for (k in 2:length(dats)) {
+ if (all.equal(colnames(dats[[k]]),pname)!=TRUE) {
+ cat(sprintf("Patient order doesn't match for %s\n",
+ names(dats)[k]))
+ browser()
+ }
+}
+
+# input nets for each category
+netSets <- lapply(dats, function(x) rownames(x))
+
+# compile data
+alldat <- do.call("rbind",dats)
+pheno_all <- pheno
+
+combList <- list(
+ clinical="clinical_cont",
+ mir="mir.profile",
+ rna="rna.profile",
+ prot="prot.profile",
+ cnv="cnv.profile",
+ dnam="dnam.profile",
+ clinicalArna=c("clinical_cont","rna.profile"),
+ clinicalAmir=c("clinical_cont","mir.profile"),
+ clinicalAprot=c("clinical_cont","prot.profile"),
+ clinicalAdnam=c("clinical_cont","dnam.profile"),
+ clinicalAcnv=c("clinical_cont","cnv.profile"),
+ all="all")
+
+# now add pathways
+pathFile <- sprintf("%s/extdata/Human_160124_AllPathways.gmt",
+ path.package("netDx.examples"))
+pathwayList <- readPathways(pathFile)
+
+rm(pheno,pheno_nosurv)
+
+# ----------------------------------------------------------
+# Build classifier
+# ----------------------------------------------------------
+if (file.exists(megaDir)) unlink(megaDir,recursive=TRUE)
+dir.create(megaDir)
+
+logFile <- sprintf("%s/log.txt",megaDir)
+sink(logFile,split=TRUE)
+tryCatch({
+## Create the mega database with all patients and all nets.
+## This will be used to predict test samples by subsetting just for feature
+## selected nets in a given round
+## Note that this is useful for all train/test splits because we can always
+## change which samples are query and can always subset based on which nets
+## are feature selected in a given round.
+netDir <- sprintf("%s/networks",megaDir)
+nonclin <- setdiff(names(netSets),"clinical")
+
+# Here we make three calls to makePSN_NamedMatrix.
+# 1) Pearson corr (writeProfiles=TRUE) for non-clinical data
+# 2) normDiff2 (writeProfiles=FALSE) for other types of data
+# 3) Pearson corr (writeProfiles=TRUE) for pathway data
+# Note
+# 1) append=TRUE for second and third.
+# 2) when classifying test you need to repeat the exercise of creating
+# three types of input networks, depending on whether networks from each
+# class scored well.
+
+netList <- makePSN_NamedMatrix(alldat,
+ rownames(alldat),netSets[nonclin],netDir,
+ verbose=FALSE,numCores=numCores,writeProfiles=TRUE)
+
+# notice append=TRUE
+netList2 <- makePSN_NamedMatrix(alldat,
+ rownames(alldat),netSets["clinical"],
+ netDir,simMetric="custom",customFunc=normDiff2,
+ verbose=FALSE,numCores=numCores,
+ sparsify=TRUE,append=TRUE)
+
+# group by pathway
+netList3 <- makePSN_NamedMatrix(dats_train$rna, rownames(dats_train$rna),
+ pathwayList,netDir,verbose=FALSE,
+ numCores=numCores,writeProfiles=TRUE,append=TRUE)
+cat(sprintf("Made %i RNA pathway nets\n", length(netList)))
+
+netList <- c(netList,netList2,netList3)
+cat(sprintf("Total of %i nets\n", length(netList)))
+
+# now create database
+megadbDir <- GM_createDB(netDir, pheno_all$ID, megaDir,numCores=numCores)
+
+# -------------------------------------------------------
+# Outer loop of nested CV: Different train/test splits
+# -------------------------------------------------------
+for (rngNum in 1:numSplits) {
+
+ rng_t0 <- Sys.time()
+ cat(sprintf("-------------------------------\n"))
+ cat(sprintf("RNG seed = %i\n", rngNum))
+ cat(sprintf("-------------------------------\n"))
+ outDir <- sprintf("%s/rng%i",megaDir,rngNum)
+ dir.create(outDir)
+
+ pheno_all$TT_STATUS <- splitTestTrain(pheno_all,pctT=trainProp,
+ setSeed=rngNum*5)
+ write.table(pheno_all,file=sprintf("%s/tt_split.txt",outDir),sep="\t",
+ col=T,row=F,quote=F)
+
+ # --------------------------------------------
+ # feature selection - train only
+ pheno <- subset(pheno_all, TT_STATUS %in% "TRAIN")
+ alldat_train <- alldat[,which(colnames(alldat) %in% pheno$ID)]
+
+ netDir <- sprintf("%s/networks",outDir)
+ nonclin <- setdiff(names(netSets),"clinical")
+
+ # Here we limit ourselves to making nets with training samples
+ netList <- makePSN_NamedMatrix(alldat_train,
+ rownames(alldat_train),netSets[nonclin],
+ netDir,verbose=FALSE,numCores=numCores,
+ writeProfiles=TRUE)
+ netList2 <- makePSN_NamedMatrix(alldat_train,
+ rownames(alldat_train),netSets["clinical"],
+ netDir,simMetric="custom",customFunc=normDiff2,
+ verbose=FALSE,numCores=numCores,
+ sparsify=TRUE,append=TRUE)
+ netList <- c(netList,netList2)
+ cat(sprintf("Total of %i nets\n", length(netList)))
+
+ # now create database
+ dbDir <- GM_createDB(netDir, pheno$ID, outDir,numCores=numCores)
+
+ # Ignore this loop -- it's example-specific and you will likely not need
+ # this. It is only to see the effect of starting with different combinations
+ # of input data
+ for (cur in names(combList)) {
+ t0 <- Sys.time()
+ cat(sprintf("%s\n",cur))
+ pDir <- sprintf("%s/%s",outDir, cur)
+ dir.create(pDir)
+
+ # run featsel once per subtype
+ subtypes <- unique(pheno$STATUS)
+
+ # -------------------------------------------------------
+ # Inner loop of nested CV: 10-fold CV
+ # -------------------------------------------------------
+ for (g in subtypes) {
+ pDir2 <- sprintf("%s/%s",pDir,g)
+ if (file.exists(pDir2)) unlink(pDir2,recursive=TRUE)
+ dir.create(pDir2)
+
+ cat(sprintf("\n******\nSubtype %s\n",g))
+ pheno_subtype <- pheno
+ ## label patients not in the current class as residual
+ nong <- which(!pheno_subtype$STATUS %in% g)
+ pheno_subtype$STATUS[nong] <- "nonpred"
+ ## sanity check
+ print(table(pheno_subtype$STATUS,useNA="always"))
+ resDir <- sprintf("%s/GM_results",pDir2)
+ ## query for feature selection comprises of training
+ ## samples from the class of interest
+ trainPred <- pheno_subtype$ID[
+ which(pheno_subtype$STATUS %in% g)]
+
+ # Cross validation
+ GM_runCV_featureSet(trainPred, resDir, dbDir$dbDir,
+ nrow(pheno_subtype),incNets=combList[[cur]],
+ verbose=T, numCores=numCores,
+ GMmemory=GMmemory)
+
+ # patient similarity ranks
+ prank <- dir(path=resDir,pattern="PRANK$")
+ # network ranks
+ nrank <- dir(path=resDir,pattern="NRANK$")
+ cat(sprintf("Got %i prank files\n",length(prank)))
+
+ # Compute network score
+ pTally <- GM_networkTally(paste(resDir,nrank,sep="/"))
+ head(pTally)
+ # write to file
+ tallyFile <- sprintf("%s/%s_pathway_CV_score.txt",resDir,g)
+ write.table(pTally,file=tallyFile,sep="\t",col=T,row=F,quote=F)
+ }
+
+ predRes <- list()
+ for (g in subtypes) {
+ pDir2 <- sprintf("%s/%s",pDir,g)
+ # get feature selected net names
+ pTally <- read.delim(
+ sprintf("%s/GM_results/%s_pathway_CV_score.txt",pDir2,g),
+ sep="\t",h=T,as.is=T)
+
+ # feature selected nets pass cutoff threshold
+ pTally <- pTally[which(pTally[,2]>=cutoff),1]
+ cat(sprintf("%s: %i pathways\n",g,length(pTally)))
+
+ # query of all training samples for this class
+ qSamps <- pheno_all$ID[which(pheno_all$STATUS %in% g &
+ pheno_all$TT_STATUS%in%"TRAIN")]
+
+ qFile <- sprintf("%s/%s_query",pDir2,g)
+ GM_writeQueryFile(qSamps,incNets=pTally,nrow(pheno_all),qFile)
+ resFile <- runGeneMANIA(megadbDir$dbDir,qFile,resDir=pDir2)
+ predRes[[g]] <- GM_getQueryROC(sprintf("%s.PRANK",resFile),
+ pheno_all,g)
+ }
+
+ # -------------------------------------------------------
+ # Predict blind test for the current split
+ # -------------------------------------------------------
+ predClass <- GM_OneVAll_getClass(predRes)
+ out <- merge(x=pheno_all,y=predClass,by="ID")
+ outFile <- sprintf("%s/predictionResults.txt",pDir)
+ write.table(out,file=outFile,sep="\t",col=T,row=F,quote=F)
+
+ acc <- sum(out$STATUS==out$PRED_CLASS)/nrow(out)
+ cat(sprintf("Accuracy on %i blind test subjects = %2.1f%%\n",
+ nrow(out), acc*100))
+
+ require(ROCR)
+ ROCR_pred <- prediction(out$SURVIVEYES_SCORE-out$SURVIVENO,
+ out$STATUS=="SURVIVEYES")
+ save(predRes,ROCR_pred,file=sprintf("%s/predRes.Rdata",pDir))
+ }
+
+ #cleanup to save disk space
+ system(sprintf("rm -r %s/dataset %s/tmp %s/networks",
+ outDir,outDir,outDir))
+ system(sprintf("rm -r %s/dataset %s/networks",
+ outDir,outDir))
+
+}
+ pheno_all$TT_STATUS <- NA
+
+ rng_t1 <- Sys.time()
+ cat(sprintf("Time for one train/test split:"))
+ print(rng_t1-rng_t0)
+
+}, error=function(ex){
+ print(ex)
+}, finally={
+ sink(NULL)
+})
diff --git a/examples/plotPredictorResults.Rmd b/examples/plotPredictorResults.Rmd
new file mode 100644
index 00000000..4c52e3b4
--- /dev/null
+++ b/examples/plotPredictorResults.Rmd
@@ -0,0 +1,315 @@
+---
+title: "Plot netDx results"
+output:
+ html_document:
+ toc: true
+ toc_depth: 2
+ html_notebook:
+ toc: true
+ toc_depth: 2
+ pdf_document: default
+ word_document: default
+---
+
+# Introduction
+After building a netDx predictor, you will want to evaluate the performance of the predictor and examine the nature of feature-selected variables. In this example, we use built-in functionality to generate the following output:
+
+* **Summary predictor performance:** Plot of average performance (AUROC/AUPR)
+* **Detailed predictor performance:** AUROC/AUPR curves for all data splits
+* **Detailed statistics:** Network scores and patient predictions in all splits
+* **Themes of selected features:** A network visualization of major themes in predictive variables (an EnrichmentMap)
+* A **network visualization of patient similarity** using only feature-selected variables
+
+# Software Requirements
+For this example to work, you will require the following components to be installed. You may have already installed these when you first setup netDx but if not, make sure these are set up now:
+
+* Cytoscape must be installed and running.
+* Cytoscape apps EnrichmentMap (v3.0.0) and AutoAnnotate (v1.2) must be installed.
+* R packages: `httr`, `RJSONIO`
+* R packages r2cytoscape and EasycyRest.
+
+**Note: If Cytoscape is not running, this example will not work!**
+
+Let us see if the required dependencies can be installed and/or loaded:
+
+```{r}
+#httr
+tryCatch(expr = { library(httr)},
+ error = function(e) { install.packages("httr")}, finally = library(httr))
+
+#RJSONIO
+tryCatch(expr = { library(RJSONIO)},
+ error = function(e) { install.packages("RJSONIO")}, finally = library(RJSONIO))
+
+#r2cytoscape
+tryCatch(expr = { library(r2cytoscape)},
+ error = function(e) { devtools::install_github('cytoscape/cytoscape-automation/for-scripters/R/r2cytoscape')}, finally = library(r2cytoscape))
+
+# EasycyRest
+tryCatch(expr = { library(EasycyRest); detach(package:EascyRest,unload=TRUE)},
+ error = function(e) { devtools::install_github('BaderLab/Easycyrest/EasycyRest@0.1')}, finally = {})
+```
+
+Now let us test if Cytoscape is open:
+
+```{r}
+# Basic settings
+version.url <-"http://localhost:1234/v1/version"
+cytoscape.open <- TRUE
+
+tryCatch(expr = { GET(version.url)},
+ error = function(e) { return (cytoscape.open = FALSE)}, finally =function(r){ return(cytoscape.open = TRUE)})
+
+if(!cytoscape.open){
+ #try and launch cytoscape
+ stop("Cytoscape is not open. Please launch cytoscape.")
+} else{
+ cytoscape.version <- GET(version.url)
+ cy.version <- fromJSON(rawToChar(cytoscape.version$content))
+ print(cy.version)
+}
+```
+
+# Predictor Design
+
+For this example, the predictor has been run using a **nested cross-validation design**. In nested cross-validation, cross-validation is performed repeatedly over multiple random splits of the data into train and blind test partitions. Feature selected networks are those that consistently score highly across the multiple splits.
+Conceptually, this is what the higher-level logic looks like for a nested cross-validation design with 10-fold CV in the inner loop, and 100 splits in the outer loop:
+
+*(Note: these aren't real function calls; this block just serves to illustrate the concept of the nested CV design for our purposes)*
+```{r,eval=FALSE}
+outerLoop <- 100 # num times to split data into train/blind test samples
+innerLoop <- 10 # num folds for cross-validation, also max score for a network
+netScores <- list() # collect set of netScores
+perf <- list() # collect set of test evaluations
+
+for k in 1:outerLoop
+ [train, test] <- splitData(80:20) # split data using RNG seed
+ netScores[[k]] <- runCV(train)
+ perf[[k]] <- collectPerformance(netScores[[k]], test)
+end
+```
+
+# Output directory structure
+netDx expects a nested directory structure with the predictor results. The top level should contain one directory for each train/test split. Within each of these directories are the predictor results for the corresponding cross-validation. Here is the directory structure for a dataset with rootDirectory `dataset_yymmdd`.
+
+```
+dataset_yymmdd/
+ + rng1/
+ + tmp/ # directory created by netDx, containing input data for GeneMANIA database
+ + networks/ # PSN created by calls to makePSN_NamedMatrix()
+ +-- Class1
+ + tmp/
+ + networks/ # networks for test classification for this split
+ + GM_results/ # results of inner loop (10-fold CV)
+ + Class1_pathway_CV_score.txt # network scores for inner CV fold
+ +--- CV_1.query # query for CV fold
+ +--- CV_1.query-results.report.txt.NRANK # network weights for CV fold
+ ...
+ +--- CV_10.query
+ +--- CV_10.query-results.report.txt.NRANK
+ +-- Class2
+ + predictionResults.txt # test predictions for this train/test split
+ + rng2/
+ + rng3/
+ + rng4/
+ ...
+ + rng100/
+```
+# Set up
+
+```{r,eval=TRUE}
+suppressMessages(require(netDx))
+suppressMessages(require(netDx.examples))
+```
+
+## Load data for plotting
+In this example, we use data from The Cancer Genome Atlas (http://cancergenome.nih.gov/), downloaded from the PanCancer Survival project (https://www.synapse.org/#!Synapse:syn1710282). We use gene expression profiles from renal clear cell carcinoma tumours to predict poor and good survival after Yuan et al. (2014) (Refs 1-2). The data consists of 150 tumours. Here we work only with the gene expression profiles generated.
+
+```{r,eval=TRUE}
+phenoFile <- sprintf("%s/extdata/KIRC_pheno.rda",path.package("netDx.examples"))
+lnames <- load(phenoFile)
+head(pheno)
+```
+
+Create a directory to store output in:
+```{r}
+outDir <- paste(getwd(),"plots",sep="/")
+if (!file.exists(outDir)) dir.create(outDir)
+setwd(outDir)
+```
+
+Now compile the paths to output files. First get the list of all directories corresponding to the outer loops:
+
+```{r}
+inDir <- sprintf("%s/extdata/KIRC_output",
+ path.package("netDx.examples"))
+all_rngs <- list.dirs(inDir, recursive = FALSE)
+print(head(basename(all_rngs)))
+```
+
+
+Each `rngX/` directory contains the results of a particular train/test split. The two classes in question are `SURVIVEYES` (good survival) and `SURVIVENO` (poor survival).
+
+```{r,eval=TRUE}
+dir(all_rngs[1])
+```
+
+# Plot overall predictor performance
+First we evaluate the average performance of the predictor over the 100 train/blind test splits. Performance is defined as the Area under the Receiver Operator Characteristic curve (AUROC) or Area under the Precision-Recall curve (AUPR).
+
+First compile prediction files for all the 100 splits in the example data. These paths are stored in `predFiles`. Then call the performance-plotting function, which is `plotPerf()`.
+
+```{r,eval=TRUE,fig.width=6,fig.height=6}
+predClasses <- c("SURVIVEYES","SURVIVENO")
+predFiles <- unlist(lapply(all_rngs, function(x)
+ paste(x, "predictionResults.txt", sep = "/")))
+predPerf <- plotPerf(inDir, predClasses=predClasses)
+```
+
+
+# Compile feature-level scores
+Feature-level scores are one of the outputs of feature-selection. These indicate the predictive strength of a feature, so that networks with higher scores indicate those with greater predictive power for the class in question. Note that in netDx, feature selection occurs once per patient class (i.e. here, once for `SURVIVEYES` and once for `SURVIVENO`), so each feature gets a predictive score for each class.
+
+This nested design has an outer loop of 100 splits with 10-fold cross-validation in the inner loop. Each network can therefore score between 0 and 10 (max number of cross-validation folds) for a given split. For the 100 splits, feature scores can be viewed as a matrix *X* which is *N*-by-100, where `N` is the number of networks. *X[i,j]* is the score for network *i* for split *j*.
+
+
+***Note:***In practice, we do not include networks that never score >0, so the number of rows is < *N*.
+
+```{r,eval=TRUE}
+featScores <- getFeatureScores(inDir,predClasses=c("SURVIVEYES","SURVIVENO"))
+```
+
+The size of `featScores` should correspond to number of networks that score >0
+at least once (rows) and the number of train/test splits (here, 100):
+
+```{r,eval=TRUE}
+dim(featScores[[1]])
+```
+The first column shows feature names:
+```{r,eval=TRUE}
+head(featScores[[1]][,1:10])
+```
+
+
+Let us define **feature-selected networks** as those that score 10 out of 10 in
+at least 70% of the train/test splits.
+```{r,eval=TRUE}
+featSelNet <- lapply(featScores, function(x) {
+ callFeatSel(x, fsCutoff=10, fsPctPass=0.7)
+})
+```
+
+Let's take a look at the top networks for each class:
+```{r,eval=TRUE}
+tmp <- lapply(featSelNet,print)
+```
+
+# Visualize functional themes in selected features (EnrichmentMap)
+An EnrichmentMap is a network-based visualization for a group of overlapping sets (Ref 3). Here, we use this visualization to examine which pathways and functional themes are feature-selected for each patient class. So, this is a network where nodes are **features** (here, pathways), and edges indicate **shared members between the features** (here, genes common to the pathway).
+
+By clustering and annotating the EnrichmentMap using the AutoAnnotate app (Ref 4), we can identify functional themes that are common among high-scoring features.
+
+Get the set of valid gene sets:
+```{r}
+pathFile <- sprintf("%s/extdata/Human_160124_AllPathways.gmt",
+ path.package("netDx.examples"))
+pathwayList <- readPathways(pathFile)
+```
+Filter for the genes measured in this dataset:
+```{r}
+xpr_genes <- sprintf("%s/extdata/EMap_input/genenames.txt",
+ path.package("netDx.examples"))
+xpr_genes <- read.delim(xpr_genes,h=FALSE,as.is=TRUE)[,1]
+head(xpr_genes)
+```
+Filter:
+```{r}
+pathwayList <- lapply(pathwayList, function(x) x[which(x %in% xpr_genes)])
+```
+
+To create an enrichment map in netDx, you need two files:
+
+1. **A .gmt file:** A file with top-scoring genesets in the GMT format, similar to the example pathway file.
+2. **Node attribute table:** A table containing the names of features and the maximum score each achieves across cross-validation.
+
+Here we generate the above two files for each patient class. We have two sets of files because each class has its own set of predictive features and therefore, its own enrichment map (EM).
+
+In addition to objects we have seen before, this step requires a table indicating what type of data each network represents. This may be used to assign visual features (e.g. node shapes) to distinguish different data sources in the EM.
+
+```{r}
+netInfoFile <- sprintf("%s/extdata/KIRC_output/inputNets.txt",
+ path.package("netDx.examples"))
+netInfo <- read.delim(netInfoFile,sep="\t",h=FALSE,as.is=TRUE)
+head(netInfo)
+```
+
+Create the EnrichmentMap input:
+```{r}
+EMap_input <- writeEMapInput_many(featScores,pathwayList,
+ netInfo,outDir=outDir)
+```
+
+Finally, plot the EnrichmentMap:
+```{r}
+pngFiles <- list()
+for (curGroup in names(EMap_input)[1:2]) {
+ pngFiles[[curGroup]] <- plotEmap(gmtFile=EMap_input[[curGroup]][1],
+ nodeAttrFile=EMap_input[[curGroup]][2],
+ netName=curGroup,outDir=outDir)
+}
+```
+
+The enrichment map showing the top-scoring networks for `SURVIVEYES` should now be active in Cytoscape. It should look like this:
+
+
+
+Nodes are coloured by the highest consistent score the particular feature achieved; consistency is defined as having achieved the score for >70% of the splits.
+
+**Note:** Node labels for individual features have been turned off to remove visual clutter. These can be enabled interactively in Cytoscape.
+
+There should also be a second EnrichmentMap for the `SURVIVENO` class:
+
+{width=70%}
+
+
+### Compute overall patient similarities
+An integrated patient network is constructed by aggregating feature-selected networks for all classes. It represents the overall view of patient similarity that results from identifying predictive networks through the netDx predictor.
+
+For the visualization we actually prefer to magnify **dissimilarity** so that similar patients are closer (i.e. have smaller edge weight) and dissimilar patients are farther apart (i.e. have larger edge weight). So the network generated in this step is really a **dissimilarity network**. The steps for computing this network are:
+
+1. Concatenate all feature-selected networks
+2. For all patient pairs, collapse edges by taking the mean similarity.
+3. Convert to dissimilarity (`1-mean_similarity`).
+4. For the Cytoscape visualization, keep a certain fraction of the top edges (e.g. top 20% disimilarities).
+
+The function `plotIntegratedPSN()` runs these computations. It therefore needs access to the folder that contains all the feature-selected networks. Here we provide `rng1/` for that purpose using the `baseDir` argument.
+
+```{r,eval=TRUE}
+netInfo <- plotIntegratedPSN(pheno=pheno,baseDir=sprintf("%s/rng1",inDir),
+ netNames=featSelNet,outDir=outDir)
+```
+The Cytoscape session should now show a network that looks like this:
+
+
+
+The returned `netInfo` object contains several pieces of information related to the integrated network shown above, including the path to the full similarity network (`aggPSN_FULL`), the pruned dissimilarity network (`aggPDN_pruned`), and information about the resulting view in Cytoscape (`network_suid` and `netView`).
+
+
+```{r,eval=TRUE}
+summary(netInfo)
+```
+
+# sessionInfo
+
+```{r}
+sessionInfo()
+```
+
+# References
+1. Yuan, Y. *et al.* (2014) Assessing the clinical utility of cancer genomic and proteomic data across tumor types. *Nat Biotechnol* **32**, 644-52.
+2. The Cancer Genome Atlas Research Network (2013). Comprehensive molecular characterization of clear cell renal cell carcinoma. *Nature* **499**, 43-9.
+3. Merico, D., Isserlin, R. & Bader, G.D. (2011). Visualizing gene-set enrichment results using the Cytoscape plug-in enrichment map. *Methods Mol Biol* **781**, 257-77.
+4. Kucera, M., Isserlin, R., Arkhangorodsky, A. & Bader, G.D. (2016). AutoAnnotate: A Cytoscape app for summarizing networks with semantic annotations. *F1000Res* **5**, 1717.
+
+
+
diff --git a/examples/plotPredictorResultsInOneShot.Rmd b/examples/plotPredictorResultsInOneShot.Rmd
new file mode 100644
index 00000000..27db11e6
--- /dev/null
+++ b/examples/plotPredictorResultsInOneShot.Rmd
@@ -0,0 +1,130 @@
+---
+title: "Plot netDx results: One function"
+output:
+ html_document:
+ toc: true
+ toc_depth: 2
+ html_notebook:
+ toc: true
+ toc_depth: 2
+ pdf_document: default
+ word_document: default
+---
+
+# Introduction
+This notebook shows you how to plot all netDx predictor results in one function call. The output generated includes:
+
+* **Summary predictor performance:** Plot of average performance (AUROC/AUPR)
+* **Detailed predictor performance:** AUROC/AUPR curves for all data splits
+* **Detailed statistics:** Network scores and patient predictions in all splits
+* **Themes of selected features:** A network visualization of major themes in predictive variables (an EnrichmentMap)
+* A **network visualization of patient similarity** using only feature-selected variables
+
+# Software Requirements
+Assumes you have all Cytoscape- and CyRest-related dependencies required to install netDx v1.0.
+
+**Note:If Cytoscape is not running, this example will not work!**
+
+Let us see if the required dependencies can be installed and/or loaded:
+
+```{r}
+#httr
+tryCatch(expr = { library(httr)},
+ error = function(e) { install.packages("httr")}, finally = library(httr))
+
+#RJSONIO
+tryCatch(expr = { library(RJSONIO)},
+ error = function(e) { install.packages("RJSONIO")}, finally = library(RJSONIO))
+
+#r2cytoscape
+tryCatch(expr = { library(r2cytoscape)},
+ error = function(e) { devtools::install_github('cytoscape/cytoscape-automation/for-scripters/R/r2cytoscape')}, finally = library(r2cytoscape))
+
+# EasycyRest
+tryCatch(expr = { library(EasycyRest); detach(package:EascyRest,unload=TRUE)},
+ error = function(e) { devtools::install_github('BaderLab/Easycyrest/EasycyRest@0.1')}, finally = {})
+```
+
+# Set up
+
+```{r,eval=TRUE}
+suppressWarnings(suppressMessages(require(netDx)))
+suppressWarnings(suppressMessages(require(netDx.examples)))
+```
+
+## Load data for plotting
+In this example, we use data from The Cancer Genome Atlas (http://cancergenome.nih.gov/), downloaded from the PanCancer Survival project (https://www.synapse.org/#!Synapse:syn1710282). We use gene expression profiles from renal clear cell carcinoma tumours to predict poor and good survival after Yuan et al. (2014) (Refs 1-2). The data consists of 150 tumours. Here we work only with the gene expression profiles generated.
+
+```{r,eval=TRUE}
+phenoFile <- sprintf("%s/extdata/KIRC_pheno.rda",path.package("netDx.examples"))
+lnames <- load(phenoFile)
+head(pheno)
+```
+
+```{r}
+pathFile <- sprintf("%s/extdata/Human_160124_AllPathways.gmt",
+ path.package("netDx.examples"))
+pathwayList <- readPathways(pathFile)
+```
+Filter for the genes measured in this dataset. For this example we have stored the names of genes measured in this dataset. In practice you would get this information from the corresponding input table.
+
+```{r}
+xpr_genes <- sprintf("%s/extdata/EMap_input/genenames.txt",
+ path.package("netDx.examples"))
+xpr_genes <- read.delim(xpr_genes,h=FALSE,as.is=TRUE)[,1]
+head(xpr_genes)
+```
+```{r}
+pathwayList <- lapply(pathwayList, function(x) x[which(x %in% xpr_genes)])
+```
+
+# Generate results
+
+```{r}
+inDir <- sprintf("%s/extdata/KIRC_output",
+ path.package("netDx.examples"))
+out <- plotAllResults(pheno, inDir,outDir=sprintf("%s/plots",getwd()),
+ fsCutoff=10,fsPctPass=0.7,pathwaySet=pathwayList)
+```
+
+The EnrichmentMaps and integrated PSN should now be present in Cytoscape, as well as PNG files in the output directory.
+
+
+# Examine results
+
+**Feature scores:**
+
+* `featScores/_featScores.txt`: Table of feature scores for all splits of nested CV.
+* `featScores/_FeatSel_cutoff_pct.png`: View of EnrichmentMap for each class
+* `EMap/_.gmt`: GMT file that serves as input for EnrichmentMap in Cytoscape
+* `EMap/_nodeAttrs_.txt`: Node attribute table with "maxScore" column that is mapped to node fill in Cytoscape.
+
+```{r}
+dir(sprintf("%s/plots",getwd()),recursive=TRUE)
+```
+
+# sessionInfo
+
+```{r}
+sessionInfo()
+```
+
+# References
+1. Yuan, Y. *et al.* (2014) Assessing the clinical utility of cancer genomic and proteomic data across tumor types. *Nat Biotechnol* **32**, 644-52.
+2. The Cancer Genome Atlas Research Network (2013). Comprehensive molecular characterization of clear cell renal cell carcinoma. *Nature* **499**, 43-9.
+3. Merico, D., Isserlin, R. & Bader, G.D. (2011). Visualizing gene-set enrichment results using the Cytoscape plug-in enrichment map. *Methods Mol Biol* **781**, 257-77.
+4. Kucera, M., Isserlin, R., Arkhangorodsky, A. & Bader, G.D. (2016). AutoAnnotate: A Cytoscape app for summarizing networks with semantic annotations. *F1000Res* **5**, 1717.
+
+
+
diff --git a/examples/plotPredictorResults_images/PSN.png b/examples/plotPredictorResults_images/PSN.png
new file mode 100644
index 00000000..a8ba509a
Binary files /dev/null and b/examples/plotPredictorResults_images/PSN.png differ
diff --git a/examples/plotPredictorResults_images/SURVIVENO.png b/examples/plotPredictorResults_images/SURVIVENO.png
new file mode 100644
index 00000000..8ddd4482
Binary files /dev/null and b/examples/plotPredictorResults_images/SURVIVENO.png differ
diff --git a/examples/plotPredictorResults_images/SURVIVEYES.png b/examples/plotPredictorResults_images/SURVIVEYES.png
new file mode 100644
index 00000000..ac4b057c
Binary files /dev/null and b/examples/plotPredictorResults_images/SURVIVEYES.png differ
diff --git a/misc/PanCancer/featSel_pathways/collectRes_batch.R b/misc/PanCancer/featSel_pathways/collectRes_batch.R
index 087a29f3..8aaa98d4 100644
--- a/misc/PanCancer/featSel_pathways/collectRes_batch.R
+++ b/misc/PanCancer/featSel_pathways/collectRes_batch.R
@@ -66,7 +66,7 @@ for (ctr in 1:nrow(setInfo)) {
outFile <- sprintf("%s/%s/KIRC_results_%s.Rdata",
outRoot,setInfo$outdir[ctr],cur)
- kset <- c(1:40,61:82)#c(1:40,41:53,61:73) #c(1:7,21:27,41:43,61:63,81:83)#c(1:49,51:68,76:92)#1:100#c(1:15,26:37,51:58,76:82) #1:maxk
+ kset <- 1:100 #c(1:40,61:82)#c(1:40,41:53,61:73) #c(1:7,21:27,41:43,61:63,81:83)#c(1:49,51:68,76:92)#1:100#c(1:15,26:37,51:58,76:82) #1:maxk
cat(sprintf("%s:Num runs=%i\n", cur,maxk))
val <