测序的PCR duplicates - I（转载）

在做NGS data处理的时候，常常会遇到需要使用picard和samtools 进行remove PCR duplicates的步骤。
在网上查了一下 PCR duplicates的来历，有两个资料讲的比较详细。
但是这两个资料侧重的点不太一样，这里节选翻译+总结一下这两份资料，供将来研究作参考。
这里是第一份资料，另一份资料在下一篇文章中。

第一篇

一般来说，建库测序的过程分为下面六步：

1， 用超声波等方法将基因组DNA打碎。
2，对打碎得到的DNA片段加接头。
3，PCR扩增加了接头的DNA片段。
4，将扩增后的产物“洒”在flowcell上，目的是一个bead上抓取到一个DNA片段分子。
5，在bead上进行PCR扩增反应，保证后面测序时可以读取到足够强的信号。
6，测序（各测序平台有自己的技术：pyrosequencing (Roche), reversible terminator chemistry (Illumina), or ion semiconductor (Ion Torrent)).

PCR duplicates是在第3步产生的。

理想情况下，对打碎的基因组DNA，每个DNA片段测且仅测到一次。

而第三步扩增了6个cycle，那么每个DNA片段有了64份拷贝。将扩增后所有产物“洒”到flowcell，来自一个DNA片段的两个拷贝，可能会锚定在两个bead上，经过测序得到的这两条read，就是PCR duplication。

一般来说，如果PCR duplication rate过高，那么同样总数目的reads，所提供的关于基因组的信息就大大减少了。

那么PCR duplication rate跟什么因素相关呢？

下面通过一个简单的计算来说明这个问题。

假设扩增后（第3步后）准备上机测序的DNA总质量为1ug。建库时候扩增了6个cycle，也就是每个DNA片段有2^6=64份拷贝。原始的（扩增前）DNA片段总质量=1ug/64=15.6ng.

一个碱基对的分子量为 660g/mol, 即 660 g per mol per bp. 假设我们打碎的DNA片段长度都为200bp，那么一个片段的分子量 = 200bp*660g/mol = 220 bp * 660 * 10^9 ng/mol

那么 原始的（扩增前）DNA片段总数目=15.6 ng / (220 bp * 660 * 10^9 ng/mol ) = 1.07438e-13 mol

由于1mol 单位中包含 6.022×10^23 个粒子，那么1.07438e-13 mol对应 1.07438e-13 *6.022×10^23 = 64699163600 ~ 7e10 个分子。

也就是说，基因组被打碎后得到了** 7e10 **独特的DNA片段。

注：阿伏伽德罗常数用于代表一摩尔物质所含的基本单元（如分子或原子）之数量，在一般计算时，常取6.02×10^{23或6.022×10}23为近似值（来源，维基百科：https://zh.wikipedia.org/wiki/阿伏伽德罗常数）

如果我们对人的基因组进行测序（3G bp），测序深度为20x，read长度为100bp，那么对应了大约 3e9 bp /100bp * 20x = 600M reads. 一个read来自一个bead，也就是有600M个beads。

到这里，我们有7e10 独特的DNA片段，每个片段在上机测序前有64份拷贝，它们会 600M beads随机进行杂交。

平均来看，一个DNA片段被一个bead“抓住”的几率 = 600M/7e10=6e8/7e10=0.0085.

为了了解PCR duplicates rate，我们需要知道这7e10个独特的DNA片段，被0个、1个、2个... bead抓住的几率，体现在最后read中就是0个read，1个read，2个reads...

由于一个DNA片段被一个bead“抓住”的几率很低，我们可以用Poission distribution来刻画这个过程。

Poisson distribution的概率密度函数 = λke-λ/k! , k=0,1,2,...n.

这里， λ 就是一个DNA片段被一个bead“抓住”的几率=0.0085. 而 k 就是一个独特的DNA片段被几个bead抓住。

k=seq(0,10,1)
names=as.character(k)
barplot(dpois(k,lambda=0.0085),names.arg=names,xlab="Number of times a given molecule is represented in reads",
ylab="Probability")

image

由于beads的数目远远小于DNA片段总数，因此大部分DNA片段根本木有被测到。

而我们关注的是“1” bar和它右边的所有bar。因此我们把“0”bar去除掉，再看：

k=seq(1,10,1)
names=as.character(k)
barplot(dpois(k,lambda=0.0085),names.arg=names,xlab="Number of times a given molecule is represented in reads",
ylab="Probability")

image

从这个图上来看，如果一个DNA片段被测到了（因为我们把k=0的都扔掉了），那么它很大机会是只被测到了一次。而如果它被测到了2次或更多次，这就造成了PCR duplicates.

我们算一下被测到的DNA片段（无论被测到几次）的比例:

sum(dpois(k,lambda=.0085)/k)/(1-dpois(0,lambda=.0085))

sum(dpois(k,lambda=.0085)/k) 就是在数 count(1)/1+count(2)/2 + ..., 也就是总计测到的、独特的DNA分子有多少个。被除数(1-dpois(0,lambda=.0085)) 是所有测出来的reads数目。

结果为0.9979， 即 0.21% PCR duplicates.

现在，我们假设，我们打碎基因组DNA后，得到了9e9个DNA片段，在第三步扩增时候，进行了9个cycle，也就是每个DNA片段被扩增了512倍。那么 Poisson distribution中的lambda=6e8/9e9 = 0.067. 即一个DNA片段被一个bead抓住的概率为0.067.

sum(dpois(k,lambda=.067)/k)/(1-dpois(0,lambda=.067))

结果为0.983，也就是有1.7% PCR duplicate rate。

k=seq(1,10,1)
names=as.character(k)
barplot(dpois(k,lambda=0.067),names.arg=names,xlab="Number of times a given molecule is represented in reads",
        ylab="Probability")

image

可以看出，相比较lambda=0.0085的柱状图，k=2的柱子明显变高了，也就是说更多的DNA片段被测到了两次，PCR duplicates rate上升了。

更极端的情况，如果我们打碎基因组DNA后，得到了1e9个DNA片段，在第三步扩增时候，进行了12个cycle，也就是每个DNA片段被扩增了4096倍。那么 Poisson distribution中的lambda=6e8/1e9 = 0.6. 即一个DNA片段被一个bead抓住的概率为0.6.

sum(dpois(k,lambda=.6)/k)/(1-dpois(0,lambda=.6))

=0.85， 即有15%的reads都是PCR duplicates。

k=seq(1,10,1)
names=as.character(k)
barplot(dpois(k,lambda=0.6),names.arg=names,xlab="Number of times a given molecule is represented in reads",
        ylab="Probability")

image

从这幅图中可以看的更明显。

到这里，我们可以看到：DNA片段/bead数目 这个比例越小，PCR duplicate rate越大。

原文作者下面还回答了一些读者提出的问题：

“为什么不准备一些DNA原材料，不做PCR，直接测序”

作者引用了Mike Talkowski （People） 的答案：

在第二步加接头的时候，只有一小部分DNA片段可以成功加上接头，第三步会把成功连上接头的DNA片段进行扩增，此时待上机的产物大部分都是有接头的DNA片段了，而上机后，只有有接头的DNA片段才可以被bead抓取到。

如果不做PCR，直接从第二步到上机，那上机产物中大部分都是没有接头的分子，上机后bead也无法抓取这些分子。