InferCNV可以用于肿瘤单细胞RNA-Seq数据中鉴定大规模染色体拷贝数变异，例如整个染色体或大片段染色体的扩增或缺失。基本思路是在整个基因组范围内，将每个肿瘤细胞基因表达与平均表达或“正常”参考细胞基因表达对比，确定其表达强度。
这段话来自InferCNV官方文档：https://github.com/broadinstitute/inferCNV/wiki
实际分析中，我们经常用inferCNV来判断肿瘤细胞。当然还可以做肿瘤异质性、克隆进化方面的探索，这部分我会在下一期的推文中介绍，本期推文介绍基本使用。公众号后台回复20210404获取今天的代码和测试数据。
1. 安装
初次安装加载时，可能会提示你安装JAGS，根据提示安装JAGS后，再加载InferCNV就没问题
BiocManager::install("infercnv")
library("infercnv")
错误: package or namespace load failed for ‘infercnv’:
loadNamespace()里算'rjags'时.onLoad失败了，详细内容：
调用: fun(libname, pkgname)
错误: Failed to locate any version of JAGS version 4
The rjags package is just an interface to the JAGS library
Make sure you have installed JAGS-4.x.y.exe (for any x >=0, y>=0) from
http://www.sourceforge.net/projects/mcmc-jags/files
2. 输入
在进行分析之前，需要准备三个文件：
- Raw Counts Matrix for Genes x Cells
- 细胞的注释文件
共两列，一列CB，一列细胞来源，tab分割，无列名；
如果肿瘤细胞的注释类似malignant_{patient}，则在后续聚类画图时，设置cluster_by_groups=T会根据patient的来源进行区分
该文件的CB可以小于count矩阵，此时inferCNV只会对这部分细胞进行分析； - 基因排序文件
tab分割，无列名
只有该文件和count矩阵共有的基因才会被分析（可以去掉不想画图的基因，比如线粒体基因），且该文件的染色体顺序决定了最终热图的染色体顺序（有些文章的图，inferCNV热图的染色体顺序是颠倒的，问题就出在这儿）。
我一般是根据cellranger提供的参考基因组gtf注释文件，得到基因排序文件
格式如下
WASH7P chr1 14363 29806
LINC00115 chr1 761586 762902
NOC2L chr1 879584 894689
MIR200A chr1 1103243 1103332
3. 流程
看一下这张流程图
我比较关注最后两步，去噪和CNV预测是分开的，一些已发表的文章这部分的分析都是基于去噪之后的结果（比如那张热图），CNV预测结果用的比较少，个人觉得使用analysis_mode = "subclusters"模式后的热图更好看，而且只有6个值
这是官网的一张对比图
4. 使用
主要步骤只有两步，如下
library(infercnv)
#使用inferCNV之前，最好过滤掉低质量的细胞
infercnv_obj = CreateInfercnvObject(raw_counts_matrix="oligodendroglioma_expression_downsampled.counts.matrix",
annotations_file="oligodendroglioma_annotations_downsampled.txt",
delim="\t",
gene_order_file="gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt",
ref_group_names=c("Microglia/Macrophage","Oligodendrocytes (non-malignant)")
)
#ref_group_names参数根据细胞注释文件填写，在示例中，这两种细胞是非恶性细胞，所以作为参照；
#ref_group_names=NULL，则会选用所有细胞的平均表达作为参照，这种模式下，最好能确保细胞内在有足够的差异
infercnv_obj = infercnv::run(infercnv_obj,
cutoff=1,
out_dir="try",
cluster_by_groups=TRUE,
#analysis_mode="subclusters", #默认是"samples"
denoise=TRUE,
HMM=TRUE,
num_threads=4
)
#去噪，HMM预测CNV这两项我一般都选上
关于cutoff参数，官网是这样说的：cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics，表示基因在所有参照细胞中，表达count的平均值的最小阈值，10X数据更稀疏，所以这个值小
cluster_by_groups：先区分细胞来源，再做层次聚类
out_dir表示输出文件夹的名字，没有会自动创建
5. 输出
我认为重要的输出文件
infercnv.png
降噪之后生成的最终热图，图中的数值是"Residual expression values"，可以简单理解为基因表达值的另一种形式；
infercnv.references.txt : the 'normal' cell matrix data values
对应热图的参照panel
infercnv.observations.txt : the tumor cell matrix data values
对应热图的观测panel
infercnv.observation_groupings.txt
group memberships for the tumor cells as clustered，对应热图的观测panel的两个bar
"Dendrogram Group" "Dendrogram Color" "Annotation Group" "Annotation Color"
"MGH36_P9_B01" "1" "#8DD3C7" "1" "#8DD3C7"
"MGH36_P3_E06" "1" "#8DD3C7" "1" "#8DD3C7"
#2 3列只有一类值
#4 5列有四类值
infercnv.observations_dendrogram.txt
树状图的newick格式，对应热图的观测panel的“树”
run.final.infercnv_obj这个rds文件，以及包含在其中的三个子文件
如下
@expr.data：对应最终热图的表达文件
@reference_grouped_cell_indices：对应最终热图的reference细胞名称
@observation_grouped_cell_indices：对应最终热图的observation细胞名称
HMM_CNV_predictions.*.pred_cnv_regions.dat:
第二列是CNV的name，唯一；
第一列是CNV所属的group，inferCNV默认的模式(analysis_mode = "samples")是将一个patient作为一个整体找CNV，所以示例会有四个group；
4 5 6列包含CNV的坐标；
第三列表示状态：
State 1: 0x: complete loss
State 2: 0.5x: loss of one copy
State 3: 1x: neutral
State 4: 1.5x: addition of one copy
State 5: 2x: addition of two copies
State 6: 3x: essentially a placeholder for >2x copies but modeled as 3x
$ less -S HMM_CNV_predictions.HMMi6.hmm_mode-samples.Pnorm_0.5.pred_cnv_regions.dat | head -n 3
cell_group_name cnv_name state chr start end
malignant_MGH36.malignant_MGH36_s1 chr1-region_2 2 chr1 6281253 144341756
malignant_MGH36.malignant_MGH36_s1 chr1-region_4 3 chr1 151336778 156213123
HMM_CNV_predictions.*.cell_groupings
tumor subclusters和cell的对应关系
只有analysis_mode = "subclusters"模式下才会生成，这个模式挺慢，但肿瘤异质性和克隆进化都是在这种模式下做的
cell_group_name cell
malignant_MGH36.malignant_MGH36_s1 MGH36_P3_E06
malignant_MGH36.malignant_MGH36_s1 MGH36_P10_E07
HMM_CNV_predictions.*.pred_cnv_genes.dat
每一个group（第一列）, 每一个CNV片段（第二列）上面每一个基因（第四列）的CNV状态（第三列），文件中基因这一列是唯一的
$ less -S HMM_CNV_predictions.*.pred_cnv_genes.dat | head -n 3
cell_group_name gene_region_name state gene chr start end
malignant_MGH36.malignant_MGH36_s1 chr1-region_2 2 ACOT7 chr1 6281253 6296032
malignant_MGH36.malignant_MGH36_s1 chr1-region_2 2 NOL9 chr1 6324329 6454451
该文件可以查询感兴趣基因的CNV状态，如下：
$ grep EGFR HMM_CNV_predictions.HMMi6.hmm_mode-samples.Pnorm_0.5.pred_cnv_genes.dat | head -n 3
malignant_MGH36.malignant_MGH36_s1 chr7-region_16 5 EGFR chr7 54819943 54827667
malignant_93.malignant_93_s1 chr7-region_98 4 EGFR chr7 54819943 54827667
malignant_97.malignant_97_s1 chr7-region_142 4 EGFR chr7 54819943 54827667
这种方法看基因的CNV状态，也有局限性，只有当这个基因与它所处的CNV region的CNV状态比较一致才会被找出来，如果这个基因的CNV状态比周围一些基因的CNV状态强很多，则找不出来。
个人猜测是因为inferCNV是画窗口针对相邻的多个基因做推断，如果只有一个基因表达很高或很低，则被平均掉
上述的这些输出文件，对于深入分析CNV很有用，但是如果只区分肿瘤细胞，则用不到这么多文件。
6. 画图
infercnv包也包含了画图函数plot_cnv，这个知道的人不多，
library(RColorBrewer)
infercnv::plot_cnv(infercnv_obj, #上两步得到的infercnv对象
plot_chr_scale = T, #画染色体全长，默认只画出（分析用到的）基因
output_filename = "better_plot",output_format = "pdf", #保存为pdf文件
custom_color_pal = color.palette(c("#8DD3C7","white","#BC80BD"), c(2, 2))) #改颜色
后续应该还有两篇inferCNV的内容，敬请期待~
因水平有限，有错误的地方，欢迎批评指正！
网友评论