inferCNV

inferCNV用与探索肿瘤单细胞RNA-seq数据，分析其中的体细胞大规模染色体拷贝数变化(copy number alterations, CNA), 例如整条染色体或大片段染色体的增加或丢失(gain or deletions)。工作原理是，以一组"正常"细胞作为参考，分析肿瘤基因组上各个位置的基因表达量强度变化. 通过热图的形式展示每条染色体上的基因相对表达量，相对于正常细胞，肿瘤基因组总会过表达或者低表达。

基本思路是在整个基因组范围内，通过计算肿瘤细胞关于参考的“正常”细胞的相对表达，利用滑窗思想对邻近的基因相对表达，计算拷贝数。

软件安装

尽管inferCNV是一个R包，但是在安装inferCNV之前还需要先下载安装JAGS ，好在它有Windows，MacOS和Linux版本，所以inferCNV在各个平台都能用。

Windows和MacOS的JAGS容易安装，而Linux的JAGS需要编译

# 手动安装BLAS和LAPACK不推荐
# yum install blas-devel lapack-devel
tar xf JAGS-4.3.0.tar.gz 
cd JAGS-4.3.0
./configure --libdir=/usr/local/lib64
make -j 20 && make install 

安装R包

install.packages("rjags")
if (!requireNamespace("BiocManager", quietly = TRUE))
     install.packages("BiocManager")
BiocManager::install("infercnv")

测试安装

library(infercnv)

infercnv_obj = CreateInfercnvObject(raw_counts_matrix=system.file("extdata", "oligodendroglioma_expression_downsampled.counts.matrix.gz", package = "infercnv"),
                                    annotations_file=system.file("extdata", "oligodendroglioma_annotations_downsampled.txt", package = "infercnv"),
                                    delim="\t",
                                    gene_order_file=system.file("extdata", "gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt", package = "infercnv"),
                                    ref_group_names=c("Microglia/Macrophage","Oligodendrocytes (non-malignant)")) 

infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics
                             out_dir=tempfile(), 
                             cluster_by_groups=TRUE, 
                             denoise=TRUE,
                             HMM=TRUE)

如果没有报错，就说明安装成功。

软件使用

准备输入文件

需要准备3个输入数据

1. 单细胞RNA-seq表达量的原始矩阵
2. 注释文件，记录肿瘤和正常细胞
3. 基因或染色体位置文件

第一个是Genes x Cells的表达矩阵(matrix)，行名是基因，列名是细胞编号。

MGH54_P16_F12	MGH54_P12_C10	MGH54_P11_C11	MGH54_P15_D06	MGH54_P16_A03	...
A2M	0	0	0	0	0	...
A4GALT	0	0	0	0	0	...
AAAS	0	37	30	21	0	...
AACS	0	0	0	0	2	...
AADAT	0	0	0	0	0	...
...	...	...	...	...	...	...

第二个是样本注释信息文件，命名为"cellAnnotations.txt"。一共两列，第一列是对应第一个文件的列名，第二列是细胞的分组

MGH54_P2_C12    Microglia/Macrophage
MGH36_P6_F03    Microglia/Macrophage
MGH54_P16_F12   Oligodendrocytes (non-malignant)
MGH54_P12_C10   Oligodendrocytes (non-malignant)
MGH36_P1_B02    malignant_MGH36
MGH36_P1_H10    malignant_MGH36

第三个是基因位置信息文件，命名为"geneOrderingFile.txt"。一共四列，第一列对应第一个文件的行名，其余三列则是基因的位置。注：基因名不能有重复,
该文件的染色体顺序决定了最终热图的染色体顺序；如果输出的图，染色体位置是颠倒的，一般是这个文件中，染色体位置排序有问题。

WASH7P  chr1    14363   29806
LINC00115       chr1    761586  762902
NOC2L   chr1    879584  894689
MIR200A chr1    1103243 1103332
SDF4    chr1    1152288 1167411
UBE2J2  chr1    1189289 1209265

两步法

最复杂的工作就是准备输入文件，而一旦上述三个文件已经创建完成，那么分析只要两步以及根据结果对参数进行调整。

第一步，根据上述的三个文件创建inferCNV对象

infercnv_obj = CreateInfercnvObject(raw_counts_matrix=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是排除了低质量细胞的count矩阵。

第二步，运行标准的inferCNV流程。

# perform infercnv operations to reveal cnv signal
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1,  # use 1 for smart-seq, 0.1 for 10x-genomics
                             out_dir="output_dir",  # 输出文件夹
                             cluster_by_groups=T,   # 聚类
                             denoise=T, #去噪
                             HMM=T) # 是否基于HMM预测CNV

关键参数是cutoff, 用于选择哪些基因会被用于分析（在所有细胞的平均表达量需要大于某个阈值）。这个需要根据具体的测序深度来算，官方教程建议10X设置为0.1，smart-seq设置为1。你可以先评估下不同阈值下的保留基因数，决定具体值。cluster_by_groups用于声明是否根据细胞注释文件的分组对肿瘤细胞进行分群。

最终会输出很多文件在out_dir的目录下，而实际有用的是下面几个

• infercnv.preliminary.png : 初步的inferCNV展示结果（未经去噪或HMM预测）
• infercnv.png : 最终inferCNV产生的去噪后的热图.
• infercnv.references.txt : 正常细胞矩阵.
• infercnv.observations.txt : 肿瘤细胞矩阵.
• infercnv.observation_groupings.txt : 肿瘤细胞聚类后的分组关系.
• infercnv.observations_dendrogram.txt : NEWICK格式，展示细胞间的层次关系.

所需要关注的参数：CreateInfercnvObject

细胞表达谱

@param raw_counts_matrix    
#基因位点参考文件：
@param gene_order_file  
包含基因组中每个染色体上每个基因的位置的数据文件。
#细胞注释文件
@param annotations_file 
细胞的描述，指示细胞的类型分类

#细胞注释文件中正常参考细胞的类型names
@param ref_group_names  

@param delim    
输入文件中使用的分隔符；默认为"\t"

@param max_cells_per_group  
每组要使用的最大单元格数。Default=NULL，使用在annotations_file中定义的所有单元格
随机测试代码的时候可以设置，如设置max_cells_per_group=50，以获得更快速的计算。

@param min_max_counts_per_cell  
参考默认即可，不需要额外设置

@param chr_exclude  
Default = c('chrX', 'chrY', 'chrM')
参考基因组注释中应排除在分析之外的染色体列表。

所需要关注的参数：infercnv::run

输入输出

@param infercnv_obj 
使用原始count数据填充的infercnv对象
@param cutoff   
对于Smart-seq最优cutoff=1对于10X数据，最优 cutoff=0.1        
@param min_cells_per_gene   
最低数量的参考细胞需要表达测量包括相应的基因。默认为3
@param out_dir  
输出文件的地址,必须提供

@param no_plot  
是否输出图像。默认为输出，相反，生成所有非图像输出作为运行的一部分(default: FALSE)

@param no_prelim_plot   
不输出初始的图像(default: FALSE)

@param output_format    
输出图像的格式。
选择："png", "pdf" and NA.  (default: "png")
NA意味着只写文本输出而不生成图形本身。

@param useRaster    
是否使用栅格化绘制热图。只有当它产生错误时才禁用它，因为它比不使用它快得多 (default: TRUE)

@param up_to_step   
运行步骤的数目；(default: 100 >> 23 steps currently in the process)

滑窗参数

@param window_length    
滑窗的长度，必须事奇数，默认为101（最好遵循默认）

@param smooth_method    
可以使用的滑窗方法: c(runmeans,pyramidinal,coordinates) 默认: pyramidinal

滑窗前表达谱的处理

@param max_centered_threshold   
**限定相对表达值分布范围**，确定肿瘤细胞基于参考细胞计算完相对表达值之后的最大值、最小值界限，(default: 3)。
 - 传入参数为“数字型”，则-1乘以该值，如果传入3，则范围被设定为[-3, 3]。
 - 传入参数为”auto” 根据细胞间的平均值分布来自动识别边界。
 - 如果为NA 则关闭该功能。
@param scale_data   
默认为FALSE，如果参考细胞和需要推测的细胞不是一种类型的细胞则需要scale

聚类设置

@param num_ref_groups   
是否对参考细胞进行分组的组别数目，默认为NULL，建议不设置组别数，在你具有确定的正常参考细胞的情况下

@param ref_subtract_use_mean_bounds 
分别确定每个ref组的平均值，然后移除平均值范围内的强度(默认值:TRUE)否则，使用组间平均值。

@param cluster_by_groups    
如果观察结果是根据组来定义的(例如。患者)，每一组细胞将单独聚在一起。(default=FALSE，改为使用k_obs_groups设置)

@param cluster_references   
是否在参考细胞中也引用聚类

@param k_obs_groups 
默认为1，根据拷贝数重新聚类的组数

@param hclust_method    
用于细胞分层聚类的方法 ，可选"ward.D", "ward.D2", "single", "complete", "average", "mcquitty", "median", "centroid". default("ward.D2")

基于肿瘤亚簇的下游分析

@param HMM  
当设置为True时，运行HMM预测CNV级别(默认为FALSE)
@param HMM_transition_prob  
HMM中的跃迁概率，默认为1e-6
@param HMM_report_by    
选择：c("subcluster", "consensus", "cell") default: subcluster
注意，report 的执行完全独立于HMM预测。因此，您可以预测子集群，但获得每个cell级别的报告(更大的输出)。

@param HMM_type 
HMM的模型类型；default：i6
选择: (i6 or i3): 
i6: infercnv 6个状态模型(0, 0.5, 1, 1.5, 2, >2) 
i3: infercnv3个状态的模型(del, neutral, amp) 基于正常细胞和HMM_i3_pval配置

@param HMM_i3_pval  
i3状态的p值overlop (default: 0.05)

@param HMM_i3_use_KS    
boolean: 使用KS统计检验估计扩增/缺失分布的均值 (default=TRUE)

肿瘤亚聚类

@param analysis_mode    
选择：("samples","subclusters","cells")
用于图像过滤或HMM预测的分组水平。
default: samples (更快，但是subclusters是更理想的)

@param tumor_subcluster_partition_method    
肿瘤亚簇的定义方法；
选择：('random_trees', 'qnorm')
default：random_trees，慢但最好
qnorm: 根据tumor_subcluster_pval定义的分位数定义树的高度

@param tumor_subcluster_pval    
定义显著肿瘤亚簇的最大p值(默认值:0.1)

去噪参数

@param denoise  
如果为真，根据下面的选项开始去噪
@param noise_filter 
Values +- from the reference cell mean will be set to zero (whitening effect) default(NA, instead will use sd_amplifier below.

sd_amplifier    
Noise is defined as mean(reference_cells) +- sdev(reference_cells) * sd_amplifier default: 1.5

noise_logistic  
use the noise_filter or sd_amplifier based threshold (whichever is invoked) as the midpoint in a logistic model for downscaling values close to the mean. (default: FALSE)

离群修剪

@param outlier_method_bound 
用于边界离群值的方法。
default: "average_bound"
如果设置，将优先使用outlier_lower_bounda和outlier_upper_bound。
@param outlier_lower_bound  
下界以下的异常值将被设置为这个值。
@param outlier_upper_bound  
这个上界上的异常值将被设置为这个值。

其他：大多数忽略即可，采用默认

@param final_scale_limits   
输出的最终热图的比例限制。
default： "auto".
或者选择:"low","high"
@param final_center_val 
输出的最终热图的中心值。
@param debug    
如果为true，则输出调试级别的日志。default： FALSE.
@param num_threads  
并行阶段采用的核数。允许的情况下可以尽量采用多个核，使计算速度更快 (default: 4)
@param lot_steps    
如果为真，保存infercnv objects并在中间步骤绘制数据。
@param resume_mode  
尽可能利用预计算和存储的infercnv对象。(default = TRUE)
@param png_res  
图像的分辨率
@param plot_probabilities   
绘制后验概率的选项(default: TRUE)
@param save_rds 
是否将当前步骤对象结果保存为.rds文件 (default: TRUE)
@param save_final_rds   
是否将最终对象结果保存为.rds文件(default: TRUE)
@param diagnostics  
选择在运行贝叶斯模型后创建诊断图 (default: FALSE)

结果图的解释

得到的quickstart.pdf文件：
图中的数值是"Residual expression values"，可以简单理解为基因表达值的另一种形式；

image.png

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

<meta charset="utf-8">

画图

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包也包含了画图函数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))) #改颜色

image

转载来自：
作者：TOP生物信息
链接：https://www.jianshu.com/p/c01f56bf031f
作者：Geekero
链接：https://www.jianshu.com/p/82fdb69d6083
作者：被格格巫抓到的蓝精灵
链接：https://www.jianshu.com/p/6b44e511f641
作者：xuzhougeng
链接：https://www.jianshu.com/p/cf15feeb4870

参考资料

• 软件安装: https://github.com/broadinstitute/inferCNV/wiki/Installing-infercnv
• 文件定义: https://github.com/broadinstitute/inferCNV/wiki/File-Definitions
• 运行inferCNV: https://github.com/broadinstitute/inferCNV/wiki/Running-InferCNV
• 引用: inferCNV of the Trinity CTAT Project. https://github.com/broadinstitute/inferCNV