“PCR放大成百上千倍,为什么NGS的Dup rate只是个位数?


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

其实测序产生的Dup reads来源很广泛,主要有以下几个方面:

1、样本本身的Dup;

2、文库构建中扩增引入的Dup;

3、测序前信号放大(荧光信号采集单元生成过程)引入的Dup;

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

很多人认为NGS数据的Dup主要来自于上述第2种,其实大错特错,真实情况应该是3.荧光信号采集单元生成的过程中引入的Dup和4.芯片测序过程中引入的光学Dup占据了主要部分。而第1种样本本身的Dup恰恰是我们生物学中应该真实关注并需要和测序数据的Dup相区分的一种“有用的Dup”

补充说明1:荧光信号采集单元的定义我在上一篇文章中已经给出:序列信息采集时基本信号单元的载体。1个荧光信号采集单元的制备或者生成其实对应的是1个Cluster(illumina平台)或者1个DNA纳米球的生成(DNBSEQ平台),之所以需要把待测的单条DNA分子通过BridgePCR或者滚环扩增(Rolling circle amplification)的方式放大,主要是因为单单1个拷贝的DNA分子发出的荧光实在太微弱,相机无法捕捉这么微弱的荧光信号。

补充说明2:PacBio的SMRT测序技术和瀚海基因(现为GeneMind真迈生物)的CenoCare测序系统是也基于荧光标记的测序平台,为什么他们就可以仅使用1个拷贝的DNA分子(单分子)就实现了荧光信号的检出呢?

PacBio单分子检测的实现主要依托于其ZMW(零模波导孔),具体的原理大家参考我关于PacBio测序原理的文章。GenoCare实现单分子荧光信号检测的原理我在此不想展开讲(或许等它热度哪天再上来的时候,我为了蹭热度会写篇它测序原理和应用场景分析的文章),但可以概括为主要依托于TIFR(Total internal Reflection Fluorescence,全内反射荧光显微镜)来实现的。Pacbio测序原理及SMRT bell文库构建流程简述(二)

我个人认为,不管是ZMW还是TIFR,其在实现单分子识别中起到的核心作用是一致的——将激光的激发区域控制在一个很小的范围内,从而大大降低了测序背景噪音。

我们先说第1种Dup——样本本身的Dup:

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

补充说明3:一个人的基因组的重量大约3.3pg,虽然我们常常直接使用这个数字来进行一系列的计算,但可能还是有很多筒子们不知道3.3pg是怎么计算出来的,这里我们也来梳理一下3.3pg是如何计算出来的。

阿伏伽德罗常数是6.02×1023,4种脱氧核糖核苷一磷酸(dNMP)的平均分子质量约为330g,即6.02×1023个dNMP是330g,那么人的基因组有3×109(3Gb)个碱基对,对应的是多少g?

因为正常的人源细胞是二倍体,所以6.6pg/Cell。

补充说明4:既然上面公式得出是3.3pg/Copy,文章为何是按照6.6pg/Copy来计算,而不是3.3pg呢?这里我举例说明,例如共有100个人源细胞,理论上这100个细胞应该有100×2条1号染色体,但需要注意的是同一个细胞中2条同源染色体(1A和1B)之间碱基序列是存在差异的(即杂合),即使断点相同的1A和1B,测序所得的Reads也可能不是Duplicates Read,所以理论上,100个细胞中序列完全一致的染色体Copy数是100条,而不是100×2.

下面我们主要说一下第2种Dup,也就是大家平时讨论最多的PCR Dup。这里还是以人的基因组DNA文库制备为例,我们通常认为理论情况下PCR扩增6个循环已经把样本分子数量放大64倍(26=64),最起码也是20-30倍(如按照1.6-1.7的扩增效率计算)。这些PCR扩增产生了Dup reads理论上应该体现在最终的测序Reads中,约占(19/20)×100% ~(29/30) ×100%,这样算来应该是高达96%以上Dup rate?但我们日常实测数据来看并没有这么高,而且华大智造的DNBSEQ平台的Dup rate稳定的在个位数波动。

那么问题来了,“PCR放大成百上千倍,为什么NGS的Dup rate只有十位数甚至是个位数呢?”

这里我们采用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循环的扩增,要求最终得到约1ug PCR扩增产物,那这也就意味着我们需要1μg/(26)=15.6ng的DNA,假设待测分子片段长度平均为200bp,根据{(15.6×103pg)/(6.6pg/细胞基因组)}×2倍体×(30亿对碱基/基因组)≈14×1012bp的公式计算出PCR前样本本身的Unique Molecules(此处不考虑样本本身Dup)数量如下:

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

经过6轮PCR扩增后,7×1010个Unique分子,每一个都有64个Copy,即有6×108个孔,而有7×1010个Unique分子簇(每簇又有64个Copy)想公平的、完全随机的进入这些孔。

前方高能!忘记了什么是二项分布、泊松分布、正态分布、超几何分布… …的筒子们先去自行百度。

或者直接滑到文末看总结。

就此,这完全转化成了一个数学问题,一个概率(p)=1/(7×1010),n=6×108的二项分布模型,一个符合二项分布模型的实验应具有以下基本属性

1.需要进行n次实验(需进行6×108次,每次从7×1010个分子簇中随机选1个)

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

3.事件发生的概率p是固定的(每次实验,每一个分子簇被选中的概率为1/(7×1010%

4.每次实验的结果只有两种,只有成功选中(p)或者未选中(q=1-p)

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

其中λ为期望值,λ=np≈6×108/7×1010≈0.0085

根据以上参数,使用R来做7×1010个Unique分子簇中同一个分子簇被随机抽取0次、1次、2次… …n次而获得对应呈现机会的概率直方图。(获得呈现>2次的,即为Dup)

图1.使用R调用dpois(x,lambda)作图

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

图2.λ=0.0085时的呈现次数(0~10)概率直方图

显而易见的是7×1010个Unique分子簇绝大部分的都没有幸运的获得呈现的机会(即“0”)我们睁大眼睛能看到呈现1次的幸运儿,但为了更好的看看更幸运的“>2”的部分(即PCR Dup),我们将“0”值去掉再重新做直方图

图3.去除“0”值,重新计算呈现概率

图4.λ=0.0085呈现次数(1~10)的概率直方图

我们利用上述数据调用dpois(x,Lambda)来计算PCR Dup的结果,计算公式如下

(7×1010个Unique分子簇呈现1次的概率)/(7×1010个Unique分子簇呈现不小于1次的概率)约为99.79%,即PCR Dup rate为0.21%。(这与实际实验操作中大概4%的PCR Dup rate存在一定差距,主要原因是由于此算法忽略PCR Bias因素)
当我们将PCR循环数增加到9个时(PCR产物总量仍为1μg,起始样本的Unique Molecules数量降低为9×109个,每个分子拥有29个Copy,λ=np≈0.067)同样可以得到概率直方图(图5)
图5.λ=0.067的呈现次数(1~10)概率直方图

经运算,PCR Dup rate为1.7%,Eric还计算了当PCR循环数为15时,PCR Dup rate的理论值已经达到了15%(由于文章篇幅原因不再展开,有兴趣可查看原文),文章题目:How PCR duplicates arise in next-generation sequencing。

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

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


好了,上面讲的这些没看懂?这都不重要,我们来类比一个例子吧,我们把待测的7×1010个Unique 分子比作700道菜,而Flowcell上的6×108个Nanowell或者Spot位点数比作你能夹菜的次数。

如果假设

1.你对700道菜的喜爱程度一样

2.每道菜你夹到的难易程度一样

3.每道菜的初始菜量一样

但是这700道菜却只允许你夹6次,那么请问你这6次中夹到同一个菜两次以上(>2)的概率是多少呢?但如果逐渐减少菜品种类,增加夹菜次数,那么同一道菜夹到多次的概率又会怎么变化呢?


(0)

相关推荐