自己动手，丰衣足食~改造inferCNV的热图
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
豆豆写于2020.1.3
关于自己改造函数，之前介绍过
如何改造Seurat包的DoHeatmap函数？
(https://www.jianshu.com/p/b2816a72e79e)
这次问题又来了：
首先花半小时来看什么是inferCNV：使用inferCNV分析单细胞转录组中拷贝数变异
划重点：
安装说明：https://github.com/broadinstitute/inferCNV/wiki/Installing-infercnv
需要依赖
jags程序mac安装：
brew install jagswindows安装：
https://sourceforge.net/projects/mcmc-jags/files/JAGS/4.x/Windows/JAGS-4.3.0.exelinux安装：直接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结果，也是可以做的
点击底部的“阅读原文”，获得更好的阅读体验哦😻
初学生信，很荣幸带你迈出第一步
🤓生信星球 🌎~ 一个不拽术语、通俗易懂的生信知识平台