基于多平臺基因表達數(shù)據(jù)的水稻干旱和鹽脅迫相關(guān)基因預(yù)測

劉亞文 張紅燕,,* 曹 丹 李蘭芝

1 湖南農(nóng)業(yè)大學(xué)信息與智能科學(xué)技術(shù)學(xué)院, 湖南長沙 410128; 2 湖南農(nóng)業(yè)大學(xué) / 湖南省農(nóng)業(yè)大數(shù)據(jù)分析與決策工程技術(shù)研究中心, 湖南長沙 410128

水稻在生長發(fā)育的過程中受到干旱、高鹽等非生物逆境因素脅迫時, 易導(dǎo)致大面積的減產(chǎn)、品質(zhì)下降甚至壞死[1], 提高其逆境抗性將增加農(nóng)業(yè)產(chǎn)量并擴大適宜耕種面積, 緩解人口壓力。水稻的逆境抗性受多基因控制, 基于基因組學(xué)數(shù)據(jù)挖掘水稻非生物脅迫相關(guān)基因, 對培育抗逆水稻新品種具有重大意義。近年來, 隨著大規(guī)?；虮磉_水平測量技術(shù)的發(fā)展, 基于雜交原理的基因芯片技術(shù)[2]和基于高通量測序技術(shù)的RNA-seq[3], 被眾多學(xué)者用于植物脅迫響應(yīng)基因的挖掘研究中[4]。

然而, 大多實驗組測序樣本有限, 僅單一地從與水稻某非生物脅迫相關(guān)的單個實驗組測序數(shù)據(jù)來挖掘脅迫相關(guān)基因, 結(jié)果不穩(wěn)定, 也很難讓人信服[5]。當(dāng)前公共數(shù)據(jù)庫中積累的大量水稻脅迫相關(guān)基因芯片和RNA-seq表達數(shù)據(jù), 為多平臺數(shù)據(jù)分析提供了研究空間。研究表明, 融合多平臺數(shù)據(jù)能夠提高基因表達分析的準(zhǔn)確性和可靠性, 多平臺表達數(shù)據(jù)的整合分析成為水稻非生物脅迫相關(guān)基因預(yù)測研究的趨勢[6]。當(dāng)前, 多平臺基因表達數(shù)據(jù)的融合通?？煞譃閮深? (1) 基于輸出層面融合的元分析法。它通過對多個研究結(jié)果進行合并匯總, 增大樣本總量, 提高檢測準(zhǔn)確率和統(tǒng)計分析結(jié)果的一致性[7]。(2) 基于原始數(shù)據(jù)融合的數(shù)據(jù)轉(zhuǎn)換法。它先通過把不同平臺基因表達數(shù)據(jù)按一定規(guī)則轉(zhuǎn)換到同一個數(shù)據(jù)范圍內(nèi), 再將轉(zhuǎn)換后的多個平臺實驗數(shù)據(jù)直接合并成一個表達數(shù)據(jù)矩陣, 以此來增加樣本數(shù)目緩解“高維數(shù)、小樣本”維數(shù)災(zāi)難問題[8]。綜合考慮單個實驗組水稻測序的小樣本, 及水稻基因芯片數(shù)據(jù)與RNA-seq數(shù)據(jù)之間尺度與維度的差異, 本研究首先對同一脅迫相關(guān)的多個基因芯片數(shù)據(jù)或多個RNA-seq數(shù)據(jù)分別采取數(shù)據(jù)轉(zhuǎn)換法融合, 再分別基于融合后的基因芯片表達數(shù)據(jù)集和RNA-seq表達數(shù)據(jù)集進行脅迫相關(guān)基因挖掘, 最后將二者的結(jié)果實施元分析, 獲取最終的脅迫響應(yīng)基因。

為了有效利用多平臺基因表達數(shù)據(jù), 本文選用加權(quán)基因共表達網(wǎng)絡(luò)分析(Weighted Gene Co-expression Network Analysis, WGCNA)法來挖掘關(guān)鍵基因。WGCNA利用基因表達數(shù)據(jù)構(gòu)造協(xié)同表達的基因模塊, 并根據(jù)基因模塊與表型的關(guān)聯(lián)性以及基因模塊的內(nèi)連性來鑒定關(guān)鍵基因[9], 其基本假定是“表達模式相似的基因功能相似”。它可將表達模式相似的基因進行聚類, 并分析模塊與特定性狀或表型之間的關(guān)聯(lián)關(guān)系, 因此在作物的干旱脅迫、鹽脅迫等非生物脅迫相關(guān)基因的挖掘研究中被廣泛應(yīng)用。例如李旭凱等[10]利用WGCNA挖掘到2599個與水稻冷脅迫、干旱脅迫和鹽脅迫都相關(guān)的基因, 并預(yù)測出25個抗逆關(guān)鍵基因; Zhu等[11]通過對轉(zhuǎn)錄組數(shù)據(jù)進行WGCNA分析, 確定了水稻鹽脅迫響應(yīng)核心差異基因和模塊; Lv等[12]以轉(zhuǎn)錄組數(shù)據(jù)為基礎(chǔ)進行WGCNA分析, 預(yù)測了各模塊重要的Hub差異基因和調(diào)控水稻干旱應(yīng)答基因表達的主要轉(zhuǎn)錄因子;Hopper等[13]使用時間序列轉(zhuǎn)錄方法結(jié)合WGCNA網(wǎng)絡(luò)分析, 為葡萄耐旱性研究提供了候選基因; 秦天元等[14]使用WGCNA挖掘馬鈴薯根系抗旱核心基因,并進一步利用RT-qPCR驗證出挖掘到的核心基因確實響應(yīng)干旱脅迫。經(jīng)典的WGCNA 以Pearson相關(guān)系數(shù)度量2個基因表達量間的線性相似性(記為WGCNA-P), 但無法捕獲基因間可能廣泛存在的非線性關(guān)聯(lián)。Reshef等[15]學(xué)者基于信息論中的互信息理論提出了一種可度量兩變量非線性相關(guān)性的普適性測度最大信息系數(shù)(Maximal Information Coefficient, MIC), 論文提出以MIC作為相似性度量替代WGCNA中的Pearson相關(guān)系數(shù)來構(gòu)建基因共表達網(wǎng)絡(luò)(記為WGCNA-MIC), 以捕捉基因間的非線性關(guān)聯(lián)。同時, 考慮到特定線性情形下MIC的統(tǒng)計功效[16]不如Pearson相關(guān)系數(shù), 所以本研究對同一數(shù)據(jù)集分別基于WGCNA-P和WGCNA-MIC兩種方法來構(gòu)建基因共表達網(wǎng)絡(luò), 并對各自獲取的Hub基因集進行整合分析。

綜上, 本研究以多平臺水稻非生物脅迫(以干旱和鹽脅迫為代表)相關(guān)的基因芯片數(shù)據(jù)和RNA-seq數(shù)據(jù)為研究對象, 分別以WGCNA-P和WGCNAMIC挖掘脅迫相關(guān)Hub基因, 進而對同一脅迫不同平臺數(shù)據(jù)使用以上2種網(wǎng)絡(luò)分析法得到的Hub基因進行整合分析, 得到最終的脅迫相關(guān)Hub基因集。最后, 從預(yù)測性能、基因功能富集分析、文獻報道和互作網(wǎng)絡(luò)分析等多角度解析了Hub基因的生物學(xué)意義。

1 材料與方法

1.1 數(shù)據(jù)的獲取及預(yù)處理

1.1.1 水稻基因芯片數(shù)據(jù)的獲取及預(yù)處理 水稻的基因芯片數(shù)據(jù)來源于NCBI的GEO (gene expression omnibus)數(shù)據(jù)庫(GPL2025平臺)。芯片數(shù)據(jù)的預(yù)處理利用R (v3.5.1)軟件完成, 其過程如圖1所示。首先利用arrayQualityMetrics包對數(shù)據(jù)進行質(zhì)量控制; 然后利用affy包的RMA算法(背景處理、歸一化處理、匯總)計算芯片表達水平; 隨后再利用biomaRt包[17-18]進行探針號注釋, 當(dāng)多個探針注釋到同一基因時, 取多探針表達量的平均值作為該基因表達量。分別合并與干旱脅迫相關(guān)的4個數(shù)據(jù)集GSE6901、GSE21651、GSE23211、GSE26280獲62個樣本, 與鹽脅迫相關(guān)的3個數(shù)據(jù)集GSE6901、GSE14403、GSE16108獲32個樣本(詳見附表1)。用limma包的removeBatchEffect函數(shù)去除批次效應(yīng),且對低表達基因進行了過濾用于后續(xù)分析。

附表1 來源于NCBI的Affymetrix基因芯片數(shù)據(jù)數(shù)據(jù)Table S1 Affymetrix microarray data from NCBI

1.1.2 水稻轉(zhuǎn)錄組數(shù)據(jù)的獲取及預(yù)處理 水稻轉(zhuǎn)錄組RNA-seq數(shù)據(jù)來源于NCBI的SRA (sequence read archive)數(shù)據(jù)庫(Illumina平臺), 干旱脅迫相關(guān)有SRR7054176-83、SRR3051740-45、SRR3051752-57共20個runs, 鹽脅迫相關(guān)有ERR266221-38、SRR3647326-31共24個runs (選用李旭凱等[10]所用數(shù)據(jù), 詳見附表2)。數(shù)據(jù)預(yù)處理過程如圖2所示。首先利用fasterq-dump (v2.10.7)工具將下載的SRA格式數(shù)據(jù)轉(zhuǎn)換為fastq格式序列文件, 并利用FastQC(v0.11.9)[19]軟件對原始測序數(shù)據(jù)進行質(zhì)量評估; 接著利用fastp (0.20.1)[20]軟件做質(zhì)量控制, 得到clean data;然后根據(jù)MSU Rice Genome Annotation Project數(shù)據(jù)庫(http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/)的水稻參考基因組和注釋信息,利用Hisat2 (v2.2.10)軟件對clean data進行序列比對;隨后利用Samtools (v0.1.19)軟件將SAM文件轉(zhuǎn)換為BAM文件并重新排序后, 用featureCounts (v2.0.1)[21]軟件得到每個基因在各個樣本中的原始reads計數(shù); 本研究使用R包DESeq2[22]獲取RNA-seq數(shù)據(jù)標(biāo)準(zhǔn)化后的基因表達量用于后續(xù)分析。

附表2 來源于NCBI的RNA-seq數(shù)據(jù)(RNA-seq)Table S2 RNA-seq data from NCBI

經(jīng)過上述對同一平臺同一脅迫相關(guān)的多個數(shù)據(jù)集的數(shù)據(jù)融合, 共獲4個水稻數(shù)據(jù)集: 干旱脅迫相關(guān)的基因芯片數(shù)據(jù)集D_affy和RNA-seq數(shù)據(jù)集D_rnaseq, 鹽脅迫相關(guān)的基因芯片數(shù)據(jù)集S_affy和RNA-seq數(shù)據(jù)集S_rnaseq, 數(shù)據(jù)詳見表1。

表1 水稻數(shù)據(jù)集Table 1 Data set of rice

1.2 基于WGCNA-P和WGCNA-MIC的共表達網(wǎng)絡(luò)分析

數(shù)據(jù)經(jīng)預(yù)處理后仍然包含2萬多個基因(表1),考慮到直接進行共表達網(wǎng)絡(luò)分析計算量過大, 本研究采用前文提及的最大信息系數(shù)MIC進行基因初篩。分別計算各數(shù)據(jù)中基因與表型之間的MIC值,MIC值越高, 意味著該基因與表型相關(guān)性越大, 我們選取MIC值較高的前30%基因用于后續(xù)的共表達網(wǎng)絡(luò)分析。

本研究中, 經(jīng)典的加權(quán)基因共表達網(wǎng)絡(luò)WGCNA-P構(gòu)建直接利用R語言中的WGCNA包提供的一系列函數(shù)實現(xiàn), 而改進的WGCNA-MIC法則基于WGCNA包中的相關(guān)函數(shù)自編代碼實現(xiàn)(代碼見附件)。二者構(gòu)建的主要步驟如下:

(3)構(gòu)建拓?fù)渲丿B矩陣TOMij;

(4)計算距離矩陣disTOM =1 - TOMij, 構(gòu)建層次聚類樹, 并利用動態(tài)剪枝算法獲得基因模塊, 模塊最小基因數(shù)設(shè)為30。接著對相似模塊進行合并,合并閾值為0.2 (cutHeight=0.2)。

1.3 表型相關(guān)基因模塊識別

為識別網(wǎng)絡(luò)中的與表型相關(guān)的顯著模塊, 通常有以下2種方法:

(1)計算基因模塊特征基因(module eigengenes,MEs)與表型的相關(guān)系數(shù), 設(shè)為ME, 其中某一模塊的第一主成分被定義為該模塊的特征基因。

(2)計算模塊的顯著性系數(shù)(module significance,MS), 模塊顯著性MS是該模塊內(nèi)所有基因的顯著性(Gene Significance, GS)的均值[9], GS為基因與表型性狀的相關(guān)系數(shù)絕對值。某模塊的ME和MS值越大, 與表型越相關(guān)。本研究中, WGCNA-P方法所有涉及相關(guān)系數(shù)的計算均采用皮爾遜相關(guān)系數(shù), 而WGCNA-MIC方法中則均用最大信息系數(shù)MIC。綜合考慮ME、MS值, 模塊數(shù)及所選模塊的代表性, 本文對模塊數(shù)小于10的選擇1個顯著模塊、大于等于10且小于15的選擇2個顯著模塊、大于15的選擇4個顯著模塊。

1.4 Hub基因選擇

利用網(wǎng)絡(luò)中連接度高的樞紐節(jié)點來確定基因的優(yōu)先級, 是一種理解和解釋網(wǎng)絡(luò)和整體生物復(fù)雜性的簡便方法[23]。Hub基因是依據(jù)基因與表型性狀之間的相關(guān)性GS值、基因與其所在模塊特征基因間的相關(guān)性MM值來選取。對同一脅迫的2個不同平臺數(shù)據(jù)分別基于WGCNA-P和WGCNA-MIC可獲得該脅迫相關(guān)的4個Hub基因子數(shù)據(jù)集, 對其進行元分析, 取并集, 可獲得該脅迫相關(guān)的Hub基因總集。

1.5 Hub基因的預(yù)測性能

支持向量機(support vector machine, SVM)提供了一種高效分兩類或兩類以上數(shù)據(jù)的方法[24], 為驗證Hub基因選擇的合理性, 本研究基于干旱脅迫和鹽脅迫的8個Hub基因子集及最終的2個Hub基因總集依次構(gòu)建SVM模型對表型進行分類預(yù)測。通過5次5折交叉驗證進行測試, 最終以平均精度作為最后的預(yù)測結(jié)果。

1.6 Hub基因功能分析

1.6.1 GO富集及文獻報道分析 利用AgriGo(http://systemsbiology.cau.edu.cn/agriGOv2/index.p hp)[25]富集分析工具對Hub基因進行GO富集分析。從國家水稻數(shù)據(jù)中心(http://www.ricedata.cn/)的ontology系統(tǒng)分別以檢索條件“干旱”、“鹽”進行檢索,并獲取到文獻已報道的250個干旱脅迫相關(guān)基因和363個鹽脅迫相關(guān)基因。隨后分別分析Hub基因總集中已報道基因情況, 并結(jié)合結(jié)果進一步挖掘可能的相關(guān)基因。

1.6.2 蛋白質(zhì)互作網(wǎng)絡(luò)構(gòu)建與分析 利用STRING和Cytoscape工具構(gòu)建Hub基因的蛋白互作網(wǎng)絡(luò)。將Hub基因?qū)隨TRING (v11.0)[26]蛋白互作在線分析工具(https://string-db.org/)構(gòu)建蛋白質(zhì)互作網(wǎng)絡(luò), 采用默認(rèn)設(shè)置, 獲得并導(dǎo)出蛋白互作數(shù)據(jù)。利用Cytoscape (v3.7.1)[27]工具提取已報道脅迫相關(guān)Hub基因及其相關(guān)基因的子網(wǎng)絡(luò)進行可視化分析,每一個基因由網(wǎng)絡(luò)中的一個節(jié)點表示, 相互連接的2個基因之間存在著某種關(guān)系。

2 結(jié)果與分析

2.1 WGCNA-P、WGCNA-MIC及顯著模塊識別

如圖3所示, 干旱脅迫基因芯片數(shù)據(jù)D_affy基于WGCNA-P進行基因共表達網(wǎng)絡(luò)分析時, 動態(tài)剪切得到35個基因模塊, 合并后得到23個模塊; 基于WGCNA-MIC方法分析, 30個模塊合并后得到10個模塊?？v坐標(biāo)的不同顏色代表不同的模塊, 各模塊與干旱脅迫之間的相關(guān)性及模塊顯著性詳見圖3,基于WGCNA-P識別的顯著模塊及模塊內(nèi)基因數(shù)分別為brown (833)、red (383)、darkgrey (106)、purple(260) 四個模塊共1582個基因, 而基于WGCNAMIC識別的模塊為darkturquoise (1114)和midnightblue (265)兩個模塊共1379個基因。

如圖4所示, 干旱脅迫RNA_seq數(shù)據(jù)D_rnaseq基于WGCNA-P和WGCNA-MIC方法分別得到13個和7個模塊, 且前者分別選取了saddlebrown(983)、darkorange (395)兩個模塊共1378個基因, 而后者則選取magenta (4089)模塊以用于后續(xù)分析。

鹽脅迫基因芯片數(shù)據(jù)(附圖1), 使用WGCNA-P方法時, 動態(tài)剪切得到20個模塊, 經(jīng)合并后得到17個模塊; 使用WGCNA-MIC方法時, 43個模塊合后得到40個模塊?；赪GCNA-P方法最終選取的模塊及模塊內(nèi)基因數(shù)分別magenta (213)、red (331)、purple (199)、pink (803)四個模塊共1546個基因, 而基于WGCNA-MIC方法則分別選取了turquoise(2818)、darkturquoise (74)、red (210)和brown (409)四個模塊共3511個基因以用于后續(xù)分析。鹽脅迫的RNA-seq數(shù)據(jù)(附圖2), 基于WGCNA-P和WGCNAMIC網(wǎng)絡(luò)最終分別得到27個和19個模塊, 且前者分別選取了brown (1148)、plum1 (73)、darkgreen(183)、magenta (351)四個模塊共1755個基因, 而后者則分別選取了pink (308)、lightcyan (42)、darkgreen(130)和darkturquoise (193)四個模塊共673個基因以用于后續(xù)分析。

2.2 Hub基因的預(yù)測性能

本研究中, 干旱脅迫相關(guān)Hub基因挑選閾值設(shè)為GS>0.4且MM>0.83, 鹽脅迫相關(guān)Hub基因篩選閾值設(shè)為GS>0.3且MM>0.75。4個數(shù)據(jù)集分別基于2種網(wǎng)絡(luò)分析方法共獲得8個Hub基因子集(D_affy_P、D_affy_MIC、D_rnaseq_P、D_rnaseq_MIC和S_affy_P、S_affy_MIC、S_rnaseq_P、S_rnaseq_MIC),對基因子集元分析后得到干旱脅迫相關(guān)Hub基因總集D_meta_hub和鹽脅迫相關(guān)Hub基因總集S_meta_hub?；诟鱄ub基因集對表型的SVM分類精度如表2所示, Hub基因的預(yù)測性能整體表現(xiàn)優(yōu)異, 其中,基于WGCNA-MIC方法獲取的Hub基因, 較之基于WGCNA-P方法獲取的Hub基因預(yù)測精度略高, 元分析后的Hub基因總集D_meta_hub和S_meta_hub,在各數(shù)據(jù)集上的平均預(yù)測精度比各Hub基因子集的精度略高。結(jié)果表明, Hub基因與表型性狀相關(guān)性強,WGCNA-MIC方法和元分析有效。

表2 Hub基因的分類精度Table 2 Classification accuracy of Hub genes

2.3 GO功能富集分析

利用AgriGO在線功能富集分析工具, 分別對干旱/鹽脅迫相關(guān)Hub基因集進行基因功能富集分析, 在生物學(xué)過程(biological process, BP)、分子功能(molecular function, MF)和細(xì)胞組分(cellular component, CC)三大分類中都顯著富集到了多個相關(guān)GO通路。具體富集結(jié)果如表3所示。干旱脅迫相關(guān)富集結(jié)果顯示, 生物學(xué)過程中, 顯著富集到的通路,包括應(yīng)對刺激的通路: 內(nèi)源性刺激響應(yīng)(GO:0009719)、激素刺激響應(yīng)(GO:0009725)和非生物刺激響應(yīng)(GO:0009628)等; 參與特殊代謝物代謝過程的通路: 萜類化合物代謝過程(GO:0006721)等;與干旱脅迫較為直接相關(guān)的通路: 對水的響應(yīng)(GO:0009415)和滲透脅迫響應(yīng)(GO:0006970)等。分子功能中, 顯著富集到了與信號傳導(dǎo)相關(guān)的通路:受體活性(GO:0004872)、翻譯因子活性與核酸結(jié)合(GO:0008135)等; 一些參與調(diào)控某些蛋白質(zhì)酶相關(guān)的通路: 蛋白質(zhì)酪氨酸激酶活性(GO:0004713)等。另外, 還有不少顯著富集到細(xì)胞組分相關(guān)的通路: 薄膜(GO:0016020)等。鹽脅迫相關(guān)富集結(jié)果顯示, 生物學(xué)過程中, 顯著富集的通路, 包括參與各種物質(zhì)代謝過程的通路: 草酸代謝過程(GO:0043436)和有機酸代謝過程(GO:0006082)等; 響應(yīng)脅迫相關(guān)的功能:內(nèi)源性激素的響應(yīng)(GO:0009719)等; 參與光合作用:光刺激響應(yīng)(GO:0009416)等。分子功能中, 最顯著富集到的通路是受體活性(GO:0004872)。細(xì)胞組分中富集到了很多與膜組分相關(guān)參與滲透作用的通路,如薄膜(GO:0016020)和細(xì)胞質(zhì)膜等(GO:0005886);參與光合作用的組件: 葉綠體(GO:0009507)等。

表3 Hub基因的GO富集部分分析結(jié)果Table 3 GO enrichment of Hub partial genes

綜上, 基于元分析獲取的2種脅迫的Hub基因,均富集到了內(nèi)源性刺激響應(yīng)(GO:0009719)、激素刺激響應(yīng)(GO:0009725)和非生物刺激響應(yīng)(GO:0009628)等脅迫響應(yīng)的相關(guān)通路上。

2.4 文獻報道情況分析

為了驗證研究結(jié)果的可靠性, 根據(jù)從國家水稻數(shù)據(jù)中心獲取到的已報道干旱和鹽脅迫相關(guān)基因分析所選Hub基因的文獻報道情況。本研究所選Hub基因中有已報道干旱脅迫相關(guān)基因31個和鹽脅迫相關(guān)基因22個, 如表4所示。

表4 已報道與脅迫相關(guān)的Hub基因Table 4 Hub genes related to stress have been reported

2.5 Hub基因互作網(wǎng)絡(luò)構(gòu)建

利用在線分析工具STRING和Cytoscape軟件挖掘Hub基因總集中蛋白互作關(guān)系, 重點關(guān)注與已報道脅迫相關(guān)的Hub基因互作情況。Hub基因總集中, 與2個及2個以上已報道Hub基因有較強關(guān)系,即網(wǎng)絡(luò)中節(jié)點度≥2且STRING中的combined_score≥0.9的Hub基因考慮作為脅迫候選基因可被進一步挖掘。如圖5和圖6, 圖中紅色節(jié)點表示前文得到的已報道脅迫相關(guān)Hub基因(不包括STRING庫中未匹配到蛋白質(zhì)的基因和無相關(guān)蛋白的基因), 節(jié)點越大, 表示與之相關(guān)的基因越多, 線條越粗且顏色越暗, 表示基因之間關(guān)系越強。最終找到了與已報道Hub基因存在蛋白互作關(guān)系的干旱脅迫候選基因11個(圖5中橙色節(jié)點), 鹽脅迫候選基因5個(圖6中橙色節(jié)點), 詳見表5。

表5 候選基因在STRING中的注釋Table 5 Candidate genes annotation in STRING

3 討論

植物為應(yīng)對干旱脅迫環(huán)境, 在生化、細(xì)胞和分子等水平上進化出了很多機制[28], 需要改變基因表達來激活促進耐旱性的代謝過程, 這包括特殊代謝物的合成與積累, 并涉及到物種和基因型特異性的酚類化合物、類黃酮、萜類化合物和含氮化合物的產(chǎn)生[29]。鹽脅迫威脅作物生長主要體現(xiàn)在滲透和氧化2個方面, 這不僅會導(dǎo)致葉片脫落、根芽壞死等不良癥狀的發(fā)生, 而且潛在地延遲了光合作用、植物激素功能、代謝途徑和基因/蛋白質(zhì)功能等生理活動[30]。

(續(xù)表5)

本研究以水稻干旱和鹽脅迫相關(guān)的Affymetrix基因芯片和RNA-seq兩種不同平臺的數(shù)據(jù)為研究對象, 基于WGCNA-P和WGCNA-MIC對其進行了脅迫相關(guān)Hub基因的挖掘。從Hub基因預(yù)測性能來看,各Hub基因集的預(yù)測精度均達80%以上, 預(yù)測性能整體較好。從GO富集分析和文獻報道來看, 一方面2種脅迫的Hub基因集都富集到了內(nèi)源性刺激響應(yīng)(GO:0009719)、激素刺激響應(yīng)(GO:0009725)和非生物刺激響應(yīng)(GO:0009628)等干旱和鹽響應(yīng)相關(guān)通路;另一方面也找到了一些已報道的干旱/鹽脅迫相關(guān)Hub基因,EDT1(Os05g0437700)和OsMYB6(Os04g0 676700)既是已報道干旱脅迫相關(guān)Hub基因, 又是已報道鹽脅迫相關(guān)Hub基因。通過總結(jié)部分已報道基因的文獻發(fā)現(xiàn), 某個脅迫可能由多個基因協(xié)同互作參與調(diào)控, 某個基因也可能參與了多個非生物脅迫。比如, 31個已報道干旱脅迫響應(yīng)基因中, 包括OsP5CS(Os05g0455500)[31]的表達受高鹽、干旱、冷脅迫和ABA處理的誘導(dǎo);OsMADS26(Os08g0112 700)[32]是水稻響應(yīng)多種脅迫的調(diào)控中心;OsbZIP23(Os02g0766700)[33]增強了水稻的抗旱耐鹽性和對ABA的敏感性;OsSIK1(Os06g0130100)[34]在水稻耐鹽和耐旱過程中起重要作用;OsRPK1(Os09g0552 300)[35]在鹽脅迫下表達水平增加, 其表達也受寒冷、干旱與脫落酸等因素的誘導(dǎo)。22個已報道鹽脅迫響應(yīng)基因中, 包括OsSRFP1(Os03g0348 900)[36]負(fù)向調(diào)控水稻的耐鹽性和耐低溫性; 干旱和鹽處理誘導(dǎo)OsPTR8(Os03g0719900)[37]的表達上調(diào);OsSCP(Os07g0129200)[38]在非生物脅迫應(yīng)答中通過調(diào)控脅迫應(yīng)答基因而發(fā)揮作用; 冷脅迫和鹽脅迫會導(dǎo)致OsPIMT1(Os08g0557000)[39]表達量增加2倍;OsLEA5(Os05g0584200)[40]與多種非生物脅迫抗性相關(guān)等。同時, 與李旭凱等[10]、Zhu等[11]和Lv等[12]的研究結(jié)果相比, 我們挖掘到的Hub基因中包含了部分已被廣泛報道的水稻干旱脅迫和鹽脅迫響應(yīng)相關(guān)的轉(zhuǎn)錄因子, 如干旱脅迫相關(guān)的bZIP轉(zhuǎn)錄因子家族(Os02g0766700、Os06g0211200、Os05g0569300)、MYB轉(zhuǎn)錄因子家族(Os04g0676700)、NAC轉(zhuǎn)錄因子家族(Os11g0126900)和HSF轉(zhuǎn)錄因子家族(Os03g0745000)等; 鹽脅迫相關(guān)的bZIP轉(zhuǎn)錄因子家族(Os05g0437700)和MYB轉(zhuǎn)錄因子家族(Os04g06 76700)等。最后, 通過分析Hub基因總集中已報道與脅迫相關(guān)的Hub基因及其相關(guān)基因之間的互作網(wǎng)絡(luò), 進一步挖掘到了與干旱或鹽脅迫相關(guān)較為緊密的候選基因。

綜上, 利用元分析的思路對水稻多平臺基因表達數(shù)據(jù)進行整合分析, 可挖掘到水稻干旱和鹽脅迫的關(guān)鍵基因, 對農(nóng)作物非生物脅迫響應(yīng)的基因挖掘具有一定的參考價值。STRING分析時, 參數(shù)閾值的設(shè)置不同, 所獲候選基因的數(shù)量也會有所不同, 本研究中利用combined_score≥0.9獲得的候選基因,可根據(jù)實際情況適當(dāng)調(diào)整閾值, 并有待進一步利用實時熒光定量PCR (RT-qPCR)驗證。

4 結(jié)論

對多平臺數(shù)據(jù), 通過加權(quán)基因共表達網(wǎng)絡(luò)分析、元分析和蛋白互作網(wǎng)絡(luò)分析, 最終獲得水稻干旱脅迫和鹽脅迫相關(guān)的Hub基因分別為1936個和1504個, 其中文獻已報道的干旱脅迫和鹽脅迫相關(guān)Hub基因分別是31個和22個, 預(yù)測得到的干旱脅迫和鹽脅迫候選基因分別是11個和5個。水稻其他非生物脅迫(如冷脅迫、高溫脅迫等)多平臺數(shù)據(jù)數(shù)據(jù)結(jié)構(gòu)及其實驗原理與干旱和鹽脅迫類似, 故此方法可推廣至其他非生物脅迫相關(guān)基因挖掘。本研究為充分利用多平臺數(shù)據(jù)挖掘水稻非生物脅迫相關(guān)基因提供了新的思路, 也為進一步研究抗逆性水稻品種提供了參考。