WGCNA新手入门笔记2(含代码和数据)

WGCNA新手入门笔记2（含代码和数据）上次我们介绍了WGCNA的入门（WGCNA新手入门笔记（含代码和数据）），大家在安装WGCNA包的时候，可能会遇到GO.db这个包安装不了的问题。主要问题应该是出在电脑的防火墙，安装时请关闭防火墙。如果还有问题，请先单独安装AnnotationDbi这个包，biocLite(AnnotationDbi)再安装GO.db，并尝试从本地文件安装该包。如果还有问题，请使用管理员身份运行R语言，尝试上述步骤。另外如果大家问题解决了请在留言处留个言，告知大家是在哪一步解决了问题，谢谢！因为本人没有进行单因素实验，不知道到底是哪个因素改变了实验结果。。。今天给大家过一遍代码。网盘中有代码和数据。链接：密码：w7g4##导入数据##library(WGCNA)options(stringsAsFactors=FALSE)enableWGCNAThreads()enableWGCNAThreads()是指允许R语言程序最大线程运行，像我这电脑是4核CPU的，那么就能用上3核：当然如果当前电脑没别的事，也可以满负荷运作samples=read.csv('Sam_info.txt',sep='t',row.names=1)expro=read.csv('ExpData.txt',sep='t',row.names=1)dim(expro)这部分代码是为了让R语言读取外部数据。当然了在读取数据之前首先改变一下工作目录，这一点在周二的文章中提过了。R语言读取外部数据的方式常用的有read.table和read.csv，这里用的是read.csv，想要查看某一函数的具体参数，可以用？函数名查看，比如：大家可以注意到read.table和read.csv中header参数的默认值是不同的，header=true表示第一行是标题，第二行才是数据，header=false则表示第一行就是数据，没有标题。##筛选方差前25%的基因##m.vars=apply(expro,1,var)expro.upper=expro[which(m.varsquantile(m.vars,probs=seq(0,1,0.25))[4]),]dim(expro.upper)datExpr=as.data.frame(t(expro.upper));nGenes=ncol(datExpr)nSamples=nrow(datExpr)这一步是为了减少运算量，因为一个测序数据可能会有好几万个探针，而可能其中很多基因在各个样本中的表达情况并没有什么太大变化，为了减少运算量，这里我们筛选方差前25%的基因。当然了，以下几种情况，你可以忽略此步，转而执行下列代码datExpr=as.data.frame(t(expro));nGenes=ncol(datExpr)nSamples=nrow(datExpr)情况1：所用的数据是比较老的芯片数据，探针数量较少；情况2：你的电脑足够强大，不必减少运算；情况3：你第一步导入的数据是差异基因的数据，已经作过初筛，比如这一文章的套路（PMID：27133569）（个人不建议这么做）。##样本聚类检查离群值##gsg=goodSamplesGenes(datExpr,verbose=3);gsg$allOKsampleTree=hclust(dist(datExpr),method=average)plot(sampleTree,main=Sampleclusteringtodetectoutliers,sub=,xlab=)save(datExpr,file=FPKM-01-dataInput.RData)执行这一段代码我们会得到下面这张图：从结果上来，我们分析的样本没啥离群值，所以代码里说不作处理。在网上的一个案例中，离群的样本就比较明显了。（）如果需要去除离群样本，则执行下列代码，其中cutHeight=多少就看你自己了。clust=cutreeStatic(sampleTree,cutHeight=20000,minSize=10)table(clust)keepSamples=(clust==1)datExpr=datExpr[keepSamples,]nGenes=ncol(datExpr)nSamples=nrow(datExpr)save(datExpr,file=FPKM-01-dataInput.RData)执行上述代码的话，就会去掉8个样本##软阈值筛选##powers=c(c(1:10),seq(from=12,to=20,by=2))sft=pickSoftThreshold(datExpr,powerVector=powers,verbose=5)par(mfrow=c(1,2));cex1=0.9;plot(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab=SoftThreshold(power),ylab=ScaleFreeTopologyModelFit,signedR^2,type=n,main=paste(Scaleindependence));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=SoftThreshold(power),ylab=MeanConnectivity,type=n,main=paste(Meanconnectivity))text(sft$fitIndices[,1],sft$fitIndices[,5],labels=powers,cex=cex1,col=red)软阈值是WGCNA的算法中非常重要的一个环节，简单的说硬阈值是一种一刀切的算法，比如高考分数500分能上一本，低于500就不行，软阈值的话切起来比较柔和一些，会考虑你学校怎么样，平时成绩怎么样之类。网盘中的数据跑起来其实是不太好的，没有合适的软阈值，这根线是要划到0.9的。或者右边这张图趋近于0像下面这个就是比较正常的。。。（）这一步是为了算powers的值。一般来说，powers取红线（0.9）左右的数字都是可以的。如果你天秤座特征比较明显，你也可以运行下列代码，让程序推荐你一个值（本案例中返回值是NA，所以后面为了让程序能够进行下去，选了powers=14）。sft$powerEstimate这个powers的值在后面的代码中会一直用到，所以你在跑别的数据的时候一定要更改powers的数值。##一步法网络构建：One-stepnetworkconstructionandmoduledetection##net=blockwiseModules(datExpr,power=14,maxBlockSize=6000,TOMType=unsigned,minModuleSize=30,reassignThreshold=0,mergeCutHeight=0.25,numericLabels=TRUE,pamRespectsDendro=FALSE,saveTOMs=TRUE,saveTOMFileBase=AS-green-FPKM-TOM,verbose=3)table(net$colors)这一步就比较烧电脑了，关于网络构建的方法，分为一步法和多步法，一步法比较无脑，多步法能设置更多参数。想要了解多步法的请参考此文继续说上面的代码，结果跑出来如下图结果是每个模块中包含的基因数量。一般来说，结果包含十几个到二十几个模块是比较正常的，此外一个模块中的基因数量不宜过多，像我们这个结果里模块1的基因数量达到了2651，这个就有点太多了，主要是因为我们powers=14，软阈值太低了导致的。所以说上述软阈值的筛选可以对我们的模块分析起到微调的作用。##绘画结果展示###openagraphicswindow#sizeGrWindow(12,9)#ConvertlabelstocolorsforplottingmergedColors=labels2colors(net$colors)#PlotthedendrogramandthemodulecolorsunderneathplotDendroAndColors(net$dendrograms[[1]],mergedColors[net$blockGenes[[1]]],Modulecolors,dendroLabels=FALSE,hang=0.03,addGuide=TRUE,guideHang=0.05)由于我们的软阈值比较低，所以这一结果中几乎没有grey模块，grey模块中的基因是共表达分析时没有被接受的基因，可以理解为一群散兵游勇。当然如果分析结果中grey模块中的基因数量比较多也是不太好的，表示样本中的基因共表达趋势不明显，不同特征的样本之间差异性不大，或者组内基因表达一致性比较差。##结果保存###moduleLabels=net$colorsmoduleColors=labels2colors(net$colors)table(moduleColors)MEs=net$MEs;geneTree=net$dendrograms[[1]];save(MEs,moduleLabels,moduleColors,geneTree,file=AS-green-FPKM-02-networkConstruction-auto.RData)这一步就是保存上面跑出来的结果了，同时哪个模块有多少基因一目了然。##表型与模块相关性##moduleLabelsAutomatic=net$colorsmoduleColorsAutomatic=labels2colors(moduleLabelsAutomatic)moduleColorsWW=moduleColorsAutomaticMEs0=moduleEigengenes(datExpr,moduleColorsWW)$eigengenesMEsWW=orderMEs(MEs0)modTraitCor=cor(MEsWW,samples,use=p)colnames(MEsWW)modlues=MEsWWmodTraitP=corPvalueStudent(modTraitCor,nSamples)textMatrix=paste(signif(modTraitCor,2),n(,signif(modTraitP,1),),sep=)dim(textMatrix)=dim(modTraitCor)labeledHeatmap(Matrix=modTraitCor,xLabels=colnames(samples),yLabels=names(MEsWW),cex.lab=0.9,yColorWidth=0.01,xColorWidth=0.03,ySymbols=colnames(modlues),colorLabels=FALSE,colors=blueWhiteRed(50),textMatrix=textMatrix,setStdMargins=FALSE,cex.text=0.5,zlim=c(-1,1),main=paste(Module-traitrelationships))这一