二代测序中Duplication Rate杂谈 (一)
2018/4/27 14:44:01 功能医学与荷尔蒙
什么是Duplication？
DNA分子多样性是影响Duplication rate的本质因素。减少PCR循环是降低Duplication rate的有效措施。由于随机事件的存在，二代测序各个环节都存在影响Duplication rate的因素。本文将深入探讨Duplication相关内容。
Duplication rate的计算公式如下：
Duplication Rate=1 - Unique Reads / Total Reads
本文中，相关的概念解释说明如下：
相关概念
Reads：二代测序平台对芯片上的Cluster连续顺序读取75、150或300个碱基信息的过程。Reads除记录碱基信息外，还包含每个碱基的质量标签及其它附加内容。
Cluster：在二代测序芯片表面或单个磁珠表面，由单个DNA模板生成的数百至数千个DNA分子的集合，类似单个细菌在LB培养基表面生成的单克隆菌落。
Duplication Reads：多个完全相同的DNA片段形成多个Cluster，产生的Reads信息也完全相同，称为Duplication Reads。
完全相同的DNA片段：指这些片段在参考基因组中的起始、终止位置以及序列信息都一样，被称为“三一样”。
“三一样”片段来源
如果没有完全相同的DNA片段，也就不会出现Duplication Rate，“三一样”片段源于以下几个地方。
01
Taq酶类介导的PCR
PCR过程本身就是对DNA分子的复制，是完全相同的DNA片段的主要来源。在靶向测序过程中，无论是建库时还是捕获后都不可避免地有PCR过程。
02
Phi29，Bst酶类介导的等温扩增
滚环扩增或链置换扩增可较为均一地扩增全基因组。某品牌单细胞DNA扩增就是用的等温扩增。等温扩增的产物均为结构复杂的大分子，进入NGS流程时仍需要被打断。由于等温扩增效率较高，可省去后续的PCR过程。
03
片段化过程
无论用超声波打断等物理手段还是酶切等生物手段，都有可能获得“三一样”片段。
片段化过程与PCR过程产生的“三一样”片段有着本质的区别，前者源于同源染色体，后者源于同一染色体上的同一条单链分子。分子标签可以给“三一样”片段加上“条形码(Barcode)”，甚至能利用Duplex原理给“三一样”片段的正负链加上不同的分子标签，人为产生不一样的片段(详情请参考NGS接头暗语)。
04
基因组的重复序列
基因组中的SINE、LINE、Alu等重复序列含量不低，它们在片段化过程中很容易形成“三一样”片段。
05
同一个片段的正负链
建库时如果未利用“成对的”分子标签对文库分子正负链进行标记，正负链就是一对天然的“三一样”片段。
影响Duplication rate的因素
样本来源、样本片段化、文库构建以及PCR扩增都会影响二代测序结果中的Duplication rate。下面我们将对引起Duplication的因素进行深入探讨。
样本来源
01
FFPE
FFPE样本由于保存条件的特殊性，存在交联、胞嘧啶脱氨成尿嘧啶、物理损伤等问题，其DNA的扩增较其它样本(血液、组织、体液，常温条件或者冷冻条件)更困难)。这意味着在起始量(Input)和连接效率相同的情况下，FFPE样本需要发生更多的PCR循环，这也导致了其Duplication rate的升高。
市面上有修复FFPE DNA的产品。推荐FFPE样本在进入NGS流程前进行质检，常见的参数有DIN(DNA Integrity Number)和Q-score(Quality score)。
02
cfDNA / ctDNA
cfDNA的断裂位置据说具有一定的模式，意味着断裂不是随机的，具有偏向性。自然产生的cfDNA与人工随机打断形成的cfDNA标准品差别较大，自然cfDNA的分子多样性可能较差，但连接效率更胜一筹。人工cfDNA由于存在物理损伤，即使进行末端修复也会影响其连接效率(具体可以参考Swift公司的一些说明[1])。需要注意的是，这里强调的是“酶切损伤”比“物理损伤”更容易修复，并不是比较Blunt-end或者Nick-end的连接效率。
另外，自然cfDNA集中在165bp左右，可能较为集中的片段分布比物理打断的片段分布更不容易丢失片段，自然cfDNA的分子多样性可能稍差但连接效率较高，而且起始量较低，理论上可能引起Duplication rate进一步升高。这时便体现出cfDNA建库好搭档“分子标签”的优势。建库过程中分子标签对文库分子“单链”或“姐妹链”的标记使得在后续数据分析中，相同的分子标签形成“家族”(Family)，对其进一步分析，Duplication reads能作为被分析的对象，共同判定稀有突变的存在。
03
甲基化DNA
经过亚硫酸氢盐处理的DNA，可以认为碱基类型都少了一种，分子多样性不但下降，更是引入了尿嘧啶，需要用尿嘧啶耐受的扩增酶，外加一些建库方式有着明显的GC偏好性，最后就是x-Ten平台测出来的Q30可能在75%以下，经过重重险阻之后，Pair-end mapping rate低得要靠Single end mapping rate去撑门面，Duplication rate已经不算什么大事了。
有经验表明使用Illumina EpiGnome / TruSeq DNA Methylation kit和HiSeq sequencing，40M reads的duplication rate在10%左右。Hiseq3000/4000在2x75bp模式下Duplication rate在20~25%之间。DNA甲基化测序质量可能受测序平台，建库试剂盒、读长，甚至和“spike in”的内容等多方面的影响。
04
RNA
RNA可以认为是全外显子组，只占全基因组的2%不到，少了这么多内含子等非编码区域的搀和，以及有些基因高表达，有些基因不表达，可想而知，分子多样性肯定就弱了很多。有研究表明，RNA-seq的Duplication rate在40~50%之间，能发文章的数据尚且如此，意味着什么情况都可能发生。
05
体液样本DNA
比如唾液样本DNA，个体差异很大，值得注意的是凋亡的情况，简单的凝胶电泳可以对凋亡的程度进行简单的评估，凋亡严重的，不仅影响Duplication rate，据说还影响数据分析，值得注意。
06
PCR产物
主要指多重PCR产物。值得思考的是针对cfDNA/ctDNA的多重PCR产物，建议在适当的引物位置引入分子标签，万分之一的灵敏度不是梦。
07
单细胞
单细胞DNA的有些区域在测序结果中压根就测不到。单细胞全基因测序的覆盖率是80%~90%，而多细胞全基因测序的覆盖率能达到99%以上，原因就是分子多样性不同。
对于多细胞测序，相同区域的分子，这个细胞没有扩增到，另一个细胞可能就扩增出来了；而单细胞只有两次机会，没了就是没了。
08
植物DNA / 微生物DNA
微生物16s区域测序，建议搭配分子标签使用，并且使用长读长测序模式。其他更深的理解有待经验积累。
片段化
01
片段化方法
样本之后，就是片段化，无论是超声波打断、红宝石孔拉断、高压气体喷断，还是酶切切断，这里强调的是随机打断，利用内切酶混合物进行酶切一定要测试其分子多样性即Duplication rate是否符合自己的要求。
02
片段化参数
片段化环节需要保证获得适当的长度，片段长度越小，导致扩增越容易，加剧了PCR Bias，最后引起PCR产物复杂度降低，Duplication rate升高。
末端修复连接
01
修复连接效率
对末端修复连接的效率的考量应该根据样本类型来考虑，比如ctDNA，单细胞样本，连接效率越高越好，不然低频的目标片段就会消失。某种程度上，连接效率越高，分子多样性越好，Duplication rate也就越低。
02
Adapter接头
当Adapter只有一种时，无论如何连接，你需要的分子数等于两头都接上Adapter的分子数，比如Illumina平台，但有些产品是两种Adapter，这种情况下，所有两头都接上Adapter的分子中，只有50%是你要的，这就需要根据具体需求进行考量。
PCR扩增
01
6个PCR循环!6个PCR循环!!6个PCR循环!!
无论是捕获前的文库扩增，还是捕获后的文库扩增，两者相加的PCR循环数能控制多低就控制多低。
02
NGS在操作的各个环节，都有各种随机事件发生，我们这里尝试套用泊松分布的公式，P= λke-λ/k! (k=0,1,2,3....n), 来计算PCR Duplication rate，其中参数λ是单位时间(或单位面积)内随机事件的平均发生率。
假设：NGS建库中，模板DNA打断、末端连接、加A和PCR扩增之后，获得500 ng PCR产物，这个产物在2100上的峰值为500 bp，并且知道1 ng 500 bp dsDNA约为3.29 fmol，即1.98×109个分子，那么可以推算出PCR循环数为n时需要的起始模板及其分子数如下：
其中，PCR起始量(ng) = 500ng/2n unique分子数=PCR起始量×1.98×109，这是为了方便推算而认为这些分子数都是unique的。
我们知道xTen每条lane可以产出375M reads；假设，λ=375M/unique分子数，此时λ的意义是单个reads在单位时间内碰到unique分子的平均发生率(这里设想的是碰撞后就分开，不占位，随机碰撞的过程)，那么不同λ值获得的P值如下表：
我们可以看到k>=2时，PCR循环数越多P值越大；其中k=0表示空值Blank，k=1时表示单个reads在单位时间内只碰到一个片段，即unique/monoclonal，k>=2时，表示多克隆(polyclonal)。
然后我们在R中，用以下公式来计算PCR Duplication rate：
sum(dpois(x1,lambda)/x1)/(1-dpois(0,lambda))，x1=seq(1,10,1)
Duplication rate=1-sum
结果如下：
这样，我们就推算出了PCR过程中的Duplication rate，6个循环的时候，Dup rate = 0.60524%，12个循环的时候，Duplication rate = 34.38%。很明显，为了获得500ng的PCR产物，如果起始分子数目越低，会需要更多的循环，这个时候就会带来更大的Duplication rate。
03
我们用不同起始量推算了不同PCR循环数时的Duplication rate，但实际上大家更关心相同起始量但不同PCR循环下的情况。考虑到PCR bias，这个过程很难用数学公式推算，那么我们设想，50%GC含量的片段扩增了1024倍的时候，而80%GC含量的片段只扩增了64倍，可想而知，最后50%GC含量的片段更有可能或者更多的进入到下一个环节，80%GC含量的片段甚至在纯化等各种步骤之后就丢失了，因此，PCR循环越多，扩增困难和扩增容易的片段之间相差就会越大。
有文章评估了不同GC含量片段在PCR前后的相对丰度[2]，如下图所示，可以看到PCR之后，GC或者AT含量高的片段丰度降低了10~100倍，毫无疑问，这些丰度低的片段在后续的环节中的存在感堪忧。解决的办法可能是选择适当的扩增酶以及扩增体系，提高扩增的均一性以及降低PCR循环数。
图1. 不同GC含量的片段在PCR前后丰度变化[2]。
04
在对PCR Duplication rate的推算中，我们假设了参加PCR过程的起始分子都是unique的，那么实际上这部分起始分子有多少比例是unique的？我们也尝试来推算一下。
我们假设: 一个300bp片段，当第300个位置的G必须要包括时，最多只有300种断裂的可能性；如下图所示：
图2. 300 bp片段断裂可能性示意图。
推而言之，需要300套基因组才能让这300种断裂方式同时存在。
再推而言之，人的hg19基因组有3137161264个碱基(以下计算按照单倍型来考虑)，当按照均一长度断裂，并穷尽各种断裂方式时，其片段数是3.14×109个，比如: 有300种方式断裂成300 bp的片段，总的断裂分子数=(3.14×109/300)×300=3.14×109，因此在300 ~ 500 bp范围(Size selection的范围)，Unique的断裂分子总数是201×3.14×109=6.31×1011，需要80400个基因组拷贝或者40200个细胞才能穷尽所有可能性，一个HeLa细胞基因组DNA=15 pg，这么多细胞DNA的质量是40200×15 pg=603 ng。
我们继续假设，当随机打断的长度为400 bp，在不同起始量DNA条件下会得到不同的断裂分子数，如果将6.31×1011视为总Unique片段数，那么随机打断过程就像在这么多unique片段中进行随机抽取，因此，我们依然套用泊松分布的公式，λ=打断后的片段数/总的Unique片段数(我们没有办法获得准确的打断后的弥散产物中的分子数，只能用400 bp进行估算)
不同Input质量(200 ng~1000 ng)DNA获得400 bp片段的分子数如下所示：
再根据λ=打断后的片段数/总的unique片段数(6.31×1011)，获得P值：
可以发现，k>=2的几率随着input质量的提高而升高，我们推算打断后的Duplication rate。
*** 这里提示我们随机打断的片段中，有较大部分是“三一样”片段。
经过纯化、末端修复，加A，连接，再纯化等5个步骤，假设每个步骤的得率都是60%，则最后只有7.776%(这是一个经验值)留下：
此时，λ=5步骤后片段数/总的Unique片段数(6.31×1011)，获得P值
此时，也就是PCR前，起始分子中存在的Duplication rate为：
***这里提示我们，经过5个步骤的损失，PCR之前获得DNA分子数可能不到起始DNA的10%，但依然存在着“三一样”片段。
那么我们现在有了上表中的Duplicatioin rate，就可以计算出5步骤后片段中含有的Unique片段数=5步骤后片段数×Duplication rate：
***这里提示我们，200 ~ 1000 ng的起始DNA，最后能获得1010~1011级别的Unique片段数，这些片段将进入到PCR过程，由于PCR过程具有偏好性，并不是一个很好的随机体系，暂时不能进行PCR过程的推算。
Cluster
01
Cluster 的生成
Cluster在Flowcell上的生成也是一个PCR过程。这个PCR比较容易被人遗忘。Cluster变少会影响Duplication rate。原因是比例少的分子可能不能产生cluster，唯一性分子数减少，进而影响Duplication rate。
适当的Cluster生成密度，不仅能够获得最佳的数据产量，也能够获得较低的Duplication rate。无论Illumina平台还是PGM平台，我们都希望Cluster是单克隆(Monoclonal)的，多克隆(Polyclonal)的Cluster甚至是因空间距离过近而导致图像识别时相互Overlap的Cluster都会被测序识别程序过滤掉，造成的直接影响就是Cluster密度过高，数据产量降低，整张芯片的Cluster多样性降低，造成Duplication rate升高。
02
Cluster模版分子的多样性
Cluster生成过程的最开始，是DNA模板与芯片上引物结合的过程，这个过程是随机事件，各种分子的比例越均一，模板分子的多样性和复杂度越高，Duplication rate就会越低，极端的情况就是每个Cluster只是一种模板分子的单克隆，这个时候Duplication rate是最理想的。我们尝试推算Cluster的Duplication rate。
xTen上机每条lane是3 nM样本取5 μL，即15 fmol=9.03×109个分子，xTen每条lane能产生375M reads，相当于在9.03×109个分子中随机抓取375M个，也即模板分子与引物结合的过程是一个随机碰撞的过程，我们利用泊松分布公式进行推算，将9.03×109个分子看成是unique的，λ=3.75×108/9.03×109=0.0417。
进而，Sum=0.9896，Duplication rate=1-0.9896=1.04%，1.04%只是表示碰撞的重复几率，假设模板引物与模板碰撞后不结合。
Illumina测序平台
图3. Illumina平台引起Duplication的因素[3]。
01
光学分辨
光学分辨引起的Duplication，如上图左上，一个不规则的Cluster被识别成了两个，这种情况会发生在Hiseq 2500/Miseq/Nextseq数据中。被识别成两个的Cluster通常比较大，或者不规则，或者是由于生成Cluster的模板DNA分子被释放后在很近的位置又被再次捕获，就近形成第二个克隆[3]。
02
排他性扩增
Exclusion amplification duplicates (ExAmp，排他性扩增) 。Hiseq X或者Hiseq 4000采用Patterned flowcell，它密布着人工的规则的Nanowells，完全不同于最初的Solexa的随机Cluster generation方式。当一个文库模板分子进入Nanowell时，这个Nanowell已经有另一个文库分子在进行极早期Cluster generation，那么还是可能会生成多克隆的Cluster，这个多克隆Cluster会被程序识别并去掉。
另外，还存在一种情况就是，一个文库模板分子进入Nanowell开始桥式扩增，PCR开始后这个文库模板分子从Nanowell中游离出来，重新进入到大的溶液体系中，并进入另一个Nanowell，生成了第二个Cluster，这种情况也就是上文分析的只是发生碰撞的状况，约有1.04%的Duplication rate。ExAmp duplicates是由于Cluster generation方式引起的，一定几率发生，减少的办法是给予适当浓度的模板分子并保证模板分子中的分子多样性尽可能的高。
03
Cluster
Cluster在Flowcell上的生成也是一个PCR过程。这个PCR比较容易被人遗忘。Cluster变少会影响Duplication rate。原因是比例少的分子可能不能产生cluster，唯一性分子数减少，进而影响Duplication rate。
适当的Cluster生成密度，不仅能够获得最佳的数据产量，也能够获得较低的Duplication rate。无论Illumina平台还是PGM平台，我们都希望Cluster是单克隆(Monoclonal)的，多克隆(Polyclonal)的Cluster甚至是因空间距离过近而导致图像识别时相互Overlap的Cluster都会被测序识别程序过滤掉，造成的直接影响就是Cluster密度过高，数据产量降低，整张芯片的Cluster多样性降低，造成Duplication rate升高。
Cluster生成过程的最开始，是DNA模板与芯片上引物结合的过程，这个过程是随机事件，各种分子的比例越均一，模板分子的多样性和复杂度越高，Duplication rate就会越低，极端的情况就是每个Cluster只是一种模板分子的单克隆，这个时候Duplication rate是最理想的。我们尝试推算Cluster的Duplication rate。
xTen上机每条lane是3 nM样本取5 μL，即15 fmol=9.03×109个分子，xTen每条lane能产生375M reads，相当于在9.03×109个分子中随机抓取375M个，也即模板分子与引物结合的过程是一个随机碰撞的过程，我们利用泊松分布公式进行推算，将9.03×109个分子看成是unique的，λ=3.75×108/9.03×109=0.0417。
进而，Sum=0.9896，Duplication rate=1-0.9896=1.04%，1.04%只是表示碰撞的重复几率，假设模板引物与模板碰撞后不结合。
04
测序错误
一般认为，测序错误在千分之一左右，意味着测1000个碱基，就有一个碱基被测错的，后续的分析步骤是有能力识别这种错误的，从某一个方面来说测序错误增加了reads的多样性，也算有点好处，但这种好处对ctDNA的测序来说可能就让情况变得更为复杂了。
其他因素
01
磁珠纯化
磁珠纯化可能会影响Duplication rate。磁珠吸附本质上是电荷吸附。DNA序列不同，电荷密度不同，洗涤条件如果不严谨，就会造成磁珠吸附有偏好性，也就意味着分子多样性降低，进而影响Duplication rate。
02
杂交过程
不考虑模板本身的多样性，如果杂交时重复序列(比如LINE、Alu序列)未能有效封闭，必然造成Duplication rate升高，有效数据降低。
有效封闭还包括Adapter的封闭。如果Adapter不能有效封闭，造成脱靶率升高。单位质量的DNA文库中，中靶率降低，势必造成中靶部分中比例低的分子减少或者缺失，容易扩增的分子占优势，即分子中靶部分的分子多样性减少，进一步影响Duplication rate。
03
试剂质量
比如某个批次的SBS测序试剂出了问题，WES的Duplication rate有可能从正常的<10%升高到30%。
04
探针因素
Duplication rate跟探针关系最小。极端的例子就是扩增子(PCR产物)测序，Duplication rate可以很高，但完全没有探针存在。探针杂交过程中最影响分子多样性的因素是探针分子与目标分子的比例，其次是杂交的时间。探针与目标分子的比例最低应该在100：1(一定的体积内，有体积的要求)，高于这个比例，探针是能够将目标分子捕获的，目前的探针法NGS中，这个比例是很高的，探针数量是高度冗余的。杂交法的目标是获得最多的分子包容性，获得最多的与参考序列不一样的序列。只有增加分子包容性，才能最大可能地包含变异的类型，因此探针分子要长，杂交时间要长。探针杂交不是为了特异性，而是为了更高的产量，获得更多的不准确。因为越准确，就意味着测序数据与参考序列(也就是探针序列)是一样的，这样就没有测序的必要了，也就没有发现新的信息的功能了。因此，杂交时间短，不能容纳最大序列可能性，造成分子多样性降低，进而影响Duplication rate。
05
建库的次数
目前用到深度测序的应用越来越多，最明显的就是ctDNA，动辄10000x的深度，一些全基因组测序也需要高深度测序，从正常的30x变成300x，当需要这种应用时，我们建议在条件允许的情况下，多建几次文库，比如全基因组300x测序，建一个文库测序300x和建3个文库各测100x，其分子的多样性是完全不一样，值得我们去尝试。
06
其他因素
还包括DNA提取质量、操作人员的操作习惯等。比如同样取200 ng基因组DNA做NGS，一个取0.1 μL获得200 ng，一个取5 μL获得200 ng，也影响分子多样性，后面这个例子可能也比较极端，比如提取DNA的时候，在乙醇沉淀的那一步，如果DNA量很大，我们甚至能看到白色絮状物，这个就是大片段的DNA，这种样本在用TE溶解之后，如果溶解不完全，我们用移液枪吸取的时候，甚至能看到粘稠的胶状物，这个时候如果你取0.1 μL，那很有可能什么DNA都没有取到，需要进一步溶解或者稀释后再吸取。
Duplication Rate经验值
WES：数据量约10G时，Duplication rate在10%左右；
WGS：数据量约90G时，Duplication rate在10%左右；
RNA-seq：Duplication rate在40% ~ 50%左右；
WGBS：>10G, Duplication rate > 10%。
多重PCR试剂盒和靶向捕获Panel，差异很大，跟需要测序的区域以及测序量有关，通常情况下只要中靶部分数据质量足够好，Duplication rate不是一个重要的考虑指标。
结论
二代测序过程中的Duplication rate受很多因素的影响，其本质是(各环节)起始时的分子多样性。因此对Duplication rate的质控是一个系统工程，体现了一个实验室的管理水平与能力。
降低Duplication rate的简单方法就是控制PCR循环数，如果在合适的PCR循环数下得不到你需要的DNA量，建议向前追溯原因，而不是简单的增加PCR循环。
免责申明：涉及公式和数字的部分不保证正确，只想作为一种提示，欢迎给予指导。
来源：那昂达科技
参考文献
[1] https://swiftbiosci.com/applications/circulating-cell-free-dna/#exp-performance-data
[2]Aird D, Ross M G, Chen W S, et al. Analyzing and minimizing PCR amplification bias in Illumina sequencing libraries[J]. Genome Biology, 2011, 12(2):1-14.
[3] http://core-genomics.blogspot.hk/2016/05/increased-read-duplication-on-patterned.html
http://weixin.100md.com
返回 功能医学与荷尔蒙 返回首页 返回百拇医药