Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
143 lines (106 sloc) 4.38 KB
title author date output
Run WGCNA
Russ Poldrack
November 1, 2014
html_document

Gene network analysis for MyConnectome data

Code available at: https://github.com/poldrack/myconnectome/blob/master/myconnectome/rnaseq/Run_WGCNA.Rmd

This code loads the variance stabilized data (residualized against RIN and the first 3 PCs across subjects) and performs weighted gene coexpression network analysis (WGCNA). See the WGCNA page at UCLA for more info on this technique.

library(WGCNA)
allowWGCNAThreads(2)
options(stringsAsFactors = FALSE);
library(knitr)

basedir=Sys.getenv('MYCONNECTOME_DIR')
save_data=TRUE

outdir=sprintf("%s/rna-seq/WGCNA",basedir)

####Load the data
vsdata=read.table(sprintf('%s/rna-seq/varstab_data_prefiltered_rin_3PC_regressed.txt',basedir) )
datExpr=t(vsdata)
#### Call the network topology analysis function

powers = c(c(1:10), seq(from = 12, to=20, by=2))
if (file.exists(sprintf('%s/sft.Rdata',outdir))) {
  load(sprintf('%s/sft.Rdata',outdir))
} else {
  sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
  if (save_data) {
    save(sft,file=sprintf('%s/sft.Rdata',outdir))
  }
}

Plot of network topology results to examine the scale-free topology metric:

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],
labels=powers,cex=cex1,col="red");

abline(h=0.90,col="red")

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


Based on these results, we choose power=8 to obtain fit of r^2~0.9

##Run WGCNA estimation
if (!file.exists(sprintf('%s/net-thr8-prefilt-rinPCreg-48sess.Rdata',outdir))) {
  net = blockwiseModules(datExpr, power = 8,corType='bicor',
  TOMType = "unsigned", minModuleSize = 20,
  reassignThreshold = 1e-6, mergeCutHeight = 0.15,
  numericLabels = TRUE, pamRespectsDendro = FALSE,
  saveTOMs = FALSE,
  verbose = 3,maxBlockSize=4000)
  if (save_data) {
    save(net,file=sprintf('%s/net-thr8-prefilt-rinPCreg-48sess.Rdata',outdir))
    write.table(net$MEs, file=sprintf('%s/MEs-thr8-prefilt-rinPCreg-48sess.txt',outdir),row.names=FALSE,col.names=TRUE)
  }
} else {
  load(sprintf('%s/net-thr8-prefilt-rinPCreg-48sess.Rdata',outdir))
}


Plot WGCNA results


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

pdf(sprintf('%s/wgcna_tree.pdf',outdir))
mergedColors = labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
dev.off()

#### Save module assignments per gene
modules=as.data.frame(net$colors)
row.names(modules)=row.names(vsdata)
if (save_data) {
  write.table(modules,file=sprintf('%s/module_assignments_thr8_prefilt_rinPCreg.txt',outdir),col.names=FALSE,row.names=TRUE,quote=FALSE)
}

get hub genes using module membership (KME)

skme=signedKME(datExpr,net$MEs)
hubgenes=c()
nhubs=5
for (i in 1:(ncol(skme)-1)) {
	kme=sprintf('kME%d',i)
	menum=which(names(skme)==kme)
	idx=order(skme[,menum],decreasing=TRUE)
	hubs=row.names(skme)[idx[1:nhubs]]
	hubgenes=rbind(hubgenes,c(names(skme)[menum],hubs))
	}
kable(as.data.frame(hubgenes))
if (save_data) {
  write.table(hubgenes,file=sprintf('%s/hubgenes_thr8_prefilt_rinPCreg.txt',outdir),col.names=FALSE,row.names=TRUE,quote=FALSE)
}