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

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

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

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

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

1 材料與方法

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

試驗(yàn)選取轉(zhuǎn)基因大豆分別為阿根廷轉(zhuǎn)基因大豆MON89788品系(“A1”)、巴西轉(zhuǎn)基因大豆MON89788品系(“B1”)、美國轉(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光譜儀，檢測器為Si和InGaAs，光譜掃描范圍400～2 600 nm，北京偉創(chuàng)英圖科技有限公司；ME204E電子天平，梅特勒-托利多有限公司；臺式真空干燥箱DZF-6050，上海捷呈實(shí)驗(yàn)儀器有限公司。

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

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

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

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

1.3.2 樣品的劃分和編號

Kennard-Stone(KS)算法，是在總樣本中選出訓(xùn)練集樣品，首先選擇歐氏距離最遠(yuǎn)的2個樣品進(jìn)入訓(xùn)練集，其后通過計算剩下的每1個樣品到訓(xùn)練集內(nèi)每1個已知樣品的歐式距離，找到擁有最大最小距離的待選樣品放入訓(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ù)測集60份樣品。分別對訓(xùn)練集、預(yù)測集和驗(yàn)證集樣品編號，訓(xùn)練集中“A1”編號為X1～X34，“B1”編號為X35～X90，“M1”編號為X91～X124，“J1”編號為X125～X180。預(yù)測集中“A1”編號為Y1～Y11，“B1”編號為Y12～Y30，“M1”編號為Y31～Y41，“J1”編號為Y42～Y60。驗(yàn)證集中“A1”編號為Z1～Z5，“B1”編號為Z6～Z10，“M1”編號為Z11～Z15，“J1”編號為Z16～Z20。

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

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

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

PLS-DA本質(zhì)上是一種基于特征變量的回歸方法，當(dāng)構(gòu)建分類模型區(qū)分基于同一訓(xùn)練集的不同樣品時，訓(xùn)練集中的每份樣品會被分配1個虛擬變量作為期望值，預(yù)測集的樣品分類取決于模型中的預(yù)測值Yi[18-19]。本試驗(yàn)是對4種轉(zhuǎn)基因大豆的判別，人為的把模型響應(yīng)變量期望值分別設(shè)定為：“A1”為“-1.5”；“B1”為“-0.5”；“M1”為“0.5”；“J1”為“1.5”。模型的判別閾值誤差設(shè)置為±0.5，當(dāng)模型對4種轉(zhuǎn)基因大豆預(yù)測時，由預(yù)測值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)作為評價標(biāo)準(zhǔn)。表達(dá)式如公式(1)所示：

(1)

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

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)的動態(tài)系統(tǒng)。其通過正向和反向的學(xué)習(xí)和校正，實(shí)現(xiàn)輸出和輸入之間的高度的非線性映射。BP-ANN一般包含3個結(jié)構(gòu)，輸入層、隱含層和輸出層。其包含2個過程：信號的正向傳播和誤差的反向傳播[22-24]。

本試驗(yàn)為判別4種轉(zhuǎn)基因大豆，由表1可知，4種轉(zhuǎn)基因大豆NIR光譜數(shù)據(jù)經(jīng)PCA后，前7個主成分包含了原始光譜的絕大多數(shù)信息，累計貢獻(xiàn)率達(dá)到99.1%，所以可由前面7個主成分作為BP-ANN的輸入。4種轉(zhuǎn)基因大豆期望值輸出可以設(shè)為：“A1”為“-3”，“B1”為“-1”，“M1”為“1”，“J1”為“3”，模型的判別閾值誤差設(shè)置為±1。經(jīng)過對BP-ANN的多次訓(xùn)練，建立了1個輸入層(輸入節(jié)點(diǎn)為7)，2個隱含層(隱含層節(jié)點(diǎn)數(shù)分別為5和8)和1個輸出層(輸出節(jié)點(diǎn)為1)的ANN。當(dāng)模型對4種轉(zhuǎn)基因大豆預(yù)測時，由預(yù)測值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ù)。過多的冗余信息，不僅計算量大，而且還會降低模型的精度。利用PCA方法可對數(shù)據(jù)進(jìn)行降維，得到各主成分的得分矩陣。

表1為NIR光譜數(shù)據(jù)中前7個主成分累計貢獻(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 主成分累計貢獻(xiàn)率
Table 1 Cumulative contribution rate of principal components

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

2.3 PLS-DA

采用PLS-DA方法對4種轉(zhuǎn)基因大豆判別，訓(xùn)練集包含180份樣品和預(yù)測集包含60份樣品。圖5為RMSECV與選擇的主成分?jǐn)?shù)的關(guān)系，設(shè)置最大主成分?jǐn)?shù)為20，當(dāng)主成分?jǐn)?shù)為4時，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ù)測集樣品的預(yù)測結(jié)果如圖6和圖7所示，訓(xùn)練集和預(yù)測集的預(yù)測結(jié)果大多數(shù)在相應(yīng)的區(qū)間內(nèi)。表2為4種轉(zhuǎn)基因大豆的判別結(jié)果，由表2可知，訓(xùn)練集中“A1”的識別率為88.2%，“B1”的識別率為96.40%，“M1”的識別率為91.1%，“J1”的識別率為96.4%。預(yù)測集中“A1”的識別率為72.7%，“B1”的識別率為94.7%，“M1”的識別率為90.9%，“J1”的識別率為89.5%。4種轉(zhuǎn)基因進(jìn)口大豆的總識別率為92.5%。使用PLS-DA方法可以較好的識別這“B1”、“M1”和“J1”這3種轉(zhuǎn)基因大豆，“A1”的識別率偏低，不能滿足現(xiàn)實(shí)要求，需要找到更適合的方法提高“A1”的識別率。

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

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

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

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

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

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

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

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

3 結(jié)論

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