heatmap_cell.jpg
今天我们复现一幅2020年发表在Cell上的热图。
DOI:10.1016/j.cell.2020.07.009
22
读图
Snipaste_2022-04-21_09-49-27.png
将不同样本的基因表达情况用热图展示并将倍数变化条形图展示在右侧，除此之外我们还可以把P值映射在条形图的颜色上，因为没有下载到原始数据我们将会用示例数据进行教学。
示例数据和代码领取
详见：https://mp.weixin.qq.com/s/tZ10XYjNCCNbygCPsG5bUQ
绘制
数据准备
数据文件有三个table，一个表达矩阵文件expr.csv，一个差异检验结果文件t_test.csv，这里用的是t检验，一个分组信息文件group.csv。表达矩阵中有一些空值，我们直接用0把空值替换掉，group表中将group转化为因子变量。接下来我们把倍数变化和P值信息都汇入到表达矩阵后方便后面绘图使用。
library(ComplexHeatmap)
library(circlize)
rm(list=ls())
expr=read.csv('expr.csv',header = T,row.names = 1)
t_test=read.csv('t_test.csv')
group=read.csv('group.csv',header = T,row.names = 1)
group$group=as.factor(group$group)
expr[is.na(expr)]=0
expr$FC=t_test$FC[match(rownames(expr),t_test$geneID)]#增加倍数变化列
expr$P.value=t_test$P.Value[match(rownames(expr),t_test$geneID)]#增加P值列
expr=expr[order(expr$FC,decreasing = T),]#按照倍数变化降序排列
heatmapdata=log(expr[1:30,1:70]+1)#取倍数变化前30的基因取对数作为热图数据
bardata=data.frame(FoldChange=expr$FC[1:30],row.names = rownames(expr)[1:30])#取FC数据做条形图
热图绘制
复杂热图需要用到Heatmap函数，其中的行注释和列注释都是需要用HeatmapAnnotation函数来生成，右侧的条形图实际上也属于行注释的一种，而注释信息的颜色映射则需要用到colorRamp2来生成颜色函数，我们首先生成行注释和列注释信息。
col_fun = colorRamp2(c( 0,0.5*max(heatmapdata),max(heatmapdata)), c( "white",'#FF99A3', "red"))#根据表达量的对数范围生成热图的颜色范围函数
col_fun_row=colorRamp2(c(min(expr$P.value[1:30]),0.05,max(expr$P.value[1:30])),c('#EC211C','#316CD7','blue'))#根据P值范围生成P值映射颜色范围
barcol=col_fun_row(expr$P.value[1:30])#生成所需基因的P值对应的颜色
names(barcol)=rownames(bardata[1:30,])
right_anno = rowAnnotation('FC(T/P)'=
row_anno_barplot(bardata[1:30,],
width = unit(3, "cm"),#条形图高度
bar_width = 0.8, #条形图条的宽度
outline = FALSE,
gp = gpar(fill = barcol,#条形图颜色
lwd=1)))#条形图边框宽度
topanno=HeatmapAnnotation(df=group,#列注释
border = T,
show_annotation_name = F,
col = list(group=c('Tumor'='#5CB85C',
'Normal'='#337AB7' )))
紧接着就是热图的绘制和图例的生成，在R中图例调整稍微有些复杂，建议生成图例后期在AI中进行调整。
Heatmap(heatmapdata,
col = col_fun,
cluster_columns = F,
cluster_rows = F,
row_names_side = 'left',#行注释置于左侧
top_annotation = topanno,
right_annotation = right_anno,
show_heatmap_legend = T,
show_column_names = F,
border = T,
rect_gp = gpar(col = "white", lwd = 1),#热图单元边框颜色和宽度
row_names_gp = gpar(fontsize = 12),
border_gp = gpar(lwd = 2),
column_split = group,#按照分组来分隔热图
column_gap = unit(2, "mm")#分隔间距
)
[图片上传失败...(image-80da51-1650511079983)]
热图基本的元素已经具备了，我们可以再为P值添加一个标签，再在AI中进行后期处理。
lgd = Legend(col_fun = col_fun_row,
title = "P.Value",
direction = "horizontal",
legend_width=unit(0.165,'npc'),
at = c(0,0.01,0.05,0.1))
draw(lgd,x = unit(0.715, "npc"), y = unit(0.925, "npc"))
Snipaste_2022-04-20_21-29-57.png
往期内容
[图片上传失败...(image-88523d-1650511079983)]
