本发明属于植物分子生物技术和基因工程领域，具体涉及一种基于基因表达量与性状动态相关性预测玉米未知基因功能的方法。
背景技术：
玉米是世界三大主要农作物之一，20世纪90年代以来，世界玉米总产量首次超越水稻和小麦，成为第一位的粮食作物。玉米籽粒中累积了大量的储存物质包括淀粉、油份和蛋白。随着人们生活水平的提高和膳食结构的变化，以及淀粉和油脂加工业的发展，玉米品种由产量型逐渐向质量型转变，玉米品质及其专用性变得越来越重要。
复杂的性状受到多个基因位点的调控，而基因之间相互作用又形成复杂的基因调控网络，来控制细胞中各种生物学反应的进行。高通量测序技术的发展使我们获得了大规模海量的组学数据，如基因型数据、基因表达量数据、蛋白相互作用数据等。研究表明，功能相近基因，其表达模式相关。因此构建共表达网络为预测未知基因功能提供了思路。然而，在构建共表达网络过程中，我们发现很多功能相近基因，其表达模式并不相关。因此，利用共表达分析预测未知基因功能具有一定的局限性。研究表明单个基因/蛋白本身对复杂数量性状的影响有限，往往需要通过高阶的细胞组织形式来行使功能，很多功能相关基因的表达量并不相关，挖掘控制表型性状的遗传位点，而这些遗传位点之间相对独立，它们之间的调控关系未知，并且传统分析方法需要多年多点的表型鉴定，费时费力。
技术实现要素：
针对现有技术存在的问题，本发明提供了一种基于基因表达量与性状动态相关性预测玉米未知基因功能的方法，该方法通过动态关联分析，鉴定玉米基因组基因共表达模式的动态关联，根据la结果，预测未知基因功能。
本发明是通过以下技术方案来实现的：
本发明提供了一种基于基因表达量与性状动态相关性预测玉米未知基因功能的方法，包括以下步骤：
（1）收集玉米自交系授粉后15天的籽粒转录本测序获得基因表达量数据；
（2）动态关联分析la模型的建立；
（3）la显著性评估；
（4）挖掘玉米全基因组基因共表达模式的动态关联；
（5）对显著la结果的基因进行功能注释，预测未知基因的功能。
进一步的，所述玉米自交系分成了2组：热带和亚热带、温带，小组内采用完全随机区组法，设2个重复，每个自交系每个重复播种1行，所有材料均进行自交，收获未成熟的授粉后15天的籽粒，每个自交系的两个重复各取3-4穗，每穗取1-2粒籽粒，混合提取籽粒总rna，随机选择自交系数量的样品用于rna-seq。
上述rna-seq包括以下步骤：首先，用ploy（t）寡聚核苷酸从总rna中抽取全部带ploy（a）尾的rna，主要为mrna，然后将截获的mrna随机打断成片段，用六碱基随机引物合成cdna第一链，并加入逆转录酶合成cdna第二链，经过试剂盒纯化并对cdna片段进行末端修饰，连接测序接头，再经琼脂糖凝胶电泳回收目的大小片段，进行pcr扩增，用illuminagaⅱ基因分析系统进行序列测定及分析，获得基因的表达量数据。
进一步的，所述动态关联分析la模型具体采用以下方法建立：la的数学定义如下：
la(x,y|z)=eg'(z)公式1
其中，所述x、y和z为玉米籽粒中基因表达量数据；假设x，y，z是均值为0，方差为1的连续随机变量，则x，y的相关性表示为e(xy)；当z=z时，g(z)=e(xy|z=z)，g(z)检测的是当z=z时，xy基因对的共表达模式。g(z)的导数表示为g'(z)，该值可用于共表达模式变化的期望测定，
当z符合标准正态分布时，la值可简单的表示为la(x,y|z)=e(xyz)。
x，y，z代表具有正态分布表达谱的三个基因，则la(x,y|z)表示为：e(xyz)=(x1y1z1+x2y2z2+...+xmymzm)/m公式2
la用来反应基因对共表达模式的动态变化，即当z基因表达量较高时，xy基因对的表达量呈正相关（co-regulated），e(xy|z=1)为正数；当z基因表达量较低时，xy基因对的表达量呈负相关（contra-regulated，），e(xy|z=0)为负数，因此基因对的表达调控模式由正相关（co-regulated）转变为负相关（contra-regulated），la值记为正；相反，基因对的表达调控模式由负相关（contra-regulated）转变为正相关（co-regulated），la值记为负。
本发明所建立的动态关联分析模型的评估步骤如下：混合所有基因的表达量值；在每次模拟中，用放回随机抽样法随机抽取一对基因（x，y）的表达量值,z基因取全基因组所有基因，计算xy基因对在全基因组的la值，可分别得到la的正极大值和负极小值；重复模拟一百万次，分别得到la的正值参考分布和负值参考分布，用la正负参考分布的99%分位数作为la正负显著性阈值。
进一步的，所述全基因组动态关联分析的结果按照la值的大小进行过滤，对显著la的基因进行功能注释，预测未知基因功能。
研究表明功能相近基因表达模式不相关的原因主要包含以下两个假设，一是这些功能相近基因的表达调控不在mrna水平上，二是功能相近基因的表达模式只在特定的细胞环境中才相关，即共表达模式的动态关联，动态关联分析（liquidassociation,la）为验证第二种假设提供了理论支持。本发明基于功能相近基因，表达模式相关的科学假设，采用la方法鉴定玉米全基因组基因共表达模式的动态关联，根据显著la结果中基因的功能注释，预测未知基因功能，并根据未知基因在拟南芥中的同源基因功能，验证la预测结果，思路创新，在植物学领域该项研究尚无报道。
本发明的有益效果为：
（1）本发明以玉米籽粒中基因对共表达模式动态关联这一现象为突破口，预测未知基因功能。相比较于传统的共表达网络构建，动态关联分析可以快速找到调控共表达模式的调控基因；
（2）本发明通过对显著la结果的基因进行功能注释，推测未知基因功能，并通过同源基因的功能验证预测结果，是预测未知基因功能的有效方法。
附图说明
图1为随机模拟生成la值评估la分析的显著性。
具体实施方式
下面结合附图和具体实施例对本发明作进一步说明，下述说明仅是实例性的，不限定本发明的保护范围。
实施例1
一种本发明所述基于动态关联分析预测玉米未知基因功能的方法，主要包括三步，基因表达量数据的收集、la显著性评估和全基因组la分析。
（1）基因表达量数据的收集：
368份自交系（本发明所使用的玉米品种可为任意品种，包括中国农业大学宋同明教授培育的35份高油玉米自交系（yang等，2010b））于2010年在湖北荆州种植，根据系谱信息分成了2组（热带和亚热带、温带），小组内采用完全随机区组法，设2个重复，每个自交系每个重复播种1行。所有材料均进行自交，收获未成熟的授粉后15天（15dap）的籽粒，每个自交系的两个重复各取3-4穗，每穗取1-2粒籽粒，混合提取籽粒总rna，随机选择368个样品用于rna-seq。样品的rna-seq工作是由深圳华大基因研究院（beijinggenomicsinstitute，bgi）完成，测序方法简要描述如下：首先，用ploy（t）寡聚核苷酸从总rna中抽取全部带ploy（a）尾的rna，主要为mrna，然后将截获的mrna随机打断成片段，用六碱基随机引物（randomhexamers）合成cdna第一链，并加入逆转录酶等合成cdna第二链，经过试剂盒（ampurexpbeads）纯化并对cdna片段进行末端修饰，连接测序接头，再经琼脂糖凝胶电泳回收目的大小片段，进行pcr扩增，从而完成整个文库构建工作，构建好的文库用illuminagaⅱ基因分析系统进行序列测定及分析。转录本测序获得的368个玉米自交系中28769个基因的表达量数据，对基因表达量数据集进行的缺失值预处理如下：基因表达数据因为实验中的噪声、检测技术等原因而存在缺失。对于数据集中的每个基因，如果其表达值在高于30%的样本中缺失，则在后续的分析中舍弃该基因。共获得24,907个基因表达量数据（该部分数据根据需要可以直接从数据库下载）用于后续的全基因组la分析；
（2）动态关联分析la模型的建立：
所述动态关联分析la模型具体采用以下方法建立：la的数学定义如下：
la(x,y|z)=eg'(z)公式1
其中，所述x、y和z为玉米籽粒中基因表达量数据；假设x，y，z是均值为0，方差为1的连续随机变量，则x，y的相关性表示为e(xy)；当z=z时，g(z)=e(xy|z=z)，g(z)检测的是当z=z时，xy基因对的共表达模式。g(z)的导数表示为g'(z)，该值可用于共表达模式变化的期望测定，
当z符合标准正态分布时，la值可简单的表示为la(x,y|z)=e(xyz)。
x，y，z代表具有正态分布表达谱的三个基因，则la(x,y|z)表示为：e(xyz)=(x1y1z1+x2y2z2+...+xmymzm)/m公式2
la用来反应基因对共表达模式的动态变化，即当z基因表达量较高时，xy基因对的表达量呈正相关（co-regulated），e(xy|z=1)为正数；当z基因表达量较低时，xy基因对的表达量呈负相关（contra-regulated，），e(xy|z=0)为负数，因此基因对的表达调控模式由正相关（co-regulated）转变为负相关（contra-regulated），la值记为正；相反，基因对的表达调控模式由负相关（contra-regulated）转变为正相关（co-regulated），la值记为负。
（3）la显著性评估
混合所有基因的表达量值；在每次模拟中，用放回随机抽样法随机抽取一对基因（x，y）的表达量值,z基因取全基因组所有基因，计算xy基因对在全基因组的la值，可分别得到la的正极大值和负极小值；重复模拟一百万次，分别得到la的正值参考分布和负值参考分布，用la正负参考分布的99%分位数作为la正负显著性阈值，具体结果见图1。
（4）全基因组la分析
以x=全基因组基因，y=全基因组基因，z=全基因组基因进行la分析，重点关注la绝对值最大的前50个共表达基因对（liquidassociationpari,lap）列表。对x、y和z进行功能注释，表1为grmzm5g858880调控的参与蛋白翻译过程的基因列表。基因grmzm5g858880调控多对共表达基因对（liquidassociationpair,lap），玉米基因组数据库（maizegenomedatabase）对该基因的功能注释为“编码包含ww功能域的蛋白”。在grmzm5g858880调控的lap列表中，发现部分基因参与蛋白质翻译过程，包括核糖体蛋白合成、蛋白翻译起始以及蛋白磷酸化，并出现多次，grmzm2g092663(编码核糖体s5蛋白家族，4次)，grmzm2g099352（编码核糖体s3蛋白家族），grmzm2g168149（编码核糖体s5蛋白家族），grmzm2g129015（编码核糖体s26e蛋白家族，2次），grmzm2g164352（编码蛋白磷酸酶2a亚基a2，4次），grmzm2g122135（编码蛋白磷酸酶2a亚基a2，2次），grmzm2g064133（编码真核生物翻译起始因子3g1），因此推测调控基因grmzm5g858880也参与蛋白翻译过程。研究表明，grmzm5g858880在拟南芥中的同源基因（at3g13225）通过核糖体减速和降低再起始效率来调控蛋白翻译过程（tran等，bmcgenomics，2008）。
表1grmzm5g858880调控的参与蛋白翻译过程的基因列表
以上的结果证明了本发明的有效性，通过全基因组基因对共表达模式的动态关联分析，并结合功能注释，预测未知基因功能，为玉米功能基因组学研究提供了新的思路和方法。