《NGS攻略》之Dup传

导语

高通量测序技术现已被广泛应用于科研和临床等多个领域。近十年来各种测序技术和平台百花齐放，其中有一些已经退出了历史舞台，而另外一些在数千个实验室广泛使用，这些技术之间的区别在哪里？能够被大范围推广和延续的技术奥秘是什么？在NGS领域哪些指标是关键性技术指标？针对这些疑问，从今天起我们将推出《NGS攻略》，与大家共同探讨其中的技术细节，包括测序质量、测序读长、数据产量、重复序列等关键参数。这一期我们从重复序列Duplicate reads（以下简称Dup）讲起。全文约5000字，推荐阅读时长8分钟。

阅读索引

1、Dup的背景

2、Dup的来源

3、Dup的分类解读

4、实际数据测评

5、总结

阅读重点

• 造成Dup rate偏高的真实原因是什么？

• PCR放大成百上千倍，为什么NGS的Dup rate只有十位数甚至是个位数？

• 有没有办法降低测序过程中的Dup？

1、Dup的背景

所谓Dup，即重复序列Duplicate reads，指通过测序得到的两对或两对以上的Pair-End Reads同时比对到参考基因组上相同的起始和结束位置的序列（如图1所示）。这些重复序列在总测序序列中占比简称为Dup rate。由于这些重复序列并不能带来额外信息，相反会影响变异检测结果准确性，因此下游生信分析中这些重复序列是需要去除的去掉，这也就意味着Dup rate越高，数据利用率越低，测序成本浪费的也就越多。此外，对于不适于去除Dup的基因表达定量测序中，高Dup rate也很大程度上影响了中低表达量转录本定量的准确性。

图1.Duplicate reads展示

（点击查看大图）

2、Dup的来源

对于Dup的来源，在NGS的“江湖”中有各种各样的“传说”，很多大侠认为主要是从样本建库端引入的（或者说认为建库PCR是造成Dup rate偏高的罪魁祸首），但其实测序产生的Dup reads来源很广泛，客观来讲主要有以下几个方面：

1、样本本身的Dup；

2、文库构建中扩增引入的Dup（即PCR Dup）；

3、测序前信号放大（荧光信号采集单元生成过程）引入的Dup；

4、芯片测序过程中引入的光学Dup。

大多数人认为NGS数据的Dup主要来自于上述第二种。但据真实案例显示，第三种荧光信号采集单元生成的过程中引入的Dup和第四种芯片测序过程中引入的光学Dup居多。而第一种样本本身的Dup反倒是我们生物学中应真正关注，并与测序数据的Dup加以区分的一种“有用的”Dup。

3、Dup的分类解读

第一种Dup——样本本身的Dup

大多数学者主要以组织样本为研究对象，这里以人基因组(大概6.6pg/copy)为例，如果获得100ng的gDNA就意味着大概有100ng/3.3pg≈30000个copy的人基因组。不同细胞相同编号染色体在基因组片段化过程中是有可能产生一些起始位置和终止位置相同的分子片段的，这些分子产生的几率实际上与样本DNA的总量有关，而不是建库和测序过程中引入的Dup。这就是为什么在做PCR-free文库时大家还是会发现存在1.5-3%比例Dup rate的原因之一。当然，由于这部分比例非常之低，对于结果影响基本上可以忽略不计。

第二种Dup——PCR Dup

同样以人的基因组DNA文库制备为例，我们通常认为理论情况下PCR扩增6个循环已经把样本分子数量放大64倍（2⁶=64），最起码也是20-30倍（如按照1.6-1.7的扩增效率计算）。这些PCR扩增产生了Dup reads理论上应该体现在最终的测序Reads中，约占20/21*100% ~ 30/31*100%,这样算来应该是高达96%以上Dup rate？但一组真实数据为例，如表1所示，相同文库基于不同平台的数据表现，实际的Dup rate并没有这么高，而且其中BGISEQ平台的Dup rate更是稳定的在个位数波动。这是怎么做到的呢？PCR放大成百上千倍，为什么NGS的Dup rate只有十位数甚至是个位数呢？

表1  相同文库在不同平台测序前后结果对比

（点击查看大图）

这里我们采用Broad研究所的大神Eric Vallabh（著名的遗传性阮病毒学家）对此问题的专业解释（高能预警！此处有大量公式出没，可直接跳至模型推导结果），基于一个数学模型大体估算PCR Dup rate:

假设1.每条待测DNA分子链/球进入纳米孔或与测序芯片Spot结合是完全随机的，不具有任何的Loading Bias。

假设2.待测DNA分子链/球Loading到测序芯片纳米孔/点中，只考虑单克隆的情况（即每个孔/点只进入或结合的1条分子链或1个DNB，其他多克隆或无结合的情况无法产生有效荧光信号，所以也不会对后续总reads数中的Dup rate计算造成影响）。

假设3.不考虑PCR扩增的Bias(即认为每一条分子在每一轮PCR中都获得了扩增机会)。

我们还是以人基因测序为例并基于以上假设，对于接头连接产物进行6个PCR循环的扩增，要求最终得到约1μg PCR扩增产物，那这也就意味这我们需要1μg/（2⁶）=15.6ng的DNA，假设待测分子片段长度平均为200bp，根据{（15.6×10³pg）/（6.6pg/细胞基因组）}×2倍体×（30亿对碱基/基因组）≈14×10¹²bp的公式计算出PCR前样本本身的Unique Molecules（此处不考虑样本本身Dup）数量如下：

1个标准的人WGS按40X测深120G数据量计算，在100PE读长条件下需要600M reads，即需要6×10⁸个单克隆纳米孔或点产生600M的PE reads。

经过6轮PCR扩增后，7×10¹⁰个Unique分子，每一个都有64个Copy，即有6×10⁸个孔，而有7×10¹⁰个Unique分子簇（每簇又有64个Copy）想公平的、完全随机的进入这些孔（前方高能！忘记了什么是二项分布、泊松分布、正态分布、超几何分布的童鞋们请自行百度）。

就此，这完全转化成了一个数学问题，一个概率的二项分布模型，一个符合二项分布模型的实验应具有一下基本属性：

    1.需要进行n次实验（需进行6×10⁸次，每次从7×10¹⁰个分子簇中随机选1个）；

    2.每次实验都是相互独立的；

    3.事件发生的概率p是固定的（每次实验，每一个分子簇被选中的概率为    ）；

    4.每次实验的结果只有两种，只有成功选中（p）或者未选中（q=1-p）。

又因为当p足够小，n足够大时，泊松分布（Poisson distribution）是二项分布模型的近似，因此我们直接使用泊松分布公式对PCR Dup rate进行估算，公式如下：

其中λ为期望值：

根据以上参数，使用R来做7×10¹⁰个Unique分子簇中同一个分子簇被随机抽取获得0次、1次、2次… …n次呈现机会的概率直方图（获得呈现>2次的，即为Dup）。

 图2：使用R调用dpois（x，lambda）作图

由于直方图展示的内容有限，因此设置只展示从呈现0到10次的概率，步长为1，直方图结果见图3 ：

图3.λ=0.0085时的呈现次数（0～10）概率直方图

显而易见的是7×10¹⁰个Unique分子簇绝大部分的都没有幸运的获得呈现的机会（即“0”）我们睁大眼睛能看到呈现1次的幸运儿，但为了更好的看看更幸运的“>2”的部分（即PCR Dup），我们将“0”值去掉（操作步骤如图4所示），再重新做直方图，如图5所示。

   图4：去除“0”值，重新计算呈现概率

图5.λ=0.0085时呈现次数（1～10）的概率直方图

我们利用上述数据调用dpois（x，Lambda）来计算PCR Dup的结果，计算公式如下：

---

sum(dpois(×,lambda=0.0085)/×1)/(1-dpois(0, lambda=0.0085))

按此公式计算所得，（7×10¹⁰个Unique分子簇呈现1次的概率）/（7×10¹⁰个Unique分子簇呈现不小于1次的概率）约为99.79%，即PCR Dup rate为0.21%（这与实际实验操作中大概4%的PCR Dup rate存在一定差距，主要原因是由于此算法忽略PCR Bias因素）。

当我们将PCR循环数增加到9个时（PCR产物总量仍为1μg，起始样本的Unique Molecules数量降低为9×10⁹个，每个分子拥有2⁹个Copy，λ=np≈0.067）同样可以得到概率直方图（图6）。

图6.λ=0.067时呈现的次数（1～10）概率直方图

经运算，PCR Dup rate为1.7%，Eric还计算了PCR循环数为15时，PCR Dup rate的理论值已经达到了15%（由于文章篇幅原因不再赘述，有兴趣的童鞋可查阅参考文献原文数据）。^【1】

综上，该简化的泊松分布模型有助于我们理解关于PCR Dup rate的问题：虽然PCR将待测分子放大了成百上千倍，但是相对于数量远远多于纳米孔/点数的Unique分子来说，能在茫茫人海中被1个孔随机选中已是万万幸，更何况是再次随机选中同一个Unique分子簇中的Copy形成Dup呢？

另外PCR前Unique分子数量也至关重要，当分子的数目（）很少，或者抽样次数（n↑）大幅增加时，多次抽到同一个Unique分子簇的概率也就会大大增加了，这就是为什么我们做降解样本或者做小基因组物种，当数据量很高时得到Dup rate相对正常完整的大基因组高的原因。

那么，我们如何能够避免这些分子被误判为Dup reads？其实现阶段一个通用的方案就是分子序列标签（Unique Identifier，UID或者Unique Molecular Identifer，UMI），全基因组片段化后产生的DNA片段随机的被UMI标记，在后续生信分析中，利用UMI就可以将样本本身来源的Dup与PCR引入的Dup区分开来（原理见图7）。

另外，还有一些采用ExAmp（Exclusion Amplification）+NanoWell技术的测序平台上突出的“标签跳跃”（index hoppig）问题也正是采用UDI（Unique Dual index）+UMI技术来“缓解”，否则index hopping很大程度上降低了低频突变检测的精确性。^【2，3】

图7.UMI原理示意图

第三种Dup——测序前信号放大（荧光信号采集单元生成过程）过程中产生的 Dup

如图8所示，DNA双链加上接头（A B接头）后经过后续的PCR放大，其正链和负链会分别形成双链PCR产物文库（两种PCR产物文库插入片段相同，但接头连接方式是180°反转）。

图8 分子信号放大Dup来源示意图

其他平台的文库形态即为上图所示的2种双链PCR产物，该双链产物随后在NaOH碱溶液中变性为2条游离态单链，如果在文库加载过程中，这2条序列互补的单链分别掉到不同的纳米孔中，则会形成2个一模一样的Cluster，从而产生一对Duplication reads。

华大智造MGI测序平台得到PCR双链文库后并不直接变性为单链，而是在后续环化实验中只环化双链中的1条，另外一条则被核酸外切酶消化，不用于形成DNA Nanoball。因此，在理论上，MGI平台的DUP rate要远低于其他平台。

此外，由于部分其他平台全部采用NanoWell结构＋排他性扩增技术（Exclusion Amplification），双链PCR产物变性解链后直接加入排他性扩增试剂Mix（DNA聚合酶、dNTP、Mg²⁺、Capture Primer… …）。该Mix中含有可以跟单链分子接头互补的Capture Primer（3’端带有双脱氧BLOCK），已经被BLOCK的单链模板即使已经加载到Flowcell并进入Nanowell中，也不能跟孔里面的草坪接头（P5、P7）互补进行Bridge PCR（互补序列被BLOCK），这种设计使得如果有两条单链模板掉到同一个纳米孔里也不会形成两个Cluster导致Polyclone，这种专利技术只会让每个孔中有且仅有一条的单链分子的BLOCK的3'端恢复为-OH状态并快速在聚合酶的作用下形成互补双链，随后又在解旋酶和单链结合蛋白的作用下解开双链，完全恢复自由状态的2条单链跟草坪接头快速结合进行Bridge PCR形成Cluster占满整个孔。但双链解链后产生了两条单链，如果其中一条漂到了其他的孔里，则又会形成Duplication reads增加Dup rate。^【4，5】

第四种Dup——测序过程中流体可能给不同类型芯片带来的光学 Dup

在信号放大过程中，无论是指数扩增而来的线性分子还是线性滚环扩增而来的DNA Nanoball，待测分子在测序芯片上加载和固定的过程中都是靠流体来保证芯片表面利用率的，高芯片利用率是数据高产出的基础，此乃流体的作用之一。另外，当前市面上应用最为广泛的待测分子结合芯片的方式有两类，一类是通过特殊蛋白质介导分子（如DNB）结合到“平坦的”芯片表面上，另外一类则是把单链分子退火到底部附着有特异序列的纳米孔内。这两类方式的主要区别源自于同属Patterned Array芯片的芯片表面是否采用凹陷的纳米孔，流体保证了文库分子能够最大限度铺遍整个芯片之外，其作用之二就是将测序生化反应所需的酶、dNTP等物料“运输”到每个信号点进行生化反应，这也不可避免的导致了不同芯片类型之间、同芯片上不同信号点之间生化反应速率的差异，流体在“一马平川”的芯片表面更新酶和dNTP等物料的效率理论上大于经过一个个小孔“减速带”并将孔内的物料更新的效率，因此纳米孔中的酶反应有更大几率的出现反应不充分的情况，这也就可能导致反应不充分的信号点因为信号强度显著弱于反应充分的“邻居”，从而被映射成两个孔表达出一样的信号，也就是一种光学上的Dup。

4、实际数据测评

基于上述方法学差异，可知：

1、基于矩阵芯片的NGS平台的Dup rate要明显低于非矩阵芯片的；

2、规则光信号点大小的NGS平台的Dup rate要低于不规则平台的。

咳咳咳，敲黑板，划重点，基于矩阵芯片和规则光信号大小的BGISEQ-500和MGISEQ-2000平台的Dup rate就稳定的在个位数波动，而其他平台文库的Dup则高达20%（图9），而基于其现成文库再扩增几个循环后使用BGISEQ平台测序，Dup rate则获得了显著的降低（见表1）。因此即使同一个文库，分别使用基于不同技术原理的测序平台进行测序，得到的Dup rate也有很大的差别。

图9. 270份BGISEQ-500平台商业数据与150份其他平台商业数据的Duplicate数据统计比较

那没有没有办法降低测序过程中的Dup呢？当然是可以的，由于此类Dup特点是在相邻位置出现，所以在图片转换为碱基信号时可以通过去除相邻碱基中同样信号的DNA分子来除此类Dup，但这导致另一个问题——芯片数据产量降低，例如去除此类原因导致的15% 重复序列，数据产量就降低了15%，这种情况尚需进一步讨论。

5、总结

综上，第一种样本本身的Dup出现几率最低，但也有可能给生信分析带来困扰；第二种PCR Dup在扩增循环数目不高的情况下并不会带来“令人窒息”的Dup rate，而只是在个位数水平波动。但是从文章前面两个平台上百个商业样本的测试数据（如图9）来看，Dup rate的差距还是很明显的（MGI平台<5% ，其他平台>20%），这个差距应该主要来自于我们前面提到的第三类测序前信号放大带来的Dup（线性扩增DNB VS 指数扩增Cluster）以及第四类不同芯片上的光学Dup了。

纸上得来终觉浅，绝知此事要躬行。《NGS攻略》之Dup传到此可以先暂告一个章节，但究竟如何，大家可以拿实际数据测评练练手再说，扫描下面二维码即可获得下载链接。欢迎大家分享实测结果或不同看法呦！

WGS数据—MGISEQ-2000_PE100

Samples: (Genome in a Bottle) Human DNA-NA12878

参考文献：

【1】Eric Vallabh Minikel. How PCR duplicates arise in next-generation sequencing[Z].2012,12.

【2】illumina. Effects of Patterned and Nonpatterned Flow Cells [Z].

【3】Steven Wingett. Illumina Patterned Flow Cells Generate Duplicated Sequences[Z].2017,03.

【4】Qiaoling Li, Xia Zhao, Wenwei Zhang, etc. Reliable Multiplex Sequencing with Rare Index Mis-Assignment on DNB-Based NGS Platform[J]. Biorxiv, 2018,06.

【5】Rahul Sinha, Geoff Stanley, Gunsagar Singh Gulati, etc. Index Switching Causes “Spreading-Of-Signal” Among Multiplexed Samples In Illumina HiSeq 4000 DNA Sequencing[J]. Biorxiv, 2017,04.

注1：部分图片来自网络，侵删。

注2：文章仅代表作者个人观点。

END

往期精彩内容推荐