前言

PCR duplicate在NGS文库中是再常见不过的名词，对于需要自己建库的小伙伴来说可能会被折磨的“痛不欲生”，这篇帖子，我希望能用自己弱的不能再弱的数学知识加上一些Google到的内容来分析一下为什么我们的文库总是充满PCR duplicate？

二项分布和泊松分布

我们以ChIP-seq为例，来简单介绍一下NGS建库过程中需要用到的简单数学概率分布：二项分布和泊松分布。

• 二项分布

实际上，我们送测的文库中只有一部分会被测序，也就是说，一个有个DNA fragments的文库最终只会被测个DNA fragments，这里。那么这个被测序的DNA fragments中到底会有多少是来自于我们的某个位点附近呢（记作），显然这和原始文库中位点附近的DNA fragments所占的比例有关，所以这个将会服从以下二项分布。

这里表示非负整数，也就是能够测到多少个位点附近的DNA fragments。为什么可以这样做模拟呢？实际上二项分布的最原始情形是这样的：n次独立重复实验当中，事件m出现的次数。对应到我们的测序事件当中：从个DNA中抽取个进行测序，其中有多少个DNA来自于位点附近。把这个过程简化一下，我们每次从个DNA中抽1个，重复这个过程次，由于我们的，而且本身也很小，所以这个不放回抽样的过程我们可以近似理解为放回抽样。简单画个图理解一下这个过程，当然这个例子和我们测序的例子是不符的，只是帮助大家理解什么是二项分布而已：

x <- 1:30
plot(x, dbinom(x, size = 100, prob = 0.1))

• 泊松分布

先给结论：当n很大，p很小的时候，二项分布可以近似地用泊松分布，，来进行表示，提供简单证明：

进一步的：

所以：

上下同时约掉：

这里：

根据极限：



所以我们最终就证明了：

这正好符合我们NGS测序的特点：高通量大数据（很大，而一般来说每个位点的DNA fragments占比也很小）。

一样的，简单绘个图：

plot(x, dpois(x, lambda = 10))

data <- data.frame(x = 1:30, 
                   y1 = dbinom(x, size = 10000, prob = 0.001),
                   y2 = dpois(x, lambda = 10))
ggplot(data = data) + 
  geom_point(aes(x = x, y = y1), color = 'red', size = 2, alpha = .5, shape = 6) +
  geom_point(aes(x = x, y = y2), color = 'blue', size = 2, alpha = .5, shape = 1) + 
  ylab(label = 'y') +
  theme_classic()

可以看到几乎完全重叠，足以证明我们的这种“近似”是靠谱的。

为什么建库要合适的PCR轮数？

现在测序公司要我们送样的时候都会要求一个最低浓度和最低体积，假如一个公司要求你最少送样15ng，你的理论最适PCR循环数为6轮，那么显然你的起始量DNA为：。注意这个0.234375ng里面的每个分子都是独特的，也就是大家经常说的unique。

同时，DNA分子的平均质量为660g/mol/bp，所以15ng的长300bp的文库里面大概有个DNA，原始的“独特”DNA大概有个，如果你要产生6G的测序数据（PE150），就相当于你测了其中的个DNA。对应到前面的泊松分布：，，所以。所以原始DNA中一个“独特”的DNA分子被测到的次数所服从的泊松分布就为：

data <- data.frame(x = 1:10,
                   y = dpois(1:10, lambda = 0.02807018))
ggplot(data = data) + 
  geom_bar(aes(x = as.factor(x), y = y), stat = 'identity') + 
  xlab(label = 'Sequencing time of an unique DNA') + 
  ylab(label = 'Frequency') + 
  theme_classic()

lambda=0.028

你可能会注意到，不管我PCR多少个循环，似乎最终测6G的数据的话，为一个定值，难道PCR循环与PCR duplicate无关？实际上不是这样，现在很多测序公司的测序都是包G的，也就是说你只需要告诉他你需要多少G数据的产出即可。但一个测序的lane能产出很多数据，所以一般公司会同时在一个lane上测很多用户的数据，这个时候就需要控制每个用户的上样量，防止由于某个用户的上样量太多，竞争性结合更多的cluster，导致其他用户的测序量不够。

什么意思呢？举个例子，现在一样是满足公司要求，15ng的DNA送测，最理想的手段是前面所讲的0.234375ng的DNA扩增6轮，但是如果你不知道PCR duplicate的概念，你扩增了8个循环，那么你最终将得到DNA，而为了保证用户之间的上样量平衡，最终公司仍然只会取其中的15ng去测序，这也就相当于只取了你起始0.234375ng中的ng去测序，显然这个时候，，再去看泊松分布：

data <- data.frame(x = 1:10,
                   y = dpois(1:10, lambda = 0.112))
ggplot(data = data) + 
  geom_bar(aes(x = as.factor(x), y = y), stat = 'identity') + 
  xlab(label = 'Sequencing time of an unique DNA') + 
  ylab(label = 'Frequency') + 
  theme_classic()

lambda=0.112

什么情况下会产生PCR duplicate？

根据前面的信息，能够增加PCR duplicate的情况最终都是导致的增大：

• 本身建库起始DNA量很少，为了满足送测要求而PCR了很多轮；

• 盲目追求多的DNA扩增后的量，对起始DNA进行过度的PCR，即使早就达到送测量的要求了，还在进行额外的PCR。

如何确定合适的PCR轮数？

简单来说，目前很多人借鉴的都是ATAC-seq建库的PCR扩增轮数确定方法，这是一个基于经验的方法，至少我目前不知道怎么去解释它……大概的思路如下：

• 按照一定的引物浓度，聚合酶浓度配制50uL PCR体系；

• 对上述体系进行PCR扩增，到第5个循环的时候暂停；

• 取出5uL，配成最终体积为20uL的qPCR新体系，需要注意的是，你需要补入额外的、一定量的引物和聚合酶，直至反应体系中这两个组分的浓度和起始的50uL体系一样，相当于让DNA在相同的环境下进行扩增，这才是有意义的，所以你的酶也必须要是一样的，qPCR就单独加SYBR，而不要用市场上一般的qPCR的mix了，因为你们的酶都不一样了，没有意义；

• 进行qCPR，观测荧光信号值与扩增轮数之间的关系，取达到最高荧光信号值1/4~1/3处对应的PCR循环数作为剩余的45uL反应体系还需要扩增的循环数对其进行进一步的PCR，最终纯化产物即可。

需要详细protocol的筒子们可以在评论区留言~

最后祝愿自己和大家，实验顺利，数据多多！！！