基于C語(yǔ)言模擬計(jì)算的cdPCR最佳反應(yīng)通孔數(shù)分析

周德良

（北京旌準(zhǔn)醫(yī)療科技有限公司，北京 100190）

C語(yǔ)言是面向過(guò)程的高級(jí)程序設(shè)計(jì)語(yǔ)言[1-6]，具有功能強(qiáng)大、應(yīng)用面廣、使用簡(jiǎn)單方便、可移植性好等優(yōu)點(diǎn)[2-4]。C語(yǔ)言不僅可以用于系統(tǒng)軟件編寫、單片機(jī)系統(tǒng)開(kāi)發(fā)，還可以用于大數(shù)據(jù)統(tǒng)計(jì)與分析[3-4，7-8]。數(shù)字PCR（dgital PCR，dPCR）是一種核酸分子絕對(duì)定量技術(shù)[9-12]，其不需要采用看家基因或標(biāo)準(zhǔn)曲線，具有高靈敏度、高準(zhǔn)確性、高耐受性、絕對(duì)定量等優(yōu)點(diǎn)[13-16，22]。dPCR主要包括液滴數(shù)字PCR（ddPCR）和芯片數(shù)字PCR（cdPCR）兩種，其中cdPCR通過(guò)把反應(yīng)液均勻分到芯片上的反應(yīng)通孔中進(jìn)行PCR反應(yīng)，然后根據(jù)泊松分布和熒光信號(hào)陽(yáng)性比例來(lái)計(jì)算目的核酸的拷貝數(shù)[13]。

利用C語(yǔ)言程序來(lái)模擬芯片數(shù)字PCR中目的核酸進(jìn)入反應(yīng)通孔的過(guò)程，從而可以對(duì)干擾粒子對(duì)目的核酸分布的影響[17-20]、目的核酸最佳反應(yīng)通孔數(shù)、目的核酸分布均勻性等問(wèn)題進(jìn)行分析與計(jì)算，旨在為數(shù)字PCR的芯片設(shè)計(jì)提供一定的參數(shù)依據(jù)與理論指導(dǎo)。

1 模擬分析算法

通過(guò)C語(yǔ)言模擬“干擾作用完全隨機(jī)的情況下目標(biāo)粒子進(jìn)入反應(yīng)通孔的過(guò)程”的模擬算法如圖1所示，具體過(guò)程如下：

（1）在[0，h+n*h/m）區(qū)間內(nèi)生成m+n個(gè)隨機(jī)數(shù)，其中[0，h-1]區(qū)間內(nèi)的隨機(jī)數(shù)代表目標(biāo)粒子，總數(shù)為m；[h，h+n*h/m）區(qū)間內(nèi)的隨機(jī)數(shù)代表干擾粒子，總數(shù)為n；h代表反應(yīng)通孔數(shù)；

（2）統(tǒng)計(jì)[0，h-1]內(nèi)非重復(fù)隨機(jī)數(shù)的總數(shù)N，N表示m個(gè)目標(biāo)粒子在h個(gè)反應(yīng)通孔中的入孔總數(shù)；

（3）將上述過(guò)程重復(fù)K次，并統(tǒng)計(jì)不同N對(duì)應(yīng)的出現(xiàn)次數(shù)Y；

（4）根據(jù)不同h對(duì)應(yīng)的N-Y曲線的重疊度或N-Y主峰N值隨h的變化率來(lái)判斷最佳反應(yīng)通孔數(shù)。

2 無(wú)干擾的最佳反應(yīng)通孔數(shù)

在沒(méi)有干擾粒子的情況下（n=0），為了分析目的核酸（或目標(biāo)粒子）在檢測(cè)芯片反應(yīng)通孔內(nèi)的分布情況，根據(jù)圖1所示采用如下統(tǒng)計(jì)過(guò)程：

圖1 模擬算法程序圖Fig.1 The program graph of simulation algorithm

設(shè)目標(biāo)粒子個(gè)數(shù)m為1 000，反應(yīng)通孔個(gè)數(shù)h為1 000 a（a取1到10的所有整數(shù)）。當(dāng)反應(yīng)通孔數(shù)為1 000時(shí)（a=1），利用C算法程序隨機(jī)生成1 000個(gè)介于[0，999]的隨機(jī)數(shù)，代表1 000個(gè)目標(biāo)粒子。當(dāng)生成隨機(jī)數(shù)x（0≤x≤999）時(shí)表示目標(biāo)粒子落入標(biāo)號(hào)為x的反應(yīng)通孔內(nèi)，統(tǒng)計(jì)1 000個(gè)隨機(jī)數(shù)的非重復(fù)數(shù)值的總數(shù)N，N就表示1 000個(gè)目標(biāo)粒子在1 000個(gè)反應(yīng)通孔中的入洞總數(shù)。將上述過(guò)程重復(fù)20萬(wàn)次，統(tǒng)計(jì)不同入洞總數(shù)N對(duì)應(yīng)的出現(xiàn)次數(shù)Y，從而得到1 000個(gè)目標(biāo)粒子在1 000個(gè)反應(yīng)通孔中的分布情況，如圖2所示，其入洞總數(shù)非零區(qū)間為[589，674]，其中當(dāng)入洞總數(shù)N=632時(shí)出現(xiàn)次數(shù)最大。

圖2 1 000個(gè)目標(biāo)粒子在1 000個(gè)反應(yīng)通孔中的分布Fig.2 The distribution of 1 000 target particles in 1 000 reaction holes

當(dāng)反應(yīng)通孔數(shù)為2 000時(shí)（a=2），利用C算法程序隨機(jī)生成1 000個(gè)介于[0，1 999]的隨機(jī)數(shù)，當(dāng)生成隨機(jī)數(shù)x（0≤x≤999）時(shí)表示目標(biāo)粒子落入標(biāo)號(hào)為x的反應(yīng)通孔內(nèi)，并統(tǒng)計(jì)1 000個(gè)隨機(jī)數(shù)的非重復(fù)數(shù)值的總數(shù)N。將上述過(guò)程重復(fù)20萬(wàn)次，統(tǒng)計(jì)不同入洞總數(shù)N對(duì)應(yīng)的出現(xiàn)次數(shù)Y，從而得到1 000個(gè)目標(biāo)粒子在2 000個(gè)反應(yīng)通孔中的分布情況，如圖3所示，其入洞總數(shù)非零區(qū)間為[745，829]，其中當(dāng)入洞總數(shù)N=787時(shí)出現(xiàn)次數(shù)最大。

依此類推，可以得到當(dāng)反應(yīng)通孔數(shù)為1 000 a時(shí)（a分別為3到10的整數(shù)）1 000個(gè)目標(biāo)粒子的分布情況，見(jiàn)圖4。

對(duì)應(yīng)圖4，各曲線波峰對(duì)應(yīng)的入洞總數(shù)N見(jiàn)表1，而由表1得到的曲線點(diǎn)圖見(jiàn)圖5。

圖3 1 000個(gè)目標(biāo)粒子在2 000個(gè)反應(yīng)通孔中的分布Fig.3 The distribution of 1 000 target particles in 2 000 reaction holes

圖4 1 000個(gè)目標(biāo)粒子在1 000 a（a為從1到10的全部整數(shù)）個(gè)反應(yīng)通孔中的分布Fig.4 The distribution of 1 000 target particles in 1 000 a（a is integer from 1 to 10）reaction holes

表1 N-Y曲線波峰對(duì)應(yīng)的入洞總數(shù)NTable 1 The value N corresponding to peak of curve

從圖4發(fā)現(xiàn)，隨著芯片反應(yīng)通孔數(shù)的增多，1 000個(gè)目標(biāo)粒子的入洞總數(shù)N也越來(lái)越大，即1 000個(gè)目標(biāo)粒子的分布均勻度越來(lái)越高，但隨著反應(yīng)通孔數(shù)的進(jìn)一步增多，不同反應(yīng)通孔數(shù)對(duì)應(yīng)的N-Y曲線的毗鄰程度越來(lái)越高，重疊程度也越來(lái)越大；而表1和圖5同樣表明，雖然隨著反應(yīng)通孔數(shù)的增加目標(biāo)粒子的分布均勻度會(huì)越來(lái)越高，但當(dāng)反應(yīng)通孔數(shù)大于某一值時(shí)分布均勻度的提高已開(kāi)始趨于飽和（增加非常平緩）。因此，為提高目標(biāo)粒子在反應(yīng)通孔中的分布均勻度，直接有效的方法是增加芯片上的反應(yīng)通孔數(shù)，但當(dāng)芯片反應(yīng)通孔數(shù)大于等于8 000或9 000時(shí)反應(yīng)通孔數(shù)增加對(duì)提高目標(biāo)粒子分布均勻度的作用將趨于飽和，如果繼續(xù)增加反應(yīng)通孔數(shù)，則將帶來(lái)工藝難度的提升與設(shè)計(jì)成本的增加，特對(duì)是對(duì)于目的核酸片段較少、檢測(cè)靈敏度要求不高的實(shí)驗(yàn)，高昂PCR儀顯然與實(shí)際需求相違背。所以，在1 000個(gè)目標(biāo)粒子的基礎(chǔ)上，cdPCR檢測(cè)芯片的最佳反應(yīng)通孔數(shù)應(yīng)為8 000或9 000，并由此推斷，在無(wú)干擾粒子的情況下，不論目標(biāo)粒子個(gè)數(shù)m取多少，檢測(cè)目標(biāo)粒子的最佳反應(yīng)通孔數(shù)與目標(biāo)粒子個(gè)數(shù)之比應(yīng)為 8∶1 或 9∶1。

圖5 表1數(shù)據(jù)對(duì)應(yīng)的點(diǎn)線圖Fig.5 The dot-line graph corresponding to the datum of table 1

3 含干擾的最佳反應(yīng)通孔數(shù)

3.1 干擾粒子對(duì)目標(biāo)粒子分布的影響

文獻(xiàn)[21]分析了干擾作用完全隨機(jī)的情況下，干擾粒子有無(wú)及其多少對(duì)目標(biāo)粒子分布的影響，并得到了1 000個(gè)目標(biāo)粒子在1 000 b個(gè)（b分別取0、10、100、1 000、10 000） 干擾粒子作用下分別在 1 000、3 000、7 000個(gè)反應(yīng)通孔中的分布情況[21]，分別見(jiàn)圖6、圖 7與圖 8，其中圖 6～圖8的（a）表示無(wú)干擾粒子對(duì)應(yīng)的目標(biāo)粒子分布；圖 6～圖 8 的（b）、（c）、（d）、（e）分別對(duì)應(yīng)干擾粒子個(gè)數(shù)是目標(biāo)粒子的10、100、1 000、10 000 倍的目標(biāo)粒子分布；圖 6～圖 8 的（e）為前述5類情況的合成圖。文獻(xiàn)[21]得到結(jié)論如下：

不論反應(yīng)通孔數(shù)h取多少，含干擾粒子與不含干擾粒子的5類N-Y曲線分布基本重合或相同，且干擾粒子個(gè)數(shù)越小，含干擾粒子的N-Y曲線分布與不含干擾粒子的N-Y曲線分布就具有越高的重合度或相似度。上述結(jié)果表明，在干擾作用完全隨機(jī)的情況下，對(duì)于任意取值的反應(yīng)通孔數(shù)h，干擾粒子有無(wú)及其多少均不會(huì)破壞目標(biāo)粒子的隨機(jī)分布程度，因此，無(wú)論有無(wú)干擾粒子或其存在多少，目標(biāo)粒子的N-Y曲線分布都基本相同。

圖6 目標(biāo)粒子在1 000個(gè)反應(yīng)通孔中的分布Fig.6 The distribution of target particles in 1 000 reaction holes

3.2 分析結(jié)果

分析表明，在干擾作用完全隨機(jī)的情況下，對(duì)于任意取值的反應(yīng)通孔數(shù)h，干擾粒子有無(wú)或其多少均不會(huì)破壞目標(biāo)粒子的隨機(jī)分布程度，因此，無(wú)論有無(wú)干擾粒子或其存在多少，目標(biāo)粒子的N-Y曲線分布[21]都將基本相同。因此，“含干擾的最佳反應(yīng)通孔數(shù)”與“無(wú)干擾的最佳反應(yīng)通孔數(shù)”應(yīng)該幾乎一致，故仍可將含干擾粒子的最佳反應(yīng)通孔數(shù)定為8 000或9 000（目標(biāo)粒子個(gè)數(shù)為1 000的情況下），并由此推斷，在含干擾粒子的情況下，不論目標(biāo)粒子個(gè)數(shù)m取多少，檢測(cè)目標(biāo)粒子的最佳反應(yīng)通孔數(shù)與目標(biāo)粒子個(gè)數(shù)之比仍為 8∶1 或 9∶1。

圖7 目標(biāo)粒子在3 000個(gè)反應(yīng)通孔中的分布Fig.7 The distribution of target particles in 3 000 reaction holes

圖8 目標(biāo)粒子在7 000個(gè)反應(yīng)通孔中的分布Fig.8 The distribution of target particles in 7 000 reaction holes

4 目標(biāo)粒子分布的均勻性

目標(biāo)粒子分布的均勻性定義為“含目標(biāo)粒子的反應(yīng)通孔數(shù)”與“目標(biāo)粒子總數(shù)”的百分比。為了分析目標(biāo)粒子在不同反應(yīng)通孔數(shù)下的分布均勻性，采用如下統(tǒng)計(jì)方法：

設(shè)目標(biāo)粒子個(gè)數(shù)為1 000，反應(yīng)通孔數(shù)h為1 000 a（a取1到200的任一整數(shù)）。利用C算法程序隨機(jī)生成1 000個(gè)介于[0，1 000 a-1]的隨機(jī)數(shù)，代表1 000個(gè)目標(biāo)粒子。當(dāng)生成隨機(jī)數(shù)x（0≤x≤1 000 a-1）時(shí)表示目標(biāo)粒子落入標(biāo)號(hào)為x的反應(yīng)通孔內(nèi)，統(tǒng)計(jì)1 000個(gè)隨機(jī)數(shù)的非重復(fù)數(shù)值的總數(shù)N，N就表示1 000個(gè)目標(biāo)粒子在1 000 a個(gè)反應(yīng)通孔中的入洞總數(shù)。將上述模擬過(guò)程重復(fù)5萬(wàn)次，并計(jì)算入洞總數(shù)N的平均值，該平均值即“含目標(biāo)粒子的反應(yīng)通孔數(shù)”，再用其除以“目標(biāo)粒子總數(shù)”，從而得到1 000個(gè)目標(biāo)粒子在1 000 a個(gè)反應(yīng)通孔中的分布均勻性，具體模擬計(jì)算結(jié)果見(jiàn)表2。

表2 目標(biāo)粒子分布均勻性Table 2 The distribution uniformity of target particles

依據(jù)上述結(jié)果，可以得到最佳反應(yīng)通孔數(shù)8 000或9 000對(duì)應(yīng)的目標(biāo)粒子分布均勻性為94.0%或94.6%；同時(shí)根據(jù)上述方法，可以得到滿足指定要求的分布均勻性對(duì)應(yīng)的最低反應(yīng)通孔數(shù)，從而為cd－PCR的芯片設(shè)計(jì)提供了理論參數(shù)指導(dǎo)。另外，由表2同樣可以發(fā)現(xiàn)，隨反應(yīng)通孔數(shù)的增多目標(biāo)粒子分布均勻性逐漸增大，但當(dāng)反應(yīng)通孔數(shù)大于某值時(shí)，分布均勻性的增加將變得愈發(fā)緩慢（增速低于1%）并逐漸趨于飽和，由此同樣可以得出最佳反應(yīng)通孔數(shù)應(yīng)為8 000或9 000的結(jié)論。

5 總結(jié)

利用C語(yǔ)言模擬了芯片數(shù)字PCR中目的核酸進(jìn)入反應(yīng)通孔的計(jì)算過(guò)程，并對(duì)cdPCR中干擾粒子對(duì)目的核酸分布的影響、檢測(cè)目的核酸的最佳反應(yīng)通孔數(shù)、目的核酸分布均勻性等問(wèn)題進(jìn)行了分析與計(jì)算，得到如下結(jié)論：

（1）在干擾作用完全隨機(jī)的情況下，對(duì)于任意取值的反應(yīng)通孔數(shù)，干擾粒子有無(wú)或其存在多少均不會(huì)破壞目標(biāo)粒子的隨機(jī)分布程度，因此，無(wú)論有無(wú)干擾粒子或其存在多少，目標(biāo)粒子的N-Y曲線分布都將基本相同。

（2）在干擾作用完全隨機(jī)的情況下，對(duì)于任意取值的反應(yīng)通孔數(shù)，無(wú)論有無(wú)干擾粒子或其存在多少，檢測(cè)目標(biāo)粒子的最佳反應(yīng)通孔數(shù)與目標(biāo)粒子個(gè)數(shù)之比都應(yīng)為 8∶1 或 9∶1。

（3）當(dāng)目標(biāo)粒子個(gè)數(shù)為1 000時(shí)檢測(cè)目標(biāo)粒子的最佳反應(yīng)通孔數(shù)為8 000或9 000，而對(duì)應(yīng)的目標(biāo)粒子分布均勻性為94.0%或94.6%。