PCR duplicates in NGS - I

在做NGS data處理的時候，常常會遇到需要remove PCR duplicates的步驟。

在網上查了一下 PCR duplicates的來歷，有兩個資料講的比較詳細。

但是這兩個資料側重的點不太一樣，這裡節選翻譯+總結一下這兩份資料，供將來研究作參考。

這裡是第一份資料，另一份資料在下一篇文章中。

原文：How PCR duplicates arise in next-generation sequencing

Dec 11, 2012 by ericminikel.

http://www.cureffi.org/2012/12/11/how-pcr-duplicates-arise-in-next-generation-sequencing/How PCR duplicates arise in next-generation sequencinghttp://www.cureffi.org/2012/12/11/how-pcr-duplicates-arise-in-next-generation-sequencing/?

www.cureffi.org

一般來說，建庫測序的過程分為下面六步：

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")

由於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")

從這個圖上來看，如果一個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")

可以看出，相比較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")

從這幅圖中可以看的更明顯。

到這裡，我們可以看到：DNA片段/bead數目 這個比例越小，PCR duplicate rate越大。

原文作者下面還回答了一些讀者提出的問題：

「為什麼不準備一些DNA原材料，不做PCR，直接測序」

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

在第二步加接頭的時候，只有一小部分DNA片段可以成功加上接頭，第三步會把成功連上接頭的DNA片段進行擴增，此時待上機的產物大部分都是有接頭的DNA片段了，而上機後，只有有接頭的DNA片段才可以被bead抓取到。

如果不做PCR，直接從第二步到上機，那上機產物中大部分都是沒有接頭的分子，上機後bead也無法抓取這些分子。