龔 浩，張莉莉，陳富榮，林麗霞，陳意君，張 樂，孫春 蓮，孫 鍵

（1.惠州學(xué)院生命科學(xué)學(xué)院，廣東 惠州 516007；2.惠州學(xué)院經(jīng)濟(jì)管理學(xué)院，廣東 惠州 516007）

【研究意義】茶樹〔Camelliasinensis(L.)O.Kuntze〕屬山茶科山茶屬多年生常綠木本植物，原產(chǎn)于熱帶及亞熱帶，是一種喜暖喜濕的葉用植物，其嫩葉經(jīng)過加工后即為茶葉。茶葉具有防輻射、提神醒腦、利尿、助消化、減肥和預(yù)防疾病的作用，因此茶葉的飲用及流傳從古至今都極受重視，是中華民族的舉國之飲、世界三大飲品之首［1］。但是茶樹異花傳粉和長期自交不育的特性，使茶樹高度雜合、親緣關(guān)系復(fù)雜，茶葉品種難以輕易分辨、分類標(biāo)準(zhǔn)難以統(tǒng)一、鑒別結(jié)果有誤差等，這就需要對(duì)茶葉不同品種進(jìn)行區(qū)分和產(chǎn)地溯源。SSR 廣泛應(yīng)用于植物基因定位和QTL 分析、DNA 指紋和品種鑒定［2］、種質(zhì)資源保存和利用、系譜分析以及標(biāo)記輔助育種，通常呈共顯性遺傳，其多態(tài)位點(diǎn)豐富，實(shí)驗(yàn)操作簡單易行［3］?；谏疃壬窠?jīng)網(wǎng)絡(luò)的簡單重復(fù)序列標(biāo)記對(duì)茶葉產(chǎn)地的溯源研究不僅有利于茶葉的分類和產(chǎn)地溯源［4］，還能為其他植物分類提供參考。

【前人研究進(jìn)展】目前已發(fā)表相關(guān)論文的茶樹測序群體一般為100～200 個(gè)樣本，過于零散且群體覆蓋性和代表性較弱，無法用于深度的群體遺傳分析［5］。目前，國內(nèi)對(duì)茶葉品種、產(chǎn)地、產(chǎn)季、年份和等級(jí)等真實(shí)屬性的鑒別還主要停留在傳統(tǒng)的理化分析與感官評(píng)定相結(jié)合的水平上，例如GB/T 19598-2006《地理標(biāo)志產(chǎn)品 安溪鐵觀音》中評(píng)價(jià)標(biāo)準(zhǔn)是以感官為主，輔以部分理化檢測。一方面，具有感官評(píng)定能力的專家非常少，特別是面對(duì)我國品種繁多的茶葉，具有特定品種茶葉感官評(píng)定能力的專家更為稀缺；另一方面，人的感官靈敏度容易受到外界因素的干擾而改變，采用感官評(píng)價(jià)方法受人為主觀影響很大，可操作性較差，且目前還沒有明確且易于實(shí)現(xiàn)的評(píng)定指標(biāo)或參數(shù)，易造成判定結(jié)果的偏差［6］?；诖水a(chǎn)生電子鼻來分類茶葉，利用氣敏傳感器陣列對(duì)揮發(fā)性氣味物質(zhì)響應(yīng)，使氣味成為量化指標(biāo)的新技術(shù)手段，具有檢測時(shí)間短、樣品預(yù)處理簡單、檢測結(jié)果可靠等優(yōu)點(diǎn)［7］，可以高效、快速、無損檢測不同種類的食品，可應(yīng)用于茶葉貯藏時(shí)間［8］、加工方式、品質(zhì)［9］和等級(jí)［10］等檢測。但這種方法會(huì)因?yàn)闄z驗(yàn)材料部位不同而出現(xiàn)較大的結(jié)果誤差。

【本研究切入點(diǎn)】近年來，分子生物學(xué)技術(shù)和生物信息學(xué)的發(fā)展有力地推動(dòng)了DNA 分子標(biāo)記的研究。與形態(tài)學(xué)標(biāo)記、細(xì)胞標(biāo)記、生化標(biāo)記等相比，DNA 分子標(biāo)記技術(shù)不易受外界環(huán)境及個(gè)體本身的影響，具有結(jié)果準(zhǔn)確、信息量大、檢測簡單、重復(fù)性及穩(wěn)定性較好等優(yōu)點(diǎn)［11］。DNA 分子標(biāo)記技術(shù)在植物分類學(xué)［12］、遺傳多樣性分析［13］、遺傳圖譜構(gòu)建［14］和輔助育種［15］等方面的研究廣為應(yīng)用，但在茶樹種質(zhì)資源方面的研究應(yīng)用較少，主要集中在遺傳多樣性及特異標(biāo)記方面?！緮M解決的關(guān)鍵問題】本研究旨在解決茶葉主成分分析、茶葉產(chǎn)地溯源及品種鑒定、DNN 模型構(gòu)建等問題。

1 材料與方法

1.1 試驗(yàn)材料

本研究根據(jù)Accession No.PRJNA595795 和PRJNA562973，從NCBI database 中下載323 份茶葉的轉(zhuǎn)錄組數(shù)據(jù)，其中來自福建、云南、浙江、湖南省的茶葉分別有130、96、54、33 份，其余10 份屬于外類群樣本即研究類群之外親緣關(guān)系最近的物種，這10 個(gè)外類群為全國收集的茶梅CamelliasasanquaThunb（表1）。

表1 茶葉樣本來源屬地統(tǒng)計(jì)Table 1 Statistics on the origin and locality of tea samples

1.2 試驗(yàn)方法

1.2.1 鑒定SSR 標(biāo)記位點(diǎn) 本研究先從323 份樣本中獲得樣本數(shù)據(jù)，使用PSR 軟件（Polymorphic SSR retrieval,PMID:26428628），鑒定茶葉參考基因組（Tea treeCamelliasinensis,舒茶早）中所有可能的SSR 標(biāo)記位點(diǎn)。首先利用PSR 軟件，設(shè)置參數(shù)支持的reads 總數(shù)大于5，同時(shí)支持reads 的比例大于10%，其他參數(shù)均為默認(rèn)參數(shù)；再過濾單個(gè)位點(diǎn)缺失率較大的位點(diǎn)；最后進(jìn)行線性回歸分析，并結(jié)合不同SSR 位點(diǎn)的相關(guān)性，保留最終的SSR 位點(diǎn)［16］。

1.2.2 主成分分析 本研究將323 個(gè)茶葉樣本進(jìn)行樣本間SSR 序列的相互比對(duì)，計(jì)算每個(gè)樣本與其他樣本的差異度，再基于樣本間的差異度計(jì)算323 個(gè)樣本的基因差異矩陣［17］。利用PCA 對(duì)323 個(gè)樣本基因差異矩陣進(jìn)行分析，然后使用R語言中的read.table 函數(shù)讀入數(shù)據(jù)、ovun.sample 函數(shù)清理處理數(shù)據(jù)，最后利用內(nèi)置函數(shù)princomp 進(jìn)行PCA。

1.2.3 構(gòu)建整體樣本的進(jìn)化樹 先用perl 的自編腳本獲得所有個(gè)體的明氏距離矩陣，然后用PHYLIP 的neighbour 模塊構(gòu)建原始進(jìn)化樹，再用dendroscope 對(duì)進(jìn)化樹進(jìn)行展示和修飾。

1.2.4 建立及優(yōu)化模型 本研究使用Matlab 軟件的神經(jīng)網(wǎng)絡(luò)工具，建立線性回歸模型、隨機(jī)森林模型和DNN 模型。在建立線性回歸模型的過程中，根據(jù)樣本數(shù)據(jù)集，分別生成x 矩陣和y 矩陣，利用線性回歸代碼建模，并使用其模型進(jìn)行預(yù)測。在建立隨機(jī)森林模型的過程中，先將整體樣本讀到內(nèi)存中，按照8∶2 的比例分為80%的訓(xùn)練集、20%的測試集；然后將訓(xùn)練集的樣本先分詞，再轉(zhuǎn)換為詞向量；接著將訓(xùn)練集的樣本和標(biāo)簽統(tǒng)一傳入算法中，得到擬合后的模型；繼而將測試集的樣本先分詞，再得到詞向量；最終把測試集得出的詞向量添加到擬合后的模型中，得出結(jié)果并將結(jié)果轉(zhuǎn)換為準(zhǔn)確率的形式。在建立DNN 模型的過程中，本研究通過下載WeightWatcher 安裝包，導(dǎo)入樣本數(shù)據(jù)，利用神經(jīng)網(wǎng)絡(luò)代碼直接預(yù)測準(zhǔn)確率。選取準(zhǔn)確率最高的深度神經(jīng)網(wǎng)絡(luò)模型作進(jìn)一步優(yōu)化。

2 結(jié)果與分析

2.1 SSR 標(biāo)記位點(diǎn)鑒定結(jié)果

首先需要初步鑒定SSR 位點(diǎn)，通過PSR 軟件，從茶葉參考基因組數(shù)據(jù)庫中獲得所有可能的SSR 標(biāo)記，最終得到3 668 個(gè)標(biāo)記位點(diǎn)，其中，比對(duì)到染色體上的位點(diǎn)有3 304 個(gè)（表2）［18］。SSR 標(biāo)記位點(diǎn)的鑒定：利用PSR 軟件，經(jīng)過篩選后得到2 924 個(gè)位點(diǎn)；過濾單個(gè)位點(diǎn)缺失率大于20%的位點(diǎn)后，獲得2 155 個(gè)多態(tài)性位點(diǎn)；通過線性回歸分析，篩選在不同省份特異性存在的位點(diǎn)（P＜0.001），獲得700 個(gè)位點(diǎn)；結(jié)合不同SSR位點(diǎn)的相關(guān)性，在兩個(gè)及其相關(guān)的位點(diǎn)中只保留差異性較大的位點(diǎn)，最終獲得54 個(gè)SSR 位點(diǎn)。

表2 各染色體中含有SSR 位點(diǎn)數(shù)目的統(tǒng)計(jì)Table 2 Statistics on the number of SSR loci contained in each chromosome

2.2 不同來源茶葉樣本PCA 結(jié)果

如圖1 所示，圖中每個(gè)點(diǎn)代表1 個(gè)樣本，兩點(diǎn)距離代表茶葉樣品受主成分影響下的相似性距離。全部樣本的PCA 結(jié)果表明，Dim1（7.6%）表示第一主成分貢獻(xiàn)率為7.6%，Dim2（4.3%）表示第二主成分貢獻(xiàn)率為4.3%，即前兩個(gè)主成分的累計(jì)貢獻(xiàn)率為 11.9%（圖1A）；外類群與福建、湖南、云南、浙江4 省茶葉樣本差異顯著，部分與云南省樣品個(gè)體相近。本研究通過對(duì)4 省份數(shù)據(jù)進(jìn)行PCA 來做進(jìn)一步判斷。根據(jù)4 個(gè)省份間的PCA 結(jié)果（圖 1B），并排除外類群的影響，可以發(fā)現(xiàn)不同省份間的整體差異較明顯，而4 個(gè)省份內(nèi)個(gè)體相對(duì)聚集。其中，云南省內(nèi)的個(gè)體較其他省份差異大；福建、浙江、湖南的樣本分別聚集，這表明福建、浙江、湖南3 個(gè)省份間茶葉差異顯著，但有少量交叉，具有一定相似的遺傳結(jié)構(gòu)特性，3 個(gè)省份間的親緣關(guān)系較近。其親緣關(guān)系遠(yuǎn)近與地理來源并不呈現(xiàn)一致性，原因可能與茶葉人工馴化程度有關(guān)。PCA 也存在一定的不足之處，簡單的PCA 只能解釋部分個(gè)體的產(chǎn)地溯源問題，若要進(jìn)一步研究溯源問題，則還需要用其他方法，如構(gòu)建進(jìn)化樹、神經(jīng)網(wǎng)絡(luò)模型等方法，來進(jìn)一步解釋和驗(yàn)證交叉?zhèn)€體的溯源問題。

圖1 主成分分析結(jié)果Fig.1 Principal component analysis results

2.3 進(jìn)化樹構(gòu)建結(jié)果

從以上PCA 分析結(jié)果可以看出，不同省份的個(gè)體分別聚集，差異較為顯著，但也有少量的交叉，其中福建主要與浙江、云南鄰近，湖南與云南較近，外類群主要分布在云南附近，而云南個(gè)體分類較其他省份分散，由此構(gòu)建不同省份茶葉的進(jìn)化樹（圖2），其結(jié)果與PCA 結(jié)果相似。

2.4 模型構(gòu)建與優(yōu)化結(jié)果

2.4.1 不同模型預(yù)測結(jié)果 本研究利用3 種不同的模型對(duì)54 個(gè)SSR 分子標(biāo)記矩陣構(gòu)建模型，再初步鑒定不同模型的差異。通過線性回歸模型（81%）、隨機(jī)森林模型（77%）及DNN 模型（86%）對(duì)54 個(gè)SSR marker 矩形構(gòu)建模型，發(fā)現(xiàn)深度神經(jīng)網(wǎng)絡(luò)模型準(zhǔn)確率最高、為86%，故選擇DNN 模型進(jìn)行預(yù)測［19］。

2.4.2 DNN 模型的優(yōu)化結(jié)果 本研究利用Matlab軟件的神經(jīng)網(wǎng)絡(luò)工具對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行建模。使用54 個(gè)SSR 和323 個(gè)樣本，構(gòu)建預(yù)測模型，再用Tensorflow2.0 優(yōu)化DNN 模型的一次訓(xùn)練樣本個(gè)數(shù)（Batch size)、訓(xùn)練次數(shù)（Step size)、隱藏層層數(shù)和每層節(jié)點(diǎn)數(shù)4 個(gè)參數(shù)。

將323 份樣本中除了外類群以外的數(shù)據(jù)分成訓(xùn)練集、測試集和驗(yàn)證集3 個(gè)部分，其中訓(xùn)練集、測試集、驗(yàn)證集的測試比例分別為0.8、0.1、0.1，即訓(xùn)練集273 份，測試集20 份，驗(yàn)證集20 份。先用訓(xùn)練集訓(xùn)練模型，再用測試集進(jìn)行最后優(yōu)化，并使用驗(yàn)證集對(duì)優(yōu)化后的模型進(jìn)行驗(yàn)證。

2.4.3 參數(shù)Batch size 和Step 的優(yōu)化 本研究通過對(duì)參數(shù)Batch size 和Step 進(jìn)行優(yōu)化，測試不同參數(shù)對(duì)準(zhǔn)確率的影響。對(duì)每次訓(xùn)練選取的Batch size 分別設(shè)為150、200、250、300，而迭代的次數(shù)Step 分別為5 000、10 000、15 000、20 000、25 000、30 000。理論上Step 越高模型準(zhǔn)確率就越高，但Step 過高會(huì)導(dǎo)致模型過度擬合。通過對(duì)測試集10 次重復(fù)驗(yàn)證，發(fā)現(xiàn)參數(shù)Batch size 為150 和Step 為20 000 綜合起來表現(xiàn)效果最好（表3、表4）。

表3 測試集和驗(yàn)證集最優(yōu)準(zhǔn)確率Table 3 Optimal accuracy of the test and validation set

表4 測試集和驗(yàn)證集平均準(zhǔn)確率Table 4 Average accuracy of the test set and validation set

2.4.4 隱藏層層數(shù)和每層節(jié)點(diǎn)數(shù)的優(yōu)化（1）隱藏層層數(shù)的優(yōu)化：利用不同的隨機(jī)參數(shù)模擬2～4 層神經(jīng)網(wǎng)絡(luò)的測試集和驗(yàn)證集的準(zhǔn)確率。經(jīng)對(duì)比，發(fā)現(xiàn)神經(jīng)網(wǎng)絡(luò)為2 層時(shí)驗(yàn)證集和測試集的準(zhǔn)確率最高，約95%（圖 2 A）。（2）每層節(jié)點(diǎn)數(shù)的優(yōu)化：確定隱藏層為2 層后，分別產(chǎn)生25～150間隔為5 的26 個(gè)可能節(jié)點(diǎn)數(shù)，隱藏層的兩層網(wǎng)絡(luò)組合一起是26×26 共676 個(gè)組合的矩陣，檢測不同參數(shù)對(duì)應(yīng)的準(zhǔn)確率，每個(gè)組合進(jìn)行10 次重復(fù)。

然后按以下打分規(guī)則對(duì)最優(yōu)準(zhǔn)確率進(jìn)行確定，通過統(tǒng)計(jì)不同指標(biāo)對(duì)所有組合進(jìn)行打分，每一種指標(biāo)都能進(jìn)10%得1 分：測試集和驗(yàn)證集準(zhǔn)確率的平均值；驗(yàn)證集準(zhǔn)確率的平均值；最優(yōu)驗(yàn)證準(zhǔn)確率。最后統(tǒng)計(jì)2 分以上的次數(shù)（圖3），圖3A 為在最優(yōu)Batch size 和Step size 時(shí)不同神經(jīng)層數(shù)的柱狀圖，對(duì)比發(fā)現(xiàn)神經(jīng)網(wǎng)絡(luò)為2 層時(shí)驗(yàn)證集和測試集的準(zhǔn)確率最高；圖3B、C、D、E 為2 層隱藏層神經(jīng)網(wǎng)絡(luò)參數(shù)的優(yōu)化，其中B 為不同維度模擬的自測數(shù)據(jù)的平均準(zhǔn)確率，C 為不同維度模擬的驗(yàn)證數(shù)據(jù)的平均準(zhǔn)確率，D 為不同維度模擬的自測數(shù)據(jù)的最優(yōu)準(zhǔn)確率。綜合準(zhǔn)確率方差等因素，本研究選擇隱藏層第一層95、第二層40 的模型為最優(yōu)模型，其中自測集的平均準(zhǔn)確率95%，自測集合驗(yàn)證平均值準(zhǔn)確率89%，驗(yàn)證集的平均準(zhǔn)確率75%以上，最優(yōu)準(zhǔn)確率為100%。

圖3 深度神經(jīng)網(wǎng)絡(luò)層數(shù)和節(jié)點(diǎn)數(shù)的參數(shù)優(yōu)化結(jié)果Fig.3 Oprtimization results of layer number and node for each layer for the Deep Neural Network

3 討論

我國是茶葉消費(fèi)大國，隨著人們生活水平的提高，消費(fèi)者對(duì)茶品質(zhì)的要求也日益提高。為了確定茶葉的真實(shí)產(chǎn)地，研究者運(yùn)用各種方法進(jìn)行研究。目前，生物信息學(xué)在基因測序分析中發(fā)揮著舉足輕重的作用，國內(nèi)主要是以實(shí)驗(yàn)為基礎(chǔ)，通過測定農(nóng)產(chǎn)物及其土壤中的礦質(zhì)元素，再進(jìn)行相關(guān)性分析、聚類分析、主成分分析等多種統(tǒng)計(jì)分析方法，進(jìn)而對(duì)農(nóng)產(chǎn)品進(jìn)行溯源分析［20］。本研究主要以生物信息學(xué)為基礎(chǔ)，通過分析相關(guān)的基因位點(diǎn)，構(gòu)建模型并進(jìn)行優(yōu)化，最終對(duì)茶葉溯源進(jìn)行分析。

本研究通過基因組的SSR 位點(diǎn)進(jìn)行基因數(shù)據(jù)分析。SSR 作為第二代分子標(biāo)記，具有重復(fù)性好、多態(tài)性高、變異豐富、呈共顯性且廣泛分布于植物基因組等優(yōu)點(diǎn)，已被廣泛應(yīng)用于高粱、大麥、小麥、青稞等作物遺傳多樣性分析和基因研究［21］。與SNP 標(biāo)記相比，SSR 標(biāo)記的優(yōu)勢(shì)是成本低、試驗(yàn)技術(shù)簡單［22］。本研究先利用PSR 軟件從茶葉參考基因組中鑒定所有可能的SSR 標(biāo)記位點(diǎn)，再比對(duì)到染色體上，利用PSR 設(shè)置參數(shù)支持的reads 總數(shù)大于5 同時(shí)支持reads 的比例大于10%［23］，得到樣本后進(jìn)行位點(diǎn)篩選，最終獲得54 個(gè)SSR 位點(diǎn)；再利用3 種不同的模型對(duì)54 個(gè)SSR 分子標(biāo)記矩陣構(gòu)建模型，初步鑒定不同模型的差異［24］；選擇準(zhǔn)確率最高的神經(jīng)網(wǎng)絡(luò)模型，進(jìn)行人工神經(jīng)網(wǎng)絡(luò)模型的優(yōu)化和參數(shù)選擇、Batch size 和Step size 的優(yōu)化、隱藏層數(shù)目和每層節(jié)點(diǎn)數(shù)優(yōu)化、2 層隱藏層神經(jīng)網(wǎng)絡(luò)參數(shù)的優(yōu)化，最后選擇準(zhǔn)確率在95%左右最優(yōu)的2 層神經(jīng)網(wǎng)絡(luò)模型［25］。

在研究地理溯源領(lǐng)域中，大部分研究都是利用分子標(biāo)記或化學(xué)標(biāo)記構(gòu)建變異圖譜，然后查看變異圖譜的相似性來進(jìn)行溯源。本研究使用深度學(xué)習(xí)預(yù)測方法，在研究產(chǎn)地溯源領(lǐng)域使用量較少，主要通過建立樣本的基因差異矩陣，使用PCA 分析323 個(gè)樣本間的差異度，結(jié)果非常直觀。通過分析圖片，發(fā)現(xiàn)外類群與福建、湖南、云南、浙江4 省份間的差異顯著，而各省份內(nèi)個(gè)體相對(duì)聚集，其中云南省內(nèi)的個(gè)體差異較其他省份大；4省份有部分材料重疊在一起，表明不同省份的部分茶葉也具有一定的遺傳相似性［26］。構(gòu)建整體樣本的進(jìn)化樹，結(jié)果表明不同省份的茶葉個(gè)體分別聚集，差異顯著，但也有少量交叉，此結(jié)果與PCA 結(jié)果相似［27］。

本研究只研究福建、湖南、云南、浙江4 省份和10 個(gè)外類群共323 個(gè)樣本，茶葉轉(zhuǎn)錄組數(shù)據(jù)存在樣本量少的局限性，后續(xù)需要增加樣本容量，對(duì)茶葉溯源作進(jìn)一步研究。

4 結(jié)論

本研究對(duì)來自湖南、云南、福建和浙江省的313 個(gè)茶葉樣本的來源屬地及10 個(gè)外類群關(guān)系進(jìn)行研究，以篩選出的54 個(gè)高質(zhì)量的SSR 位點(diǎn)為基礎(chǔ)，對(duì)樣本進(jìn)行主成分分析，并通過3 種不同的分類模型比對(duì)及優(yōu)化，得出2 層神經(jīng)網(wǎng)絡(luò)模型對(duì)茶葉分析效果最佳，準(zhǔn)確率約95%。本研究構(gòu)建的分類模型也可以用于其他物種重測序數(shù)據(jù)的屬地來源鑒定。