10X单细胞个性化分析之CNV篇
最近父亲住院，需要人照顾，需要资金，所以不要脸的求一发大家的支持。
目前单细胞分析而言，分析方向大致包括以下几个方面，1）器官发育（这个用空间转录组更为合适）；2）疾病样本，尤其是肿瘤样本的分析研究；3）其他非模式物种的细胞图谱。其中对于肿瘤样本的分析，在基因组研究中CNV的分析占了很重要的一部分，CNV(Copy number variation, 拷贝数变异)是由基因组发生重排而导致的, 一般指长度为1 kb 以上的基因组大片段的拷贝数增加或者减少, 主要表现为亚显微水平的缺失和重复。CNV 是基因组结构变异(Structural variation, SV) 的重要组成部分。CNV位点的突变率远高于SNP(Single nucleotide polymorphism), 是人类疾病的重要致病因素之一。而对于单细胞转录组，识别肿瘤细胞和发生的CNV事件同样重要，实际分析中，也经常用软件来判断肿瘤细胞。当然还可以做肿瘤异质性、克隆进化方面的探索，而本篇来介绍单细胞数据的CNV分析。关于单细胞CNV分析，目前主流的分析软件为inferCNV和后起的“新秀”copyCAT，本篇就从这两个软件着手，体现CNV分析在单细胞研究中的重要作用。
InferCNV章节
InferCNV用于探索肿瘤单细胞 RNA-Seq 数据，以确定体细胞大规模染色体拷贝数改变的证据，例如整个染色体或大段染色体的gain或loss。这是通过与一组参考“正常”细胞相比，探索肿瘤基因组位置上基因的表达强度来完成的。 生成的热图说明了每个染色体上的相对表达强度，并且与正常细胞相比，肿瘤基因组的哪些区域过度表达或更少表达通常变得很明显。
用inferCNV判断肿瘤细胞的CNV事件通常包括以下几个步骤(如下图)：
1）样本的基础质控和注释；
2）选择合适的reference；
3）依据基因在染色体上的位置对基因进行排序；
4）数据处理，包括肿瘤细胞与ref的信号比较去除、数据均一化处理、降低噪音等过程；
5）CNV最终的预测。
从分析过程中来讲，inferCNV需要的输入文件包括：表达矩阵、细胞注释信息、基因在染色体上的位置信息。
inferCNV分析过程示意图
Reference的选择
使用inferCNV分析单细胞转录组，确定reference是最关键、也是最开始需要考虑的内容，如果不指定reference，那么软件默认会把样本中所有细胞的基因平均表达值作为“基线”来识别肿瘤细胞，这种方法目前没有文章引用，原因也很简单，混淆所有细胞作为reference，其中也包括了肿瘤细胞，无法确定分析结果的准确性。所以做inferCNV最为基础和关键的地方，还是前期对样本的质控和细胞注释，选择合适的reference，在此基础上才可以合理地进行inferCNV分析。
最佳的reference选择是对应肿瘤细胞类型的正常细胞类型，也是高分文章通常的做法，例如上皮细胞癌变，那么就以正常的上皮细胞作为reference来分析肿瘤细胞的CNV事件，这样分析的结果可靠，但是有一个问题，尤其对于人的肿瘤样本，往往取不到正常的组织区域，就会给CNV分析带来不小的麻烦，有些肿瘤样本会带有癌旁区域，癌旁部分含有正常的细胞类型，但是在细胞解离的过程中跟肿瘤细胞混淆在一起，后续的分析无法很好的区分，这种情况下，只能退而求其次，选择免疫细胞（T、NK等）作为reference，同时遵守一个原则，尽量多的选择reference细胞类型，最大限度保证结果可信，在文章A single-cell and spatially resolved atlas of human breast cancers中，就将免疫和内皮细胞最为reference来推断肿瘤细胞的CNV事件，下图为中设置E8细胞作为reference分析得到的CNV结果,可见选择合适的reference会得到良好的分析结果，不仅可以判断细胞类型发生的CNV事件，也可以分析肿瘤细胞内部的异质性。
拷贝数变异分析肿瘤异质性
InferCNV的算法原理
InferCNV算法的详细步骤涉及以下内容：
1）过滤基因：从计数矩阵中删除那些在少于“min_cells_per_gene”中表达的基因，这一步类似于样本质控过程中的基因去除。
2）测序深度的归一化（总和归一化）：read counts per cell are scaled to sum to the median total read count across cells。 值不是每百万计数 (cpm) 等指标，而是每中位数总和的计数（这一点区别于Seurat分析单细胞的均一化）。
3）对数转换：单个矩阵值 (x) 转换为log(x+1)，这里对数转换的作用与Seurat分析中的相同。
4）center by normal gene expression：从对应基因的所有细胞中减去正常（参考）细胞中每个基因的平均值。 由于此减法是在对数空间中执行的，因此这有效地导致了相对于正常细胞平均值的对数倍变化值。
5）对数倍数变化值的阈值动态范围。 abs(log(x+1)) 超过'max_centered_threshold' (default=3) 的任何值都被设定为该值（设置了最高上限）。
6）chromosome-level smoothing：对于每个细胞，沿每个染色体排序的基因具有使用加权运行平均值拟合的表达强度。 默认情况下，这是一个包含 101 个基因的窗口，具有pyramidinal weighting scheme。
7）centering cells：如果大多数基因不在 CNV 区域中，每个细胞的中心表达强度中值设定为零。
8）相对于正常细胞的调整：再次从肿瘤细胞中减去正常值的平均值。 这进一步补偿了拟合处理后产生的差异。
9）log转换被还原，这使得amplification 或 deletion的证据在平均值周围更加对称。
上述就是推断CNV分析的基本过程，但是通常为了更加准确的推断CNV事件，往往还要添加两个步骤de-noising filters和HMMs算法。
inferCNV de-noising filters
降噪的目的是降低噪音（正常细胞中的残余信号），同时保留肿瘤细胞中可被解释为 CNV 的信号。
基础分析结束后的正常信号保存在初步的 inferCNV 对象，该对象已被smoothed、centered，并减去了正常（参考）细胞的平均值，如下图：
reference的信号表达特征
为了确定分析得到的是真正的CNV事件，需要对肿瘤细胞的CNV信号进行检验，也就是降噪，inferCNV通常有三种方法处理这一过程。
1）可以使用“noise_filter”属性设置与平均值的特定阈值偏差，如下图：
设定阈值降噪
如上图，设置0.1为过滤阈值，也就是在这种情况下，reference基因表达0.9~1.1以外的基因表达被判定为CNV事件，高于1.1为gain，低于0.9为loss，这也是inferCNV默认的过滤方法。
2）动态阈值设置：可以使用“sd_amplifier”设置调整阈值。 可以使用 1.5 * reference基因表达的标准差进行过滤，如下图：
动态阈值
如前所述，低于最小阈值为loss，高于最大阈值为gain。
3）通过 sigmoidal（逻辑）函数调整强度（软阈值）：可以通过应用 sigmoidal 函数来应用过滤梯度，而不是应用严格的阈值，该函数可以减少接近均值的强度，而不是更远离均值的强度，如下图：
sigmoidal软阈值分析
HMMs预测模型
目前inferCNV支持两种基于 HMM 的 CNV 预测模型，称之为 i3 和 i6 模型。每种方法都对已通过标准 inferCNV 处理的对象操作，包括减去与“正常（参考）”细胞对应的信号和smoothed操作。
1）i3模型：loss、normal、gain三种状态，如前所述，大多数信号对应于normal而异常信号强度对应于CNV。
2）i6 模型：一种六态 CNV 模型，可预测以下 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.
此外，预测的 CNV 区域使用贝叶斯网络进一步分析，以计算每个细胞属于给定状态的 CNV 区域的后验概率。 具有高于最大阈值的平均后验概率正常（无 CNV）的 CNV 区域作为可能的假阳性预测被移除。
InferCNV的实际运用
在肿瘤研究中，可以通过CNV预测分析区分肿瘤细胞和非恶性细胞。2018年纽约大学计算医学研究所等单位的研究人员在Nature Biotechnology发表了利用单细胞和空间转录组研究胰腺导管癌（PDAC）异质性的文章。为了区分癌细胞和非恶性导管细胞，该研究对PDAC-A和PDAC-B 2例单细胞数据进行了CNV预测分析。发现PDAC-A中高表达TM4SF1（簇1）和S100A4（簇2）的两个细胞群及PDAC-B中高表达TM4SF1的一个细胞群表现出拷贝数变异特征。通过免疫荧光验证发现PDAC-A中TM4SF1和S100A4在恶性导管细胞中表达，PDAC-B中TM4SF1与恶性细胞标志物KRT19共定位，结合CNV预测结果证实了PDAC样本存在转录不同的肿瘤细胞群。
拷贝数变异分析区分胰腺导管癌肿瘤细胞与非恶性细胞
在肿瘤研究中，可以通过CNV预测分析探索肿瘤的克隆进化。2020年美国迈阿密大学等单位的研究人员在Nature Communications发表了利用单细胞测序研究葡萄膜黑色素瘤进化复杂性的文章。该研究对8例原发癌和3例转移癌进行单细胞CNV预测分析，发现不同样本间存在显著的拷贝数变异差异，揭示了葡萄膜黑色素瘤潜在的肿瘤间异质性。进一步根据某个CNV在细胞中的占比构建进化树，发现驱动葡萄膜黑色素瘤突变的3条进化轨迹—低度转移肿瘤中的EIF1AX突变、中度转移肿瘤中的SF3B1突变及高度转移肿瘤中的BAP1突变，绘制了葡萄膜黑色素瘤的进化轨迹及发展机制。
拷贝数变异揭示葡萄膜黑色素瘤的克隆进化
InferCNV分析的不足
1）为分析来自第一代 scRNA-seq 技术的数据而设计的，技术具有较低的细胞通量和较高的覆盖深度。
2）不适用于分析来自新开发的高通量 scRNA-seq 平台（微滴和纳米孔平台）的数据，这些平台执行全转录组扩增和仅在非常稀疏的覆盖深度下对 mRNA 的 3' 或 5' 端进行测序（10X的单细胞技术具有这个特点）。
3）不能准确地解决特定染色体断点的基因组位置或从非整倍体拷贝数谱中对肿瘤和正常细胞进行分类。
CNV分析的“新贵”---copyCAT
CopyKAT 的工作流程将贝叶斯方法与层次聚类相结合（inferCNV其实也用到了层次聚类），以计算单个细胞的基因组拷贝图谱，并从高通量 3' scRNA-seq 数据中定义克隆亚型。分析流程将唯一分子标识符 (UMI) 计数的基因表达矩阵作为计算的输入。分析从每行的基因注释开始，按照它们的基因组坐标对它们进行排序（跟inferCNV的原理一致）。执行 Freeman-Tukey 变换以稳定方差，然后执行多项式动态线性建模 (DLM) 以smoothed单细胞 UMI 计数中的异常值。下一步是检测具有高置信度的正常细胞（reference），以推断正常 2N 细胞的拷贝数基线值（软件CopyCAT自动检测）。为此，将细胞细分为几个小的聚类（层次聚类），并使用高斯混合模型 (GMM) 估计每个聚类的方差。通过遵循严格的分类标准，具有最小估计方差的cluster被定义为“reference”。当数据只有少数正常细胞或肿瘤细胞具有接近二倍体基因组且拷贝数畸变 (CNA) 事件有限时，可能会发生潜在的错误分类。在这种情况下，CopyKAT 提供了一种“GMM 定义”模式来逐个识别二倍体正常细胞，其中假设单个细胞中基因表达的三种高斯模型的混合代表基因组gain、loss和中性状态。当处于中性状态的基因占表达基因的至少 99% 时，细胞被定义为“normal”细胞。
为了检测染色体断点（chromosome breakpoints），整合了泊松伽马模型和马尔可夫链蒙特卡洛 (MCMC) 迭代来生成每个基因窗口的后验均值，然后应用 Kolmogorov-Smirnov (KS) 检验来加入在它们之间没有显著差异的相邻窗口方法。为了加快计算速度，将数千个单细胞分成clusters，找到一致的染色体断点并将它们合并在一起，形成样本中整个细胞群的基因组断点的联合。然后将每个窗口的最终拷贝数值计算为跨越每个细胞中相邻染色体断点的所有基因的后验平均值。通过将基因重新排列到 220-kb 可变基因组bin中，进一步将得到的拷贝数值从基因空间转换为基因组位置，从而以大约 5 Mb 的分辨率获得每个单细胞的全基因组拷贝数谱。基因组分辨率是根据整个基因组的中位相邻基因距离（~20 kb）乘以基因窗口的大小（25个基因）来估计的（精度高于inferCNV）。然后对单细胞拷贝数数据进行层次聚类，以确定非整倍体肿瘤细胞和二倍体基质细胞之间的最大距离；但是，如果基因组距离不显著，切换到 GMM 定义模型来逐个预测单个肿瘤细胞。最后，对单细胞拷贝数数据进行聚类以识别克隆亚群并计算代表亚克隆基因型的共有谱，以进一步分析它们的基因表达差异，流程图如下：
copyCAT分析流程示意图
Copycat基因组拷贝数分辨率的理论估计
为了估计从单细胞 RNA 数据推断出的拷贝数谱的预期分辨率，需要GRCh38 (v28) 中所有基因的 BED 文件。因为染色体 Y 不包括在拷贝数计算中，只考虑了位于染色体 1-22 和染色体 X 上的基因，它们共有 56,051 个基因。通过取基因起始位置和基因结束位置的平均值来估计单个基因的基因组中心位置。接下来，根据基因组位置对所有基因进行排序，并通过计算两个基因中心之间的距离来估计两个相邻基因之间的距离。总的来说，在整个基因组中定义了 56,028 个基因区间。从染色体 1-22 和染色体 X 中，基因区间的数量如下：5,127, 3,872, 2,925, 2,430, 2,779, 2,802, 2,292, 2,189, 2,137, 3,189, 2,857, 1,279, 2,152, 2,081, 2,440, 1,133、2,917、1,350、795、1,300 和 2,281。整个基因组中基因间隔的第一四分位数、中位数、平均值、第三四分位数和最大值如下：9,430 bp、24,532 bp、52,806 bp、58,485 bp 和 21,765,992 bp。因为基因区间的大小分布严重向右倾斜，计算了中值来估计拷贝数分辨率。因为需要在pipeline中的整个单细胞群中检测到至少 7,000 个基因，所以这个数字相当于基因检测率的中位数 7,000/56,051 ≈ 12.5%。最后，将分析中的最小基因间隔计算为每个基因间隔 24,532 bp ÷ 12.5% ≈ 200 kb。使用 25 个基因窗口启动拷贝数分析；因此，估计片段的最小大小为 200 kb × 25 = 5 Mb，用于检测每个细胞基因组中的拷贝数事件的基因组分辨率。
Copycat分析细节
同样地，copycat的输入文件也需要三个：表达矩阵、注释信息和基因位置文件，窗口的设置在25~200之间（inferCNV默认是50），在数据处理和分析结果方面大多借鉴了inferCNV，下图是copycat和inferCNV的分析结果比较。
图片.png
从结果来看，copycat检测的CNV与inferCNV基本一致，在细节方面copycat表现更好一点，尤其在断点处基因的分析，分析更加精细化。
Copycat的局限性
并非所有癌症类型都具有可用于区分正常细胞和肿瘤细胞的非整倍体拷贝数事件。特别是，小儿癌症和造血系统癌症（例如AML和CLL）的拷贝数变化很少，因此可能不适合CopyKAT分析。另一个限制是，CopyKAT主要限于基于整个基因组读取深度的变化来检测CNA事件，而不能用于检测其他有助于基因组多样性的基因组事件，包括染色体结构重排、插入、缺失和体细胞突变。此外，由于3''scRNA-seq数据的技术差异，CopyKAT无法在具有独特基因型的单个细胞的基因组上提供可靠的拷贝数信息。这使得CopyKAT更适合于分析许多细胞已扩增并具有相似基因型的肿瘤中亚克隆，而不是分析复杂细胞或极为罕见的亚群。CopyKAT一个潜在问题是，当scRNA-seq数据集没有任何肿瘤细胞时，CopyKAT可能会尝试错误地检测具有最高基因表达水平的簇中的CNA事件。在这种情况下，推断的CNA事件将与这些癌症中已知的细胞遗传事件不一致，具体需要忽略。
写在后面
CNV分析在单细胞肿瘤样本中占据了重要的分析篇幅，在预测基因发生的CNV事件中即表征了肿瘤内的关键变化，也体现了瘤内的异质性，对于我们认识肿瘤起到了非常关键的作用；同时也要认识到，单细胞肿瘤样本中的CNV推断对于样本前期的质控处理有很高的要求，同时也要添加注释信息，以此为基础来判断CNV事件，这就要求再分析的过程一定要做好基础分析，个性化的分析才足够的可靠、可信。
文献
[1] Anoop P. Patel, Itay Tirosh, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014 Jun 20: 1396-1401.
[2] Gao R , Bai S , Ying C H , et al. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes[J]. Nature Biotechnology, 2021:1-10.
[3] Moncada R, Barkley D, Wagner F, et al. Integrating Microarray-based Spatial Transcriptomics and Single-cell RNA-seq Reveals Tissue Architecture in Pancreatic Ductal Adenocarcinomas[J]. Nature Biotechnology, 2018, 38(3):333-342.
[4] Durante M A, Rodriguez D A, Kurtenbach S, et al. Single-cell Analysis Reveals New Evolutionary Complexity in Uveal Melanoma[J]. Nature Communications, 2020, 11(1):496.
