「转录组」WGCNA实战原理两不误

这篇具有很好参考价值的文章主要介绍了「转录组」WGCNA实战原理两不误。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

表达矩阵的处理

后续分析所用到的数据,均为FPKM标准化后的表达矩阵。

从流程上对WGCNA进行解读

1)当对芯片数据或者RNA-Seq等数据完成分析之后,我们可以得到一张表达矩阵

一行为一个样本,一列为一个gene。

2)读取表达矩阵之后,对其进行adjacency矩阵的计算

adjacency矩阵,基于gene之间相关性的矩阵。每一个单元代表了2个gene之间的关联性(similarity)。

adjacency矩阵的构建,涉及到软阈值的选择,即构建一个幂函数,选择一个指数(power),强化强相关性,弱化弱相关性。

「adjacency矩阵的拓展」

adjacency矩阵有2种计算方式,

1)unsigned: a i j = ∣ c o r ( x x , x j ) ∣ β a_{ij} = |cor(x_{x}, x_{j})|^{β} aij=cor(xx,xj)β

2)signed: a i j = ( 0.5 ∗ ( 1 + c o r ) ) β a_{ij}=(0.5*(1+cor))^{β} aij=(0.5(1+cor))β

第一种情况,认为具有高负相关的gene之间,是有联系的。

第二种情况,认为具有高负相关的gene之间,是没有联系的。

举个例子,

假设cor = -1,β=2,unsigned情况下,计算得到的adjacency为1,即gene之间高度关联。signed情况下,计算得到的adjacency为0,即gene之间无关联。

3)TOM矩阵的构建

引入一个指标,即topological overlaps的计算,用于定义该gene是否有多个高关联性的gene(用于gene module的构建
C o r ( x i , x j ) − > T O M ( x i , x j ) Cor(x_{i},x_{j}) -> TOM(x_{i}, x_{j}) Cor(xi,xj)>TOM(xi,xj)

「TOM矩阵的拓展」

但是针对unsigned类型的adjacency矩阵,可能会出现无法判断几个基因之间的关联关系。

比如现有i,j,k 3个gene,计算得到它们之间的关联关系为 ( + , + , − ) (+,+,-) (+,+,)

1)i和j之间为正相关;2)i和k之间为正相关;3)j和k之间为负相关。

这就矛盾了,说明上面的关联关系可能受到了一些噪音的影响等等。

此时,引入signed TOM,其相较于unsigned TOM能够更好的解决gene之间的冲突情况。

Note:需要注意的是,当使用signed adjacency进行分析时,不会出现上述矛盾。

同时,这些矛盾作者都已经考虑到了,但还是有一定区别,

  • TOMsimilarity()中,默认设置为TOMType = "unsigned"
  • TOMsimilarityFromExpr()中,默认设置为TOMType = "signed"

4)TOM dissimilarity矩阵的构建

TOM dissimilarity矩阵代表了某个gene与其他gene之间的距离。

用上述矩阵进行聚类分析,得到gene module。

5)初步构建gene module

以TOM dissimilarity矩阵作为输入,进行聚类分析。

代码如下,

# Turn data expression into topological overlap matrix
power=sft$powerEstimate   # 使用sft$powerEstimate调用预估出的软阈值
TOM = TOMsimilarityFromExpr(datExpr, power = power)
dissTOM = 1-TOM
# Plot gene tree
geneTree = hclust(as.dist(dissTOM), method = "average");   # 用于后续cutreeDynamic(),对gene tree进行裁剪
# pdf(file = "3-gene_cluster.pdf", width = 12, height = 9);
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
     labels = FALSE, hang = 0.04);
# dev.off()

6)使用cutreeDynamic()进行聚类分析的优化

使用pairwise eigengene,进一步计算得到eigengene dissimilarity,用于聚类分析,筛选指标,最终合并gene module。

Note:eigengene为gene表达模式的指标(PCA降维得到的第一个主成分)

cutreeDynamic()代码如下,

  • cutreeDynamic
  • labels2colors
  • plotDendroAndColors
# Module identification using dynamic tree cut
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2, 
                            pamRespectsDendro = FALSE, minClusterSize = 30);
table(dynamicMods)
length(table(dynamicMods)) 
# Convert numeric labels into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
# Plot the dendrogram and colors underneath
# pdf(file = "4-module_tree.pdf", width = 8, height = 6);
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",dendroLabels = FALSE,
                    hang = 0.03,addGuide = TRUE, guideHang = 0.05,main = "Gene dendrogram and module colors")
# dev.off()

绘制图如下,

「转录组」WGCNA实战原理两不误

gene module合并代码如下 -> 基于cutreeDynamic()的分析结果进行聚类分析的gene module的合并

  • mergeCloseModules:合并gene module
  • plotDendroAndColors
# Merge close modules
MEDissThres=0.25
abline(h=MEDissThres, col = "red")
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3) 
mergedColors = merge$colors  
mergedMEs = merge$newMEs  
# Plot merged module tree
# pdf(file = "5-merged_Module_Tree.pdf", width = 12, height = 9)  
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), 
                    c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = FALSE, 
                    hang = 0.03, addGuide = TRUE, guideHang = 0.05)  
# dev.off()
# merge$oldMEs,为数据框,行为样本,列为对应的gene module,其中的数值代表了它们之间的关联度
write.table(merge$oldMEs,file="oldMEs.txt");
write.table(merge$newMEs,file="newMEs.txt");

绘制图如下,可以看到有一些gene module被合并了,

「转录组」WGCNA实战原理两不误

7)将gene module与表型特征相联系

使用标准化的eigengene进行计算。
C o r ( M E s , T r a i t s ) Cor(MEs, Traits) Cor(MEs,Traits)
代码如下,

  • moduleEigengenes,计算每个gene module的特征值
  • moduleTraitCor = cor(MEs, datTraits, use = "p")
  • moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
# Define numbers of genes and samples
nGenes = ncol(datExpr);
nSamples = nrow(datExpr);
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(datExpr, mergedColors)$eigengenes
MEs = orderMEs(MEs0)
# r$> MEs[1:5, 1:5]
#           MElightcyan1 MEdarkolivegreen MEgreenyellow       MEcyan     MEwhite
# LH38_Z1_1    0.3093871      -0.01665288     0.1583588 -0.017809267 -0.01305706
# LH38_Z1_2    0.3318247      -0.01679628     0.1629585 -0.010450719  0.04613664
# LH38_Z1_3    0.3087032      -0.01762908     0.1527085 -0.018456470 -0.03661408
# LH38_Z2_1    0.3050361      -0.01431782     0.1589569 -0.020199995  0.06161998
# LH38_Z2_2    0.3422982      -0.01466196     0.1722949 -0.006663298  0.08192540

# 一些数据处理部分
# Read microbial data as traits
bac_traits = read.table("traits_file/b_order_234.txt", header = T, sep = "\t")
rownames(bac_traits) = bac_traits[, 1]
bac_traits = bac_traits[, -1]
# r$> bac_traits[1:5, 1:5]
#           Pseudomonadales Enterobacteriales Xanthomonadales Burkholderiales Verrucomicrobiales
# LH38_Z1_1      0.02120943       0.006338742      0.07261663      0.05920385         0.02674949
# LH38_Z1_2      0.04192444       0.009089757      0.06880071      0.05583164         0.02156440
# LH38_Z1_3      0.01393256       0.004525862      0.06961207      0.05100152         0.02189402
# LH39_Z1_1      0.11288033       0.013045132      0.07138692      0.04186106         0.01743154
# LH39_Z1_2      0.01214503       0.003410243      0.08236562      0.03733519         0.01953600


rownames(MEs) = paste(substr(rownames(MEs), 1, nchar(rownames(MEs))-1), rep(c("1", "2", "3"), 60), sep = "")
# sample names should be consistent in eigen genes and traits !!!!
bac_traits = bac_traits[match(rownames(MEs), rownames(bac_traits)), ]
table(rownames(MEs) == rownames(bac_traits))

# Calculate pearson correlation coefficients between module eigen-genes and traits
moduleTraitCor = cor(MEs, bac_traits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
write.table(moduleTraitCor,file="moduleTrait_correlation.txt");
write.table(moduleTraitPvalue,file="moduleTrait_pValue.txt");

结果图可视化代码如下,

sizeGrWindow(10,6)
# Will display correlations and their p-values
# 合并
textMatrix =  paste(signif(moduleTraitCor, 2), "\n(",
                    signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
# pdf("module-traits-bacteria-order.pdf", width = 100, height = 30)
par(mar = c(15, 12, 5, 5));
# Display the correlation values within a heatmap plot
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(bac_traits),
               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"))
# dev.off()

结果图如下,

Note:右边的图例,代表相关性程度

「转录组」WGCNA实战原理两不误

8)Hub gene的鉴定

2种方法,

  • 与目标module有强关联性的gene(gene significance)
  • 使用module membership来鉴定关键gene, C o r ( i , M E ) Cor(i,ME) Cor(i,ME)

一般先绘制出以module membership为x,gene significance为y的散点图。

代码如下,

# 以相关性最高的模块为例
which(moduleTraitCor == max(moduleTraitCor, na.rm = TRUE), arr.ind = TRUE)
#            row col
# MEskyblue3  25  30

Nitrosomonadales <- as.data.frame(bac_traits[, 30])
names(Nitrosomonadales) = "Nitrosomonadales"

modNames = substring(names(MEs), 3) # 去除ME前缀

# 计算module membership
# 使用的是WGCNA自带cor函数,使用皮尔逊计算相关性
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
geneModuleMembership[1:5, 1:5]
#                MElightcyan1 MEdarkolivegreen MEgreenyellow      MEcyan       MEwhite
# Zm00001d001763 -0.045547074      -0.10334193   -0.20030129  0.22671189  0.1501121833
# Zm00001d001766  0.004752275      -0.01548396   -0.11130522 -0.22346048  0.0178421845
# Zm00001d001770 -0.286230379      -0.23340743   -0.39930773  0.22791746  0.0925346831
# Zm00001d001774  0.029505956      -0.06714112   -0.06323784 -0.34302212  0.1677206174
# Zm00001d001775 -0.029767437      -0.07642270   -0.13732818 -0.07876773 -0.0008678361

MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
MMPvalue[1:5, 1:5]

names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");

# 计算gene significance
geneTraitSignificance = as.data.frame(cor(datExpr, Nitrosomonadales, use = "p"));
# head(geneTraitSignificance)
# nrow(geneTraitSignificance)
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(Nitrosomonadales), sep="");
names(GSPvalue) = paste("p.GS.", names(Nitrosomonadales), sep="");


module = "skyblue3"
column = match(module, modNames);  # match(x, y),找到y中x的索引位置
moduleGenes = mergedColors==module;
# length(mergedColors)
# nrow(geneModuleMembership)
# nrow(geneTraitSignificance)
# table(mergedColors)
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 Nitrosomonadales",
                   main = paste("Module membership vs. gene significance\n"),
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

散点图结果如下,可以看到module membership值越高的gene,其gene significance也越高。

「转录组」WGCNA实战原理两不误

其他

绘制gene module eigengenes与samples之间的热图

library("pheatmap")

# Heatmap of old module eigen-genes and samples
pdf(file="oldMEs.pdf",heigh=80,width=20)
row.names(merge$oldMEs)=names(data0)  # oldMEs,是一个矩阵
pheatmap(merge$oldMEs,cluster_col=T,cluster_row=T,show_rownames=T,show_colnames=T,fontsize=6)
dev.off()

# Heatmap of new module eigen-genes and samples
# pdf(file="newMEs.pdf",heigh=60,width=20)
row.names(merge$newMEs)=names(data0)
pheatmap(merge$newMEs,cluster_col=T,cluster_row=T,show_rownames=T,show_colnames=T,fontsize=6)
# dev.off()

参考资料

[1] https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/

[2] https://github.com/PengYuMaize/Yu2021NaturePlants

[3] 跟着Nature Plants学数据分析:R语言WGCNA分析完整示例

[4] https://www.youtube.com/watch?v=BzYfg1lO3jw

[5] https://www.biostars.org/p/288153/

[6] https://peterlangfelder.com/2018/11/25/signed-or-unsigned-which-network-type-is-preferable/

[7] https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/TechnicalReports/signedTOM.pdf文章来源地址https://www.toymoban.com/news/detail-423969.html

到了这里,关于「转录组」WGCNA实战原理两不误的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包