本帖最后由 基迪奥-Jt桃 于 2021.3.22 10:01 编辑
热图是用来展示基因表达模式最直观的工具。早在2002年,就有大佬使用热图展示复杂网络潜在的模块性(modularity),如下,图A为假设的仅有11个结点(node)分成3个模块(红、绿、蓝)的小网络,结点之间连线上的数值称为topological overlap,用于描述两个结点的关联程度;而图B为使用图A相应的TOM(topological overlap matrix)矩阵绘制的聚类热图,聚类结果可直观展示出A网络中的3个模块。
(science, 2002) ,
基于以上思路,Steve Horvath等开发出了一套大家熟知的加权基因共表达网络分析方法(WGCNA, weighted gene co-expression network analysis)。WGCNA的核心步骤其实就是通过对TOM矩阵进行分层聚类,将表达模式相似的基因聚类到相同的模块,如下。图中的聚类树分支对应不同的基因模块(深色区块),模块内连通性较高的基因位于聚类树分支的尖端。
,, (BMC bioinformatics, 2008)
本文的主要内容是对WGCNA的分析结果进行可视化,包括基因共表达热图和模块表型相关性热图的个性化绘制。
WGCNA包是存放在CRAN的,可直接安装。但部分依赖的包放在Bioconductor,需单独安装。
[AppleScript] 纯文本查看 复制代码 #WGCNA包安装方法;
install.packages("WGCNA")
#Warning:
#dependencies ‘impute’, ‘preprocessCore’, ‘GO.db’ are not available
#我的电脑这里需要安装这3个包;
BiocManager::install("GO.db")
BiocManager::install("impute")
BiocManager::install("preprocessCore")
注意:大家可能之前已经安装过某些依赖的R包,具体需要安装那些包,可以根据安装时的警告信息逐一安装,原则是“缺什么,装什么”!
范例数据脚本下载
什么?到这里就卡住了?如果想学习R语言,或者想对WGCNA有更深入的了解,推荐你参加即将于2021年3月29日-4月2日举办的线上转录组培训班,课表如下:
报名方式:
方式1: 扫描下方二维码,填写信息报名
方式2:
客服:020-39341079或18054271626 邓小姐
1.计算TOM矩阵
[AppleScript] 纯文本查看 复制代码 # 导入WGCNA分析的结果数据;
load(file = "WGCNA_01_data.RData")
load(file = "WGCNA_02_networkConstruction-auto.RData")
ls()
#明确基因数和样本数;
nGenes <- ncol(datExpr)
#[1] 3600
nSamples <- nrow(datExpr)
#[1] 134
#重新计算TOM矩阵;
TOM <- TOMsimilarityFromExpr(datExpr, power = 6)
dissTOM <- 1-TOM
#对dissTOM矩阵进行指数转换,使中等强度的关系更容易在热图上展示出来;
plotTOM <- dissTOM^7
#将对角线数据设为NA;
diag(plotTOM) <- NA;
到这里,不知道你会不会有疑问,WGCNA包的作者为什么不直接用TOM矩阵而是用(1﹣x)7 画热图呢?
我这里用x=seq(from=0,to=1,by=0.02)生成0~1之间共51个数据视作相关系数,然后进行幂运算,指数取不同的整数,如下图,可以发现取指数后结果要小于原始值(0,1除外),靠近0端数据偏低,靠近1端数据急速升高。通过这种数据处理方式,使整体数值降低,但靠近1端的数据降幅较小。这样更符合无尺度网络的规律(幂律分布),少数核心结点关联性较强,大多数结点关联性较低,也就是大家常听说的80/20法则。
如果是1﹣x,这里相当于对x进行倒序排列,我们再看下(1﹣x)取指数后的结果,如下图,我们发现(1﹣x)6 相较于x6(x6 信息丢失较多)“温和”多了,甚至和(1﹣x)7 相差不大。因此,如果使用(1﹣x)6 作图,相对于直接使用TOM值画热图,并且可适当“压低”中等数值,突出展示接近于1的数据。
这里直接用热图直观看一下效果,如下图,我们可以发现,使用(1﹣x)6 或(1﹣x)7 更能“温和地”突出展示0.8以上的数据,在共表达网络中更容易展示强相关的核心基因。
2.共表达网络热图绘制
[AppleScript] 纯文本查看 复制代码 #使用所有基因绘制TOM热图,耗时较久,个人电脑约30min,不推荐运行;
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
[AppleScript] 纯文本查看 复制代码 #由于TOM值的范围为0~1,且TOM热图无颜色条图例,故可以通过以下方式反转颜色映射;
TOMplot((1-plotTOM), geneTree, moduleColors, main = "Network heatmap plot, all genes")
由于使用所有的基因绘图太慢了,可能需要10~30min无法执行其他命令,如下图。
因此,我们需要通过随机取样的方式降低热图的分辨率,也就是通过减少基因数量进行提速,通过这种方式可用较低的分辨率展示基因互作网路。
[AppleScript] 纯文本查看 复制代码 #设置抽取的基因数,建议设置为400~800之间,数值越小运行越快;
nSelect <- 400
# 为了反复重现结果, 这里设置了随机种子;
set.seed(10)
#使用sample()函数从3600个基因中随机抽取400个;
select <- sample(nGenes, size = nSelect)
head(select)
#[1] 2207 2239 3358 3230 2507 2227
#提取抽取基因相应的TOM矩阵;
selectTOM <- dissTOM[select, select]
selectTOM[1:6,1:6]
[AppleScript] 纯文本查看 复制代码 # 根据TOM矩阵重新对基因聚类;
selectTree <- hclust(as.dist(selectTOM), method = "average")
#提取选择基因的模块颜色;
selectColors <- moduleColors[select]
#[1] "grey" "green" "yellow" "black" "turquoise" "blue"
# 对TOM矩阵进行指数转化,使热图能展示更丰富的信息;
plotDiss <- selectTOM^7
diag(plotDiss) <- NA
#绘制热图,默认数值越大(越接近1)颜色越深(底色);
TOMplot(plotDiss, selectTree, selectColors,
main = "Network heatmap plot, selected genes")
[AppleScript] 纯文本查看 复制代码 #如果觉得颜色不好看可反转映射颜色;
TOMplot((1-plotDiss), selectTree, selectColors,
main = "Network heatmap plot, selected genes")
[AppleScript] 纯文本查看 复制代码 #自定义渐变色,一般从左到右,对应的数值越大;
mycol = colorRampPalette(c("Firebrick","orange","lemonchiffon"))(256)
mycol2 = colorRampPalette(c("lemonchiffon","orange","Firebrick"))(256)
#蓝色渐变;
mycol2 = colorRampPalette(c("darkblue","lightskyblue","azure"))(256)
#绿色渐变;
mycol2 = colorRampPalette(c("forestgreen","greenyellow","mintcream"))(256)
#使用自定义颜色;
TOMplot(plotDiss, selectTree, selectColors, col=mycol2)
当然如果对绘图结果不满意还可以使用Ai(Adobe illustrator)等矢量图编辑软件进一步美化,比如调整聚类树分支的粗细;颜色等,如下。
如果不用颜色标签表示特定模块,还可以使用重新着色图稿功能进一步调整模块的配色方案, 如下。
最终的调整效果如下:
当然还可以是这样的,也挺好看!
[AppleScript] 纯文本查看 复制代码 #如果直接使用tom矩阵作图呢?
tom <- 1-selectTOM
tomv <- tom^7
diag(tom) = NA;
diag(tomv) = NA;
#直接绘制tom热图,其实绘图结果非常相似,只不过模块不如上文的明显;
TOMplot(tom, selectTree, selectColors,col=mycol2)
[AppleScript] 纯文本查看 复制代码 #如果直接绘制指数化的tom矩阵,图中的信息较少,不推荐;
TOMplot(tomv, selectTree, selectColors,col=mycol2)
3.特征基因与表型性状的关系
[AppleScript] 纯文本查看 复制代码 #重新计算模块特征基因;
MEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
#提取表型数据框的weight列;
weight <- as.data.frame(datTraits$weight_g)
names(weight) <- "weight"
#将weight列添加到特征基因表格中;
MET = orderMEs(cbind(MEs, weight))
#单独绘制聚类树(dendrogram);
par(cex = 1.0)
plotEigengeneNetworks(MET, "Eigengene dendrogram", marDendro = c(0,4,2,0),plotHeatmaps = FALSE)
[AppleScript] 纯文本查看 复制代码 #单独绘制热图;
par(cex = 1.0,mar=c(5,5,2,2))
plotEigengeneNetworks(MET, "Eigengene adjacency heatmap", marHeatmap = c(5,5,2,2),plotDendrograms = FALSE, xLabelsAngle = 90)
[AppleScript] 纯文本查看 复制代码 #绘制特征基因与表型性状的关系组合图;
par(cex = 0.9)
par(mai=c(1,1,1,1))
plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2),
marHeatmap = c(5,5,2,2), cex.lab = 0.8, xLabelsAngle= 90)
好啦,今天的内容就分享到这里啦~
参考文献 Ravasz E, Somera A L, Mongru D A, et al. Hierarchical organization of modularity in metabolic networks[J]. science, 2002, 297(5586): 1551-1555. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis[J]. BMC bioinformatics, 2008, 9(1): 1-13.
本文作者:基迪奥-莫北
|