转录组数据分析笔记（9）——pheatmap绘制差异表达基因热图

前面说到怎么用ggplot做一个火山图来查看各个基因的表达情况，火山图是以log2FC值为横坐标，以-log10（FDR）值作为纵坐标，将所有的基因都做了点状图。虽然能比较直观地看到所有基因表达情况，但我们真正感兴趣的是处理后差异表达的基因。因此，我们也可以通过前面得到的表达矩阵获得差异表达的基因名，对raw count数据进行提取和均一化，然后做一个差异基因的热图，能更直观地看到差异基因在各个样本中的上调下调情况。

做热图我们用的最多的R包是pheatmap，可以直接用biocmanager下载。

后面所有Rstudio操作都是在同一个Rproject中进行，引用的变量如果不理解就翻前面的笔记，最后我会把一整个转录组下游分析的R流程代码写成一个文件上传到github备份。

1. 热图介绍

废话不多说，先上一段热图的定义介绍：热图是用来对采集的因子响应强度或其他的一些因素进行均一化，从而利用颜色条的变化来直观地表示不同样本之间的含量变化情况的图。

定义很简单，这里我们的因子响应强度就是每个基因的raw count值，但是raw count值从0到几千上万差别非常之大，作图不方便。所以我们通常会用均一化的方法，使每个基因的raw count值变化程度处于同一个数量级，再通过不同颜色变化得到基因在不同样品的含量变化。

R自带的均一化函数是scale()，注意下scale默认的均一化方式是按列进行的，我们还可以通过函数 t() 进行矩阵的行列转化，只需要将差异基因挑出来按行（也就是基因名）进行均一化，导入pheatmap包即可做成一个最简单的热图。

2. 简化版代码

先上一个最最简易版的，比如我要分析前25个最可能发生差异表达的基因，代码如下：

library("pheatmap")
nDEGs <- gene[which(gene$group != "NOT_CHANGE"),]   # 筛选差异基因
sort_DEGs <- arrange(nDEGs,padj)        # 按照padj值升序排序
choose_gene <- head(sort_DEGs[,1],25)       # 取padj值最小的前25个基因
choose_matrix <- mycounts[choose_gene,]     # 从raw count矩阵中挑出这25个基因数据
heat_matrix <- t(scale(t(choose_matrix)))   # 转换了两次行列并均一化，实际就是按row进行了均一化
pheatmap(heat_matrix)   # 以默认参数做热图

在plots窗口可以预览生成的热图：

因为没有加任何参数调整，所以不好看（已经比R自带的热图函数做出来的好看了），先解释一下上面代码实现的原理。gene是我们前面做火山图的矩阵，里面已经有了我们差异基因分组的一列group

nDEGs <- gene[which(gene$group != "NOT_CHANGE"),] 这个代码是将group列中字符不等于“NOT_CHANGE”的数据挑出来赋值给nDEGs，注意下赋值后的nDEGs也是矩阵，可以直接查看。

sort_DEGs <- arrange(nDEGs,padj) arrange() 函数的功能是升序排列，这里按照padj值升序排列。

choose_gene <- head(sort_DEGs[,1],25) head()函数用法不说了，取了前25个基因。注意下我们取的第一列是基因名，如果你前面已经将基因名作为rownames导入了，那就要用rowname。

choose_matrix <- mycounts[choose_gene,]，返回到我们前面的raw count矩阵，将基因名对应的数据挑出来，可以看下这个时候的choose_matrix矩阵是怎么样的：

heat_matrix <- t(scale(t(choose_matrix))) 先进行一次行列转换，对列数据进行均一化，再进行一次行列转换，说白了就是对每行基因的raw count数据进行均一化，得到如下矩阵：

pheatmap(heat_matrix)以默认参数做热图，大功告成。

如果要对所有差异表达的基因做热图，只需要修改一下输入的矩阵就行：

1
2
3

all_matrix <- mycounts[(sort_DEGs[,1]),]
heat_matrix_all <- t(scale(t(all_matrix)))
pheatmap(heat_matrix_all)

因为这是默认参数作图，所以输出结果非常感人：

你论文敢用这种图？所以还是需要了解一下pheatmap包的各种参数，对热图进行调整和修改。

3. 参数详解

此部分内容参考CSDN博客：跳动的喵尾巴

pheatmap(mat, color = colorRampPalette(rev(brewer.pal(n = 7, name =
  "RdYlBu")))(100), kmeans_k = NA, breaks = NA, border_color = "grey60",
  cellwidth = NA, cellheight = NA, scale = "none", cluster_rows = TRUE,
  cluster_cols = TRUE, clustering_distance_rows = "euclidean",
  clustering_distance_cols = "euclidean", clustering_method = "complete",
  clustering_callback = identity2, cutree_rows = NA, cutree_cols = NA,
  treeheight_row = ifelse((class(cluster_rows) == "hclust") || cluster_rows,
  50, 0), treeheight_col = ifelse((class(cluster_cols) == "hclust") ||
  cluster_cols, 50, 0), legend = TRUE, legend_breaks = NA,
  legend_labels = NA, annotation_row = NA, annotation_col = NA,
  annotation = NA, annotation_colors = NA, annotation_legend = TRUE,
  annotation_names_row = TRUE, annotation_names_col = TRUE,
  drop_levels = TRUE, show_rownames = T, show_colnames = T, main = NA,
  fontsize = 10, fontsize_row = fontsize, fontsize_col = fontsize,
  angle_col = c("270", "0", "45", "90", "315"), display_numbers = F,
  number_format = "%.2f", number_color = "grey30", fontsize_number = 0.8
  * fontsize, gaps_row = NULL, gaps_col = NULL, labels_row = NULL,
  labels_col = NULL, filename = NA, width = NA, height = NA,
  silent = FALSE, na_col = "#DDDDDD", ...)

参数内容非常之多，我这里仅挑选一些可能用得上的做个记录：

参数	描述
color	表示热图颜色，colorRampPalette(rev(brewer.pal(n = 7, name = “RdYlBu”)))(100)表示颜色渐变调色板，“n” 的数量取决于调色板中颜色的数量，“name” 为调色板的名称，(100)表示100个等级；color = colorRampPalette(c(“blue”, “white”, “red”))(100)则是通过设置三种不同的颜色进行渐变显示
scale	表示进行均一化的方向，值为 “row”, “column” 或者”none”
kmeans_k	默认为NA，即不会对行进行聚类；如果想在进行层次聚类之前，先对行特征(因子)进行 k-means 聚类，则可在此调整热图的行聚类数
cluster_rows	表示仅对行聚类，值为TRUE或FALSE
cluster_cols	表示仅对列聚类，值为TRUE或FALSE
clustering_distance_cols	表示列聚类使用的度量方法，与行聚类的度量方法一致
clustering_method	表示聚类方法，包括：‘ward’, ‘ward.D’, ‘ward.D2’, ‘single’, ‘complete’, ‘average’, ‘mcquitty’, ‘median’, ‘centroid’
cutree_rows	若进行了行聚类，根据行聚类数量分隔热图行
cutree_cols	若进行了列聚类，根据列聚类数量分隔热图列
treeheight_row	若进行了行聚类，其热图行的聚类树高度，默认为 “50”
treeheight_col	若进行了列聚类，其热图列的聚类树高度，默认为 “50”
breaks	用来定义数值和颜色的对应关系，默认为 “NA”
border_color	表示热图每个小的单元格边框的颜色，默认为 “NA”
cellwidth	表示单个单元格的宽度，默认为 “NA”，即根据窗口自动调整
cellheight	表示单个单元格的高度，默认为 “NA”，即根据窗口自动调整
fontsize	表示热图中字体大小
fontsize_row	表示行名字体大小，默认与fontsize一致
fontsize_col	表示列名字体大小，默认与fontsize一致
fontsize_number	表示热图上显示数字的字体大小
angle_col	表示列标签的角度，可选择 “0”，“45”，“90”，“270”，“315”
display_numbers	表示是否在单元格上显示原始数值或按照特殊条件进行区分标记
number_format	表示热图单元格上显示的数据格式，如 “%.2f” 表示两位小数； “%.1e” 表示科学计数法
number_color	表示热图单元格上显示的数据字体颜色
legend	表示是否显示图例，值为TRUE或FALSE
annotation_row	表示是否对行进行注释
annotation_col	表示是否对列进行注释
annotation_colors	表示行注释及列注释的颜色
annotation_legend	表示是否显示注释的图例信息
annotation_names_row	表示是否显示行注释的名称
annotation_names_col	表示是否显示列注释的名称
show_rownames	表示是否显示行名
show_colnames	表示是否显示列名
main	表示热图的标题名字
gaps_row	仅在未进行行聚类时使用，表示在行方向上热图的隔断位置，如 gaps_row = c(２, 4)表示在第２与第4列进行隔断
gaps_col	仅在未进行列聚类时使用，表示在列方向上热图的隔断位置，同 gaps_row
labels_row	表示使用行标签代替行名
labels_col	表示使用列标签代替列名
filename	表示保存图片的位置及命名
width	表示输出绘制热图的宽度
height	表示输出绘制热图的高度
margins	表示热图距画布的空白距离

好吧，看上去基本上都能用到。制作热图应用的统计学原理就不多说了，我也没研究明白，我们用上面的这些参数能做出好看点的图就行。所以其实pheatmap中也有对应的参数scale来对行或列进行均一化，在pheatmap中设置参数或者多加一行代码做个均一化转换都是一样的。

4. 最终流程代码

知道了简化版代码的各行命令和pheatmap各参数作用，稍加一点点修改得到最终流程代码：

library("pheatmap")
nDEGs <- gene[which(gene$group != "NOT_CHANGE"),]
sort_DEGs <- arrange(nDEGs,padj)
all_matrix <- mycounts[(sort_DEGs[,1]),]
heat_matrix_all <- t(scale(t(all_matrix)))
# 要做热图的每列分组，底下两行代码是必须的
annotation_col <- data.frame(sample = factor(c(rep("SD",4),rep("LD1",4)))) # 简单粗暴地写一个分组矩阵
row.names(annotation_col) <- colnames(heat_matrix_all)
pheatmap(heat_matrix_all,
         color = colorRampPalette(c("blue", "white", "red"))(100),
         treeheight_col = 30,
         treeheight_row = 10,
         show_rownames = FALSE,
         angle_col = 45,
         annotation_col = annotation_col,
         filename = "test.pdf")