单细胞分析实录(11): inferCNV的基本用法

2021-04-04 本文已影响0人 TOP生物信息

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的内容，敬请期待~

因水平有限，有错误的地方，欢迎批评指正！