基于多平台基因表达数据的水稻干旱和盐胁迫相关基因预测

刘亚文 张红燕,,* 曹 丹 李兰芝

1 湖南农业大学信息与智能科学技术学院, 湖南长沙 410128; 2 湖南农业大学 / 湖南省农业大数据分析与决策工程技术研究中心, 湖南长沙 410128

水稻在生长发育的过程中受到干旱、高盐等非生物逆境因素胁迫时, 易导致大面积的减产、品质下降甚至坏死[1], 提高其逆境抗性将增加农业产量并扩大适宜耕种面积, 缓解人口压力。水稻的逆境抗性受多基因控制, 基于基因组学数据挖掘水稻非生物胁迫相关基因, 对培育抗逆水稻新品种具有重大意义。近年来, 随着大规模基因表达水平测量技术的发展, 基于杂交原理的基因芯片技术[2]和基于高通量测序技术的RNA-seq[3], 被众多学者用于植物胁迫响应基因的挖掘研究中[4]。

然而, 大多实验组测序样本有限, 仅单一地从与水稻某非生物胁迫相关的单个实验组测序数据来挖掘胁迫相关基因, 结果不稳定, 也很难让人信服[5]。当前公共数据库中积累的大量水稻胁迫相关基因芯片和RNA-seq表达数据, 为多平台数据分析提供了研究空间。研究表明, 融合多平台数据能够提高基因表达分析的准确性和可靠性, 多平台表达数据的整合分析成为水稻非生物胁迫相关基因预测研究的趋势[6]。当前, 多平台基因表达数据的融合通常可分为两类: (1) 基于输出层面融合的元分析法。它通过对多个研究结果进行合并汇总, 增大样本总量, 提高检测准确率和统计分析结果的一致性[7]。(2) 基于原始数据融合的数据转换法。它先通过把不同平台基因表达数据按一定规则转换到同一个数据范围内, 再将转换后的多个平台实验数据直接合并成一个表达数据矩阵, 以此来增加样本数目缓解“高维数、小样本”维数灾难问题[8]。综合考虑单个实验组水稻测序的小样本, 及水稻基因芯片数据与RNA-seq数据之间尺度与维度的差异, 本研究首先对同一胁迫相关的多个基因芯片数据或多个RNA-seq数据分别采取数据转换法融合, 再分别基于融合后的基因芯片表达数据集和RNA-seq表达数据集进行胁迫相关基因挖掘, 最后将二者的结果实施元分析, 获取最终的胁迫响应基因。

为了有效利用多平台基因表达数据, 本文选用加权基因共表达网络分析(Weighted Gene Co-expression Network Analysis, WGCNA)法来挖掘关键基因。WGCNA利用基因表达数据构造协同表达的基因模块, 并根据基因模块与表型的关联性以及基因模块的内连性来鉴定关键基因[9], 其基本假定是“表达模式相似的基因功能相似”。它可将表达模式相似的基因进行聚类, 并分析模块与特定性状或表型之间的关联关系, 因此在作物的干旱胁迫、盐胁迫等非生物胁迫相关基因的挖掘研究中被广泛应用。例如李旭凯等[10]利用WGCNA挖掘到2599个与水稻冷胁迫、干旱胁迫和盐胁迫都相关的基因, 并预测出25个抗逆关键基因; Zhu等[11]通过对转录组数据进行WGCNA分析, 确定了水稻盐胁迫响应核心差异基因和模块; Lv等[12]以转录组数据为基础进行WGCNA分析, 预测了各模块重要的Hub差异基因和调控水稻干旱应答基因表达的主要转录因子;Hopper等[13]使用时间序列转录方法结合WGCNA网络分析, 为葡萄耐旱性研究提供了候选基因; 秦天元等[14]使用WGCNA挖掘马铃薯根系抗旱核心基因,并进一步利用RT-qPCR验证出挖掘到的核心基因确实响应干旱胁迫。经典的WGCNA 以Pearson相关系数度量2个基因表达量间的线性相似性(记为WGCNA-P), 但无法捕获基因间可能广泛存在的非线性关联。Reshef等[15]学者基于信息论中的互信息理论提出了一种可度量两变量非线性相关性的普适性测度最大信息系数(Maximal Information Coefficient, MIC), 论文提出以MIC作为相似性度量替代WGCNA中的Pearson相关系数来构建基因共表达网络(记为WGCNA-MIC), 以捕捉基因间的非线性关联。同时, 考虑到特定线性情形下MIC的统计功效[16]不如Pearson相关系数, 所以本研究对同一数据集分别基于WGCNA-P和WGCNA-MIC两种方法来构建基因共表达网络, 并对各自获取的Hub基因集进行整合分析。

综上, 本研究以多平台水稻非生物胁迫(以干旱和盐胁迫为代表)相关的基因芯片数据和RNA-seq数据为研究对象, 分别以WGCNA-P和WGCNAMIC挖掘胁迫相关Hub基因, 进而对同一胁迫不同平台数据使用以上2种网络分析法得到的Hub基因进行整合分析, 得到最终的胁迫相关Hub基因集。最后, 从预测性能、基因功能富集分析、文献报道和互作网络分析等多角度解析了Hub基因的生物学意义。

1 材料与方法

1.1 数据的获取及预处理

1.1.1 水稻基因芯片数据的获取及预处理 水稻的基因芯片数据来源于NCBI的GEO (gene expression omnibus)数据库(GPL2025平台)。芯片数据的预处理利用R (v3.5.1)软件完成, 其过程如图1所示。首先利用arrayQualityMetrics包对数据进行质量控制; 然后利用affy包的RMA算法(背景处理、归一化处理、汇总)计算芯片表达水平; 随后再利用biomaRt包[17-18]进行探针号注释, 当多个探针注释到同一基因时, 取多探针表达量的平均值作为该基因表达量。分别合并与干旱胁迫相关的4个数据集GSE6901、GSE21651、GSE23211、GSE26280获62个样本, 与盐胁迫相关的3个数据集GSE6901、GSE14403、GSE16108获32个样本(详见附表1)。用limma包的removeBatchEffect函数去除批次效应,且对低表达基因进行了过滤用于后续分析。

附表1 来源于NCBI的Affymetrix基因芯片数据数据Table S1 Affymetrix microarray data from NCBI

1.1.2 水稻转录组数据的获取及预处理 水稻转录组RNA-seq数据来源于NCBI的SRA (sequence read archive)数据库(Illumina平台), 干旱胁迫相关有SRR7054176-83、SRR3051740-45、SRR3051752-57共20个runs, 盐胁迫相关有ERR266221-38、SRR3647326-31共24个runs (选用李旭凯等[10]所用数据, 详见附表2)。数据预处理过程如图2所示。首先利用fasterq-dump (v2.10.7)工具将下载的SRA格式数据转换为fastq格式序列文件, 并利用FastQC(v0.11.9)[19]软件对原始测序数据进行质量评估; 接着利用fastp (0.20.1)[20]软件做质量控制, 得到clean data;然后根据MSU Rice Genome Annotation Project数据库(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文件转换为BAM文件并重新排序后, 用featureCounts (v2.0.1)[21]软件得到每个基因在各个样本中的原始reads计数; 本研究使用R包DESeq2[22]获取RNA-seq数据标准化后的基因表达量用于后续分析。

附表2 来源于NCBI的RNA-seq数据(RNA-seq)Table S2 RNA-seq data from NCBI

经过上述对同一平台同一胁迫相关的多个数据集的数据融合, 共获4个水稻数据集: 干旱胁迫相关的基因芯片数据集D_affy和RNA-seq数据集D_rnaseq, 盐胁迫相关的基因芯片数据集S_affy和RNA-seq数据集S_rnaseq, 数据详见表1。

表1 水稻数据集Table 1 Data set of rice

1.2 基于WGCNA-P和WGCNA-MIC的共表达网络分析

数据经预处理后仍然包含2万多个基因(表1),考虑到直接进行共表达网络分析计算量过大, 本研究采用前文提及的最大信息系数MIC进行基因初筛。分别计算各数据中基因与表型之间的MIC值,MIC值越高, 意味着该基因与表型相关性越大, 我们选取MIC值较高的前30%基因用于后续的共表达网络分析。

本研究中, 经典的加权基因共表达网络WGCNA-P构建直接利用R语言中的WGCNA包提供的一系列函数实现, 而改进的WGCNA-MIC法则基于WGCNA包中的相关函数自编代码实现(代码见附件)。二者构建的主要步骤如下:

(3)构建拓扑重叠矩阵TOMij;

(4)计算距离矩阵disTOM =1 - TOMij, 构建层次聚类树, 并利用动态剪枝算法获得基因模块, 模块最小基因数设为30。接着对相似模块进行合并,合并阈值为0.2 (cutHeight=0.2)。

1.3 表型相关基因模块识别

为识别网络中的与表型相关的显著模块, 通常有以下2种方法:

(1)计算基因模块特征基因(module eigengenes,MEs)与表型的相关系数, 设为ME, 其中某一模块的第一主成分被定义为该模块的特征基因。

(2)计算模块的显著性系数(module significance,MS), 模块显著性MS是该模块内所有基因的显著性(Gene Significance, GS)的均值[9], GS为基因与表型性状的相关系数绝对值。某模块的ME和MS值越大, 与表型越相关。本研究中, WGCNA-P方法所有涉及相关系数的计算均采用皮尔逊相关系数, 而WGCNA-MIC方法中则均用最大信息系数MIC。综合考虑ME、MS值, 模块数及所选模块的代表性, 本文对模块数小于10的选择1个显著模块、大于等于10且小于15的选择2个显著模块、大于15的选择4个显著模块。

1.4 Hub基因选择

利用网络中连接度高的枢纽节点来确定基因的优先级, 是一种理解和解释网络和整体生物复杂性的简便方法[23]。Hub基因是依据基因与表型性状之间的相关性GS值、基因与其所在模块特征基因间的相关性MM值来选取。对同一胁迫的2个不同平台数据分别基于WGCNA-P和WGCNA-MIC可获得该胁迫相关的4个Hub基因子数据集, 对其进行元分析, 取并集, 可获得该胁迫相关的Hub基因总集。

1.5 Hub基因的预测性能

支持向量机(support vector machine, SVM)提供了一种高效分两类或两类以上数据的方法[24], 为验证Hub基因选择的合理性, 本研究基于干旱胁迫和盐胁迫的8个Hub基因子集及最终的2个Hub基因总集依次构建SVM模型对表型进行分类预测。通过5次5折交叉验证进行测试, 最终以平均精度作为最后的预测结果。

1.6 Hub基因功能分析

1.6.1 GO富集及文献报道分析 利用AgriGo(http://systemsbiology.cau.edu.cn/agriGOv2/index.p hp)[25]富集分析工具对Hub基因进行GO富集分析。从国家水稻数据中心(http://www.ricedata.cn/)的ontology系统分别以检索条件“干旱”、“盐”进行检索,并获取到文献已报道的250个干旱胁迫相关基因和363个盐胁迫相关基因。随后分别分析Hub基因总集中已报道基因情况, 并结合结果进一步挖掘可能的相关基因。

1.6.2 蛋白质互作网络构建与分析 利用STRING和Cytoscape工具构建Hub基因的蛋白互作网络。将Hub基因导入STRING (v11.0)[26]蛋白互作在线分析工具(https://string-db.org/)构建蛋白质互作网络, 采用默认设置, 获得并导出蛋白互作数据。利用Cytoscape (v3.7.1)[27]工具提取已报道胁迫相关Hub基因及其相关基因的子网络进行可视化分析,每一个基因由网络中的一个节点表示, 相互连接的2个基因之间存在着某种关系。

2 结果与分析

2.1 WGCNA-P、WGCNA-MIC及显著模块识别

如图3所示, 干旱胁迫基因芯片数据D_affy基于WGCNA-P进行基因共表达网络分析时, 动态剪切得到35个基因模块, 合并后得到23个模块; 基于WGCNA-MIC方法分析, 30个模块合并后得到10个模块。纵坐标的不同颜色代表不同的模块, 各模块与干旱胁迫之间的相关性及模块显著性详见图3,基于WGCNA-P识别的显著模块及模块内基因数分别为brown (833)、red (383)、darkgrey (106)、purple(260) 四个模块共1582个基因, 而基于WGCNAMIC识别的模块为darkturquoise (1114)和midnightblue (265)两个模块共1379个基因。

如图4所示, 干旱胁迫RNA_seq数据D_rnaseq基于WGCNA-P和WGCNA-MIC方法分别得到13个和7个模块, 且前者分别选取了saddlebrown(983)、darkorange (395)两个模块共1378个基因, 而后者则选取magenta (4089)模块以用于后续分析。

盐胁迫基因芯片数据(附图1), 使用WGCNA-P方法时, 动态剪切得到20个模块, 经合并后得到17个模块; 使用WGCNA-MIC方法时, 43个模块合后得到40个模块。基于WGCNA-P方法最终选取的模块及模块内基因数分别magenta (213)、red (331)、purple (199)、pink (803)四个模块共1546个基因, 而基于WGCNA-MIC方法则分别选取了turquoise(2818)、darkturquoise (74)、red (210)和brown (409)四个模块共3511个基因以用于后续分析。盐胁迫的RNA-seq数据(附图2), 基于WGCNA-P和WGCNAMIC网络最终分别得到27个和19个模块, 且前者分别选取了brown (1148)、plum1 (73)、darkgreen(183)、magenta (351)四个模块共1755个基因, 而后者则分别选取了pink (308)、lightcyan (42)、darkgreen(130)和darkturquoise (193)四个模块共673个基因以用于后续分析。

2.2 Hub基因的预测性能

本研究中, 干旱胁迫相关Hub基因挑选阈值设为GS>0.4且MM>0.83, 盐胁迫相关Hub基因筛选阈值设为GS>0.3且MM>0.75。4个数据集分别基于2种网络分析方法共获得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),对基因子集元分析后得到干旱胁迫相关Hub基因总集D_meta_hub和盐胁迫相关Hub基因总集S_meta_hub。基于各Hub基因集对表型的SVM分类精度如表2所示, Hub基因的预测性能整体表现优异, 其中,基于WGCNA-MIC方法获取的Hub基因, 较之基于WGCNA-P方法获取的Hub基因预测精度略高, 元分析后的Hub基因总集D_meta_hub和S_meta_hub,在各数据集上的平均预测精度比各Hub基因子集的精度略高。结果表明, Hub基因与表型性状相关性强,WGCNA-MIC方法和元分析有效。

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

2.3 GO功能富集分析

利用AgriGO在线功能富集分析工具, 分别对干旱/盐胁迫相关Hub基因集进行基因功能富集分析, 在生物学过程(biological process, BP)、分子功能(molecular function, MF)和细胞组分(cellular component, CC)三大分类中都显著富集到了多个相关GO通路。具体富集结果如表3所示。干旱胁迫相关富集结果显示, 生物学过程中, 显著富集到的通路,包括应对刺激的通路: 内源性刺激响应(GO:0009719)、激素刺激响应(GO:0009725)和非生物刺激响应(GO:0009628)等; 参与特殊代谢物代谢过程的通路: 萜类化合物代谢过程(GO:0006721)等;与干旱胁迫较为直接相关的通路: 对水的响应(GO:0009415)和渗透胁迫响应(GO:0006970)等。分子功能中, 显著富集到了与信号传导相关的通路:受体活性(GO:0004872)、翻译因子活性与核酸结合(GO:0008135)等; 一些参与调控某些蛋白质酶相关的通路: 蛋白质酪氨酸激酶活性(GO:0004713)等。另外, 还有不少显著富集到细胞组分相关的通路: 薄膜(GO:0016020)等。盐胁迫相关富集结果显示, 生物学过程中, 显著富集的通路, 包括参与各种物质代谢过程的通路: 草酸代谢过程(GO:0043436)和有机酸代谢过程(GO:0006082)等; 响应胁迫相关的功能:内源性激素的响应(GO:0009719)等; 参与光合作用:光刺激响应(GO:0009416)等。分子功能中, 最显著富集到的通路是受体活性(GO:0004872)。细胞组分中富集到了很多与膜组分相关参与渗透作用的通路,如薄膜(GO:0016020)和细胞质膜等(GO:0005886);参与光合作用的组件: 叶绿体(GO:0009507)等。

表3 Hub基因的GO富集部分分析结果Table 3 GO enrichment of Hub partial genes

综上, 基于元分析获取的2种胁迫的Hub基因,均富集到了内源性刺激响应(GO:0009719)、激素刺激响应(GO:0009725)和非生物刺激响应(GO:0009628)等胁迫响应的相关通路上。

2.4 文献报道情况分析

为了验证研究结果的可靠性, 根据从国家水稻数据中心获取到的已报道干旱和盐胁迫相关基因分析所选Hub基因的文献报道情况。本研究所选Hub基因中有已报道干旱胁迫相关基因31个和盐胁迫相关基因22个, 如表4所示。

表4 已报道与胁迫相关的Hub基因Table 4 Hub genes related to stress have been reported

2.5 Hub基因互作网络构建

利用在线分析工具STRING和Cytoscape软件挖掘Hub基因总集中蛋白互作关系, 重点关注与已报道胁迫相关的Hub基因互作情况。Hub基因总集中, 与2个及2个以上已报道Hub基因有较强关系,即网络中节点度≥2且STRING中的combined_score≥0.9的Hub基因考虑作为胁迫候选基因可被进一步挖掘。如图5和图6, 图中红色节点表示前文得到的已报道胁迫相关Hub基因(不包括STRING库中未匹配到蛋白质的基因和无相关蛋白的基因), 节点越大, 表示与之相关的基因越多, 线条越粗且颜色越暗, 表示基因之间关系越强。最终找到了与已报道Hub基因存在蛋白互作关系的干旱胁迫候选基因11个(图5中橙色节点), 盐胁迫候选基因5个(图6中橙色节点), 详见表5。

表5 候选基因在STRING中的注释Table 5 Candidate genes annotation in STRING

3 讨论

植物为应对干旱胁迫环境, 在生化、细胞和分子等水平上进化出了很多机制[28], 需要改变基因表达来激活促进耐旱性的代谢过程, 这包括特殊代谢物的合成与积累, 并涉及到物种和基因型特异性的酚类化合物、类黄酮、萜类化合物和含氮化合物的产生[29]。盐胁迫威胁作物生长主要体现在渗透和氧化2个方面, 这不仅会导致叶片脱落、根芽坏死等不良症状的发生, 而且潜在地延迟了光合作用、植物激素功能、代谢途径和基因/蛋白质功能等生理活动[30]。

(续表5)

本研究以水稻干旱和盐胁迫相关的Affymetrix基因芯片和RNA-seq两种不同平台的数据为研究对象, 基于WGCNA-P和WGCNA-MIC对其进行了胁迫相关Hub基因的挖掘。从Hub基因预测性能来看,各Hub基因集的预测精度均达80%以上, 预测性能整体较好。从GO富集分析和文献报道来看, 一方面2种胁迫的Hub基因集都富集到了内源性刺激响应(GO:0009719)、激素刺激响应(GO:0009725)和非生物刺激响应(GO:0009628)等干旱和盐响应相关通路;另一方面也找到了一些已报道的干旱/盐胁迫相关Hub基因,EDT1(Os05g0437700)和OsMYB6(Os04g0 676700)既是已报道干旱胁迫相关Hub基因, 又是已报道盐胁迫相关Hub基因。通过总结部分已报道基因的文献发现, 某个胁迫可能由多个基因协同互作参与调控, 某个基因也可能参与了多个非生物胁迫。比如, 31个已报道干旱胁迫响应基因中, 包括OsP5CS(Os05g0455500)[31]的表达受高盐、干旱、冷胁迫和ABA处理的诱导;OsMADS26(Os08g0112 700)[32]是水稻响应多种胁迫的调控中心;OsbZIP23(Os02g0766700)[33]增强了水稻的抗旱耐盐性和对ABA的敏感性;OsSIK1(Os06g0130100)[34]在水稻耐盐和耐旱过程中起重要作用;OsRPK1(Os09g0552 300)[35]在盐胁迫下表达水平增加, 其表达也受寒冷、干旱与脱落酸等因素的诱导。22个已报道盐胁迫响应基因中, 包括OsSRFP1(Os03g0348 900)[36]负向调控水稻的耐盐性和耐低温性; 干旱和盐处理诱导OsPTR8(Os03g0719900)[37]的表达上调;OsSCP(Os07g0129200)[38]在非生物胁迫应答中通过调控胁迫应答基因而发挥作用; 冷胁迫和盐胁迫会导致OsPIMT1(Os08g0557000)[39]表达量增加2倍;OsLEA5(Os05g0584200)[40]与多种非生物胁迫抗性相关等。同时, 与李旭凯等[10]、Zhu等[11]和Lv等[12]的研究结果相比, 我们挖掘到的Hub基因中包含了部分已被广泛报道的水稻干旱胁迫和盐胁迫响应相关的转录因子, 如干旱胁迫相关的bZIP转录因子家族(Os02g0766700、Os06g0211200、Os05g0569300)、MYB转录因子家族(Os04g0676700)、NAC转录因子家族(Os11g0126900)和HSF转录因子家族(Os03g0745000)等; 盐胁迫相关的bZIP转录因子家族(Os05g0437700)和MYB转录因子家族(Os04g06 76700)等。最后, 通过分析Hub基因总集中已报道与胁迫相关的Hub基因及其相关基因之间的互作网络, 进一步挖掘到了与干旱或盐胁迫相关较为紧密的候选基因。

综上, 利用元分析的思路对水稻多平台基因表达数据进行整合分析, 可挖掘到水稻干旱和盐胁迫的关键基因, 对农作物非生物胁迫响应的基因挖掘具有一定的参考价值。STRING分析时, 参数阈值的设置不同, 所获候选基因的数量也会有所不同, 本研究中利用combined_score≥0.9获得的候选基因,可根据实际情况适当调整阈值, 并有待进一步利用实时荧光定量PCR (RT-qPCR)验证。

4 结论

对多平台数据, 通过加权基因共表达网络分析、元分析和蛋白互作网络分析, 最终获得水稻干旱胁迫和盐胁迫相关的Hub基因分别为1936个和1504个, 其中文献已报道的干旱胁迫和盐胁迫相关Hub基因分别是31个和22个, 预测得到的干旱胁迫和盐胁迫候选基因分别是11个和5个。水稻其他非生物胁迫(如冷胁迫、高温胁迫等)多平台数据数据结构及其实验原理与干旱和盐胁迫类似, 故此方法可推广至其他非生物胁迫相关基因挖掘。本研究为充分利用多平台数据挖掘水稻非生物胁迫相关基因提供了新的思路, 也为进一步研究抗逆性水稻品种提供了参考。