如何绘制基因共表达热图？

基迪奥-Jt桃

热图是用来展示基因表达模式最直观的工具。早在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包，具体需要安装那些包，可以根据安装时的警告信息逐一安装，原则是“缺什么，装什么”！

范例数据脚本下载

https://www.omicshare.com/forum/thread-6999-1-1.html

什么？到这里就卡住了？如果想学习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.

本文作者：基迪奥-莫北           

2012csc

顺利毕业凹

谢谢分享

给祖国添麻烦了

明明昊YG

十分有用

2012csc

wangyankun

kiz2011

吴丹丹

给祖国添麻烦了

给祖国添麻烦了

给祖国添麻烦了

mengtw

顶！！！！！！！！！！！！！！！！！！！

wangyankun

wangyankun

2012csc

流澜灵泽

mengtw

顶！！！！！！！！！！！！

wangyankun

2012csc