搜索
查看: 11536|回复: 10

初学WGCNA

[复制链接]

29

主题

131

帖子

1188

积分

金牌会员

Rank: 6Rank: 6

积分
1188
发表于 2017-1-23 14:40:00 | 显示全部楼层 |阅读模式
本帖最后由 anlan 于 2017-1-24 08:39 编辑

最近看了博主的视频以及考虑一些自身需要,尝试着对WGCNA进行了一些学习,内容几乎不是原创,来自于一些公司的论坛以及网上的博客,下面是我对WGCNA的一些总结:

WGCNA分析步骤可以分为以下几步:
1. 通过数据处理,将基因间的相互作用关系强度符合无尺度分布
2. 将基因分类,并把表达模式相似的基因归为一个模块
3. 研究模块,找出跟我们研究相关的模块
4. 研究模块内基因之间的调控关系

从上四步可以看出,WGCNA应该就是:处理大样本,筛选出一定量表达模式相同的基因,进而研究基因之间的关系。感觉就是类似PPI,但是是没有数据库支持的PPI,属于预测范畴。

对于WGCNA的一些概念这里不说了,因为我也说不清楚,刚学。。。
测试数据就用博主博客里的那个56个样本的表达矩阵,http://www.bio-info-trainee.com/2297.html,我主要是对R代码进行个人理解的解释,以及一些分析策略。

R代码如下:

第一步:对表达矩阵进行处理,由于这个是测试数据,我就不进行筛选差异基因了。虽然我觉得不应该筛选差异基因再去做WGCNA,但是有些论文以及一些期刊是这么做的。。。我觉得只要进行变异系数过滤,低表达量过滤等等就够了。。。下面我就同博主视频里面的代码进行过滤,实在是原始数据基因数有点多了。。然后取3000个基因来分析
[AppleScript] 纯文本查看 复制代码
library(WGCNA)

data <- read.table(file = "tmp.txt", sep = "\t", stringsAsFactors = F)
library(reshape2)
fpkm <- dcast(data, formula = V2~V1)
rownames(fpkm) <- fpkm[,1]
fpkm <- fpkm[,-1]
names(fpkm) <- unlist(lapply(names(fpkm), function(x){
  tmp <- strsplit(x, '_')[[1]][2]
  tmp <- strsplit(tmp, '\\.')[[1]][1]
})
)
data_matrix_mv = t(fpkm[order(apply(fpkm,1,mad), decreasing = T)[1:3000],])
#data_matrix_mv <- data_matrix_mv[1:20,]


第二步:筛选合适的阈值
其实就是选一个合适的值,使得样本变为一个不尺度分布,即少部分基因处于绝对优势的位置(表达量高),大部分处于劣势(表达量低)。WGCNA用的是幂指数的方式,对基因之间的相关系数取一定的幂指数,最终确定一个合适的阈值,结果如图1所示。从左边那张图可以看出,
[AppleScript] 纯文本查看 复制代码
#2. 选择合适的阀值
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(data_matrix_mv, powerVector = powers, verbose = 5)
# Plot the results:
##sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
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");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
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")
图1左边那张图表示,红线附近的值都是可以作为阈值的,R^2已经在0.9左右了。如果想知道WGCNA建议的阈值,可以输入sft,列表的第一项值就是估计的最优值。。。也就是后续代码里面的power值。
第三步:
筛选到合适的阈值后,我们就可以构建网络图了(也就是 模块图)。官网给出了两种方法,一个是一步法:简单明了,一步出结果。另一个是分步法:可以细调一系列的参数,以达到满意的结果。官网也没说哪种适合什么条件,都可以用吧。具体参数我就不列了,可以看这个博客:http://blog.sina.com.cn/s/blog_61f013b80101lcpr.html 讲比我详细多了。
[AppleScript] 纯文本查看 复制代码
#3. 一步法网络构建:One-step network construction and module detection
net = blockwiseModules(data_matrix_mv, power = 6, maxBlockSize = 5000,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.45,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "AS-green-FPKM-TOM",
                       verbose = 3)
table(net$colors)

table(net$colors)可以查看不同模块的颜色以及每个模块里的基因数目。我为了减少模块数,所以mergeCutHeight用了0.45,这个根据自身需求来吧,一般都是0.25就行了。然后下面的代码就可以展示出网络图了,如图2所示。
[AppleScript] 纯文本查看 复制代码
#4. 绘画结果展示
# 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)

然后我们需要把一些参数保留,以用于后续的分析,如下:
[AppleScript] 纯文本查看 复制代码
#5.结果保存
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
table(moduleColors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]
第四步:
接下来的就是不同的策略的选择了。
方法1:最简单的就是对于每个模块进行GO以及KEGG富集,因为从理论上来说,每个模块都应该都有生物学意义的,因此可以看每个模块的GO以及KEGG注释的结果,看看都富集到哪个通路,再结合每个模块的特征值(查看MEs即可)来确定哪个模块的基因作为后续研究的切入点。
方法2:如果再有前人研究或者自己的调查的结果,知道某几基因是自己研究方向的重要基因,可以通过这几个基因在哪些或者某个模块,来确定接下来的选择哪个模块作为研究重点。
方法3:就像官网说的,结合模块与样本性状的相关性来分析,但一般都没有样本性状数据的,我猜。。。
方法4:虽然我觉得有点扯,但是我还是会选择用,至少在论文里看到过。。就是结合模块与样本间的相关系数。。如图3所示。。
[AppleScript] 纯文本查看 复制代码
MEs0 <- moduleEigengenes(data_matrix_mv, moduleColors)$eigengenes
MEs <- orderMEs(MEs0)
datTraits <- diag(56)
datTraits <- as.data.frame(datTraits)
rownames(datTraits) <- rownames(data_matrix_mv)
colnames(datTraits) <- rownames(data_matrix_mv)
nSamples = nrow(data_matrix_mv)
moduleTraitCor <- cor(MEs, datTraits, use = "p")
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(4.5, 8.5, 3, 1))
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"))

由于样本数量有点多,所以会不太清晰。。y轴是模块名,x轴是样本名。。每个模块都跟每个样本有一个空格,空格的上部分是模块与样本的相关系数(代码里面我构建了一个标量矩阵代表样本,然后求了模块和样本的相关系数),空格的下部分是一个p值,我猜是代表相关的显著性吧(不知道理解对不对,但是反正都是挑p值越小的越好)。。由于我不知道这个测试数据的样本信息各个代表咋样的实验,但是如果是有生物学重复实验的话,有人是这样做的:在3个生物学重复中有2个样本在某个模块是显著相关的话,就挑选作为后续研究。。。。有的人是这样的:关注与某一个或者几个样本,看看有哪个模块与这些样本是显著相关的。。。。还有的人。。我也不知道了。。文献还没怎么看。。
但是我觉得还是有点扯,有人可以能说说这样做的好处以及意义吗。。。其他策略嘛,网上也有,论文里面也有,都可以看看。。看数据选吧。。。

第五步:调选完模块后,我们可以做个该模块的特征值所对应的柱状图和聚类图,我这里就不做了(其实是我没查到WGCNA里有无该内置函数,反正别人是做了。。实在找不到函数的话,只好自己画了)
然后就可以做一张TOM图,如图4所示
[AppleScript] 纯文本查看 复制代码
#1. 可视化全部基因网络
# 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(data_matrix_mv, power = 12);
# 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")

第六步:
导出网络图到cytoscape,记得在代码里面修改自己所需要模块(颜色)。。然后再筛选下txt表格里面的weight值(应该就是TOM值?),可以选择硬阈值筛,还有种说法是用软阈值筛。。等我搞懂了再来修改。。。最后就导入cytoscape作图即可。。
[AppleScript] 纯文本查看 复制代码
#6. 导出网络到Cytoscape
# Recalculate topological overlap if needed
TOM = TOMsimilarityFromExpr(data_matrix_mv, power = 6);
# Read in the annotation file
# annot = read.csv(file = "GeneAnnotation.csv");
# Select modules需要修改,选择需要导出的模块颜色
modules = c("turquoise");
# Select module probes选择模块探测
probes = colnames(data_matrix_mv)
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]);


最后最后重点就是,请教各位做过WGCNA的大神们。。。你们有人用过WGCNA来处理质谱数据吗。。。或者看到用WGCNA处理质谱数据发表的文章吗?
无论发表时间,多多益善。。。我现在才找了没几篇。。真的好少诶。。。有的话,发个文章名即可。。多谢




本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:在R里面实现read.table(file, sep = &quot;\t&quot;, fill=TRUE)的效果
下一篇:【直播】我的基因组(一):直播的目的及意义
回复

使用道具 举报

0

主题

5

帖子

75

积分

注册会员

Rank: 2

积分
75
发表于 2017-2-16 21:58:30 | 显示全部楼层
赞一个,现在正在学,新手上路,没有老司机带
回复 支持 反对

使用道具 举报

6

主题

36

帖子

453

积分

中级会员

Rank: 3Rank: 3

积分
453
发表于 2017-3-26 21:13:39 | 显示全部楼层
给老哥赞一个,帮助很大,稳
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1188

积分

金牌会员

Rank: 6Rank: 6

积分
1188
 楼主| 发表于 2017-3-26 23:10:50 | 显示全部楼层
017中大普外 发表于 2017-3-26 21:13
给老哥赞一个,帮助很大,稳

这个是我刚开始初略写的,自己也不太懂。我最近写了个帖子,代码讲的比较详细,可以看看那个,帖子里面那个链接的内容更赞!别人总结的很好,也有实例文件。
回复 支持 反对

使用道具 举报

2

主题

17

帖子

145

积分

注册会员

Rank: 2

积分
145
发表于 2017-5-23 23:18:11 | 显示全部楼层
新手上路,没有老司机带,赞一个
回复 支持 反对

使用道具 举报

2

主题

41

帖子

365

积分

中级会员

Rank: 3Rank: 3

积分
365
发表于 2017-8-29 23:36:54 | 显示全部楼层
谢谢分享,很精彩的帖子,按照你的步骤把图做出来了~
还需要要时间理解每一步骤代码的具体用法以及含义~
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-8-30 20:19:48 | 显示全部楼层
用你的这个发公众号了哦,生信菜鸟团
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1188

积分

金牌会员

Rank: 6Rank: 6

积分
1188
 楼主| 发表于 2017-8-30 22:24:47 | 显示全部楼层
学习最快乐 发表于 2017-8-30 20:19
用你的这个发公众号了哦,生信菜鸟团

这个太老了~~!要不要过几天我重新写一份吧,我找篇文献数据来一遍
回复 支持 反对

使用道具 举报

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-8-31 13:19:54 | 显示全部楼层
anlan 发表于 2017-8-30 22:24
这个太老了~~!要不要过几天我重新写一份吧,我找篇文献数据来一遍

好呀~~waiting for u
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1188

积分

金牌会员

Rank: 6Rank: 6

积分
1188
 楼主| 发表于 2017-9-11 13:44:42 | 显示全部楼层

发你私信咯
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-7-22 11:12 , Processed in 0.117905 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.