NGS中的错误总结（三）——Illumina平台错误率评估
image.png
1 摘要
这是一篇德国波恩大学 LIMES研究所于2018年发表在《SCIENTIFIC REPORTS》杂志的一篇文章，名为《系统评估短序列样品中NGS的错误率和发生原因》。该文章在Illumina测序仪上对单条序列进行了测序，并对获得的数百万reads进行了分析。得出了以下结论：1、尽管PCR这一步骤被传统认为是NGS中主要的错误来源，但在本研究中样品制备过程的index-PCR步骤没有对错误率造成影响；2、pre-phasing 效应在测序过程中持续存在；3、测序平均错误率为每碱基0.24±0.06%，序列突变的百分比为6.4±1.24%；4、在5 ' 和3 ' 端添加固定区域，例如引物结合位点，对突变率没有影响，重测序的样品仍获得高重复性；5、由于NGS设备以及测序流程不同，也会导致phasing效应和其他测序问题各不相同，这篇文章推荐NGS使用者进行错误率以及错误类型分析评估以提高NGS的数据质量及分析结果。
2 简介
目前最常用的测序技术是边合成边测序（SBS）。目前报导的单碱基平均错误率是0.1%，大部分情况是单碱基的替换。此外，该技术会导致固有测序错误：1、颜色或激光干扰、相邻簇间的干扰、2、phasing（移相）3、dimming(调光)。颜色干扰的原因是读取合并碱基时，不同的荧光团之间的激发和发射光谱重叠。虽然颜色干扰可以被校正，但是由于相同的原因，相邻簇之间干扰的问题仍然会存在。phasing(移相)描述了两种现象，这两种现象都可以导致单碱基与其他簇的相位不一致：(1) 如果在一个循环中纳入了两个或者更多的核苷酸则发生pre-phasing(前相移)，这是因为dNTP上标记的叠氮基团(N3)不正确的掉落，当丢失了叠氮基团的dNTP加到（合成链的）3’端之后，它的聚合反应不会终止，而是会继续往前走。当再加上了一个带叠氮基团的dNTP之后，这个聚合反应才停下来。（2）post-phasing(后相移) 是由于酶活性不足或者终止子没有完全移除，导致后续的dNTP碱基无法加入，序列滞后于簇的其余部分(如图1所示)。无法移除的终止子以及激光对DNA链的损伤会导致一个簇中测序的序列减少，从而使读出的荧光变暗。
image.png
目前，碱基检测软件Bustard已经包含了相移的校正（假设相移有固定的频率），包括考虑周围碱基、根据逐个纳入的碱基调整算法。除了这些内在的测序错误，还有一些突变是由样本制备和PCR过程导致的。通过研究paired-end 或者 duplex-DNA重叠区域的互补配对关系，可以去除那些错误的碱基。以上因素导致的突变，可以通过indices 或者barcodes 的方法分析，紧密监控这些错误的发生。
所有以上这些方法，都是为了确定比单端NGS reads更长的序列错误。然而，NGS也可通过体外分析较短的序列（单端reads可以覆盖整条序列）进行参数的优化。本文目的是对体外制备的样品进行测序，获得测序错误的详细信息: 这里采用index-PCR技术在5 ' 和3 ' 末端添加barcode，使12个样品在单个flowcell中同时进行测序。
最终文章研究表明，相移效应是初始错误率的主要原因。去除短序列可以排除相移序列，并确定单个碱基0.25%的真实错误率。此外，相同样本的测序本次实验具有很好的可重复性。文章认为这些发现有助于提高人们对NGS过程中相移效应和实际错误率的认识并更好的使用NGS。
image.png
3 结果
3.1 样本制备对错误率的影响
3.1.1 实验过程
1、为了研究样品制备对错误率的影响，对C12序列进行了测序和分析。C12 是一段结合GFP配体、通过click-chemistry化学修饰的序列，序列模板都是用标准的核苷酸进行合成。
2、对于C12_T_PWO和C12_T_Taq样本，使用PWO或者Taq 聚合酶进行index-PCR。对于C12_T_w/o样本，合成模板时已经包括了index，因此不需要再进行index-PCR。index-PCR之后，所有样品进行混合，消除其他处理步骤造成的影响。
3、对产生的序列进行汇总，分析所有突变序列的占比以及每个碱基的平均突变率，即错误率。
3.1.2 实验结果
1、如表1所示，无论是否进行index-PCR步骤，样品的序列突变频率都没有发生变化。没有进行index-PCR的C12_T_w/o样品，碱基错误率略低于C12_T_Taq、C12_T_PWO样品。
3、如附表1所示，列出了C12_T_PWO 样本占比TOP25的序列信息，并列出了各种突变类型。
4、如图Figure2所示，a-c图显示了C12_T_Taq、C12_T_PWO和 C12_T_w/o样品的各个位置的突变频率，随着位置的增加，突变频率变大。d图显示了4个样品中原碱基的平均突变频率。C12_T_Taq样品所有碱基的突变频率比C12_T_PWO样本略低，没有进行index-PCR的样本C12_T_w/o，突变频率最低。e图显示了4个样品中突变为指定碱基的平均突变频率。f 图显示了C12 序列中的原始碱基分布。图e 中突变为指定碱基的突变频率与原始碱基分布一致。
image.png
注：EdU与BrdU是一种胸腺嘧啶核苷类似物，它们能够在细胞增殖时期代替胸腺嘧啶(T)渗入正在复制的DNA分子，然后分别利用免疫荧光技术、与Apollo荧光染料特异性反应检测细胞增殖情况
image.png
image.png
3.2 核苷酸修饰对错误率的影响
1、为了研究碱基修饰对NGS误差率的影响。以5’-乙炔-脱氧尿苷(EdU)代替胸腺嘧啶合成了C12_EdU 样本序列模板。在非保护条件下，约20%的EdU转化为酮副产物(KdU)，这可能对PCR的保真性造成影响。
2、与其他只含标准碱基的C12样品相比，C12_ EdU突变序列比例和错误率明显增加（表1所示），原碱基突变频率和指定碱基突变频率也显著增加（图2 d,e所示）。与其他C12样本类似，随着位置的增加，突变频率变大（图2 g 所示）。
3.3 分析重复序列
1、Table2 中描述了分析的重复序列。GATC和G4A4T4C4由于正链和负链不能退火，所以不能使用NGS进行测序。
2、所有初始重复序列使用引物结合的FT2文库进行分析。如Table3所示，重复序列的突变序列占比和碱基错误率都低于C12样品。FT2_T4G4C4A4样品虽然具有较低的碱基错误率，但是与FT2_GATC样品相比具有较高的突变序列占比。如图3 a b 所示，这是因为FT2_T4G4C4A4样品每个重复单元的前3个碱基的具有极低的突变频率，但是最后一个碱基却有较高的突变频率。FT2_GATC样品与C12样品一致，随着位置的增加，突变频率变大。
3、如图3 c d所示，原始碱基突变为哪种碱基具有特殊的偏好性。图3 e 汇总了突变碱基的关系，碱基更偏好于突变成随后的碱基。
4、Table4 汇总了各个样品中一个碱基突变成另一个碱基的比例。在完全随机状态下，一个碱基突变成另一个碱基的比例是33.3%，而表4 中的比例为FT2_GATC(64%)至FT2_T4G4C4A4(84%)。
5、如图3f 所示，随着所有测试样本中相同的连续核苷酸的数量的增加，后续核苷酸的突变频率稳步增加(从65%到85%)，并且与核苷酸的顺序无关。
image.png
image.png
image.png
image.png
3.4 测序数据的重复性和引物结合位点序列对突变率的影响
1、为了评估测序数据的重复性，FT2_GATC 和 FT2_T4G4C4A4样品进行了重新建库、测序分析，见FT2_GATC_II、FT2_T4G4C4A4_II 样品。为了评估引物结合位点对突变率的影响，用D3 文库结合TCGA和T4G4C4A4进行了重新分析。
2、图4a,b以及从Table 3和Table 4可以看出，尽管FT2-GATC获得的序列数量相差5倍，但错误率、突变频率、突变序列数量以及对后续核苷酸的突变频率的变化很小.
3、尽管D3-TGCA的突变频率和错误率略低于FT2-TGCA，但对后续核苷酸的突变频率高于FT2-TGCA，而D3-T4G4C4A4与FT2-TGCA的突变频率无差异。
image.png
3.5 去除缩短序列排除phasing效应
1、随着序列长度的增加，突变频率的增加和后续核苷酸的高突变率可以在所有样本中被识别出来，这可能是由于相位效应，我们的目的是将这些从分析中排除。
2、因此，通过评估了不同样本中26个最丰富的序列，发现包含pre-phasing 效应的序列被缩短。显然，序列的缩短也可能是由于碱基缺失导致的，pre-phasing效应和碱基缺失是无法区分的。
3、如图5 a b 和 Table 5所示，去除缩短序列之后，所有测试样本的突变序列百分比和错误率都有大幅度的下降。 如图5 cd 所示，去除缩短序列之后，之前观察到所有样本随着序列长度增加突变频率增加的趋势也完全消失。
4、如图5 e 和 Table 5 Tabel6所示，碱基更偏好于突变成随后的碱基的突变频率下降到预期的33.3%左右，现在与相同连续核苷酸的数量无关。
5、Table 6 汇总了去除缩短序列后的变化。分析序列的数量平均减少5.2%，未突变序列增加5.6%。相比之下，错误率下降了79%。所有这些都是非常明确的迹象，表明我们忽略了大部分通过pre-phasing 产生的突变序列，而没有排除高百分比的序列。
6、因此，重新分析了样本，以确定NGS中的“真实”错误率。在排除C12_EdU样品后，所有样品的每个碱基的平均错误率为0.24± 0.06%，平均突变序列百分比为6.4 ± 1.24%。
7、图6 展示了去除短序列后，所有样品中一个核苷酸替换另一个核苷酸的突变百分比。平均突变率是在排除C12_EdU样品后计算的。
image.png
image.png
image.png
image.png
3.6 去除缩短序列对SELEX样品的影响
1、为了确定在体外筛选过程中排除短序列对样品的影响，我们重新分析了选择了碱基修饰的GFP配体样品。图7显示了在不同的选择周期的四种不同模式频率。在去除缩短序列的前后略有不同，但总体趋势和绝对频率都没有改变。
image.png
4 材料和方法
4.1给读者的建议
我们的突变数据显示了phasing效应的巨大影响，我们可以通过去除所有缩短的序列来排除phasing效应。尽管文献中已经显示了问题的原因，但它并不是每一个测序仪器的所面临的突出问题。因此文章建议NGS用户可以定期的使用这篇文章发布的重复序列进行测试，以更好的了解测序仪器的错误类型和频率。关于使用NGS分析指数富集(SELEX)，缩短序列的去除也可能导致绑定序列的去除，因为缩短的序列也可能是富集文库所固有的。如果缩短的序列占富集文库的比例较大，我们的策略并不能很好的起作用。如果是这样的话，我们建议用计算方法排除phasing 效应。此外，单序列分析在测序过程中受到实际突变和测序错误的影响要比这里基于序列家族的分析影响大得多。
4.2 样本制备
样本准备和测序在不同的runs。
- C12_EdU,
- C12_T_wo, C12_T_PWO, and C12_T_Taq,
- GATC, G4A4T4C4, FT2_GATC, and FT2_G4A4T4C4,
- FT2_GATC_II, FT2_G4A4T4C4_II, FT2_G2A2T2C2, FT2_G3A3T3C3, FT2-TGCA, D3-TGCA, FT2-T4G4C4A4,and D3-T4G4C4A4,
- GFP-SELEX samples.
C12_EdU、C12_T_PWO、GFP-SELEX、C12_T_Taq 进行 index-PCR操作。其余样品进行常规的PCR操作。PCR产物进行跑胶、回收、PCR产物清洁。然后连接接头进行测序。测序采用7HiSeq1500 和 NextSeq500 平台，76 base pair 和 7 index bases进行测序。
4.3 NGS分析
分析NGS数据所用的软件是 COMPAS，先根据样本唯一barcode信息对数据进行拆分。按照序列的相似度进行聚类，并计算每种类型序列的reads支持数。为了去除phasing的影响，进行了了除缩短的序列去除，只有正确长度和较长的序列才会进行后续的分析。
4.4 突变分析
突变序列频率用于计算有突变的reads占所有reads的比例。每个核苷酸的突变率是由核苷酸分布从1减去特定位置正确核苷酸的频率计算出来的。一个特定核苷酸的每个核苷酸的突变频率的平均值和标准差被给出为“突变nt”。所有突变核苷酸的总体平均值和标准偏差就是“错误率”。
为了计算“突变成nt”的平均值和标准差，我们考虑了该核苷酸不是原始核苷酸的所有位置。
4.5 统计分析
Shapiro-Wilk 用于检测正态性。单因素方差分析、双尾t检验、非正态分布采用Kruskal-Wallis检验。
5 重要信息
5.1 文献中Illumina测序仪的错误率
image.png
5.2 pre-phasing 效应造成错误影响的原因
如下图所示，pre-phasing效应会在一个测序周期中加入互补配对的TA两个碱基，但是多插入的A碱基测序成像是不能显示出来的，虽然后续的序列会逐渐的加入，但是最终的测序序列会变短,而且多添加碱基之后的所有碱基会被认为是突变。这也是虽然pre-phasing影响reads很少，但是会对错误率造成重大影响的原因。
image.png
6 参考文献
[1] Pfeiffer, F., Gröber, C., Blank, M. et al. Systematic evaluation of error rates and causes in short samples in next-generation sequencing. Sci Rep 8, 10950 (2018). https://doi.org/10.1038/s41598-018-29325-6
[2] 测序原理参考：https://www.cnblogs.com/think-and-do/p/6638157.html
