Relating modules to external clinical traits and identifying important genes using the WGCNA package in R

Here is my code, can someone add the next step that is listed in the UCLA horvath tutorial which is: Relating modules to external clinical traits and identifying important genes. Also if there are any errors that you find could you fix them. My code is below:

Install required packages

install.packages(c('WGCNA', 'matrixStats'))
if (!require("BiocManager", quietly = TRUE))
if (!require("BiocManager", quietly = TRUE))
if (!require("BiocManager", quietly = TRUE))



Load the required packages


Display the current working directory


Set the working directory

workingDir = "/Users/ayanpatel/LumiereProject"

The following setting is important, do not omit.

options(stringsAsFactors = FALSE)

Read in the dataset

Data = read.csv("GSE194313_FPKM.matrix.csv", head = TRUE, sep = "")

Preprocessing steps

rownames(Data) = Data[,"GeneID"]
Data = Data[,-c(1:3)]
Data = Data[!rowSums(as.matrix(Data)) == 0,]

Take a quick look at the data set


Additional preprocessing steps

Data <- DGEList(Data)
Data <- calcNormFactors(Data, method = "TMM")
Data <- cpm(Data)

Continue with the rest of your code

powers <- c(c(1:10), seq(from = 12, to = 20, by = 2))
sft <- pickSoftThreshold(Data, powerVector = powers, verbose = 5)
sizeGrWindow(9, 5)
par(mfrow = c(1, 2))
cex1 <- 0.9

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],)

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")

Network Construction and Module Detection

net = blockwiseModules(Data, power = 6,
TOMType = "unsigned", minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs = TRUE,
saveTOMFileBase = "lunghemorrhageTOM",
verbose = 3)

Visualizing the network and module colors

sizeGrWindow(12, 9)

Convert labels to colors for plotting

mergedColors <- labels2colors(net$colors)

Plot the dendrogram and the module colors underneath

plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)

Saving Results

moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]
save(MEs, moduleLabels, moduleColors, geneTree, file = "Lunghemorrhage.RData")

This topic was automatically closed 21 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.