生物中的数学【二】为什么你的文库总是充满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^7p=\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的筒子们可以在评论区留言~

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

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

推荐阅读更多精彩内容