Hello,

My name is David Brohawn and I am a 4th year graduate student at Virginia Commonwealth University.

My advisor and I research ALS (Lou Gherigs disease). I have recently generated high density (50 Million 2X150 read pairs) RNA-Sequencing data for 15 human postmortem spinal cord tissue homogenates (7 ALS samples and 8 matched healthy controls). I have processed it with the Tuxedo Suite pipeline and have FPKM measurements for all known annotated genes in the hg19 build for all samples.

I find WGCNA highly intriguing, and would love to perform a consensus network analysis with the data we have prior to assessing differences between the two groups. The R code appears very manageable and does not seem unwieldy. 

However, I am unclear if we have enough samples to run this analysis.

I see the first question on the WGCNA FAQ page says:

We do not recommend attempting WGCNA on a data set consisting of fewer than 15 samples. In a typical high-throughput setting, correlations on fewer than 15 samples will simply be too noisy for the network to be biologically meaningful.

Does this mean 15 samples in each group (cases and controls, so 30 total) or 15 total?

Sifting through many publications and online forums (Seqanswers, Biostars) has not shed much light on this numbers question, and I find very few published reports using RNA-Sequencing data and WGCNA with smaller N's.

2) Presuming 15 total samples is the minimum recommended and we can run a consensus analysis, do you recommend use of the spearman correlation or the biweight midcorrelation as we are going to be susceptible to outliers? My understanding is Spearman is more robust to outliers, though it may wipe out true signals. 

Please advise, and thank you,

David Brohawn

Hi David,

I'm the WGCNA author and maintainer. 15 samples (cases plus controls combined) is indeed the minimum. 

Unless you have an independent data set, I would run a simple WGCNA on this data and see what comes out. You may very well find cell type-related modules which may or may not be interesting. Running a consensus WGCNA on the ALS and control samples viewed as two separate sets will not make much sense since you only have 7 and 8 samples, respectively.

As for Spearman vs. bicor, I suggest bicor with additional argument maxPOutliers = 0.1 or so. I personally don't have much experience with Spearman correlation but early attempts in our group to use it were not successful.

In terms of pre-processing, please read through the WGCNA FAQ at http://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/faq.html, particularly item 3.

Hope this helps.

Peter

Yes, use a binary variable coded 0 for controls and 1 for ALS as the trait. One major change against the tutorials that we recommend is to use argument "type" and/or "networkType" with value "signed hybrid". This applies to functions like pickSoftThreshold, adjacency and blockwiseModules. We have found that signed networks are usually cleaner and more easy to interpret than the default unsigned networks.

HTH,

Peter

Hello Peter,

Thanks for your quick reply. I already filtered and log transformed my FPKM values, so they are ready to go.

If we run a standard WGCNA analysis with these 15 samples' pre-processed data using the bicor and maxPOutliers flags you recommended:

Could I then set disease/control status as a binary external trait variable and relate identified modules for all 15 samples to that status (in a manner similar to what they did with weight in Tutorial 1 -http://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-03-relateModsToExt.pdf)

I understand weight is a continuous variable whereas disease status is not. However, I remain curious if our study parameters are amenable to using your software to look for disease-specific changes.

Thanks so much for your consideration,

David Brohawn

Many thanks sir! Very excited to see what comes out.

David Brohawn

Hello Peter, 

I have been practicing with the online WGCNA tutorials in preparation of running our analysis, and have a followup question regarding proper usage of bicor and maxPOutliers:

I successfully created the datExpr and datTraits variables as instructed in the data input and cleaning section (Tutorial part 1)

I opted for the step-by-step network construction choice (Tutorial part 2b)(http://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-02-networkConstr-man.pdf)

I modified the picksoftThreshold and adjacency commands (as well as the TOMsimilarity command) from what is listed in the tutorial to incorporate the bicor and signed hybrid network flags per your suggestion: 

sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5, networkType ="signed hybrid", corFnc= "bicor")

adjacencymatrix = adjacency(datExpr, power=softpower, corFnc="bicor", type="signed hybrid")

TOM = TOMsimilarity(adjacency, TOMType="signed")

However, these changes did not incorporate the maxPOutliers = 0.1 flag you recommended. I see in the WGCNA manual it is not an option for the picksoftThreshold, adjacency, or the TOM commands. 

I was able to carry out the rest of Tutorial part 2b and 3 without any issues. I see the maxPOutliers flag does appear in the blockwiseModules command for the automated network construction directions in Tutorial part 2c.

Is it appropriate to leave out the bicor command with maxPOutliers=0.1, or is there a way/need to incorporate it if I use the directives in Tutorial part2b (step-by-step network construction)?

Thank you for your help in advance!

David Brohawn

Hello, 

I should add that I was able to apply the bicor function to my datExpr variable as below:

Test <- bicor(datExpr, use="pairwise.complete.obs", maxPOutliers=0.1)

prior to using the pickSoftThreshold command, though I am not sure if that was appropriate.

Thanks again,
 

David Brohawn

Use the argument corOptions = list(maxPOutliers =0.1) to pickSoftThreshold() and to adjacency().

Peter

Hi Peter,

Thank you again for all of your help.

I was able to add corOptions=list(maxPOutliers=0.1) to the pickSoftThreshold command:

sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5, networkType ="signed hybrid", corFnc= "bicor", corOptions=list(maxPOutliers=0.1)) 

and retrieved a table of soft threshold powers with their corresponding R.sq values.

Power SFT.R.sq  slope truncated.R.sq mean.k. median.k. max.k.
1      1   0.1500  1.260          0.947 2010.00   1980.00 3220.0
2      2   0.0137 -0.235          0.943  850.00    821.00 1730.0
3      3   0.2760 -1.010          0.946  431.00    404.00 1090.0
4      4   0.4790 -1.370          0.960  244.00    221.00  740.0
5      5   0.6090 -1.650          0.969  149.00    129.00  533.0
6      6   0.6820 -1.850          0.977   96.80     80.20  403.0
7      7   0.7260 -1.990          0.981   65.50     51.70  313.0
8      8   0.7490 -2.120          0.981   45.90     34.50  247.0
9      9   0.7760 -2.160          0.987   33.10     23.70  199.0
10    10   0.7940 -2.190          0.989   24.50     16.70  162.0
11    12   0.8180 -2.220          0.991   14.20      8.82  111.0
12    14   0.8300 -2.210          0.991    8.76      4.97   78.9
13    16   0.8430 -2.170          0.992    5.68      2.94   57.5
14    18   0.8550 -2.130          0.990    3.83      1.82   42.9
15    20   0.8650 -2.060          0.984    2.67      1.16   32.6

However, when I tried to add corOptions=list(maxPOutliers=0.1) to my adjacency command:

adjacency = adjacency(datExpr, power=softpower, corFnc="bicor", type="signed hybrid", corOptions=list(maxPOutliers=0.1))

I recieved the following error message: Error in bicor(datExpr, 0.1) : 
  'x' and 'y' have incompatible dimensions (unequal numbers of rows).

If I remove the corOptions flag, adjacency will run, albeit with the following warning:

Warning message:
In bicor(datExpr, use = "p") :
  bicor: zero MAD in variable 'x'. Pearson correlation was used for individual columns with MAD=NA.

Can you advise?

Thanks!

David Brohawn

Hi Dr. Langfelder,

Have you had a chance to look over the error message I received when trying to use the maxPOutliers flag in the adjacency command?

I'm guessing my previous question got buried in the shuffle.

Thanks!

Dave Brohawn

Dear David,

I have the same problem like you, could you pass the error?

Regards

Sorry to revive a long-dead thread, I missed this message. The function adjacency takes the corOptions argument in a different format, you need to specify corOptions = "use = 'p', maxPOutliers = 0.1" (i.e. a literal string, not a list).

Hi Peter,

I had a problem similar to this discussion while trying to use the bicor function. I followed your instruction to use the corOptions as a literal string, but then the bicor was not implemented. Please see various combinations I tried. Could you please suggest the right usage? 

> adjacency = adjacency(datExpr, power = softPower, type = "signed hybrid",corFnc="bicor");
Warning message:
In bicor(datExpr, use = "p") :
  bicor: zero MAD in variable 'x'. Pearson correlation was used for individual columns with MAD=NA.

> adjacency = adjacency(datExpr, power = softPower, type = "signed hybrid",corFnc="bicor", corOptions = "use = 'p', maxPOutliers = 0.1");
Warning message:
In bicor(datExpr, use = "p", maxPOutliers = 0.1) :
  bicor: zero MAD in variable 'x'. Pearson correlation was used for individual columns with MAD=NA.

> adjacency = adjacency(datExpr, power = softPower, type = "signed",corFnc="bicor", corOptions = "use = 'p'");
Warning message:
In bicor(datExpr, use = "p") :
  bicor: zero MAD in variable 'x'. Pearson correlation was used for individual columns with MAD=NA.

> adjacency = adjacency(datExpr, power = softPower, type = "signed",corFnc="bicor", corOptions = "use = 'p', maxPOutliers = 0.1");
Warning message:
In bicor(datExpr, use = "p", maxPOutliers = 0.1) :
  bicor: zero MAD in variable 'x'. Pearson correlation was used for individual columns with MAD=NA.

> adjacency = adjacency(datExpr, power = softPower, type = "unsigned",corFnc="bicor", corOptions = "use = 'p', maxPOutliers = 0.1");
Warning message:
In bicor(datExpr, use = "p", maxPOutliers = 0.1) :
  bicor: zero MAD in variable 'x'. Pearson correlation was used for individual columns with MAD=NA.

Thanks,

Manoj.