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个输入数据
- 单细胞RNA-seq表达量的原始矩阵
- 注释文件，记录肿瘤和正常细胞
- 基因或染色体位置文件
第一个是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
网友评论