如何改造inferCNV的热图？

2020-01-03 本文已影响0人刘小泽

刘小泽写于2020.1.3
关于自己改造函数，之前介绍过如何改造Seurat包的DoHeatmap函数？

这次问题又来了：

首先花半小时来看什么是inferCNV：使用inferCNV分析单细胞转录组中拷贝数变异

划重点：

安装说明：https://github.com/broadinstitute/inferCNV/wiki/Installing-infercnv
需要依赖jags 程序
- mac安装：brew install jags
- windows安装：https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/JAGS-4.3.0.exe
- linux安装：直接conda
函数运行需要准备三个文件：
- 行为基因、列为样本的表达矩阵(matrix)
- 基因位置信息文件：第一列对应第一个文件的行名，其余三列则是基因的位置；并且表达矩阵的行名，要和基因坐标文件的基因名顺序一致
- 样本注释信息文件：基因名不能有重复，并且第二列是细胞的分组
准备好文件后可以运行”两步走“程序：CreateInfercnvObject和infercnv::run
其中在 infercnv::run中有一个参数：cutoff，官方建议10X数据设置为0.1，smart-seq数据设置为1
运行的结果会产生十多个数据和图片，注意区分！

然后找到inferCNV的官方说明文档：https://github.com/broadinstitute/inferCNV/wiki

它的结果会产生这样一张热图：

关于图片的解读在：

https://github.com/broadinstitute/inferCNV/wiki/Interpreting-the-figure

大体的意思是：

正常细胞的表达量在上半部分，肿瘤细胞的表达量在下半部分；
横坐标是各个基因在染色体的分布（染色体经过了排序），纵坐标是各个样本；
表达量差异来源于肿瘤细胞的表达量减去正常细胞的表达量，红色部分表示这部分染色体区域的基因发生扩增，蓝色表示缺失

这个图刚拿到手，的确看起来十分复杂，不知从何下手
但仔细结合官网的解释，又发现这个热图和我们平常画的图有很相似

怎么做到知识迁移？

改造图的一个好处就是锻炼自己已掌握知识的迁移能力

思考的过程

看到这个图，我首先想到如何拿到表达矩阵，然后既然是一个大热图包含了两个小热图，那么如何分别作出两个小图？

我把这整个过程都运行了一遍，速度很快（记得先在infercnv::run这一步中设置HMM=F，可以提高计算速度），然后发现得到了一堆文件

那么哪些是做热图需要的呢？

这就要去官方找答案了：官方的解答很友好，把所有可能问到的问题尽可能详细地罗列清楚：https://github.com/broadinstitute/inferCNV/wiki/Output-Files

这里就能得到一个惊喜：

有了这些储备，就可以开始探索了

step0：首先拿到必备数据

# 提取：处理后的表达矩阵
expr <- infercnv_obj@expr.data
dim(expr)
expr[1:4,1:4]
# 提取：表达矩阵样本中正常细胞的位置【这个normal_loc是一个列表】
(normal_loc <- infercnv_obj@reference_grouped_cell_indices)
normal_loc <- normal_loc$WT
# 提取：表达矩阵样本中肿瘤细胞的位置
(tumor_loc <- infercnv_obj@observation_grouped_cell_indices)
tumor_loc <- tumor_loc$Tumor

Step1: 将基因名与染色体位置对应

当看到示例图中的分隔线，我就想到了ComplexHeatmap也有这个功能。但示例图中的横坐标是染色体，而我们的表达矩阵是基因和样本名，怎么对应到染色体呢？

想到了之前准备的三个文件之一：基因位置信息文件，而且它存储的也是染色体排序后的信息

gn <- rownames(expr)
geneFile <- read.table('infercnv-heatmap/geneFile.txt')
> length(geneFile$V1);length(gn) 
[1] 19580
[1] 13337

看到我们这里的表达矩阵由于在计算过程中进行了某些处理，所以相比最初的基因数量少了一些。那么就对geneFile取子集：

sub_geneFile <-  geneFile[geneFile$V1%in%gn,]
> dim(sub_geneFile)
[1] 13337     4
# 最后可以再检查一下它们是否一致
identical(rownames(expr),sub_geneFile$V1) # TRUE

Step2: 拆分矩阵

上面说过，大的热图整体分成两部分：上面的热图是正常细胞，下面是肿瘤细胞，并且我们知道了各自的位置，就能先把各自的小表达矩阵提取出来

位置信息就存储在之前的normal_loc, tumor_loc中

norm_expr <- expr[,normal_loc]
norm_expr$chr <- as.factor(sub_geneFile$V2) #最后加上一列：对应的chr信息【保存为因子型，一会作图会用到】

> table(norm_expr$chr) #这个信息就与横坐标的间隔对应，chr间隔越大表示其中包含的基因越多
 chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19  chr2  chr3  chr4 
  830   650  1133   467   493   460   534   430   685   370   472  1113   694   884 
 chr5  chr6  chr7  chr8  chr9 
  885   693  1035   716   793 

# 同理对tumor
tumor_expr <- expr[,tumor_loc]
tumor_expr$chr <- as.factor(sub_geneFile$V2)

Step3-1: 开始画图=》第一次尝试【先画tumor的】

关于配色部分：一开始在纠结到底应该用什么颜色，自己不想去找。想到：作者是怎么定义颜色的呢？于是又继续延续之前在：狸猫换太子-"掉包"某个包装好的绘图函数中的思想，查找了作图函数 infercnv::plot_cnv

library(ComplexHeatmap)
# 调整配色
library("RColorBrewer")
# 来自函数：infercnv::plot_cnv
get_group_color_palette <- function () {
  return(colorRampPalette(RColorBrewer::brewer.pal(12, "Set3")))
}
color <- get_group_color_palette()(length(unique(tumor_expr$chr)))
# 列标题的颜色框
top_color <- HeatmapAnnotation(
  cluster = anno_block(gp = gpar(fill = color), # 设置填充色
                       labels = levels(new_cluster), 
                       labels_gp = gpar(cex = 0.9, col = "black")))

开始画图，这个画图代码是不断从complex heatmap的官方说明(https://jokergoo.github.io/ComplexHeatmap-reference/book/a-single-heatmap.html)中拷贝、粘贴、调试，得到的：

pdf("test1-heatmap.pdf",width = 15,height = 10)
n <- t(tumor_expr[,-ncol(tumor_expr)])
ht_tumor = Heatmap(as.matrix(n),
                  cluster_rows = T,
                  cluster_columns = F,
                  show_column_names = F,
                  show_row_names = F,
                  column_split = new_cluster,
                  heatmap_legend_param = list(
                    title = "Modified Expression",
                    title_position = "leftcenter-rot", # 图例标题位置
                    at=c(0.4,1.6), #图例范围
                    legend_height = unit(3, "cm") #图例长度
                  ),
                  top_annotation = top_color,
                  row_title = "Observations (Cells)",
                  row_title_side = c("right"),
                  column_title = "Genomic Region",
                  column_title_side = c("bottom"))
draw(ht_tumor, heatmap_legend_side = "left") # 图例位置

dev.off()

Step3-2: 画图=》第二次尝试【组合normal和tumor】

m <- t(norm_expr[,-ncol(norm_expr)])

pdf("test2-heatmap.pdf",width = 25,height = 30)
ht_normal = Heatmap(as.matrix(m),
                   cluster_rows = T,
                   cluster_columns = F,
                   show_column_names = F,
                   show_row_names = F,
                   column_split = new_cluster,
                   row_title = "References (Cells)",
                   row_title_side = c("right"),
                   row_title_rot = 90,
                   row_title_gp = gpar(fontsize = 25),
                   column_title = NULL, 
                   heatmap_legend_param = list(
                     title = "Modified Expression",
                     title_position = "leftcenter-rot", # 图例标题位置
                     title_gp = gpar(fontsize = 20),# 图例标题大小
                     at=c(0.4,1.6), #图例范围
                     legend_height = unit(6, "cm")),#图例长度
                   width = 20, height = 5
                   ) 

ht_tumor = Heatmap(as.matrix(n),
                   cluster_rows = T,
                   cluster_columns = F,
                   show_column_names = F,
                   show_row_names = F,
                   column_split = new_cluster,
                   show_heatmap_legend=F,
                   top_annotation = top_color,
                   row_title = "Observations (Cells)",
                   row_title_side = c("right"),
                   row_title_rot = 90,
                   row_title_gp = gpar(fontsize = 25),
                   column_title = "Genomic Region",
                   column_title_side = c("bottom"),
                   column_title_gp = gpar(fontsize = 25),
                   width = 20, height = 10,
                   heatmap_height = 15)

# 设置竖直排列
draw(ht_normal %v% ht_tumor)

dev.off()

写在最后

上面👆的图是用的preliminary.infercnv_obj做的，如果要做出其他的热图，只要能拿到其他的obj结果，也是可以做的

欢迎关注我们的公众号～_～　　
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!