用微信扫码二维码

分享至好友和朋友圈

在上一期的生信实用小技巧︱TCGA数据做表达差异分析及数据可视化，小编已介绍了如何利用TCGA数据进行表达差异分析。今天小编继续介绍用TCGA数据做WGCNA分析。

WGCNA（Weighted Gene Coexpression Network Analysis）是一个基于基因表达数据构建基因共表达网络的方法。WGCNA和差异基因分析的区别在于差异基因分析主要针对样本之间的差异，WGCNA主要针对基因之间的关系。WGCNA通过分析基因之间的关联性，将基因区分为多个模块。最后通过这些模块和样本表型之间的关联性分析，寻找特定表型的分子特征。

分析流程

简单说说分析流程：

1、从TCGA数据库下载要研究的癌种的基因、miRNA表达数据。

2、分别针对基因、miRNA进行表达差异分析。

3、分别针对基因和miRNA进行WGCNA分析，获得聚类模块。并分析与病理学分期相关联的基因模块。

4、验证基因模块与临床特征相关的基因。

5、选择模块中的hub gene和hub miRNA进行生存分析。

6、从病理分期相关联的模块中验证的基因和miRNA做代谢通路分析，构建miRNA-基因的调控网络。

下面以TCGA乳腺癌基因数据为例，运用WGCNA对乳腺癌各个亚型做共表达模块。我们首先从TCGA数据库里下载乳腺癌的RNA-seq数据和临床病理资料，有关数据挖掘的操作这里不说了，大家可浏览生信实用小技巧︱从TCGA数据库中挖掘数据。

接下来用R代码实现WGCNA:

setwd('输入TCGA data电脑路径')samples=read.csv('xxx.txt',sep = '\t',row.names = 1)dim(samples)expro=read.csv('xxx.txt.cv.txt',sep = '\t',row.names = 1)dim(expro)

数据读取完成，从结果中可知道样本数量和基因数量。然后我们选择一些具有代表性的基因做WGCNA分析，这里我们选择在100个样本中方差较大的那些基因。

m.vars=apply(expro,1,var)expro.upper=expro[which(m,vars>quantile(m,vars, probs = seq(0, 1, 0.25))[4]),]

通过上述步骤拿到的基因表达谱作为WGCNA的输入数据集，进一步看看样本之间的差异情况。

datExpr=as.data.frame(t(expro.upper));gsg = goodSamplesGenes(datExpr, verbose = 3);gsg$allOKsampleTree = hclust(dist(datExpr), method = "average")plot(sampleTree, main = "Sample clustering to detect outliers" , sub ="", xlab="")

从图中看出大部分样本表现比较相近，但是有2个离群样本，对后续分析可能造成影响，我们要把它去掉。

clust = cutreeStatic(sampleTree, cutHeight = 80000, minSize = 10)table(clust)keepSamples = (clust==1)datExpr = datExpr[keepSamples, ]nGenes = ncol(datExpr)nSamples = nrow(datExpr)save(datExpr, file = "FPKM-01-dataInput.RData")

得到最终的数据矩阵后，我们要确定软阈值

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

从图中看出软阈值7比较合适，所以选择软阈值7进行共表达模块挖掘。

pow=7net = blockwiseModules(datExpr, power = pow, maxBlockSize = 7000,TOMType = "unsigned", minModuleSize = 30,reassignThreshold = 0, mergeCutHeight = 0.25,numericLabels = TRUE, pamRespectsDendro = FALSE,saveTOMs = TRUE,saveTOMFileBase = "FPKM-TOM",verbose = 3)table(net$colors)sizeGrWindow(12, 9)mergedColors = labels2colors(net$colors)plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],groupLabels = c("Module colors","GS.weight"),dendroLabels = FALSE, hang = 0.03,addGuide = TRUE, guideHang = 0.05)

从上图可看出大部分基因在灰色区域，灰色部分一般认为没有模块接受的。

mergedColors = labels2colors(net$colors)table(mergedColors)

当我们拿到了相关性较高的模块（如Yellow模块）时，那怎样获取重要的基因呢。我们需要先计算每个基因与模块的相关关系

modTraitCor = cor(MEsFemale, datExpr, use = "p")modTraitP = corPvalueStudent(modTraitCor, nSamples)corYellow=modTraitCor[which(row.names(modTraitCor)=='MEyellow'),]head(corYellow[order(-corYellow)])

进一步分析Yellow模块的基因共表达关系，我们用cytoscape进行数据可视化。

TOM = TOMsimilarityFromExpr(datExpr, power = pow);probes = names(datExpr)mc='yellow'mcInds=which(match(moduleColorsAutomatic, gsub('^ME','',mc))==1modProbes=probes[mcInds]modTOM = TOM[mcInds, mcInds];dimnames(modTOM) = list(modProbes, modProbes)cyt = exportNetworkToCytoscape(modTOM,edgeFile = paste("edges-", mc, ".txt", sep=""),nodeFile = paste("nodes-", mc, ".txt", sep=""),weighted = TRUE,threshold = median(modTOM),nodeNames = modProbes,#altNodeNames = modGenes,nodeAttr = moduleColorsAutomatic[mcInds]);

我们来看看这篇今年3月发表在《Journal of Cancer》的文章Identification of tumor mutation burden-related hub genes and the underlying mechanism in melanoma。

这篇文章的作者从TCGA数据库中获取了472个黑素瘤样本的基因、miRNA表达数据以及临床特征数据，然后分别对基因和miRNA进行差异表达基因分析、WGCNA分析、构建蛋白互作网络和生存分析找出与肿瘤突变负荷TMB相关联的hub gene，接着作者通过构建竞争性内源RNA（ceRNA）网络去探讨hub gene功能潜在的分子机制。总的来说，这篇文章的作者通过上述方法找出三个TMB相关联的关键基因，从而发现TMB对黑素瘤具有很大的影响。TMB、与TMB有关的hub基因可被视为候选的预测性biomarkers。

分析结果

1、数据下载和筛选

从TCGA数据库中下载和筛选472例黑素瘤样本的基因表达数据和miRNA表达数据。

2、差异分析

用“edgeR”语言包分别对黑素瘤组织的低TMB和高TMB组基因和miRNA进行表达差异分析，在高TMB组中共有443个DEGs发生了显著差异，其中370个基因上调和73个基因下调。排名前10的上调基因包括FGFR3、TCHH、CNFN、TGM1、SULT2B1、ENTPD3、SLC22A1、TBX1、VIPR1、TREX2；排名前10的下调基因包括ZG16B、TG、ADAMTS8、PIGR、KLHL41、DES、CA6、MRGPRX4、RRAD、SCX。

3、对基因表达数据和TBM数据进行WGCNA分析，找出与TBM相关的hub gene，作者从56499个基因里面筛选出35312个表达极低的基因，于是获得了21183个基因，然后形成了基因模块，从而构建了WGCNA调控网络，继而找出与TMB相关度最高的模块。所以作者在这23个模块中，发现棕褐色模块与TMB相关度最高。接着，作者针对棕褐色模块中的59个基因通过PPI网络和MCODE进一步分析。最终作者筛选出与TMB相关联的hub基因，就是FLNC、NEXN和TNNT3.

好了，本期就讲到这里。如果你想要阅读上面这篇文献的全文，记得点击左下角的“阅读原文”。

欢迎分享到朋友圈，让更多人跟我们一起学习吧。

南博屹相伴，科研不孤单

生信实用小技巧︱TCGA数据做表达差异分析及数据可视化

生信实用小技巧︱从TCGA数据库中挖掘数据

单细胞转录组测序技术介绍及应用

单细胞转录组数据分析和数据可视化实战（附代码）

手把手教你做meta分析及相关软件介绍（领取免费资料）

浅谈国际学术不端案例及如何避免学术不端问题的发生

剖析常见学术论文拒稿原因及解决方法

（生信实操）用R语言处理芯片数据及结果可视化分析

生信快速入门—R语言与统计分析

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

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.