探索单细胞测序中的恶性肿瘤细胞，一定不能少了它——inferCNV

inferCNV 是一款 R 包，用来探索肿瘤单细胞转录组数据中体细胞大规模染色体拷贝数变异，例如整个染色体或大片段染色体的扩增（gain）或缺失（loss）。

inferCNV 工作原理

通过与一组参考的“正常”细胞比较，探索肿瘤基因组不同位置的基因表达强度来实现。通过热图显示每条染色体上的相对表达强度，在以正常细胞作为参考的情况下可以直观地判断肿瘤基因组上的哪些区域是过表达和低表达。为了提高推断结果的可靠性，inferCNV 采用几种 residual expression 过滤方法减少数据中噪声信号，从而发现真正的 CNV 的信号。

inferCNV 必备文件

创建 InferCNV 对象有 3 个必要文件：

1、单细胞 RNA-seq 表达量的原始矩阵；

2、细胞类型注释文件；

3、基因排序文件。

① 单细胞转录组表达量的原始矩阵

原始表达矩阵数据如下图所示，行为基因名，列为细胞名，用制表符分隔。

参考获取代码：

# seurat 3.0+counts_matrix = Seurat::GetAssayData(seurat_obj, slot="counts")

② 细胞类型注释文件

细胞类型注释文件格式只是两列，以制表符分隔，没有列标题，第一列是细胞名，第二列是定义的细胞类型名。这些细胞类型信息将被聚类并在结果热图中展示。另外，inferCNV 将只分析示例注释文件中列出的那些细胞。如果感兴趣的细胞是 total counts 矩阵的子集，则不需要创建包含感兴趣的子集的新矩阵。

参考获取代码：

cellanno = FetchData(seurat_obj, vars = "celltype" ) %>% tibble::rownames_to_column(var = "cellbarcode")write.table(cellanno, "cnv_celltype_group.xls", sep = "\\t", col.names = F,row.names = F, quote = F)

③ 基因排序文件

基因排序位置信息文件提供每个基因的染色体位置（需要按顺序）。文件以制表符分隔，没有列标题，仅提供基因名称、染色体和基因跨度。分析的计数矩阵中的每个基因都应该在这个文件中提供相应的基因名称和位置信息。只有同时存在于计数矩阵和基因排序文件中的基因才会包含在 inferCNV 分析中。另外，基因名称是唯一的。

参考获取代码：

gtf = "./refdata-gex-GRCh38-2020-A/genes/genes.gtf"gtf = plyranges::read_gff(gtf)gene.chr = gtf %>% plyranges::filter(type == "gene" & gene_name %in% rownames(seurat_obj)) %>%    as.data.frame() %>%    dplyr::select(gene_name, seqnames, start, end) %>%    dplyr::distinct(gene_name, .keep_all=T) write.table(gene.chr, "gene_order_file.xls", col.names =F, row.names =F, sep = "\\t", quote =F )

inferCNV 基本用法

第一步：创建 inferCNV 对象

# create the infercnv objectinfercnv_obj = CreateInfercnvObject(raw_counts_matrix=counts_matrix,  # 可以直接提供矩阵对象                                    annotations_file="cellAnnotations.txt",                                    delim="\\t",                                    gene_order_file="gene_ordering_file.txt",                                    ref_group_names=c("normal"))  # 用于设置参考组，正常的细胞类型组别

其中，ref_group_name 参数用于设置参考组。如果不确定哪些组是正常细胞，那么设置为 ref_group_name=NULL, inferCNV 会以所有细胞的平均信号值作为基线，这适用于当所包含的细胞之间存在显著差异时。此外，需要注意的是 raw_count_matrix 输入的是排除了低质量细胞的计数矩阵。

第二步，运行 inferCNV 流程

# perform infercnv operations to reveal cnv signalinfercnv_obj = infercnv::run(infercnv_obj,                                    cutoff=0.1,  # use 1 for smart-seq, 0.1 for 10x-genomics                                    out_dir="output_dir",  # dir is auto-created for storing outputs                                    analysis_mode="subclusters" # 更好方法，按照 cnv 亚群进行预测，而不是样本。能更好的区分肿瘤异质性。                                    cluster_by_groups=T,  # 根据细胞类型对肿瘤细胞执行单独的聚类，如 cell annotations 文件中定义的那样                                    denoise=T,  # 去噪处理                                    num_threads=10, # 设置线程数, 多线程运行，加快计算速度                                    HMM=T) # 是否基于 HMM 预测 CNV,选择 F 会加快运算速度

需要注意的是：cutoff 为每个基因在参考细胞中的最小平均 read 计数，该阈值决定了哪些基因将用于 inferCNV 分析，低于细胞平均计数的基因将被排除在外。对于 smart-seq（全长转录组测序），推荐 cutoff 为 1 ；对于 10x（以及其他可能的 3'' 端测序和液滴分析，计数矩阵往往更稀疏），推荐 cutoff 为 0.1。此外，设置 HMM=T 的计算时间会长于 HMM=F，一些已发表的文章这部分的分析都是基于去噪之后的结果，即 De-noising 热图。这一步一般比较耗时，细胞数越多，运行速度越慢。

workflow 示意图：

inferCNV 结果解读

inferCNV 结果很多，比较重要的文件如下：

infercnv.preliminary.png：初步推断 CNV 的热图 (在去噪或 HMM 预测之前生成)。

infercnv.png：利用去噪方法推断 CNV 生成的热图。

infercnv.references.txt：正常细胞的表达矩阵信息，对应 infercnv.png 顶部热图。

infercnv.observations.txt：观测细胞表达矩阵信息，对应 infercnv.png 底部热图。

infercnv.observation_groupings.txt：观测细胞表达矩阵信息分组对应的数据，对应 infercnv.png 底部热图左侧两个 bar。

infercnv.observations_dendrogram.txt：树状图的 newick 格式，对应 infercnv.png 底部热图的“树”。

run.final.infercnv_obj：最终生成的 inferCNV 对象，后续可根据该文件对 CNV 结果进行进一步处理和个性化绘图。

infercnv_obj@expr.data: 经过去噪或者 HMM 处理的表达矩阵信息。

infercnv_obj@reference_grouped_cell_indices: 包含与每个正常细胞类型对应的表达矩阵列索引列表。

infercnv_obj@observation_grouped_cell_indices: 包含与每个肿瘤细胞类型对应的表达矩阵列索引列表。

图解：正常细胞的表达值绘制在顶部热图上，肿瘤细胞绘制在底部热图上，基因在染色体上从左到右排列。从肿瘤细胞表达数据中有效地减去正常细胞表达数据以产生差异，其中染色体区域扩增显示为红色块，而染色体区域缺失显示为蓝色块。

Residual expression values（剩余表达值）图中的数值是热图的颜色强度对应于在执行一系列数据转换并且从肿瘤细胞表达数据中有效地减去正常细胞表达数据之后的残余表达值，可以简单理解为基因表达值的另一种形式。

CNV 预测

基于 HMM（隐马尔可夫模型）的 CNV 预测的两种模型，分别称为 i3 和 i6 模型，流程中 默认采用 i6 模型。六态 CNV 模型，预测以下 CNV 水平:

0x: 完全缺失， state 为 1；

0.5x: 缺失一个拷贝，state 为 2； 

1x: 正常，state 为 3； 

1.5x: 增加一个拷贝，state 为 4； 

2x: 增加两个拷贝，state 为 5； 

3x: 增加 3 个或者 3 个以上拷贝， state 为 6。

CNV region prediction reports

①HMM_CNV_predictions.*.pred_cnv_regions.dat：包含 CNV 区域坐标，state 分配和细胞分组的摘要。

less -S HMM_CNV_predictions.HMMi6.rand_trees.hmm_mode-subclusters.Pnorm_0.5.pred_cnv_regions.dat# cell_group_name    cnv_name    state   chr start   end# malignant_MGH36.malignant_MGH36.1.1.1.1    chr1-region_1   2   chr1    14363   145139569# malignant_MGH36.malignant_MGH36.1.1.1.1    chr1-region_3   4   chr1    151132224   156470620# malignant_MGH36.malignant_MGH36.1.1.1.1    chr1-region_5   2   chr1    202496438   220263804# malignant_MGH36.malignant_MGH36.1.1.1.1    chr2-region_8   4   chr2    233562009   242708231# malignant_MGH36.malignant_MGH36.1.1.1.1    chr3-region_9   4   chr3    3168600 12913415# ...

②HMM_CNV_predictions.*.pred_cnv_genes.dat：包含每个细胞、每个基因 CNV state level。

less -S HMM_CNV_predictions.HMMi6.rand_trees.hmm_mode-subclusters.Pnorm_0.5.pred_cnv_genes.dat# cell_group_name    gene_region_name    state   gene    chr start   end# malignant_MGH36.malignant_MGH36.1.1.1.1    chr1-region_1   2   WASH7P  chr1    14363   29806# malignant_MGH36.malignant_MGH36.1.1.1.1    chr1-region_1   2   LINC00115   chr1    14363   29806# malignant_MGH36.malignant_MGH36.1.1.1.1    chr1-region_1   2   NOC2L   chr1    761586  762902# malignant_MGH36.malignant_MGH36.1.1.1.1    chr1-region_1   2   SDF4    chr1    879584  894689# malignant_MGH36.malignant_MGH36.1.1.1.1    chr1-region_1   2   UBE2J2  chr1    1152288 1167411# ...

如果有感兴趣的基因，可从该文件找到该基因是否有明显的 CNV 扩增和缺失。例如，该样本对应于少突胶质细胞瘤(脑肿瘤类型)，在这种肿瘤类型中，EGFR 基因经常被发现扩增，并作为肿瘤发生的驱动因素或贡献者。

grep "EGFR" HMM_CNV_predictions.HMMi6.rand_trees.hmm_mode-subclusters.Pnorm_0.5.pred_cnv_genes.dat | head -n 2# malignant_MGH36.malignant_MGH36.1.1.1.1    chr7-region_21  5   EGFR    chr7    54819943    54827667# malignant_MGH36.malignant_MGH36.1.1.1.2    chr7-region_91  5   EGFR    chr7    54819943    54827667

Note：这种方法看基因的 CNV 状态具有局限性：只有假设感兴趣的基因的扩增或缺失水平与所处较大 CNV 区域的状态是一致的，该基因才能被找到。如果该基因的扩增或缺失状态比其他相邻基因更严重，可能找不出该基因的 CNV 状态，并且这种方法无法确定这一特征。因为 inferCNV 是基因窗口针对相邻的基因推断 CNV 状态（默认 window_length=101，即 101 个基因的窗口），如果只有 1 个基因表达过高或过低，则被平均掉。

③ HMM_CNV_predictions.*.genes_used.dat：列出了所有用于分析的基因，而不是所有的基因，因为这里只包括在推断 CNV 处理过程中在细胞中满足最低表达量并且最终保存在表达矩阵中的基因。

该部分文件对深入分析 CNV 很有用，如果只需要区分肿瘤细胞，则用不上这些文件。

结语

目前，利用 inferCNV 检测 CNV 的研究者比较多，今年 3 月份我司项目文章 《单细胞转录组测序分析子宫腺肌症患者在位内膜和异位病灶的转录特征》发表于期刊 Cell & Bioscience ，该文章应用 inferCNV 分析来研究拷贝数变异，从而证实恶性肿瘤细胞。本次对于 inferCNV 的介绍，就到此了，感谢您的阅览与学习。

参考文献

1. https://github.com/broadinstitute/infercnv/wiki （可

2. Puram S V, Tirosh I, Parikh A S, et al. Single-cell transcriptomic analysis of primary and metastatic tumor ecosystems in head and neck cancer[J]. Cell, 2017, 171(7): 1611-1624. e24.