Goals

1) Make the Sage coexpression software runnable by any data analyst in R.

Document the analysis steps, package functions, and adjustable parameters.  Write start-to-finish vignette(s) having sufficient detail to guide a new analyst through the software, including making parameter settings and other intermediate decisions.  Refactoring software to enhance readability and/or to make decision points and adjustable parameters explicit, is within scope for this goal.

2) Clearly explain the methodology underlying the coexpression algorithms.

This includes contrasting the Sage algorithms with the separately published “WGCNA” algorithms and explaining the rationale for the differences and the decision process for choosing which to use.

3) Make the Sage coexpression software publicly available.

The code must be available through a standard R distribution channel (CRAN or Bioconductor).  An option is to merge Sage-specific algorithms or steps into WGCNA (if the two are sufficiently overlapping).

4) Make the Sage coexpression software perform well, on commonly available hardware.

Currently, special high capacity hardware is required to run the Sage coexpression software on commonly encountered data set sizes.  The goal is to optimize the code to run on commonly available hardware, at a minimum the code should be runnable on the 68GB, quad core servers available through Amazon Web Services. Ideally, the code would be runnable on a “heavy” Sage laptop (8 GB total RAM).

Strategy

The steps for Sage Coexpression are: 

1. Compute correlation coefficient matrix.
2. Determine optimal value for the scale free exponent, beta, and collect regression statistics.
3. Compute the topological overlap matrix (TOM).
4. Perform hierarchical clustering of genes, based on TOM.
5. Detect and label modules in TOM, using "Dynamic Tree Cutting".
6. Merge modules based on hierarchical clustering of representative genes.
7. Cluster samples hierarchically.
8. Compute intra/inter-module network statistics, per gene.
9. Produce diagnostic plots (dendrograms, heat maps, statistical scatter plots).
10. Produce tabular output of module membership, network statistics, and scale-free regression statistics.

Elsewhere http://sagebionetworks.jira.com/wiki/display/SCICOMP/Coexpression+Evaluation we have shown that for realistic dataset sizes, the vast majority of time is spent in steps 1 and 3, and that steps 1, 3, and 4 are identical in the Sage and UCLA-WGCNA code bases, while the UCLA-WGCNA code uses compiled/optimized software for these three steps.  Further, we have seen that steps 2 and 5 (with the right parameter choices in the UCLA package) produce extremely similar results.  (Note, the observed similarity is no surprise, since the two code bases represent forks from an original set of algorithms, which have evolved separately for appx. 6 years.)

Our strategy, therefore, is:

• leverage the UCLA-WGCNA package for the "common" steps, 1->5, gaining significant performance
• provide the user a parameter choice at step 5, to do "tree cutting" in the manner of the Sage algorithm, or in that of the UCLA-WGCNA algorithm
• provide two algorithms for step 6 (module merging), allowing a user to choose the Sage or UCLA-WGCNA algorithm
• leverage the UCLA-WGCNA dendrogram/module plotting algorithm in step 9
• maintain the Sage algorithms for the Sage-specific post-processing, i.e. statistics in step 8 and the heat maps in step 9.  Accelerate step 8 using compiled code.

External dependencies

WGCNA::cor  -- the compiled/accelerated Pearson correlation computation

WGCNA::pickSoftThreshold -- optimal choice of scale free exponent

WGCNA::Tomdist -- the compiled/accelerated TOM computation

flashClust::flashClust -- the compiled/accelerated hierarchical clustering computation

WGCNA::scaleFreePlot -- scatter plot of network connectivity, with regression line

dynamicTreeCut::cutreeDynamic -- tree cutting / module determination

WGCNA::plotDendroAndColors -- graphic function to plot dendrogram with colored modules aligned underneath

Sage software dependencies

module merging, by analyzing most-highly-connected genes in each module

computation of within- and between- module per-gene connectivity statistics (but using additional compiled code to dramatically speed up analysis of large modules)

heatmap generation for correlation and TOM matrices

Reduction of scale-free regression threshold (R^2) when no solution is found.

The Package

The source code for the created package is available on our Atlassian-hosted SVN repository, under the SCICOMP project:

http://sagebionetworks.jira.com/source/browse/SCICOMP/trunk/Coexpression/SageBionetworksCoex

Package functions

The main 'points of entry' into the package are:

performCoexFromFiles - run the entire package using a file of gene expression data as input

performCoexpressionAnalysis - run the analysis portion of Coexpression, taking a data frame as input

clusterGenes - run the rote, time consuming portion of Coexpression, taking a data frame as input

modulesFromGeneTree - choose modules from a gene dendrogram, taking the output of 'clusterGenes' as input

clusterSamples, intraModularStatistics - analysis steps auxiliary to module determination

createDiagnosticPlots - create dendrogram, heatmap plots, etc., from the results of the coexpression analysis

Package features

Comparison to Sage code base

Correlation computation is faster, with identical results.
TOM computation is faster, with identical results.
Hierarchical clustering is faster, with identical results.
Scale-free exponent (beta) determination is similar, with very similar results and regression statistics.
"Dynamic tree cutting" algorithm is the same, with very similar results.

Intramodular statistics are faster, with identical results.
Diagnostic plot set is reduced from 12 to 6, omitting redundant plots.

Sample clustering was omitted, having been deemed unnecessary.

Additional Features

Option to do tree cutting and/or subsqeuent merging by UCLA-WGCNA algorithm or by 'Sage classic' algorithm.
Separation of rote, time consuming steps (correlation, TOM computation) from tree cutting.
Separation of analysis from plotting.
Separation of analysis from file system, to facilitate Synpase integration.

Performance and Limits

Coexpression has O(n^3) ("cubic") time complexity and O(n^2) ("quadratic") space complexity, where n is the number of probes in the dataset.  The time complexity is primarily due to the TOM computation.  The space complexity is due to the need to hold the nXn correlation and TOM matrices in memory.  The R language has an inherent limit on the size of a vector or matrix of about 2 billion elements (http://stat.ethz.ch/R-manual/R-devel/library/base/html/Memory-limits.html).  This means the maximum size of a square matrix is 46340 x 46340.  A square matrix in R requires 8n^2 bytes.  (Values are double precision and R has no provision for single precision representation of floating point values.)  Such a matrix requires 16.78 GB of contiguous memory.   The coexpression package uses two nXn matrices (as mentioned above) and at times creates temporary variables at sizes equal to or a large fraction of the nXn matrix size.  Therefore, to support maximum dataset sizes, machines used to run coexpression should have dozens of GB of memory.  In the empirical evaluations summarized below, we use Amazon Web Services' "High-Memory Quadruple Extra Large"  machines having 68 GB RAM, each of which costs (at the time this is written) $2.88/hour (http://aws.amazon.com/ec2/instance-types/). 

Evaluation

Dataset were run through the original Sage coexpression code and the newly created package to compare performance and similiarity of results. 

Similiarity questions:  Is the scale-free exponent 'beta' the same?  Are the dendrograms the same?  How similar are the resultant modules? 

Performance questions:  For datasets having >18,000 probes, how much time and space does each algorithm use?

 * These were run on an Amazon Elastic Compute Cloud (EC2) "High-Memory Quadruple Extra Large" unix server, having 68GB of RAM.

** Note:  UCLA-WGCNA package also runs out of memory.  (An alternative is to use the WGCNA preprocessing step of K-means decomposition, which has been shown to work with >50K genes.)

*** http://florence.acadiau.ca/collab/hugh_public/index.php?title=R:compare_partitions

Note: We skip the performance evaluation for the small data sets.

The details of the differences summarized in the table can be found here: http://sagebionetworks.jira.com/wiki/display/SCICOMP/Package+Comparison+Details

Conclusions

The new package has considerably better time performance than does the original code.  Though the two algorithms have the same approach for computing 'beta', the results can vary greatly.  When beta is the same, the dendrograms and modules are similar or identical.  However, module determination is very sensitive to beta, which can vary greatly with small changes in regression statistics, as can be seen on the 'details' page. 

Goals, Revisited

Choices for package 'publication' include:

• Making available on Sage website
• Making available in Synapse
• Publish to BioConductor
• Publish to CRAN
• Request that UCLA include as part of WGCNA package