用蒙特卡洛方法評(píng)定轉(zhuǎn)基因玉米59122的測(cè)量不確定度

宋 君， 劉月悅， 王 東， 張富麗， 李 潔

(1.四川省農(nóng)業(yè)科學(xué)院分析測(cè)試中心，四川成都 610066； 2.四川省農(nóng)業(yè)科學(xué)院農(nóng)村經(jīng)濟(jì)與農(nóng)業(yè)信息研究所，四川成都 610066

轉(zhuǎn)基因玉米59122品系是一種既能殺蟲(chóng)(鞘翅目)又能抗除草劑的“雙抗”型轉(zhuǎn)基因作物。該轉(zhuǎn)基因品系由美國(guó)陶氏農(nóng)業(yè)科學(xué)公司(DOW AgroSciences LLC)和先鋒國(guó)際公司(Pioneer Hi-Bred International Inc.)開(kāi)發(fā)，2001—2003年通過(guò)食用和環(huán)境安全評(píng)價(jià)。2004年轉(zhuǎn)基因玉米59122被投放到美國(guó)和墨西哥的食品和飼料市場(chǎng)。2005年美國(guó)和墨西哥開(kāi)始商業(yè)化種植該品系。2010年德國(guó)開(kāi)始對(duì)59122玉米品系進(jìn)行環(huán)境釋放。截至目前，59122玉米品系商業(yè)化種植時(shí)間已有12年[1]。我國(guó)于2006年批準(zhǔn)進(jìn)口轉(zhuǎn)基因59122玉米用作食品和飼料的加工原材料，但不允許商業(yè)化種植。截至目前，我國(guó)已經(jīng)連續(xù)進(jìn)口該轉(zhuǎn)基因玉米品系11年[2]。轉(zhuǎn)基因玉米59122品系包含2個(gè)殺蟲(chóng)基因Cry34Ab1、Cry35Ab1以及1個(gè)耐除草劑基因pat，共3個(gè)目的基因。抗蟲(chóng)基因分別由啟動(dòng)子Ubi Zm1、P-peroxidase和終止子PinII調(diào)控，耐除草劑基因由啟動(dòng)子CaMV35S和終止子NOS調(diào)控。為了區(qū)分轉(zhuǎn)基因產(chǎn)品與非轉(zhuǎn)基因產(chǎn)品，貫徹我國(guó)《農(nóng)業(yè)轉(zhuǎn)基因生物安全管理?xiàng)l例》《農(nóng)業(yè)轉(zhuǎn)基因生物標(biāo)識(shí)管理辦法》和《轉(zhuǎn)基因食品衛(wèi)生管理辦法》等法律、法規(guī)，宋君等針對(duì)轉(zhuǎn)基因玉米59122品系建立了轉(zhuǎn)基因玉米59122品系的實(shí)時(shí)熒光定量PCR檢測(cè)方法，并驗(yàn)證了該法的特異性和靈敏度[3]。

近年來(lái)，隨著我國(guó)經(jīng)濟(jì)和技術(shù)的快速發(fā)展，“中國(guó)制造”的各領(lǐng)域產(chǎn)品出口到世界各國(guó)，尤其是食品和農(nóng)產(chǎn)品貿(mào)易遍及全球。很多進(jìn)口國(guó)的法律和法規(guī)對(duì)含有轉(zhuǎn)基因成分的食品和農(nóng)產(chǎn)品有超過(guò)限量標(biāo)識(shí)的規(guī)定，因此，精準(zhǔn)定量分析產(chǎn)品中轉(zhuǎn)基因成分含量對(duì)我國(guó)進(jìn)出口貿(mào)易有重要意義。測(cè)量不確定度評(píng)定是繼誤差理論后發(fā)展起來(lái)的對(duì)測(cè)量質(zhì)量進(jìn)行科學(xué)評(píng)估的重要手段。目前，測(cè)量不確定度已經(jīng)被世界各國(guó)廣泛認(rèn)可和接受。產(chǎn)品進(jìn)口國(guó)往往要求出口方提供產(chǎn)品的質(zhì)量和安全檢測(cè)數(shù)據(jù)及其不確定度以評(píng)估檢測(cè)數(shù)據(jù)的可靠性。不確定度評(píng)定在我國(guó)的大量應(yīng)用始于21世紀(jì)初，但應(yīng)用范圍主要局限于機(jī)械、工程、航空、物理和化學(xué)等領(lǐng)域[4-8]，在生物領(lǐng)域的應(yīng)用相對(duì)較少[9]。由于轉(zhuǎn)基因生物安全學(xué)是近年發(fā)展起來(lái)的新學(xué)科，因此，測(cè)量不確定度評(píng)定在轉(zhuǎn)基因生物及產(chǎn)品安全檢測(cè)中的應(yīng)用更是非常有限[10]，尚未得到推廣和普及。

測(cè)量不確定度評(píng)定的蒙特卡洛法(Monte Carlo Method，MCM)是2008年由國(guó)際計(jì)量學(xué)組織(JCGM)針對(duì)《不確定度評(píng)定表達(dá)指南》(Guide to the Expression of Uncertainty in Measurement，GUM)存在線性模型近似假設(shè)、復(fù)雜模型靈敏度的偏導(dǎo)數(shù)難求、輸入量概率密度函數(shù)對(duì)稱和輸出量概率密度函數(shù)正態(tài)分布的假設(shè)等缺陷[11]，提出的基于數(shù)字計(jì)算仿真模擬實(shí)施“概率分布傳播”的測(cè)量不確定度評(píng)定的另一種方法。與GUM法[12]相比，可能是因?yàn)轫氁婕暗接?jì)算機(jī)編程等問(wèn)題，測(cè)量不確定度的蒙特卡洛評(píng)定方法在國(guó)際上應(yīng)用相對(duì)較少[13-15]。截至目前，國(guó)際上尚未有采用MCM評(píng)定轉(zhuǎn)基因成分測(cè)量不確定度的研究報(bào)道。為了明確影響轉(zhuǎn)基因成分定量檢測(cè)不確定度的因素及其貢獻(xiàn)大小，提高轉(zhuǎn)基因成分檢測(cè)質(zhì)量，維護(hù)我國(guó)農(nóng)產(chǎn)品貿(mào)易利益，本研究首次采用MCM法評(píng)定混合樣品中轉(zhuǎn)基因玉米59122的測(cè)量不確定度，以期為轉(zhuǎn)基因生物及產(chǎn)品的測(cè)量不確定度評(píng)定提供新方法探索和實(shí)踐指導(dǎo)。

1 測(cè)量過(guò)程

1.1 標(biāo)準(zhǔn)物質(zhì)及試樣

采用含量為10%的轉(zhuǎn)基因玉米59122有證標(biāo)準(zhǔn)物質(zhì)(certified reference material，CRM)作為儀器的校準(zhǔn)子(calibrator)。用含有0.5%轉(zhuǎn)基因玉米59122的混合樣品作為測(cè)試樣品。

1.2 主要儀器設(shè)備

用美國(guó)熱電公司生產(chǎn)的型號(hào)為Nanodrop-1000的分光光度計(jì)測(cè)量標(biāo)準(zhǔn)物質(zhì)和試樣的核酸濃度。轉(zhuǎn)基因玉米59122片段(載體與受體植物基因組相結(jié)合的側(cè)翼片段)的實(shí)時(shí)擴(kuò)增和59122片段的拷貝數(shù)測(cè)量在美國(guó)ABI公司生產(chǎn)的7500型實(shí)時(shí)PCR系統(tǒng)上進(jìn)行(全部試驗(yàn)工作于2014年在四川省農(nóng)業(yè)科學(xué)院分析測(cè)試中心轉(zhuǎn)基因生物安全研究室完成)。

1.3 關(guān)鍵化學(xué)試劑

本研究涉及到的主要試劑包括核酸提取和實(shí)時(shí)PCR化學(xué)試劑，購(gòu)自天根生物技術(shù)有限公司(北京)。

1.4 引物、探針合成

本研究使用的引物/探針和試驗(yàn)方法按照文獻(xiàn)[3]報(bào)道的引物/探針序列和方法開(kāi)展試驗(yàn)。引物/探針由上海生工合成和標(biāo)記。

1.5 校準(zhǔn)曲線擬合和試樣的測(cè)量

將分離得到的10%含量的轉(zhuǎn)基因玉米59122(CRM)的DNA溶液分別稀釋成100、20、4、0.8、1.6×10-4μg/μL 5個(gè)濃度。根據(jù)玉米基因組大小2 504 Mbp[16]計(jì)算，0.1 μg/μL DNA溶液相當(dāng)于1 μL DNA溶液中含36 400 copies DNA分子。用3倍體積的上述5個(gè)濃度點(diǎn)的DNA(Q)和經(jīng)PCR反應(yīng)得到的儀器響應(yīng)值(threshold of cycles，Ct)作為校準(zhǔn)子擬合校準(zhǔn)曲線。實(shí)時(shí)PCR反應(yīng)體系含探針型主要混和液(dNTPs、Mg2+、聚合酶等成分)12.5 μL，濃度為10 μmol/L的上/下游引物各1 μL，濃度為10 μmol/L的探針0.5 μL，然后添加滅菌的二次蒸餾水至25 μL。5個(gè)DNA濃度點(diǎn)分別做3次平行試驗(yàn)；測(cè)試樣品做16次重復(fù)測(cè)量。實(shí)時(shí)PCR程序?yàn)?5 ℃，10 min； 95 ℃，15 s；60 ℃，1 min(40 cycles)。

2 轉(zhuǎn)基因玉米59122分子片段的數(shù)學(xué)測(cè)量模型

59122片段絕對(duì)含量的對(duì)數(shù)值(lgQ)與儀器響應(yīng)值(Ct)之間存在函數(shù)關(guān)系，

Ct=m×lgQ+k。

(1)

式中：Q為內(nèi)源基因zSSIIb或59122片段的絕對(duì)含量，μg或copies；Ct為PCR儀對(duì)Q的響應(yīng)值；k為校準(zhǔn)曲線的截距；m為校準(zhǔn)曲線的斜率。

將試樣的zSSIIb和59122片段的Ct值分別代入各自的校準(zhǔn)曲線，計(jì)算得到試樣zSSIIb和59122片段的絕對(duì)含量Q，然后將Q代入式(2)，獲得59122片段的相對(duì)含量C。

C=Qex/Qed×100%。

(2)

式中：Qex為zSSIIb基因的絕對(duì)含量；Qed為59122片段的絕對(duì)含量；C為59122片段的相對(duì)含量，%。

3 MCM的實(shí)施

3.1 m和k的聯(lián)合分布及計(jì)算公式

按照GUIDE 98-3/SUPP.1概率分布傳遞原理[17]，使用MATLAB軟件實(shí)施MCM評(píng)定轉(zhuǎn)基因玉米59122片段的測(cè)量不確定度。由于斜率m和截距k線性相關(guān)，所以m和k屬于聯(lián)合分布(雙元變量正態(tài)分布)，用矢量E及其協(xié)方差矩陣V來(lái)表示。

(3)

(4)

式中：u(k，m)為k和m的協(xié)方差；u(m)為m的不確定度；u(k)為k的不確定度。

u(m)和u(k)分別由式(5)、式(6)計(jì)算：

(5)

(6)

而s(Ct)由式(7)計(jì)算：

(7)

式中：n為校準(zhǔn)測(cè)量次數(shù)；

Cti為PCR儀對(duì)第i個(gè)校準(zhǔn)子的響應(yīng)值。

協(xié)方差u(k，m)由式(8)計(jì)算：

u(k,m)=r(k,m)u(k)u(m)。

(8)

相關(guān)系數(shù)

(9)

3.2 MCM的輸入

5個(gè)校準(zhǔn)子的對(duì)數(shù)值(Qi和Qi#)與其對(duì)應(yīng)的響應(yīng)值Cti和Cti#(表1、表2)以均勻分布R(a，b)(a為均勻分布的下限，b為上限)作為輸入量，擬合內(nèi)源基因zSSIIb和59122片段的校準(zhǔn)曲線。在zSSIIb基因擴(kuò)增過(guò)程中，0.1 μg/μL濃度點(diǎn)的響應(yīng)值Ct1為22.44)，標(biāo)準(zhǔn)差(SD)為0.04(表1)，因此，循環(huán)閾值Ct1的分布范圍為Ct1～U[22.40，22.48]。類似地，計(jì)算獲得各濃度點(diǎn)的Ct響應(yīng)值分布范圍。均勻分布的抽樣方法為從標(biāo)準(zhǔn)均勻分布R(0，1)中抽取隨機(jī)數(shù)r，構(gòu)造ξ=a+(b-a)r。MATLAB軟件在擬合內(nèi)源基因zSSIIb和59122片段的校準(zhǔn)曲線時(shí)，矢量E、協(xié)方差矩陣V、u(m)、u(k)、u(k，m)和r(k，m)分別按照式(3)至(6)、式(8)至(9)計(jì)算。

16次重復(fù)測(cè)量試樣獲得zSSIIb基因和59122片段的Ct值(表3)概率分布遵守正態(tài)分布N[y，u2(y)][y為最佳估計(jì)值，u(y)為標(biāo)準(zhǔn)不確定度]，抽樣方法為從標(biāo)準(zhǔn)正態(tài)分布N(0，1)中抽取隨機(jī)數(shù)r，構(gòu)造ξ=y+u(y)r。

3.3 MCM的傳播

根據(jù)表1、表2、表3中的數(shù)據(jù)、“3.2”節(jié)中設(shè)定的各輸入量概率分布以及抽樣方法，編寫(xiě)MATLAB程序代碼模擬運(yùn)行抽樣，獲得106個(gè)隨機(jī)樣本以及每個(gè)輸入量輸入數(shù)學(xué)模型[式(1)至(9)]計(jì)算得到106個(gè)均勻分布或正態(tài)分布數(shù)據(jù)組，獲得輸出量的最佳估計(jì)值、標(biāo)準(zhǔn)不確定度及95%包含概率(coverage probability)下的包含區(qū)間(coverage interval)。

3.4 MCM的輸出

MATLAB程序模擬Cti、Cti#、CtSi、CtSi#、m、k、QzSSIIb、Q59122和C等9個(gè)輸入量的概率分布、抽樣及通過(guò)數(shù)學(xué)模型式(1)、式(3)至(9)獲得的平均值、標(biāo)準(zhǔn)不確定度和95%包含概率的包含區(qū)間。

表1 zSSIIb內(nèi)源基因校準(zhǔn)曲線擬合

表2 59122片段校準(zhǔn)曲線擬合

表3 試樣中內(nèi)源基因zSSIIb和59122片段的Ct值

4 結(jié)果

由表4可知，zSSIIb基因和59122片段的校準(zhǔn)曲線的斜率分別為-3.60和-3.29，截距分別為31.50和30.41，即zSSIIb基因和59122片段的校準(zhǔn)曲線分別為y=-3.60x+31.50和y=-3.29x+30.41。QzSSIIb和Q59122既作為輸入量又作為輸出量的平均值分別為178.47和1.02，不確定度分別為1.99和0.03，95%包含概率的包含區(qū)間分別為174.59～182.41和0.96～1.08。獲得的59122片段相對(duì)含量C的平均值、標(biāo)準(zhǔn)不確定度和95%包含概率的包含區(qū)間分別為 0.57%、1.82×10-4和0.54%～0.61%。此外，MATLAB程序運(yùn)行產(chǎn)生的混合樣品中轉(zhuǎn)基因玉米59122的相對(duì)含量的概率密度分布為典型的“鐘形”正態(tài)分布(圖1)。

表4 各輸入/輸出量的概率分布、均值、不確定度和95%置信度的包含區(qū)間

續(xù)表4

影響因子或輸入/輸出量概率分布平均值標(biāo)準(zhǔn)不確定度95%包含概率的包含區(qū)間CtS1均勻分布23.410.0723.29～23.53CtS2均勻分布23.450.0723.33～23.57CtS3均勻分布23.480.0723.36～23.60CtS4均勻分布23.450.0723.33～23.57CtS5均勻分布23.450.0723.33～23.57CtS6均勻分布23.490.0722.37～22.61CtS7均勻分布23.450.0723.33～23.57CtS8均勻分布23.320.0723.20～23.44CtS9均勻分布23.400.0723.28～23.52CtS10均勻分布23.450.0723.33～23.57CtS11均勻分布23.500.0723.38～23.62CtS12均勻分布23.580.0723.46～23.70CtS13均勻分布23.560.0723.44～23.68CtS14均勻分布23.490.0723.37～23.61CtS15均勻分布23.510.0723.39～23.63CtS16均勻分布23.360.0723.24～23.48CtS1#均勻分布30.250.1729.97～30.53CtS2#均勻分布30.370.1730.09～30.65CtS3#均勻分布30.380.1730.10～30.66CtS4#均勻分布29.940.1729.67～30.22CtS5#均勻分布30.520.1730.25～30.81CtS6#均勻分布30.370.1730.09～30.65CtS7#均勻分布30.350.1730.07～30.63CtS8#均勻分布30.390.1730.11～30.67CtS9#均勻分布30.230.1729.95～30.51CtS10#均勻分布30.630.1730.35～30.91CtS11#均勻分布30.340.1730.06～30.62CtS12#均勻分布30.510.1730.23～30.79CtS13#均勻分布30.570.1730.29～30.85CtS14#均勻分布30.630.1730.35～30.91CtS15#均勻分布30.320.1730.04～30.60CtS16#均勻分布30.290.1730.01～30.57QzSSIIb正態(tài)分布178.471.99174.59～182.41Q59122正態(tài)分布1.020.030.96～1.08C正態(tài)分布0.005 71.82×10～40.005 4～0.006 1

注：Ct11～Ct53(統(tǒng)稱Cti)表示zSSIIb基因測(cè)量中標(biāo)準(zhǔn)物質(zhì)各濃度點(diǎn)的響應(yīng)值；Ct11#～Ct53#(統(tǒng)稱Cti#)表示59122玉米品系測(cè)量中標(biāo)準(zhǔn)物質(zhì)各濃度點(diǎn)的響應(yīng)值；CtS1～CtS16(統(tǒng)稱CtSi)表示試樣16次重復(fù)測(cè)量zSSIIb基因獲得的Ct值；CtS1#～CtS16#(統(tǒng)稱CtSi#)表示試樣16次重復(fù)測(cè)量59122片段獲得的Ct值；QzSSIIb表示試樣中zSSIIb的絕對(duì)含量；Q59122表示試樣中59122片段的絕對(duì)含量；C表示試樣中59122片段的相對(duì)含量；mzSSIIb和m59122分別為zSSIIb基因和59122片段的校準(zhǔn)曲線的斜率；kzSSIIb和k59122分別為zSSIIb基因和59122片段的校準(zhǔn)曲線的截距。

5 討論

本研究首次報(bào)道了采用蒙特卡洛法實(shí)施“概率分布傳播”評(píng)定轉(zhuǎn)基因成分的測(cè)量不確定度。轉(zhuǎn)基因玉米59122的相對(duì)含量為0.57%，與其“理論含量”0.5%的相對(duì)偏差(bias)為14%，在國(guó)際廣泛認(rèn)可的轉(zhuǎn)基因成分定量檢測(cè)偏差范圍(-25%≤bias≤25%)內(nèi)。本次試驗(yàn)中，在95%的包含概率下轉(zhuǎn)基因玉米59122的相對(duì)含量檢測(cè)值落在0.54%～0.61%范圍內(nèi)， 不確定度極小(1.82×10-4)， 包含區(qū)間范圍較窄，顯示檢測(cè)質(zhì)量較高，數(shù)據(jù)可靠?；旌蠘悠分修D(zhuǎn)基因玉米59122的相對(duì)含量的概率密度分布呈正態(tài)分布，與用GUM法評(píng)定測(cè)量不確定的假設(shè)條件“輸出量概率密度函數(shù)為正態(tài)分布”一致，因此，揭示了轉(zhuǎn)基因成分測(cè)量過(guò)程中各個(gè)輸入量/輸出量的概率分布及數(shù)據(jù)對(duì)稱性等均滿足用GUM法評(píng)定測(cè)量不確定度的假設(shè)條件，即轉(zhuǎn)基因成分測(cè)量不確定度的評(píng)定可以采用GUM法。

在表4的9個(gè)變量的不確定度中，不確定度貢獻(xiàn)最大的是玉米內(nèi)源基因zSSIIb絕對(duì)含量的不確定度(1.99)，轉(zhuǎn)基因玉米59122的相對(duì)含量C的不確定度最小(1.82×10-4)，顯示出檢測(cè)質(zhì)量高、數(shù)據(jù)可靠。雖然玉米內(nèi)源基因zSSIIb絕對(duì)含量的不確定度最大，但其對(duì)測(cè)量結(jié)果影響非常小。因此，將玉米內(nèi)源基因zSSIIb絕對(duì)含量的概率分布輸入測(cè)量模型傳遞到轉(zhuǎn)基因玉米59122的相對(duì)含量C這個(gè)最終輸出量時(shí)，測(cè)量不確定度已經(jīng)非常小。此外，部分濃度(Ct51～Ct53、Ct11#～Ct13#、Ct41～Ct43、Ct41#～Ct43#和Ct21～Ct23)的標(biāo)準(zhǔn)物質(zhì)稀釋貢獻(xiàn)的不確定度分布在0.28～0.40之間?？傮w上，這些不確定度居于本次測(cè)量不確定度的中等水平，反映出轉(zhuǎn)基因成分測(cè)量不確定度的主要因素為標(biāo)準(zhǔn)物質(zhì)的梯度稀釋，而影響梯度稀釋的主要原因是微量移液器的使用，這與采用GUM法評(píng)定轉(zhuǎn)基因成分測(cè)量不確定度報(bào)道的主要因素為微量移液[10]一致。本研究的結(jié)果表明，影響轉(zhuǎn)基因成分測(cè)量的主要因素為微量移液，因此在轉(zhuǎn)基因成分的定量分析中應(yīng)將提高微量液體體積轉(zhuǎn)移的準(zhǔn)確度作為首要考慮對(duì)象。