Entering edit mode
9.6 years ago
madkitty
▴
690
I'm trying to compare 3 "control" datasets to 2 different samples in WGCNA. I have no clue about standards to run WGCNA nor do I understand why people use that stuff.. anyhow my boss wants to me to use it. When following the steps posted here: http://pklab.med.harvard.edu/scw2014/WGCNA.html
Apparently, we need to define a soft threshold, and I don't understand how can I determine necessary parameters. Is there any rule to help me pick up the parameters cex1, abline h and the soft power threshold?
powers = c(c(1:10), seq(from = 12, to=20, by=2)); sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,powerVector = powers,corFnc = cor,corOptions = list(use = 'p'),networkType = "unsigned")
## Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
## 1 1 0.719 -1.700 0.874 40.60 37.900 102.0
## 2 2 0.521 -0.707 0.679 15.00 12.500 40.0
## 3 3 0.814 -0.476 0.764 9.92 7.970 25.9
## 4 4 0.900 -0.671 0.887 7.73 5.790 23.6
## 5 5 0.803 -0.809 0.748 6.47 4.140 22.2
## 6 6 0.749 -0.901 0.677 5.63 3.610 21.1
## 7 7 0.880 -0.913 0.846 5.02 3.100 20.2
## 8 8 0.885 -0.928 0.853 4.55 2.610 19.4
## 9 9 0.849 -0.940 0.806 4.17 2.170 18.6
## 10 10 0.918 -0.940 0.896 3.86 1.890 18.0
## 11 12 0.897 -0.942 0.867 3.36 1.610 16.8
## 12 14 0.882 -0.949 0.849 2.99 1.360 15.8
## 13 16 0.819 -0.996 0.770 2.69 1.180 15.0
## 14 18 0.890 -0.974 0.858 2.44 0.997 14.3
## 15 20 0.885 -0.987 0.852 2.24 0.802 13.7
# Plot the results
sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit, signed R^2",type="n", main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=cex1,col="red");
# Red line corresponds to using an R^2 cut-off
abline(h=0.80,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")