...
Code Block | ||
---|---|---|
| ||
#build a adjacency "correlation" matrix enableWGCNAThreads() softPower = 18 adjacency = adjacency(datExpr, power = softPower, type = "signed") #specify network type head(adjacency) # Construct Networks- USE A SUPERCOMPUTER IRL ----------------------------- #translate the adjacency into topological overlap matrix and calculate the corresponding dissimilarity: TOM = TOMsimilarity(adjacency, TOMType="signed") # specify network type dissTOM = 1-TOM # Generate Modules -------------------------------------------------------- # Generate a clustered gene tree geneTree = flashClust(as.dist(dissTOM), method="average") plot(geneTree, xlab="", sub="", main= "Gene Clustering on TOM-based dissimilarity", labels= FALSE, hang=0.04) #This sets the minimum number of genes to cluster into a module minModuleSize = 30 dynamicMods = cutreeDynamic(dendro= geneTree, distM= dissTOM, deepSplit=2, pamRespectsDendro= FALSE, minClusterSize = minModuleSize) dynamicColors= labels2colors(dynamicMods) MEList= moduleEigengenes(datExpr, colors= dynamicColors,softPower = 18) MEs= MEList$eigengenes MEDiss= 1-cor(MEs) METree= flashClust(as.dist(MEDiss), method= "average") save(dynamicMods, MEList, MEs, MEDiss, METree, file= "Network_allSamples_signed_RLDfiltered.RData") #plots tree showing how the eigengenes cluster together png(file="clusterwithoutmodulecolors.png") plot(METree, main= "Clustering of module eigengenes", xlab= "", sub= "") #set a threhold for merging modules. In this example we are not merging so MEDissThres=0.0 MEDissThres = 0.0 merge = mergeCloseModules(datExpr, dynamicColors, cutHeight= MEDissThres, verbose =3) mergedColors = merge$colors mergedMEs = merge$newMEs dev.off() #plot dendrogram with module colors below it png(file="cluster.png") plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels= FALSE, hang=0.03, addGuide= TRUE, guideHang=0.05) moduleColors = mergedColors colorOrder = c("grey", standardColors(50)) moduleLabels = match(moduleColors, colorOrder)-1 MEs = mergedMEs dev.off() save(MEs, moduleLabels, moduleColors, geneTree, file= "Network_allSamples_signed_nomerge_RLDfiltered.RData") |
...
Code Block | ||
---|---|---|
| ||
# Correlate traits -------------------------------------------------------- #Define number of genes and samples nGenes = ncol(datExpr) nSamples = nrow(datExpr) #Recalculate MEs with color labels MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes MEs = orderMEs(MEs0) moduleTraitCor = cor(MEs, datTraits, use= "p") moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples) #Print correlation heatmap between modules and traits textMatrix= paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep= "") dim(textMatrix)= dim(moduleTraitCor) par(mar= c(6, 8.5, 3, 3)) #display the corelation values with a heatmap plot png(file="heatmap.png") labeledHeatmap(Matrix= moduleTraitCor, xLabels= names(datTraits), yLabels= names(MEs), ySymbols= names(MEs), colorLabels= FALSE, colors= blueWhiteRed(50), textMatrix= textMatrix, setStdMargins= FALSE, cex.text= 0.5, zlim= c(-1,1), main= paste("Module-trait relationships")) dev.off() |
Figure 4: Module-Trait relationships. Color scale (red-blue) represents the strength of the correlation between the module and the trait. For example, the turquoise module is highly significantly correlated with HPF, RedFluoro and GreenFluoro. Each box gives a correlation value (R^2) followed by p-value (in parenthesis).
...