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))
install.packages("BiocManager")
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("preprocessCore")
BiocManager::install("edgeR")
BiocManager::install("impute")
Load the required packages
library(WGCNA)
library(edgeR)
library(matrixStats)
library(impute)
library(preprocessCore)
Display the current working directory
getwd()
Set the working directory
workingDir = "/Users/ayanpatel/LumiereProject"
setwd(workingDir)
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
dim(Data)
names(Data)
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")