外显子拷贝数分析之cnvkit

作者，Evil Genius

跟美国人交流的时候，讨论起外显子的CNV检测，国内用的cnvkit居多，但是我翻了很多的推文，感觉都写的很肤浅，美国人问到的问题都比较深入，涉及到算法核心，所以这一篇我们好好梳理一下cnvkit的运用。

反正不详细的深入理解，糊弄不过去，这一点美国人还是很较真的,在大佬面前真的是一点侥幸心理都不能有。

cnvkit的官网网址https://cnvkit.readthedocs.io/en/stable/

文章在CNVkit: Genome-Wide Copy Number Detection and Visualization from Targeted DNA Sequencing

cnvkit方法使用targeted reads和nonspecifically captured off-target reads来均匀地推断整个基因组的拷贝数。这种组合在目标区域实现了外显子水平的分辨率，在较大的内含子和基因间区实现了足够的分辨率，以识别拷贝数的变化。

知识背景

• 拷贝数变化是包括癌症在内的许多疾病的有用诊断指标。
• 对于临床应用，通常更倾向于对基因组分区(如外显子组或一组疾病相关基因)进行测序，以富集感兴趣的区域，并在更高的覆盖范围内对其进行测序，以提高call变异的敏感性。
• 在目标区域富集过程中，通过杂交捕获目标区域;然而，文库中仍然保留了大量的脱靶DNA，这些DNA被测序，代表了相当大的一部分reads。因此，除了在目标区域获得高覆盖率测序外，脱靶reads提供了全基因组的极低覆盖率测序。虽然脱靶读取本身不能提供足够的覆盖范围来调用单核苷酸变异(snv)和其他小变异，但它们可以在更大的范围内提供有关拷贝数的有用信息。

• CNVKIT利用在靶和脱靶测序读取并应用一系列校正来提高拷贝数检测的准确性。

  
  

pipeline

CNVkit使用on-target reads和非特异性捕获的off-target reads来计算每个样本基因组的log2拷贝比。简而言之，off-target bin是从目标区域之间的基因组位置分配的，off-target bin的平均大小比平均on-target bin大得多，以匹配它们的read counts.然后分别使用目标位置和非目标位置来计算每个间隔内的平均reads深度。然后将目标reads深度和非目标reads深度结合起来，将其归一化为来自对照样本的参考，并对几个系统偏差进行校正，从而得到 log2 copy ratios的最终表。内置的segmentation algorithm可以在log2 copy ratios上运行，以推断discrete copy number segments。

Calculation of off-target intervals

用于计数脱靶reads的基因组间隔最初是从目标间隔的基因组位置计算的。CNVkit根据target regions，将每个target之间的非target区域划分为bin，通常按100 kilobases的顺序划分。作为可选的输入，在创建off-target bins时，可以使用可测序的染色体区域和低映射区域的单独列表来排除端粒、着丝粒和其他不可测序或不可映射的重复区域。
每个连续的非目标区域被划分为大小相等的bin，使得该区域内的平均bin大小尽可能接近指定大小。可以通过计算目标区域平均大小与目标区域测序reads的fold-enrichment的乘积来选择合适的off-target bin大小，使得平均大致相同数量的reads被映射到on - and - off-target bin上。为了最大限度地增加bin的数量，CNVkit将偏离指定的bin大小，将bin放入小的区域，如内含子，这些区域的大小受到限制。还可以指定bin大小的下限，以避免评估非常小的非目标区域，在这些区域中，捕获的reads太少，无法给出可靠的拷贝数估计。一旦生成了一组可靠的off-target bin并保存为BED文件，同一个BED文件可以在CNVkit中重复使用，用于使用相同panel的其他样品的拷贝数分析，并在相同的平台上测序。

Estimation of copy number by read depth

CNVkit coverage命令使用BAM格式的测序reads比对和BED或interval list格式的on或off-target bins的位置，计算样本中每个bin中的log2平均读取深度。对于每个bin，使用pysam计算和求和bin中每个碱基对的read depths，然后除以bin的大小。输出是一个表，其中显示了每个给定bin的平均reads depths，经过log2变换并以所有常染色体centered to the median read depth of all autosomes。

Construction of a copy number reference

在每个基因组bin中，提取每个给定对照样本的read depths。对每个control样本执行reads深度偏差校正。在每个bin中，计算control samples中log2 reads深度的加权平均值，以标记上具有较高或较低覆盖率的bin，the spread or statistical dispersion of log2 read depths indicates bins that have erratic coverage so that they can be de-emphasized at the segmentation step. 也可以使用单个成对的control sample，或者，在没有任何control samples的情况下，可以构建一个“通用”参考，读取深度为log2，分配给所有箱子的扩展为0。在所有情况下，都可以指定一个“男性参考”，其中X染色体bin的预期读取深度是常染色体的一半。

Bin size and resolution

• 人类基因组中的外显子平均大小约为200bp。选择目标bin大小默认值267，这样拆分较大的外显子将产生最小大小为200的bin。由于包含较少 reads的bin会导致更嘈杂的拷贝数信号，因此这种方法确保通过分裂较大外显子产生的bin的“噪声”不会比平均水平差。
• 例如，将目标bin的平均大小设置为100bp，将产生大约两倍的目标bin，这可能会导致更高分辨率的分割。然而，在每个bin中计数的读取次数将减少大约一半，从而增加了bin级覆盖率的方差或“噪声”。过多的噪声bin会使可视化变得困难，并且由于噪声可能不是正态分布的，特别是在存在许多reads为零的bin的情况下，分割算法可能会在低覆盖率样本上产生不太准确的结果。因此，建议总体靶向测序覆盖深度至少为200x至300x，reads长度为100bp，以证明将平均目标bin大小减少到100bp是合理的。
• 对于混合捕获，如果targets are not tiled with uniform density——例如，target panel is designed with a subset of targets having twice or half the usual number of tiles for a fixed number of genomic bases ——不需要做任何特别的事情来弥补这一点，as long as you are using a pooled reference。当测试样本的读深度归一化到pool引用时，log2比率将趋于平衡.

Filtering segments

• cn值，merging adjacent with the same called value.
• Keeping only high-level amplifications (5 copies or more) and homozygous deletions (0 copies) (ampdel).
• Confidence interval overlapping zero (ci).
• Standard error of the mean (sem), a parametric estimate of confidence intervals which behaves similarly.

在每种情况下，根据给定的标准将具有相同值的相邻段合并在一起，并适当地重新计算列值。即使总拷贝数相同，位于不同染色体上或具有不同等位基因特异性拷贝数值的片段也不会合并。

breaks

cnvkit.py breaks Sample.cnr Sample.cns

这有助于鉴定(a)发生不平衡融合或其他结构重排断点的基因，或(b)由于拷贝数信号不一致而难以调用CNV的基因。

Columns:

• gene, chromosome – as in .cns, the gene where the breakpoint occurs and the chromosome it lies on.
• location – the end of the segment to the left of the breakpoint, and start of the segment to the right.
• change – the difference in log2 values between the adjacent segments.
• probes_left, probes_right – the number of probes on each side of the breakpoint within the gene. (Not the same as the number of probes supporting each segment; just the portion within the gene.)

genemetrics

Identify targeted genes with copy number gain or loss above or below a threshold.

The remaining output columns have slightly different meaning depending on whether or not segments were provided. Without segments (.cnr alone):

• log2: Weighted mean of log2 ratios of all the gene’s bins, including any off-target intronic bins.
• depth: Weighted mean of un-normalized read depths across all this gene’s bins.
• weight: Sum of this gene’s bins’ weights.
• nbins: The number of bins assigned to this gene.

With segments (-s):

• log2: The log2 ratio value of the segment covering the gene, i.e. weighted mean of all bins covered by the whole segment, not just this gene.
• depth, weight, probes: As above.
• seg_weight: The sum of the weights of the bins supporting the segment.
• seg_probes: The number of probes supporting the segment.

今天就看panel部分，生活很好，有你更好