专栏首页R语言交流中心R语言实现加权共表达网络分析

R语言实现加权共表达网络分析

WGCNA(Weighted GeneCo-Expression Network Analysis,加权共表达网络分析)分析方法旨在寻找协同表达的基因模块(module),并探索基因网络与关注的表型之间的关联关系,以及网络中的核心基因。我们今天介绍下在R语言如何实现WGCNA,此包还有一个限制那就是样本总数必须大于15。

首先需要安装WGCNA包:首先需要先通过bioconductor安装相关的依赖包:‘impute’,‘preprocessCore’, ‘GO.db’, ‘AnnotationDbi’,然后install.packages(‘WGCNA’)。

包的操作实例数据来自于其官网(https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/)提供的数据:

接下来我们看下如何进行一步一步的操作。

第一部分数据的预处理

首先我们需要对输入的数据进行评估,需要用到函数goodSamplesGenes。

其函数具体的参数个人建议默认应该是最好的。除非有特殊需求。具体的参数就不一一介绍了。我们直接看实例:

library(WGCNA)
options(stringsAsFactors= FALSE)
femData= read.csv("G:/LiverFemale3600.csv")
datExpr0= as.data.frame(t(femData[, -c(1:8)]))##转置矩阵
names(datExpr0)= femData$substanceBXH;
rownames(datExpr0)= names(femData)[-c(1:8)];

gsg =goodSamplesGenes(datExpr0)

判断是否符合要求,如果不符合进行下一步的填充数据,或者去除不合格数据。

if(!gsg$allOK)
{
#Optionally, print the gene and sample names that were removed:
if(sum(!gsg$goodGenes)>0)
printFlush(paste("Removinggenes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ",")));
if(sum(!gsg$goodSamples)>0)
printFlush(paste("Removingsamples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ",")));
#Remove the offending genes and samples from the data:
datExpr0= datExpr0[gsg$goodSamples, gsg$goodGenes]
}

接下来就是对所有的样本进行聚类。查看是否存在异常值。这里用到了层析聚类算法hclust。此算法我们就不详细介绍了,网上一搜一大把。具体我们看下实例:

sampleTree= hclust(dist(datExpr0), method = "average");
# Plotthe sample tree: Open a graphic output window of size 12 by 9 inches
# Theuser should change the dimensions if the window is too large or too small.
sizeGrWindow(12,9)
#pdf(file= "Plots/sampleClustering.pdf", width = 12, height = 9);
par(cex= 0.6);
par(mar= c(0,4,2,0))
plot(sampleTree,main = "Sample clustering to detect outliers", sub="",xlab="", cex.lab = 1.5,
cex.axis= 1.5, cex.main = 2)

从上图,我们可以看出有一个异常样本。我们可以手动删除这些异常样本。同时呢,为了解决有时候不好判断的情况,我们也可以设置阈值,利用算法直接将异常值进行筛选。实例如下:

 
# Plota line to show the cut
abline(h= 15, col = "red");
#Determine cluster under the line
clust= cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
table(clust)
#clust 1 contains the samples we want to keep.
keepSamples= (clust==1)
datExpr= datExpr0[keepSamples, ]
nGenes= ncol(datExpr)
nSamples= nrow(datExpr)

接下来我们引入临床数据,可视化展示临床的量化指标,同时综合以上的聚类图。此处用到一个重要的函数numbers2colors。

此函数主要对数值化的参数进行高低的颜色标记,形成相应的热图。我们直接看下实例:

traitData= read.csv("G:/ClinicalTraits.csv");
dim(traitData)
names(traitData)
#remove columns that hold information we do not need.
allTraits= traitData[, -c(31, 16)];
allTraits= allTraits[, c(2, 11:36) ];
dim(allTraits)
names(allTraits)
# Forma data frame analogous to expression data that will hold the clinical traits.
femaleSamples= rownames(datExpr);
traitRows= match(femaleSamples, allTraits$Mice);
datTraits= allTraits[traitRows, -1];
rownames(datTraits)= allTraits[traitRows, 1]
#Re-cluster samples
sampleTree2= hclust(dist(datExpr), method = "average")
#Convert traits to a color representation: white means low, red means high, greymeans missing entry
traitColors= numbers2colors(datTraits, signed = FALSE);
# Plotthe sample dendrogram and the colors underneath.
plotDendroAndColors(sampleTree2,traitColors,
groupLabels= names(datTraits),
main ="Sample dendrogram and trait heatmap")

第二部分 网络的构建

首先需要用到一个函数pickSoftThreshold,用来计算最优软阈值的。所谓软阈值就是指的一个构建拓扑网络的阈值范围,通过设置范围评估最低能量的拓扑结构。

其中主要的参数powerVector 的设置,如果想要输出更多的运行日志信息可以增大verbose的值。接下来我们看下实例:

#Choose a set of soft-thresholding powers
powers= c(c(1:10), seq(from = 12, to=20, by=2))
# Callthe network topology analysis function
sft =pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plotthe results:
sizeGrWindow(9,5)
par(mfrow= c(1,2));
cex1 =0.9;
# Scale-freetopology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="SoftThreshold (power)",ylab="Scale Free Topology Model Fit,signedR^2",type="n",
main =paste("Scale independence"));
text(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# thisline corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Meanconnectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1],sft$fitIndices[,5],
xlab="SoftThreshold (power)",ylab="Mean Connectivity", type="n",
main =paste("Mean connectivity"))
text(sft$fitIndices[,1],sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

图中左侧的图可以看出在6的时候后面趋于稳定,那么这个最优值就是6。或者可以直接提取最优值。实例:

sft$powerEstimate

接下来就是网络结构的构建,需要用到函数blockwiseModules。

其中主要的参数:

corType网络相关性分析的算法

TOMType拓扑结构构建算法。

minModuleSize每个模块最少的基因个数

具体的实例:

net =blockwiseModules(datExpr, power = 6,
TOMType= "unsigned", minModuleSize = 30,
reassignThreshold= 0, mergeCutHeight = 0.25,
numericLabels= TRUE, pamRespectsDendro = FALSE,
saveTOMs= TRUE,
saveTOMFileBase= "femaleMouseTOM",
verbose= 3)
# opena graphics window
sizeGrWindow(12,9)
#Convert labels to colors for plotting
mergedColors= labels2colors(net$colors)
# Plotthe dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]],mergedColors[net$blockGenes[[1]]],
"Modulecolors",
dendroLabels= FALSE, hang = 0.03,
addGuide= TRUE, guideHang = 0.05)

从上图我们可以看出各个模块的分布以及聚类的情况,我们也可以将各个数据提取出来。示例如下:

moduleLabels= net$colors
moduleColors= labels2colors(net$colors)
MEs =net$MEs;
geneTree= net$dendrograms[[1]]

以上是针对数量不大的网络构建,如果基因数量太大没法一次性构建好网络那么需要函数blockwiseModules中的saveTOMFileBase= "femaleMouseTOM-blockwise"。可以把每组基因数控制在3000以内分组构建最后整合。其参数和前面一样,只是运行结果展示可以不同。实例如下:

bwnet= blockwiseModules(datExpr, maxBlockSize = 2000,
power= 6, TOMType = "unsigned", minModuleSize = 30,
reassignThreshold= 0, mergeCutHeight = 0.25,
numericLabels= TRUE,
saveTOMs= TRUE,
saveTOMFileBase= "femaleMouseTOM-blockwise",
verbose= 3)
bwLabels= matchLabels(bwnet$colors, moduleLabels);
#Convert labels to colors for plotting
bwModuleColors= labels2colors(bwLabels)
# opena graphics window
sizeGrWindow(6,6)
# Plotthe dendrogram and the module colors underneath for block 1
plotDendroAndColors(bwnet$dendrograms[[1]],bwModuleColors[bwnet$blockGenes[[1]]],
"Modulecolors", main = "Gene dendrogram and module colors in block 1",
dendroLabels= FALSE, hang = 0.03,
addGuide= TRUE, guideHang = 0.05)
# Plotthe dendrogram and the module colors underneath for block 2
plotDendroAndColors(bwnet$dendrograms[[2]],bwModuleColors[bwnet$blockGenes[[2]]],
"Modulecolors", main = "Gene dendrogram and module colors in block 2",
dendroLabels= FALSE, hang = 0.03,
addGuide= TRUE, guideHang = 0.05)

接下来我们合并两个分组与前面一步构建的网络进行对比。需要用到下面的函数plotDendroAndColors。实例如下:

sizeGrWindow(12,9)
plotDendroAndColors(geneTree,
cbind(moduleColors,bwModuleColors),
c("Singleblock", "2 blocks"),
main ="Single block gene dendrogram and module colors",
dendroLabels= FALSE, hang = 0.03,
addGuide= TRUE, guideHang = 0.05)

此包还提供了一步一步的进行网络构建以及模块的识别。我们就不展开讲了。

第三部分 寻找重要的基因

首先是前面数据的再次利用:

#Define numbers of genes and samples
nGenes= ncol(datExpr);
nSamples= nrow(datExpr);
#Recalculate MEs with color labels
MEs0 =moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs =orderMEs(MEs0)
moduleTraitCor= cor(MEs, datTraits, use = "p");
moduleTraitPvalue= corPvalueStudent(moduleTraitCor, nSamples);

然后我们对数据进行可视化的展示,需要用到函数labeledHeatmap。

具体的参数我们不一一介绍了,基本都是用的默认参数。直接看下实例:

sizeGrWindow(10,6)
# Willdisplay correlations and their p-values
textMatrix= paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue,1), ")", sep = "");
dim(textMatrix)= dim(moduleTraitCor)
par(mar= c(6, 8.5, 3, 3));
#Display the correlation values within a heatmap plot
labeledHeatmap(Matrix= moduleTraitCor,
xLabels= names(datTraits),
yLabels= names(MEs),
ySymbols= names(MEs),
colorLabels= FALSE,
colors= greenWhiteRed(50),
textMatrix= textMatrix,
setStdMargins= FALSE,
cex.text= 0.5,
zlim =c(-1,1),
main =paste("Module-trait relationships"))

上图中每行代表一个模块,每列代表一个临床指标。每个单元格代表一个相关系数和P值。

既然上面是所有的临床指标的可视化,那么如何单个指标单个模块进行相关性分析呢。我们直接上实例:

#Define variable weight containing the weight column of datTrait
weight= as.data.frame(datTraits$weight_g);
names(weight)= "weight"
#names (colors) of the modules
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="");
 
module= "brown"
column= match(module, modNames);
moduleGenes= moduleColors==module;
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)

第五部分 数据的导出

此包主要的是将网络数据导出到Cytoscape和VisANT。我们主要关注Cytoscape,下面就是导出的实例:

#Recalculate topological overlap if needed
TOM =TOMsimilarityFromExpr(datExpr, power = 6);
# Readin the annotation file
annot= read.csv(file = "G:/GeneAnnotation.csv");
#Select modules
modules= c("brown", "red");
#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 edgeand node list files Cytoscape can readcyt = exportNetworkToCytoscape(modTOM,
edgeFile = paste("CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),
nodeFile = paste("CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),
weighted = TRUE,
threshold = 0.02,
nodeNames = modProbes,
altNodeNames = modGenes,
nodeAttr = moduleColors[inModule]);

最终获得两个文件

那么这两个文件怎么用呢?接下来给大家介绍下如何使用:首先需要到官网(https://cytoscape.org/)下载相应的软件。

接下来导入我们的数据,我们只需要将edge的数据导入即可。

然后是选择对应的起始节点,目标节点以及边权重。

确定后就会形成相应的网络图,当然这个网络图还需要我们进行后期的美化。然后我们就可以将我们的节点的注释文件导入。

选择对应的网络图。

其他的筛选、局部绘图啥的在这里就不做一一介绍了。如果感兴趣可以参照官网的手册自行研究。

欢迎互相学习交流!

本文分享自微信公众号 - R语言交流中心(R_statistics),作者:one sand

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-08-07

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • R语言调用C++程序

    R语言在数据处理方面很是强大,然而也面临着很多的局限性。比如图像的分析处理,大数据的运算效率问题。今天我们介绍R语言和高效语言结合的一种方法:

    一粒沙
  • R语言利用GOplot实现功能分析可视化

    生信中大家都不陌生GO分析,然而如何将分析结果进行可视化展示是我们苦恼的问题,大部分都是画个Bar图或者列个表格啥的。今天我们给大家介绍一个可以实现功能分析可视...

    一粒沙
  • R语言实现测序原始数据的文件转化

    了解测序的同志们应该都知道有很多格式的原始文件,同时在转化过程中很是麻烦。今天我们给大家介绍一个R包,它可以进行对原始数据的读取,同时并且可以导出时进行转化格式...

    一粒沙
  • 从输入url到页面返回到底发生了什么

    1. 前言 Google应该是开发者平日里用得最多的网站之一,今早笔者在浏览器地址栏里键入www.google.com的时候,突然想了解下这背后的网络通信过程究...

    潘成涛
  • "Goole项目托管"及"CodePlex发布开源项目"要点

    一.google项目托管相对比较容易 http://code.google.com/ 先注册一个gmail邮箱,然后参考孟子的这篇文章http://blog.c...

    菩提树下的杨过
  • 计算机网络基础知识总结

     为了使不同计算机厂家生产的计算机能相互通信,在更大范围内建立计算机网络,国际标准化组织(ISO)在1978年提出了“开放系统互联参考模型”,即著名的OSI/R...

    mathor
  • Rasa 聊天机器人专栏(四):消息和语音通道

    如果您在本地计算机(即非服务器)上进行测试,则需要使用[ngrok]()。这为您的机器提供了域名,以便Facebook,Slack等知道将消息发送到本地计算机的...

    磐创AI
  • 浅析Facebook LibraBFT与比原链Bystack BBFT共识

    它是区块链的根基。无论公链或是联盟链,共识机制都从基础上限制了区块链的交易处理能力和扩展性。

    比原链Bytom
  • 浅析Facebook LibraBFT与比原链Bystack BBFT共识

    它是区块链的根基。无论公链或是联盟链,共识机制都从基础上限制了区块链的交易处理能力和扩展性。

    比原链Bytom
  • 一个小例子秒懂ThreadLocal使用及原理

    最近一个小伙伴把项目中封装的日期工具类用在多线程环境下居然出了问题,来看看怎么回事吧

    Java识堂

扫码关注云+社区

领取腾讯云代金券