File:Wgcna null dend.png

Summary

<syntaxhighlight lang='r'> library(WGCNA)

1. Set seed for reproducibility

set.seed(123)

1. 1. Create a "Scaffold" for 3 modules (patterns of expression)

nSamples = 50

1. 2. Generate 1000 genes based on these patterns + some noise
2. Module 1: 200 genes correlated with pattern 1

noise = matrix(rnorm(nSamples*1000), nrow=nSamples)

1. 3. Combine into one expression matrix (Samples as Rows, Genes as Columns)

datExpr2 = as.data.frame(noise) colnames(datExpr2) = paste0("Gene", 1:1000) rownames(datExpr2) = paste0("Sample", 1:50)

1. Check dimensions (Should be 50 x 1000)

print(dim(datExpr2))

1. Choose Soft Threshold

powers = c(c(1:10), seq(from = 12, to=30, by=2)) sft = pickSoftThreshold(datExpr2, powerVector = powers, verbose = 3)

1. Plotting the results (Scale-Free Topology)
2. you should see the $R^2$ curve rise quickly.
3. The curve should hit a plateau (usually above 0.8 or 0.9) at a relatively low power (like 6 or 8).
4. See "WGCNA: an R package for weighted gene co-expression network analysis,"
5. Official WGCNA FAQ and Documentation https://edo98811.github.io/WGCNA_official_documentation/faq.html
6. says for more than 40 samples, power=6 for Unsigned and signed hybrid networks
7. power=12 for Signed networks
8. Choose signed network, if I want to see separate modules for upregulated genes and downregulated genes.
9. For signed network, modules are easier to interpret because all genes in a cluster move together.
10. Strongly recommended in later documentation for better biological accuracy.
11. Verdict:
12. If you are studying a disease vs. control and want to know which genes are
13. "activated" as a functional unit, choose Signed. If you just want to find
14. any genes that are mathematically related regardless of direction, stay with Unsigned.

1. Check the power suggested by the package

sft$powerEstimate # NA

1. If sft$powerEstimate returns NA, it means the data never hit the internal default threshold (usually 0.90)

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 = "Scale Independence")

text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],

    labels=powers, col="red")

1. Even at power 10, your model fit is only at 0.25.
2. In a successful WGCNA run, you want that value on the y-axis to be near 0.80.

1. the "Scale-Free Topology" fit is finally behaving. It reached the magic threshold of 0.80 at power 30.

1. Typically, for 50 samples, WGCNA authors suggest a power around 6 (for unsigned networks) or 12 (for signed networks).
2. The fact that you need 30 suggests that the "Noise" we added in the simulation is very strong relative to the "Signals."

pow <- 12 net = blockwiseModules(datExpr2, power = pow, TOMType = "signed",

                      minModuleSize = 30, reassignThreshold = 0, 
                      mergeCutHeight = 0.25, numericLabels = TRUE, 
                      pamRespectsDendro = FALSE, verbose = 3)

1. Calculating module eigengenes block-wise from all genes
2. Flagging genes and samples with too many missing values...
3. ..step 1
4. ..Working on block 1 .
5. TOM calculation: adjacency..
6. ..will not use multithreading.
7. Fraction of slow calculations: 0.000000
8. ..connectivity..
9. ..matrix multiplication (system BLAS)..
10. ..normalization..
11. ..done.
12. ....clustering..
13. ....detecting modules..
14. No modules detected in block 1
15. ..merging modules that are too close..
16. mergeCloseModules: Merging modules whose distance is less than 0.25
17. Warning message:
18. In blockwiseModules(datExpr2, power = pow, TOMType = "signed", minModuleSize = 30,  :
19. blockwiseModules: mergeCloseModules failed with the following error message:
20. Error in mergeCloseModules(datExpr, colors[gsg$goodGenes], cutHeight = mergeCutHeight,  :
21. Error in moduleEigengenes(expr = exprDataset$data, colors = setColors,  :
22. Color levels are empty. Possible reason: the only color is grey and grey module is excluded from the calculation.
25. --> returning unmerged colors.

1. Visualize the clusters

moduleColors = labels2colors(net$colors) plotDendroAndColors(net$dendrograms1, moduleColors[net$blockGenes1],

                   "Module colors", dendroLabels = FALSE, hang = 0.03,
                   addGuide = TRUE, guideHang = 0.05)

1. 0 is always the unassigned/grey module.

table(net$colors)

1. 0
2. 1000

<syntaxhighlight>