查看: 1578|回复: 19

[R语言] 如何绘制基因共表达热图?

  [复制链接]

迅猛龙

Rank: 8Rank: 8

主题
214
注册时间
2020.6.16
在线时间
129 小时

发表于 2021.3.22 10:01:31 | 显示全部楼层 |阅读模式
本帖最后由 基迪奥-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:
发送姓名、学校、电话到邮箱contact@genedenovo.com,主题注明“XXX培训班”

客服: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.

本文作者:基迪奥-莫北           

本帖子中包含更多资源

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

x
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
32 小时

发表于 2021.3.22 10:49:18 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2020.12.2
在线时间
1 小时

发表于 2021.3.22 14:33:39 | 显示全部楼层
谢谢分享
又是学习的一天
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2017.6.6
在线时间
110 小时

灌水之王


发表于 2021.3.22 16:21:43 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2020.12.18
在线时间
3 小时

发表于 2021.3.22 20:01:35 | 显示全部楼层
十分有用
我开心
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
32 小时

发表于 2021.3.23 09:45:38 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
28 小时

发表于 2021.3.25 08:16:13 | 显示全部楼层
加油,加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
1
注册时间
2016.8.25
在线时间
62 小时

发表于 2021.3.25 16:05:28 | 显示全部楼层
假期愉快
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2020.12.21
在线时间
4 小时

发表于 2021.3.29 00:03:39 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2017.6.6
在线时间
110 小时

灌水之王


发表于 2021.3.30 14:14:18 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2017.6.6
在线时间
110 小时

灌水之王


发表于 2021.3.30 14:41:32 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

中华鲟

Rank: 5Rank: 5

主题
0
注册时间
2017.6.6
在线时间
110 小时

灌水之王


发表于 2021.3.30 15:21:40 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.1.28
在线时间
4 小时

发表于 2021.4.1 15:08:59 | 显示全部楼层
顶!!!!!!!!!!!!!!!!!!!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
28 小时

发表于 2021.4.2 09:10:48 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
28 小时

发表于 2021.4.4 10:44:02 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
32 小时

发表于 2021.4.5 11:40:59 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.3.29
在线时间
3 小时

发表于 2021.4.5 16:49:08 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.1.28
在线时间
4 小时

发表于 2021.4.7 09:11:14 | 显示全部楼层
顶!!!!!!!!!!!!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2017.9.8
在线时间
28 小时

发表于 2021.4.8 08:27:04 | 显示全部楼层
加油,加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
0
注册时间
2020.11.21
在线时间
32 小时

发表于 2021.4.8 10:42:20 | 显示全部楼层
坚持就是胜利!
回复

使用道具 举报

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

本版积分规则

快速回复 返回顶部 返回列表