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

1、WGCNA新手入门笔记2(含代码和数据)上次我们介绍了WGCNA的入门(WGCNA新手入门笔记(含代码和数据)，大家在安装WGCNA包的时候，可能会遇到GO.db这个包安装不了的问题。主要问题应该是出在电脑的防火墙，安装时请关闭防火墙。如果还有问题，请先单独安装AnnotationDbi这个包，biocLite(AnnotationDbi)再安装GO.db，并尝试从本地文件安装该包。如果还有问题，请使用管理员身份运行R语言，尝试上述步骤。另外如果大家问题解决了请在留言处留个言，告知大家是在哪一步解决了问题，谢谢！因为本人没有进行单因素实验，不知道到底是哪个因素改变了实验结果。今天给大家过一遍代码。网盘中有代码和数据。链接：http:/ #导入数据#library(WGCNA)options(stringsAsFactors=FALSE)enableWGCNAThreads()enableWGCNAThreads()是指允许R语言程序最大线程运行，像我这电脑是4核CPU的，那么就能用上3核：当然如果当前电脑没别的事，也可以满负荷运作samples=read.csv(Sam_info.tx

2、t,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=exprowhich(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(

3、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-

4、01-dataInput.RData)执行这一段代码我们会得到下面这张图：从结果上来，我们分析的样本没啥离群值，所以代码里说不作处理。在网上的一个案例中，离群的样本就比较明显了。(https:/www.shengxin.ren/article/88)如果需要去除离群样本，则执行下列代码，其中cutHeight=多少就看你自己了。clust=cutreeStatic(sampleTree,cutHeight=20000,minSize=10)table(clust)keepSamples=(clust=1)datExpr=datExprkeepSamples,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

5、=c(1,2);cex1=0.9;plot(sft$fitIndices,1,-sign(sft$fitIndices,3)*sft$fitIndices,2,xlab=SoftThreshold(power),ylab=ScaleFreeTopologyModelFit,signedRA2,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)软阈值是WGC

6、NA的算法中非常重要的一个环节，简单的说硬阈值是一种一刀切的算法，比如高考分数500分能上一本，低于500就不行，软阈值的话切起来比较柔和一些，会考虑你学校怎么样，平时成绩怎么样之类。网盘中的数据跑起来其实是不太好的，没有合适的软阈值，这根线是要划到0.9的。或者右边这张图趋近于0像下面这个就是比较正常的。(http:/www.bio-info- #一步法网络构建：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)这一步就比较烧电脑了，关于网络构建的方法，分

7、为一步法和多步法，一步法比较无脑，多步法能设置更多参数。想要了解多步法的请参考此文http:/tiramisutes.github.io/2016/09/14/WGCNA.html继续说上面的代码，结果跑出来如下图结果是每个模块中包含的基因数量。一般来说，结果包含十几个到二十几个模块是比较正常的，此外一个模块中的基因数量不宜过多，像我们这个结果里模块1的基因数量达到了2651，这个就有点太多了，主要是因为我们powers=14，软阈值太低了导致的。所以说上述软阈值的筛选可以对我们的模块分析起到微调的作用。# #绘画结果展示#openagraphicswindow#sizeGrWindow(12,9)# ConvertlabelstocolorsforplottingmergedColors=labels2colors(net$colors)# PlotthedendrogramandthemodulecolorsunderneathplotDendroAndColors(net$dendrograms1,mergedColorsnet$blockGenes1,Modulecolors,

8、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$dendrograms1;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=coln

《WGCNA新手入门笔记2(含代码和数据)》由会员s9****2分享，可在线阅读，更多相关《WGCNA新手入门笔记2(含代码和数据)》请在金锄头文库上搜索。