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

WGCNA新手入门笔记2（含代码和数据） 上次我们介绍了WGCNA的入门（WGCNA新手入门笔记（含代码和数据），大家在安装WGCNA包的时候，可能会遇到GO.db这个包安装不了的问题。主要问题应该是出在电脑的防火墙，安装时请关闭防火墙。如果还有问题，请先单独安装AnnotationDbi这个包，biocLite(AnnotationDbi)再安装GO.db，并尝试从本地文件安装该包。如果还有问题，请使用管理员身份运行R语言，尝试上述步骤。另外如果大家问题解决了请在留言处留个言，告知大家是在哪一步解决了问题，谢谢！因为本人没有进行单因素实验，不知道到底是哪个因素改变了实验结果。今天给大家过一遍代码。网盘中有代码和数据。链接：/s/1bpvu9Dt密码：w7g4#导入数据#library(WGCNA)options(stringsAsFactors = FALSE)enableWGCNAThreads() enableWGCNAThreads()是指允许R语言程序最大线程运行，像我这电脑是4核CPU的，那么就能用上3核：当然如果当前电脑没别的事，也可以满负荷运作samples=read.csv( Sam_info.txt,sep = t,s = 1)expro=read.csv( ExpData.txt,sep = t,s = 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(datExpr) 这一步是为了减少运算量，因为一个测序数据可能会有好几万个探针，而可能其中很多基因在各个样本中的表达情况并没有什么太大变化，为了减少运算量，这里我们筛选方差前25%的基因。当然了，以下几种情况，你可以忽略此步，转而执行下列代码datExpr= as.data.frame(t(expro);nGenes = ncol(datExpr)nSamples = nrow(datExpr) 情况1：所用的数据是比较老的芯片数据，探针数量较少；情况2：你的电脑足够强大，不必减少运算；情况3：你第一步导入的数据是差异基因的数据，已经作过初筛，比如这一文章的套路（PMID：）（个人不建议这么做）。#样本聚类检查离群值#gsg = goodSamplesGenes(datExpr, verbose = 3);gsg$allOKsampleTree = hclust(dist(datExpr), method = average)plot(sampleTree, main = Sample clustering to detect outliers, sub= , xlab= )save(datExpr, file = FPKM-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 = 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 R2,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) 软阈值是WGCNA的算法中非常重要的一个环节，简单的说硬阈值是一种一刀切的算法，比如高考分数500分能上一本，低于500就不行，软阈值的话切起来比较柔和一些，会考虑你学校怎么样，平时成绩怎么样之类。网盘中的数据跑起来其实是不太好的，没有合适的软阈值，这根线是要划到0.9的。或者右边这张图趋近于0像下面这个就是比较正常的。（/2535.html）这一步是为了算powers的值。一般来说，powers取红线（0.9）左右的数字都是可以的。如果你天秤座特征比较明显，你也可以运行下列代码，让程序推荐你一个值（本案例中返回值是NA，所以后面为了让程序能够进行下去，选了powers=14）。sft$powerEstimate 这个powers的值在后面的代码中会一直用到，所以你在跑别的数据的时候一定要更改powers的数值。#一步法网络构建：One-step network construction and module detection#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) 这一步就比较烧电脑了，关于网络构建的方法，分为一步法和多步法，一步法比较无脑，多步法能设置更多参数。想要了解多步法的请参考此文http:/tiramisutes.github.io/2016/09/14/WGCNA.html继续说上面的代码，结果跑出来如下图结果是每个模块中包含的基因数量。一般来说，结果包含十几个到二十几个模块是比较正常的，此外一个模块中的基因数量不宜过多，像我们这个结果里模块1的基因数量达到了2651，这个就有点太多了，主要是因为我们powers=14，软阈值太低了导致的。所以说上述软阈值的筛选可以对我们的模块分析起到微调的作用。#绘画结果展示# open a graphics window#sizeGrWindow(12, 9)# Convert labels to colors for plottingmergedColors = labels2colors(net$colors) # Plot the dendrogram and the module colors underneathplotDendroAndColors(net$dendrograms 1, mergedColorsnet$blockGenes 1, Module colors, 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-trait relationships) 这一步和网上的很多代码会有一些区别，有些代码会在这一环节构建样本特征的矩阵，而我们在最开始导入数据的时候，就是导入完整的样本特征矩阵，在这里直接读取就好了。（本人技术比较菜，所以在可视化的情况构建好矩阵，用代码凭空做一个矩阵，脑子不太够用）cex.lab可以更改X轴Y轴label字体的大小，cex.text可以更改热图中字体的大小，colors可以改变颜色。如果出现下述问题：请将绘图窗口最大化，然后最小化收起来：样本特征和共表达模块的相关性热图中，grey模块中的相关性应该很小，如果你与样本特征相关性最显著的模块是grey模块，那肯定是有问题的，毕竟grey模块中的基因是一群散兵游勇，它们的表达在各个样本中杂乱无章，根本说明不了问题。#导出网络到Cytoscape# Recalculate topological overlap if neededTOM = TOMsimilarityFromExpr(datExpr, power = 14); # Read in the annotation file# annot = read.csv(file = GeneAnnotation.csv);# Select modules需要修改，选择需要导出的模块颜色modules = c( lightgreen); # Select module probes选择模块探测probes = names(datExpr)inModule = is.finite(match(moduleColors, modules);modProbes = probesinModule; #modGenes = annot$gene_symbolmatch(modProbes, annot$substanceBXH);# Select the corresponding Topological OverlapmodTOM = TOMinModule, inModule;dimnames(modTOM) = list(modProbes, modProbes) # Export the network into edge and node list files Cytoscape can readcyt = exportNetworkToCytoscape(modTOM, edgeFile = paste( AS-green-FPKM-One-step-CytoscapeInput-edges-, paste(modules, collapse= -), .txt, sep= ), nodeFile = paste( AS-green-FPKM-One-step-CytoscapeInput-nodes-, paste(modules, collapse= -), .txt, sep= ), weighted = TRUE, threshold = 0.02, nodeNames = modProbes, #altNodeNames = modGenes,nodeAttr = moduleColorsinModule); 这一步就是把选定的模块中的基因导出来，结果包含edges和nodes的信息。导出不同模块的基因只需要改变modules = c(模块颜色名)即可，输出多个模块的信息时，从该行代码运行即可，前面一行的运算量很大。edges文件很大，试想一个模块中有500个基因，几乎两两基因之间都有关系，那就有上万条信息，构建出来的网络肯定密密麻麻的用不了。这里处理办法有两种：1、取Weight值前多少的作用关系；2、选定seed基因，比如某个lncRNA或者已知与表型具有密切关联的基因，构建与该基因有关的共表达网络(Time-Course Analysis of Brain Regional Expression Network Responses to Chronic Intermittent Ethanol and Withdrawal: Implications for Mechanisms Underlying Excessive Ethanol Consumption)# 可视化基因网络# # Calculate topological overlap anew: this could be done more efficiently by saving the TOM# calculated during module detection, but let us do it again here.dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 14); # Transform dissTOM with a power to make moderately strong connections more visible in the heatmapplotTOM = dissTOM 7; # Set diagonal to NA for a nicer plotdiag(plotTOM) = NA; # Call the plot function#sizeGrWindow(9,9)TOMplot(plotTOM, geneTree, moduleColors, main = Network heatmap plot, all genes) #随便选取1000个基因来可视化nSelect = 1000# For reproducibility, we set the random seedset.seed( 10);select = sample(nGenes, size = nSelect);selectTOM = dissTOMselect, select; # Theres no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.selectTree = hclust( as.dist(selectTOM), method = average)selectColors = moduleColorsselect; # Ope