生物中的数学【二】为什么你的文库总是充满PCR duplicates？

前言

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

二项分布和泊松分布

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

  • 二项分布

实际上，我们送测的文库中只有一部分会被测序，也就是说，一个有M个DNA fragments的文库最终只会被测N个DNA fragments，这里N<<M。那么这N个被测序的DNA fragments中到底会有多少是来自于我们的某个位点i附近呢（记作n_i），显然这和原始文库中i位点附近的DNA fragments所占的比例p_i有关，所以这个n_i将会服从以下二项分布。
f(k,N,p_i) = Pr(k;N,p_i) = Pr(X=k) = {N \choose k}p_i^k(1-p_i)^{N-k}
这里k表示非负整数，也就是能够测到多少个位点i附近的DNA fragments。为什么可以这样做模拟呢？实际上二项分布的最原始情形是这样的：n次独立重复实验当中，事件m出现的次数。对应到我们的测序事件当中：从M个DNA中抽取N个进行测序，其中有多少个DNA来自于位点i附近。把这个过程简化一下，我们每次从M个DNA中抽1个，重复这个过程N次，由于我们的N<<M，而且p_i本身也很小，所以这个不放回抽样的过程我们可以近似理解为放回抽样。简单画个图理解一下这个过程，当然这个例子和我们测序的例子是不符的，只是帮助大家理解什么是二项分布而已：

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


理解一下这个过程，某个实验成功的概率为10%，然后我做了100次独立重复实验，其中成功的次数为x，x应该服从一个这样的分布，例如成功10次的概率为0.1318653。这可能和你所想象的不一样：难道概率为10%，我重复100次不就应该成功10次吗？实际上不是这样，这大概就是概率的魅力吧，只能告诉你，你成功10次的概率是最高的。

  • 泊松分布

先给结论：当n很大，p很小的时候，二项分布可以近似地用泊松分布f(k;\lambda) = Pr(X=k) = \frac {\lambda^k e^{-\lambda}}{k!}，\lambda=np，来进行表示，提供简单证明：
Pr(X=k) = {n \choose k}p^k(1-p)^{n-k} = {n \choose k}{(\frac {\lambda}{n})}^k(1-{\frac {\lambda}{n}})^{n-k}
进一步的：
{n \choose k} = \frac{n!}{(n-k)!k!}
所以：
{n \choose k}{\frac {\lambda}{n}}^k(1-{\frac {\lambda}{n}})^{n-k} = \frac{n!}{(n-k)!k!} {(\frac{\lambda}{n})}^k(1-\frac{\lambda}{n})^{n-k}
上下同时约掉(n-k)!：
\frac{n!}{(n-k)!k!} {(\frac{\lambda}{n})}^k(1-\frac{\lambda}{n})^{n-k} = \frac{n(n-1)(n-2)...(n-k+1)}{n^k}\frac{\lambda^k}{k!}(1-\frac{\lambda}{n})^n(1-\frac{\lambda}{n})^{-k}
这里：
\frac{n(n-1)(n-2)...(n-k+1)}{n^k} = (1-\frac{1}{n})(1-\frac{2}{n})...(1-\frac{k-1}{n})
根据极限：
\lim\limits_{n\rightarrow\infty}(1-\frac{1}{n})(1-\frac{2}{n})...(1-\frac{k-1}{n}) = 1
\lim\limits_{n\rightarrow\infty}(1-\frac{\lambda}{n})^{-k} = 1
\lim\limits_{n\rightarrow\infty}(1-\frac{\lambda}{n})^n=e^{-\lambda}
所以我们最终就证明了：
\lim\limits_{n\rightarrow\infty} {n \choose k}p^k(1-p)^{n-k} = \frac {\lambda^k e^{-\lambda}}{k!}
这正好符合我们NGS测序的特点：高通量大数据（N很大，而一般来说每个位点i的DNA fragments占比p_i也很小）。

一样的，简单绘个图：

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


这里\lambda = 10，我们来验证一下，当n很大，p很小的时候，二项分布和泊松分布到底有没有很大的差别：

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为：15/2^6=0.234375ng。注意这个0.234375ng里面的每个分子都是独特的，也就是大家经常说的unique。

同时，DNA分子的平均质量为660g/mol/bp，所以15ng的长300bp的文库里面大概有4.56*10^{10}个DNA，原始的“独特”DNA大概有7.125*10^8个，如果你要产生6G的测序数据（PE150），就相当于你测了其中的2*10^7个DNA。对应到前面的泊松分布：n=2*10^7，p=\frac{1}{7.125*10^8}，所以\lambda=np=0.028。所以原始DNA中一个“独特”的DNA分子被测到的次数所服从的泊松分布就为X \sim P(0.028)：

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

可以看到，能测到2次以上的概率远低于测到1次的。此外还有很大一部分概率是，这个DNA根本测不到。

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

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

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

可以看到，能测到2次的概率已经上来了，也就意味着这两轮额外的PCR已经为你的数据带来了一些PCR duplicate了。

什么情况下会产生PCR duplicate？

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

  • 本身建库起始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的筒子们可以在评论区留言~

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

©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容