生信入门

您所在的位置：网站首页 › 显著性差异分析怎么算 › 生信入门

生信入门

2023-06-09 15:57| 来源: 网络整理| 查看: 265

0 分享至

用微信扫码二维码

分享至好友和朋友圈

加权基因共表达网络分析WGCNA，用于描述不同样品之间基因关联模式的系统生物学方法，鉴定高度协同变化的基因集，根据基因集的内连性和基因集与表型之间的关联鉴定出候补生物标记基因或治疗靶点。

相比只关注差异表达的基因，WGCNA利用几千甚至上万个变化最大的基因或全部基因的信息识别感兴趣的基因集，并与表型进行显著性关联分析。

WGCNA的基本分析流程就是这样的：

在开始做WGCNA分析之前，我们先理解以下这几个术语的含义。

共表达网络（加权基因网络）：点代表基因，边代表基因表达相关性。加权表示对相关性值进行幂次运算（幂次值=软阈值）。

Module（模块，高度内连的基因集）：在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。把基因聚类成模块后，可对模块进行三个层次的分析：第一层次，功能富集分析查看其功能特征是否与研究目的相符；第二层次，模块与性状进行关联分析，找出与关注性状相关度最高的模块；第三层次，模块与样本进行关联分析，找到样品特异高表达的模块。

connectivity（连接度）：类似于网络中“degree“的概念。每个基因的连接度是与其相连的基因的边属性之和。

module membership：给定基因表达谱与给定模型的eigengene的相关性。

module eigengene：给定模型的第一主成分，代表整个模型的基因表达谱。

hub gene：连接度最多或连接多个模块的基因。

adjacency matrix：邻接矩阵，基因与基因之间的加权相关性值构成的矩阵。

TOM：把邻接矩阵转换成拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可用于构建网络或绘制TOM图。

安装

source("http://bioconductor.org/biocLite.R"); biocLite("WGCNA")

这个R包是用于计算各种加权关联分析的功能集合，可用于网络构建、基因筛选、基因簇鉴定、拓扑特征计算、数据模拟和可视化等。

实战演示

# 导入数据 library(WGCNA) options(stringsAsFactors = FALSE) # 指允许R语言程序最大线程运行 allowWGCNAThreads() # 设置工作目录路径，R脚本语言和文件夹在同一个目录下 setwd("...") samples=read.csv('xxx.txt',sep = 't',row.names = 1) expro=read.csv('xxx.txt',sep = 't',row.names = 1) dim(expro) ##筛选方差前25%的基因## m.vars=apply(expro,1,var) expro.upper=expro[which(m.vars>quantile(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%的基因。

##样本聚类检查离群值## gsg = goodSamplesGenes(datExpr, verbose = 3); gsg$allOK sampleTree = hclust(dist(datExpr), method = "average") plot(sampleTree, main = "Sample clustering to detect outliers" , sub="", xlab="") save(datExpr, file = "FPKM-01-dataInput.RData")

从结果上来，我们分析的样本没啥离群值

举一个离群值的案例：

如果需要去除离群样本，则执行下列代码：

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(seq(1,10,by = 1), seq(12, 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 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")

软阈值是WGCNA的算法中非常重要的一个环节。

# 运行下列代码，让程序推荐你一个power, 数据质量太差啦，程序给了我"NA",自己设定power=14 > sft$powerEstimate [1] NA ##一步法网络构建：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) sft$powerEstimate

结果如下图：

结果是每个模块中包含的基因数量。一般来说，结果包含十几个到二十几个模块是比较正常的，此外一个模块中的基因数量不宜过多，像我们这个结果里模块1的基因数量达到了2307，这个就有点太多了，主要是因为我们powers=14，软阈值太低了导致的。所以说上述软阈值的筛选可以对我们的模块分析起到微调的作用。

##绘画结果展示### open a graphics window #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)

由于我们的软阈值比较低，所以这一结果中几乎没有grey模块，grey模块中的基因是共表达分析时没有被接受的基因，可以理解为一群散兵游勇。当然如果分析结果中grey模块中的基因数量比较多也是不太好的，表示样本中的基因共表达趋势不明显，不同特征的样本之间差异性不大，或者组内基因表达一致性比较差。

##结果保存### moduleLabels = net$colors moduleColors = 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$colors moduleColorsAutomatic = labels2colors(moduleLabelsAutomatic) moduleColorsWW = moduleColorsAutomatic MEs0 = moduleEigengenes(datExpr, moduleColorsWW)$eigengenes MEsWW = orderMEs(MEs0) modTraitCor = cor(MEsWW, samples, use = "p") colnames(MEsWW) modlues=MEsWW modTraitP = 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")) dev.off()

cex.lab可以更改X轴Y轴label字体的大小，cex.text可以更改热图中字体的大小，colors可以改变颜色。

样本特征和共表达模块的相关性热图中，grey模块中的相关性应该很小，如果你与样本特征相关性最显著的模块是grey模块，那肯定是有问题的。

###导出网络到Cytoscape#### Recalculate topological overlap if needed TOM = 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 = probes[inModule]; #modGenes = annot$gene_symbol[match(modProbes, annot$substanceBXH)]; # Select the corresponding Topological Overlap modTOM = TOM[inModule, inModule]; dimnames(modTOM) = list(modProbes, modProbes) # Export the network into edge and node list files Cytoscape can read cyt = 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 = moduleColors[inModule]);

这一步就是把选定的模块中的基因导出来，结果包含edges和nodes的信息。导出不同模块的基因只需要改变modules = c("模块颜色名")即可，输出多个模块的信息时，从该行代码运行即可，前面一行的运算量很大。

edges文件很大，试想一个模块中有500个基因，几乎两两基因之间都有关系，那就有上万条信息，构建出来的网络肯定密密麻麻的用不了。这里处理办法有两种：1、取Weight值前多少的作用关系；2、选定seed基因，比如某个lncRNA或者已知与表型具有密切关联的基因，构建与该基因有关的共表达网络

## 可视化基因网络## # 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 heatmap plotTOM = dissTOM^7; # Set diagonal to NA for a nicer plot diag(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 seed set.seed(10); select = sample(nGenes, size = nSelect); selectTOM = dissTOM[select, select]; # There's 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 = moduleColors[select]; # Open a graphical window#sizeGrWindow(9,9) # Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing# the color palette; setting the diagonal to NA also improves the clarity of the plot plotDiss = selectTOM^7; diag(plotDiss) = NA; TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

#此处根据基因间表达量进行聚类所得到的各模块间的相关性图 MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes MET = orderMEs(MEs) sizeGrWindow(7, 6) plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)

下面这个是分析共表达模块之间的相关性分析。

基因与性状和重要模块的关系：基因重要性和模块成员

通过将基因显著性GS定义为基因和性状之间的相关性的绝对值，我们可量化单个基因与我们感兴趣的性状（体重）的关联。对于每个模块，我们还定义了模块成员MM的定量度量，作为模块特征基因与基因表达谱的相关性，这样能够量化微阵列上所有基因与每个模块的相似性。

weight = as.data.frame(datTraits$weight_g) names(weight) = "weight" modNames = substring(names(MEs), 3) geneModuleMembership = as.data.frame(cor(datExpr, MEs, use= "p")) MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples)) names(geneModuleMembership) = paste("MM", modNames, sep="") names(MMPvalue) = paste("p.MM", modNames, sep="") geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = "p")) GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples)) names(geneTraitSignificance) = paste("GS.", names(weight), sep="") names(GSPvalue) = paste("p.GS.", names(weight", sep="")

模内分析：鉴定具有高GS和MM的基因

module = "brown" column = match(module, modNames) moduleGenes = moduleColors==module table(moduleGenes) sizeGrWindow(7, 7) par(mfrow = c(1,1)) verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]), abs(geneTraitSignificance[moduleGenes, 1]), xlab = paste("Module Membership in", module, "module"), ylab = "Gene significance for body weight", main = paste("Module membership vs. gene significance\n"), cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

小结

加权共表达网络主要涉及以下两个步骤，这样才能找出与目标性状显著相关的模块，对模块里面的基因进行后续分析。

共表达网络构建和模块划分

模块和性状关联分析

END

欢迎点赞、加关注

南博屹相伴，科研不孤单

往期推荐

生物信息学入门必知的知识

2021-10-15

生信入门必知︱富集分析里的统计学

2021-11-12

生物信息学分析必备︱不同物种的数据库

2021-10-29

论文投稿必备︱上传核酸序列到GenBank

2021-09-30

“上交原始数据”能解决科研诚信问题吗？

2021-09-24

OmicShare | 零代码在线生信分析工具

2021-09-17

点个在看你最好看

特别声明：以上内容(如有图片或视频亦包括在内)为自媒体平台“网易号”用户上传并发布，本平台仅提供信息存储服务。

Notice: The content above (including the pictures and videos if any) is uploaded and posted by a user of NetEase Hao, which is a social media platform and only provides information storage services.

/阅读下一篇/ 返回网易首页下载网易新闻客户端

【本文地址】

公司简介

联系我们