基于近紅外光譜結(jié)合化學(xué)計(jì)量學(xué)的轉(zhuǎn)基因大豆產(chǎn)地判別

大豆起源于中國(guó)，營(yíng)養(yǎng)豐富，可作為一種理想的食物，在世界上很多國(guó)家和地區(qū)的飲食中受到青睞，是人們不可獲缺的食物。據(jù)中國(guó)海關(guān)總署統(tǒng)計(jì)2020年我國(guó)大豆進(jìn)口總量超過1億t，已成為最大的大豆進(jìn)口國(guó)[1]。進(jìn)口大豆大多為轉(zhuǎn)基因大豆，其粗蛋白、脂肪和黃酮含量較高[2]。但是轉(zhuǎn)基因大豆的潛在風(fēng)險(xiǎn)和危害是無法預(yù)測(cè)的，如對(duì)生物多樣性的影響，對(duì)人類和動(dòng)物健康的影響。因此我國(guó)對(duì)轉(zhuǎn)基因大豆管控嚴(yán)格。目前我國(guó)的進(jìn)口大豆來源于美洲國(guó)家，主要進(jìn)口源有美國(guó)、巴西，阿根廷和加拿大。近年來我國(guó)從美國(guó)進(jìn)口轉(zhuǎn)基因大豆的數(shù)量急劇下降，調(diào)整為從美洲其他國(guó)家進(jìn)口大豆。國(guó)外有部分商販從美國(guó)進(jìn)口轉(zhuǎn)基因大豆，而后冒充本國(guó)大豆高價(jià)出口到我國(guó)，牟取不正當(dāng)利益的同時(shí)給我國(guó)海關(guān)對(duì)轉(zhuǎn)基因大豆檢測(cè)和分類造成困難。國(guó)內(nèi)有部分不法分子走私轉(zhuǎn)基因大豆而后銷售到國(guó)內(nèi)市場(chǎng)，這對(duì)我國(guó)的生物安全帶來了嚴(yán)重的威脅，因此對(duì)轉(zhuǎn)基因大豆的產(chǎn)地朔源有利于從源頭打擊不法分子的犯罪行為，保護(hù)我國(guó)的生物安全。

轉(zhuǎn)基因大豆檢測(cè)和鑒別的主要方法是蛋白質(zhì)檢測(cè)方法與核酸檢測(cè)方法[3]。蛋白質(zhì)檢測(cè)法主要包括試紙條法和酶聯(lián)免疫吸附法[4-5]；核酸檢測(cè)方法主要包括定性PCR和環(huán)介導(dǎo)等溫?cái)U(kuò)增技術(shù)[6-7]。以上檢測(cè)方法是破壞性檢測(cè)方法，需大量的實(shí)驗(yàn)試劑、試驗(yàn)過程繁瑣復(fù)雜、投入的成本較高、檢測(cè)專業(yè)性強(qiáng)、不易普及并且不能實(shí)現(xiàn)實(shí)時(shí)在線檢測(cè)。

近紅外光(near infrared，NIR)是介于可見光(visible，Vis)和中紅外(mid infrared，MIR)之間的電磁輻射波。采用NIR光譜技術(shù)分析待測(cè)樣品具有簡(jiǎn)單、高效、無損、實(shí)時(shí)、綠色環(huán)保的優(yōu)點(diǎn)[8]。但是，NIR光譜受環(huán)境和樣品影響較大，容易形成未知組分和灰色體系且有多重共線性問題[9]?；瘜W(xué)計(jì)量學(xué)具有獨(dú)特的優(yōu)勢(shì)，通過對(duì)樣品測(cè)量數(shù)據(jù)的分析，可以最大限度的呈現(xiàn)出樣品的各種化學(xué)信息。NIR光譜結(jié)合化學(xué)計(jì)量學(xué)是一種快速、準(zhǔn)確、高效，可實(shí)現(xiàn)實(shí)時(shí)在線檢測(cè)的方法[10]。

近年來，很多學(xué)者對(duì)NIR光譜結(jié)合化學(xué)計(jì)量學(xué)進(jìn)行研究，NIR光譜結(jié)合化學(xué)計(jì)量學(xué)已廣泛應(yīng)用到食品、農(nóng)業(yè)、醫(yī)藥、化工等多個(gè)領(lǐng)域[11-13]。其在鑒別轉(zhuǎn)基因大豆和非轉(zhuǎn)基因大豆上也成功應(yīng)用，但在轉(zhuǎn)基因大豆的產(chǎn)地判別上鮮有報(bào)道，針對(duì)我國(guó)進(jìn)口轉(zhuǎn)基因大豆的現(xiàn)狀，對(duì)轉(zhuǎn)基因大豆的朔源具有重要意義。

1 材料與方法

1.1 實(shí)驗(yàn)材料

試驗(yàn)選取轉(zhuǎn)基因大豆分別為阿根廷轉(zhuǎn)基因大豆MON89788品系(“A1”)、巴西轉(zhuǎn)基因大豆MON89788品系(“B1”)、美國(guó)轉(zhuǎn)基因大豆MON89788品系(“M1”)、加拿大轉(zhuǎn)基因大豆MON89788品系(“J1”)以上4種大豆為同一品系不同產(chǎn)地的轉(zhuǎn)基因大豆?！癆1”，“M1”分別取50份樣品，“B1”，“J1”分別取80份樣品。如圖1所示，從左到右分別為“A1”、“B1”、“M1”和“J1”單粒大豆特征，單粒大豆在外觀上無明顯差異，所有進(jìn)口轉(zhuǎn)基因大豆都由秦皇島海關(guān)提供。

a-“A1”；b-“B1”；c-“M1”；d-“J1”
圖1 四種轉(zhuǎn)基因大豆單粒特征
Fig.1 Single grain characteristics of four transgenic soybeans

1.2 儀器與設(shè)備

全波反射型NIR光譜儀，檢測(cè)器為Si和InGaAs，光譜掃描范圍400～2 600 nm，北京偉創(chuàng)英圖科技有限公司；ME204E電子天平，梅特勒-托利多有限公司；臺(tái)式真空干燥箱DZF-6050，上海捷呈實(shí)驗(yàn)儀器有限公司。

1.3 實(shí)驗(yàn)方法

1.3.1 NIR光譜采集與光譜預(yù)處理

NIR光譜儀開機(jī)預(yù)熱，白板校正后采集光譜。分別取4種轉(zhuǎn)基因大豆共計(jì)260份樣品，放置在干燥箱內(nèi)，干燥溫度為40 ℃，時(shí)間為36 h，光譜采集在(23±2)℃的恒溫室內(nèi)進(jìn)行，試驗(yàn)所取的光譜范圍為911～2 600 nm，每掃描10次計(jì)算1次平均光譜作為1條原始光譜。每間隔1 nm記錄1個(gè)點(diǎn)，每條原始光譜記錄了1 690個(gè)吸光度值。4種轉(zhuǎn)基因大豆共采集了260條NIR光譜。

試驗(yàn)過程中由于實(shí)驗(yàn)儀器、環(huán)境和樣品的影響，原始光譜中包含一部分的噪音。因此NIR光譜數(shù)據(jù)分析之前要對(duì)NIR原始光譜進(jìn)行預(yù)處理，以此減少或者消除噪音對(duì)試驗(yàn)的影響。平滑是一種提高光譜信噪比的方法，原始光譜經(jīng)過平滑處理后可有效減少光譜噪音。標(biāo)準(zhǔn)正態(tài)變量變換(standard normal variate transformation，SNV)主要是用來消除固體顆粒大小、表面散射以及光程變化對(duì)光譜的影響。本試驗(yàn)采取的光譜預(yù)處理方法為平滑+SNV[14]。

1.3.2 樣品的劃分和編號(hào)

Kennard-Stone(KS)算法，是在總樣本中選出訓(xùn)練集樣品，首先選擇歐氏距離最遠(yuǎn)的2個(gè)樣品進(jìn)入訓(xùn)練集，其后通過計(jì)算剩下的每1個(gè)樣品到訓(xùn)練集內(nèi)每1個(gè)已知樣品的歐式距離，找到擁有最大最小距離的待選樣品放入訓(xùn)練集，以此類推，直到達(dá)到所要求的樣品數(shù)目[15]。本試驗(yàn)共有轉(zhuǎn)基因大豆樣品260份，分別選取“A1”樣品45份、“B1”樣品75份、“M1”樣品45份、“J1”樣品75份共240份用來建立判別模型，剩余20份樣品作為模型驗(yàn)證集。采用KS算法選擇模型的訓(xùn)練集180份樣品和預(yù)測(cè)集60份樣品。分別對(duì)訓(xùn)練集、預(yù)測(cè)集和驗(yàn)證集樣品編號(hào)，訓(xùn)練集中“A1”編號(hào)為X1～X34，“B1”編號(hào)為X35～X90，“M1”編號(hào)為X91～X124，“J1”編號(hào)為X125～X180。預(yù)測(cè)集中“A1”編號(hào)為Y1～Y11，“B1”編號(hào)為Y12～Y30，“M1”編號(hào)為Y31～Y41，“J1”編號(hào)為Y42～Y60。驗(yàn)證集中“A1”編號(hào)為Z1～Z5，“B1”編號(hào)為Z6～Z10，“M1”編號(hào)為Z11～Z15，“J1”編號(hào)為Z16～Z20。

1.3.3 主成分分析(principal component analysis，PCA)

PCA方法作為化學(xué)計(jì)量學(xué)中分析NIR光譜數(shù)據(jù)的常用方法，其核心思想是利用方差最大原則，對(duì)光譜數(shù)據(jù)多個(gè)自變量進(jìn)行線性擬合。這樣就可使高維的原始光譜數(shù)據(jù)最大限度的保留有效信息投影到低維空間，從而實(shí)現(xiàn)了光譜數(shù)據(jù)的降維，實(shí)現(xiàn)數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)化。在實(shí)際應(yīng)用中取前面幾個(gè)主成分，前面幾個(gè)主成分基本包含了樣品的絕大多數(shù)信息，可計(jì)算主成分的累計(jì)貢獻(xiàn)率，當(dāng)貢獻(xiàn)率達(dá)到試驗(yàn)的要求時(shí)即可，這樣可去除多余的數(shù)據(jù)，用更少的數(shù)據(jù)表達(dá)樣品更多的信息，減少模型的計(jì)算量[16-17]。

1.3.4 偏最小二乘判別分析(partial least squares-discriminate analysis，PLS-DA)

PLS-DA本質(zhì)上是一種基于特征變量的回歸方法，當(dāng)構(gòu)建分類模型區(qū)分基于同一訓(xùn)練集的不同樣品時(shí)，訓(xùn)練集中的每份樣品會(huì)被分配1個(gè)虛擬變量作為期望值，預(yù)測(cè)集的樣品分類取決于模型中的預(yù)測(cè)值Yi[18-19]。本試驗(yàn)是對(duì)4種轉(zhuǎn)基因大豆的判別，人為的把模型響應(yīng)變量期望值分別設(shè)定為：“A1”為“-1.5”；“B1”為“-0.5”；“M1”為“0.5”；“J1”為“1.5”。模型的判別閾值誤差設(shè)置為±0.5，當(dāng)模型對(duì)4種轉(zhuǎn)基因大豆預(yù)測(cè)時(shí)，由預(yù)測(cè)值Yi大小按照Yi<-1為“A1”；-1≤Yi<0為“B1”；0≤Yi≤1為“M1”；Yi>1為“J1”的區(qū)間劃分進(jìn)行歸類。

采用留一交互驗(yàn)證法來確定模型的最佳主成分?jǐn)?shù)[20-21]。以交互驗(yàn)證均方根誤差(root mean square error of cross valdarion,RMSECV)作為評(píng)價(jià)標(biāo)準(zhǔn)。表達(dá)式如公式(1)所示：

(1)

式中：n表示交互驗(yàn)證集樣本數(shù)；ci表示第i個(gè)樣本的預(yù)測(cè)值；yi表示第i個(gè)樣本的期望值。

1.3.5 誤差反向傳播人工神經(jīng)網(wǎng)絡(luò)(back-propagation artificial neural network，BP-ANN)

ANN是通過人工建立的具有自適應(yīng)、自組織、自學(xué)習(xí)特點(diǎn)的以有向圖組成拓?fù)浣Y(jié)構(gòu)的動(dòng)態(tài)系統(tǒng)。其通過正向和反向的學(xué)習(xí)和校正，實(shí)現(xiàn)輸出和輸入之間的高度的非線性映射。BP-ANN一般包含3個(gè)結(jié)構(gòu)，輸入層、隱含層和輸出層。其包含2個(gè)過程：信號(hào)的正向傳播和誤差的反向傳播[22-24]。

本試驗(yàn)為判別4種轉(zhuǎn)基因大豆，由表1可知，4種轉(zhuǎn)基因大豆NIR光譜數(shù)據(jù)經(jīng)PCA后，前7個(gè)主成分包含了原始光譜的絕大多數(shù)信息，累計(jì)貢獻(xiàn)率達(dá)到99.1%，所以可由前面7個(gè)主成分作為BP-ANN的輸入。4種轉(zhuǎn)基因大豆期望值輸出可以設(shè)為：“A1”為“-3”，“B1”為“-1”，“M1”為“1”，“J1”為“3”，模型的判別閾值誤差設(shè)置為±1。經(jīng)過對(duì)BP-ANN的多次訓(xùn)練，建立了1個(gè)輸入層(輸入節(jié)點(diǎn)為7)，2個(gè)隱含層(隱含層節(jié)點(diǎn)數(shù)分別為5和8)和1個(gè)輸出層(輸出節(jié)點(diǎn)為1)的ANN。當(dāng)模型對(duì)4種轉(zhuǎn)基因大豆預(yù)測(cè)時(shí)，由預(yù)測(cè)值Yi大小按照Yi<-2為“A1”；-2≤Yi<0為“B1”；0≤Yi≤2為“M1”；Yi>2為“J1”的區(qū)間劃分進(jìn)行歸類。

1.3.6 數(shù)據(jù)處理

數(shù)據(jù)采用MATLAB軟件自編的PCA、PLS-DA和BP-ANN建模程序處理，Origin 2018軟件繪制圖像。

2 結(jié)果與分析

2.1 轉(zhuǎn)基因大豆原始光譜圖

圖2是240份轉(zhuǎn)基因大豆樣品NIR原始光譜圖，由圖2可知光譜帶有部分噪音肉眼無法通過光譜圖來區(qū)分4種轉(zhuǎn)基因大豆，圖3是光譜經(jīng)過平滑+SNV預(yù)處理之后的光譜，由圖3可知光譜的噪音明顯減少，但仍然無法通過肉眼區(qū)分。

圖2 轉(zhuǎn)基因大豆NIR原始光譜
Fig.2 NIR spectrum of transgenic soybean

圖3 平滑+SNV處理后的轉(zhuǎn)基因大豆NIR光譜
Fig.3 NIR spectra of transgenic soybean after smoothing +SNV treatment

2.2 PCA判別

轉(zhuǎn)基因大豆光譜經(jīng)過預(yù)處理仍然有很龐大的數(shù)據(jù)。過多的冗余信息，不僅計(jì)算量大，而且還會(huì)降低模型的精度。利用PCA方法可對(duì)數(shù)據(jù)進(jìn)行降維，得到各主成分的得分矩陣。

表1為NIR光譜數(shù)據(jù)中前7個(gè)主成分累計(jì)貢獻(xiàn)率，PC1的貢獻(xiàn)率為92.5%，PC2的貢獻(xiàn)率為3.6%。圖4為PC1和PC2的得分圖，由圖4可知，4種轉(zhuǎn)基因大豆有較好的聚類，尤其是“A1”、“B1”和“J1”。但是“A1”和“J1”分布區(qū)域比較靠近，部分樣品有覆蓋，不能區(qū)分這2種轉(zhuǎn)基因大豆。

表1 主成分累計(jì)貢獻(xiàn)率
Table 1 Cumulative contribution rate of principal components

圖4 主成分得分
Fig.4 Principal component score

2.3 PLS-DA

采用PLS-DA方法對(duì)4種轉(zhuǎn)基因大豆判別，訓(xùn)練集包含180份樣品和預(yù)測(cè)集包含60份樣品。圖5為RMSECV與選擇的主成分?jǐn)?shù)的關(guān)系，設(shè)置最大主成分?jǐn)?shù)為20，當(dāng)主成分?jǐn)?shù)為4時(shí)，RMSECV最小，因此選擇4為模型的最佳主成分?jǐn)?shù)。

圖5 RMSECV與主成分?jǐn)?shù)的關(guān)系
Fig.5 Relationship between RMSECV and number of principal component

訓(xùn)練集樣品和預(yù)測(cè)集樣品的預(yù)測(cè)結(jié)果如圖6和圖7所示，訓(xùn)練集和預(yù)測(cè)集的預(yù)測(cè)結(jié)果大多數(shù)在相應(yīng)的區(qū)間內(nèi)。表2為4種轉(zhuǎn)基因大豆的判別結(jié)果，由表2可知，訓(xùn)練集中“A1”的識(shí)別率為88.2%，“B1”的識(shí)別率為96.40%，“M1”的識(shí)別率為91.1%，“J1”的識(shí)別率為96.4%。預(yù)測(cè)集中“A1”的識(shí)別率為72.7%，“B1”的識(shí)別率為94.7%，“M1”的識(shí)別率為90.9%，“J1”的識(shí)別率為89.5%。4種轉(zhuǎn)基因進(jìn)口大豆的總識(shí)別率為92.5%。使用PLS-DA方法可以較好的識(shí)別這“B1”、“M1”和“J1”這3種轉(zhuǎn)基因大豆，“A1”的識(shí)別率偏低，不能滿足現(xiàn)實(shí)要求，需要找到更適合的方法提高“A1”的識(shí)別率。

圖6 訓(xùn)練集轉(zhuǎn)基因大豆PLS-DA方法的預(yù)測(cè)結(jié)果
Fig.6 Prediction results of PLS-DA method for transgenic soybean in training set

圖7 預(yù)測(cè)集轉(zhuǎn)基因大豆PLS-DA方法的預(yù)測(cè)結(jié)果
Fig.7 Prediction results of PLS-DA method for transgenic soybean in prediction set

表2 PLS-DA和BP-ANN方法判別結(jié)果
Table 2 Identification results for PLS-DA and BP-ANN

2.4 ANN方法判別

采用BP-ANN方法對(duì)4種轉(zhuǎn)基因大豆判別，訓(xùn)練集包含180份樣品和預(yù)測(cè)集包含60份樣品，ANN經(jīng)訓(xùn)練集訓(xùn)練優(yōu)化后確定權(quán)值和閾值，預(yù)測(cè)集的60份樣品進(jìn)行驗(yàn)證。訓(xùn)練集和預(yù)測(cè)集的預(yù)測(cè)結(jié)果如圖8和圖9所示，訓(xùn)練集和預(yù)測(cè)集中“A1”、“B1”、“M1”和“J1”的預(yù)測(cè)值與期望值高度一致。由表2可知，訓(xùn)練集和預(yù)測(cè)集識(shí)別率均為100%。

2.5 PLS-DA模型與ANN方法判別模型的驗(yàn)證

取未參與建模的驗(yàn)證集20份樣品對(duì)PLS-DA方法模型與BP-ANN方法判別模型進(jìn)行驗(yàn)證。驗(yàn)證結(jié)果如圖10、圖11所示，PLS-DA方法模型識(shí)別率為90.0%，BP-ANN方法判別模型的識(shí)別率為100%?？梢姡琍LS-DA方法模型與BP-ANN方法判別模型對(duì)轉(zhuǎn)基因大豆識(shí)別率較高。

圖8 訓(xùn)練集轉(zhuǎn)基因大豆BP-ANN方法的預(yù)測(cè)結(jié)果
Fig.8 Prediction results of BP-ANN method for transgenic soybean in training set

圖9 預(yù)測(cè)集轉(zhuǎn)基因大豆BP-ANN方法的預(yù)測(cè)結(jié)果
Fig.9 Prediction results of BP-ANN method for transgenic soybean in prediction set

圖10 驗(yàn)證集轉(zhuǎn)基因大豆PLS-DA方法的預(yù)測(cè)結(jié)果
Fig.10 Prediction results of PLS-DA method for transgenic soybean in validation set

圖11 驗(yàn)證集轉(zhuǎn)基因大豆BP-ANN方法的預(yù)測(cè)結(jié)果
Fig.11 Prediction results of BP-ANN method for transgenic soybean in validation set

3 結(jié)論

采用NIR光譜結(jié)合化學(xué)計(jì)量學(xué)對(duì)4種轉(zhuǎn)基因大豆進(jìn)行判別分析，利用平滑+SNV方法預(yù)處理原始光譜數(shù)據(jù)，KS算法劃分訓(xùn)練集樣品和預(yù)測(cè)集樣品，PCA、PLS-DA和BP-ANN方法分析預(yù)處理后的NIR光譜數(shù)據(jù)。試驗(yàn)結(jié)果顯示平滑+SNV的預(yù)處理方法能有效減少NIR光譜的噪音；PCA方法能判別出4種轉(zhuǎn)基因大豆中的3種，阿根廷轉(zhuǎn)基因大豆和加拿大轉(zhuǎn)基因大豆不能同時(shí)判別；PLS-DA方法對(duì)預(yù)測(cè)集轉(zhuǎn)基因大豆的判別正確率為88.3%；BP-ANN方法能夠準(zhǔn)確的判別4種轉(zhuǎn)基因大豆，判別正確率為100%；并用未參與建模的4種轉(zhuǎn)基因大豆對(duì)PLS-DA方法模型和BP-ANN方法模型進(jìn)行驗(yàn)證，驗(yàn)證集中PLS-DA方法模型判別正確率為90.0%，BP-ANN方法模型判別正確率為100%。本試驗(yàn)雖然只選取了4種轉(zhuǎn)基因大豆進(jìn)行建模判別，但是試驗(yàn)結(jié)果表明采用NIR光譜結(jié)合PLS-DA和BP-ANN方法對(duì)轉(zhuǎn)基因大豆產(chǎn)地朔源是可行的?？蔀槲覈?guó)相關(guān)部門對(duì)轉(zhuǎn)基因大豆的產(chǎn)地朔源提供部分方法。