CN106980777B

CN106980777B - 基于相似基因比例的基因家族关系判定方法

Info

Publication number: CN106980777B
Application number: CN201710229350.9A
Authority: CN
Inventors: 向浏欣; 吴朝锋; 蔡应繁; 汪露; 邓聿杉; 何琳
Original assignee: Chongqing University of Post and Telecommunications
Current assignee: Chongqing University of Post and Telecommunications
Priority date: 2017-04-10
Filing date: 2017-04-10
Publication date: 2019-12-31
Anticipated expiration: 2037-04-10
Also published as: CN106980777A

Abstract

本发明公开了一种基于相似基因比例的基因家族关系判定方法，包括确定一个基因家族为指定基因家族，计算两个考察物种谁与参考物的指定基因家族亲缘关系更近，计算方法为将参考物种指定基因家族的基因分别与两个考察物种指定基因家族的基因进行两两序列比对，分别获取多个阈值条件下的参考物种的考察物种非冗余基因比例，确定相同阈值条件下参考物种的考察物种非冗余基因比例间的差值的绝对值最大且非冗余基因比例较大的考察物种为与参考物种指定基因家族关系更密切，可认为该考察物种与参考物种具有更相似的相应形态发育特征而应用于优势物种或目的物种的筛选，还可反映不同物种间基因家族的进化关系。

Description

基于相似基因比例的基因家族关系判定方法

技术领域

本发明涉及生物信息学领域，具体涉及基于相似基因比例的基因家族关系判定方法。

背景技术

物种进化关系或亲缘关系研究是生物信息学领域的重要研究内容之一，研究物种进化的重要方法是对物种间重要的直系同源基因或保守基因进行氨基酸或核酸序列比对，根据序列的变化情况来推测基因间的进化关系或亲缘关系，进而推测物种间的进化关系或亲缘关系，二者进化关系或亲缘关系越近表明二者关系更密切。

基因家族是具有相同结构域序列(一段保守的氨基酸序列)的多个基因的集合，物种中有的基因家族甚至包含上百或上千个基因，属于一个大家族。同一基因家族的基因由于具有相同的结构域，因此往往具有相似的生物学功能。当然，基因之间序列相似性越高，基因的功能越相近。然而，目前的进化研究中，不管是单个或少量的直系同源基因或保守基因，还是一个基因家族内的基因的进化分析，往往都是通过序列比对后构建进化树展示“单个”基因之间的进化关系或亲缘关系。各个物种是一个独立的整体，物种之间均具有某一相同的基因家族是常见的，而把物种内的某一基因家族看做一个整体，不同物种同一基因家族间的进化研究或关系研究还未见报道。一个物种任一形态发育特征往往不是单个基因的作用，通常是基因家族或不同基因家族之间共同作用的结果，因此对物种间单个或少量直系同源基因或保守基因的进化分析或关系分析不如对基因家族之间的进化分析或关系分析能更真实反映物种进化、物种间形态发育特征的关系。由于同一基因家族的基因往往具有相似的功能，那么物种间同一基因家族相似性水平越高则表明该物种间同一基因家族的关系越密切，则相应的基因功能表现水平或物种形态发育特征具有更近的相似性，可应用于推断基因功能表现水平及相应形态发育特征的物种筛选。比如，NBS基因家族为熟知的植物抗病相关基因家族，假设已知A物种具有很强的抗病性，又假设有未知物种B、C、D，要想了解B、C、D中哪一个物种具有更强的抗病性，那么，通过本发明基于相似基因比例的基因家族关系判定方法，就可以从B、C、D中筛选出与A具有最近的NBS基因家族亲缘关系的具有较强抗病性的物种，避免了繁琐的、长时间的、大规模的、高花费的生物技术实验的筛选，节约人力物力，可以极大促进物种育种筛选效率，较快获得优秀品种、目的品种。随着基因组测序技术的进步，越来越多的物种的全基因组被测序，那么获得更多物种的完整的基因家族成为可能，则通过本发明基于相似基因比例的基因家族关系判定方法，可以更快筛选到具有特定形态发育特征的和具有特定应用价值的物种。另外，对于异源杂交物种，可反映其与2个亲本物种的某基因家族的遗传进化关系；对于近缘物种，可反映近缘物种间某基因家族的进化关系。

发明内容

为了比较不同物种同一基因家族间的关系密切程度，并了解物种间某一基因家族的功能表现水平，提高特定形态发育特征的物种筛选，提出了本发明基于相似基因比例的基因家族关系判定方法。

基于相似基因比例的基因家族关系判定方法，包括：

确定一个待分析的基因家族作为指定基因家族，通常一个基因家族所包含的基因均具有相同的结构域，结构域为一段比较保守的氨基酸序列；

选择参考物种，获取其所述指定基因家族所包含的基因作为参考基因，指定基因家族的基因的获取可以通过将物种的基因序列提交给一些在线网站或软件去分析所包含的结构域，如果包含指定基因家族的结构域，则该基因属于指定基因家族，在线网站如NCBI的结构域分析工具(https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi)、EMBL-EBI的Pfam搜索工具(http://pfam.xfam.org/search)等，软件如EMBL-EBI提供的本地HMMER安装程序等；

确定第一考察物种，获取其所述指定基因家族所包含的基因作为第一考察基因，获取方法见上述参考物种指定基因家族基因的获取；

确定第二考察物种，获取其所述指定基因家族所包含的基因作为第二考察基因，获取方法见上述参考物种指定基因家族基因的获取；

将参考基因与第一考察基因进行两两序列比对，从序列比对结果中获取大于阈值条件的第一考察物种的非冗余基因个数及该个数占第一考察基因总数的比例，将该比例作为参考物种的第一考察物种非冗余基因比例，所述阈值条件为匹配序列长度和匹配序列相似性值的组合，所述大于阈值条件为大于匹配序列长度且大于匹配序列相似性值，所述阈值条件的数量大于等于2个。序列比对可以使用ClustalX、ClustalW或者Blast软件进行，序列比对的结果通常包含两两序列间的匹配序列长度和匹配序列相似性值或匹配序列距离值，匹配序列相似性值与匹配序列距离值的关系为它们的和等于1，当比对结果数据量大时可以通过Perl程序编程获取参考物种非冗余基因数量，非冗余基因指不重复的基因；

将参考基因与第二考察基因进行两两序列比对，从序列比对结果中获取大于所述阈值条件的第二考察物种的非冗余基因个数及该个数占第二考察基因总数的比例，将该比例作为参考物种的第二考察物种非冗余基因比例；

计算相同阈值条件时参考物种的第一考察物种非冗余基因比例与参考物种的第二考察物种非冗余基因比例的差值，确定差值的绝对值最大且非冗余基因比例较大的考察物种为与参考物种指定基因家族的关系更密切。

优选地，参考物种和考察物种均已被全基因组测序，全基因组测序后得到的基因家族的基因信息更完整的，则后续分析的结果更可靠。

优选地，所述序列为氨基酸序列。

优选地，所述阈值条件中的匹配序列长度为所述指定基因家族的结构域序列长度的三分之二，匹配序列长度阈值设置太小会降低后续分析结果的可信度，而通常结构域序列长度的三分之二可以认为该结构域比较完整。

优选地，所述阈值条件中的匹配序列相似性值为至少30％以上，通常序列相似性达30％可以认为两条序列可能具有相似的功能。

优选地，所述匹配序列相似性值为至少30％以上包括30％，40％，50％，60％，70％，80％和90％。

优选地，所述阈值条件的数量大于等于2个为7个阈值条件，所述7个阈值条件的匹配序列相似性值分别为30％，40％，50％，60％，70％，80％和90％，匹配序列长度均为所述指定基因家族的结构域序列长度的三分之二。

当然考察物种可以大于2个，通过本发明的方法可以从所有考察物种中获得与参考物种指定基因家族关系最密切的物种，如先从任意两个考察物种中计算出与参考物种指定基因家族关系更密切的考察物种，然后计算此考察物种与另一考察物种谁与参考物种指定基因家族关系更密切，以此类推，获得与参考物种指定基因家族关系最密切的考察物种。

与不同物种间或同种物种内单个基因之间的亲缘关系或进化关系相比，不同物种同一基因家族之间的关系密切程度具有更重要的应用价值，前者仅仅表现基因之间的亲缘关系或进化关系，后者密程度越高表明基因家族的功能表现水平更相近，物种相应的形态发育特征更相似，从而可应用于更快更准更高效地筛选优势物种或目的物种。另外，对于异源杂交物种可反应其与2个亲本物种的某基因家族的遗传进化关系，对于近缘物种可反映近缘物种间某基因家族的进化关系。

附图说明

图1为本发明基于相似基因比例的基因家族关系判定方法优选实施例流程示意图。

具体实施方式

下面将结合实施例来详细说明本发明，这些实施例仅起说明性作用，并不局限于本发明的应用范围。本发明不限于下述实施方式或实施例，凡不违背本发明精神所做出的修改及变形，均应包括在本发明范围之内。

实验例1：海岛棉、亚洲棉与雷蒙德氏棉间NBS基因家族的亲缘关系计算

1.材料来源：亚洲棉和雷蒙德氏棉的基因组数据来自于中国农业科学院棉花研究所(http://cgp.genomics.org.cn/)，分别包含40,134和40,976个基因；海岛棉的基因组数据来自于华中农业大学(http://cotton.cropdb.org/)，包含109,918个基因。

2.方法与步骤

第一，将海岛棉、亚洲棉与雷蒙德氏棉基因组数据中的全部基因的蛋白序列(也是氨基酸序列)提交EMBL-EBI的Pfam搜索工具(http://pfam.xfam.org/search)预测每一条基因所包含的结构域信息。

第二，由于NBS基因家族为包含NB-ARC结构域(结构域编号或Pfam编号为PF00931)的基因，于是通过Perl程序编程分别从第一步中获取海岛棉、亚洲棉与雷蒙德氏棉基因中结构域信息包含PF00931或NB-ARC的基因，即为NBS基因家族成员，获取的海岛棉、亚洲棉与雷蒙德氏棉的NBS基因家族分别包含682、246和365个NBS基因。

第三，设置7个阈值条件，阈值条件的匹配序列相似性值分别为90％、80％、70％、60％、50％、40％和30％，阈值条件的匹配序列长度均为200(由于NBS基因家族的NB-ARC结构域大约有300个氨基酸，设定两两序列比对的匹配序列长度阈值为结构域的三分之二，即200)。

第四，将海岛棉(参考物种)的NBS基因与雷蒙德氏棉(第一考察物种)的NBS基因使用ClustalX软件进行两两序列比对，生成的nj文件包含两两序列间的匹配序列长度和匹配序列距离值，匹配序列相似性值即为1减去匹配序列距离值的得数，再通过Perl程序编程分别获取大于第三步所述的7个阈值条件的雷蒙德氏棉的非冗余基因个数及其占雷蒙德氏棉NBS基因家族基因总数的比例，并将此比例叫做海岛棉的雷蒙德氏棉非冗余基因比例，由于第三步所述阈值条件为7个，因此获得相应的7个海岛棉的雷蒙德氏棉非冗余基因个数及其比例，结果如表1所示，从相似性90％至30％，海岛棉的雷蒙德氏棉非冗余基因个数分别为216、310、339、345、348、350和352，占雷蒙德氏棉NBS总数的比例分别为59％、85％、93％、95％、95％、96％和96％。

第五，将海岛棉(参考物种)的NBS基因与亚洲棉(第二考察物种)的NBS基因进行两两序列比对，并获取海岛棉的亚洲棉非冗余基因个数及其占亚洲棉NBS基因家族基因总数的比例，将此比例叫做海岛棉的亚洲棉非冗余基因比例，比对方法和获取方法与第四步所述的方法相同，获得相应的7个海岛棉的亚洲棉非冗余基因个数及其比例，结果如表1所示，从相似性90％至30％，海岛棉的亚洲棉非冗余基因个数分别为115、185、212、219、222、225和226，占亚洲棉NBS基因总数的比例分别为47％、75％、86％、89％、90％、91％和92％。

第六，计算相同阈值条件时海岛棉的亚洲棉非冗余基因比例与雷蒙德氏棉非冗余基因比例的差值，如表1所示，差值绝对值最大的是相似性大于90％条件下海岛棉的雷蒙德氏棉非冗余基因比例比亚洲棉非冗余基因比例大12％，并且，在各阈值条件下海岛棉的雷蒙德氏棉非冗余基因比例均比亚洲棉非冗余基因比例大，表明雷蒙德氏棉的NBS基因家族与海岛棉的NBS基因家族的关系更密切。

表1不同阈值条件下海岛棉的雷蒙德氏棉非冗余基因与亚洲棉非冗余基因数量

由于NBS基因为抗病相关基因，而黄萎病是影响棉花产量最严重的植物疾病，可称得上是棉花的“癌症”。事实上，研究表明，雷蒙德氏棉和海岛棉均具有较强的抗棉花黄萎病的能力，而亚洲棉却易受黄萎病菌的侵染，本发明方法的分析结果显示海岛棉的NBS基因家族与雷蒙德氏棉的NBS基因家族的关系更密切，表明海岛棉与雷蒙德氏棉具有更相近的NBS基因家族对应的形态发育特征，即抗病水平，与棉花抗病现实情况一致。另外，由于海岛棉为雷蒙德氏棉和亚洲棉间异源杂交形成的新物种，根据本发明的结果表明，海岛棉从雷蒙德氏棉中遗传到了更多的NBS基因家族基因，推测海岛棉的抗病能力与雷蒙德氏棉更相近，与棉花抗病的现实情况一致。可见，不同物种间同一基因家族的密切关系程度跟基因家族功能表现水平和物种相应形态发育特征水平具有紧密的关系，表明本发明方法具有较高的应用可信度。

实验例2：陆地棉、亚洲棉与雷蒙德氏棉间NBS基因家族的亲缘关系计算

1材料来源：陆地棉、亚洲棉与雷蒙德氏棉的基因组数据均来自于中国农业科学院棉花研究所(http://cgp.genomics.org.cn/)，分别包含76,943、40,134和40,976个基因。

2方法与步骤

第一，将陆地棉、亚洲棉与雷蒙德氏棉基因组数据中的全部基因的蛋白序列(也是氨基酸序列)提交EMBL-EBI的Pfam搜索工具(http://pfam.xfam.org/search)预测每一条基因所包含的结构域信息。

第二，由于NBS基因家族为包含NB-ARC结构域(结构域编号或Pfam编号为PF00931)的基因，于是通过Perl程序编程分别从第一步中获取陆地棉、亚洲棉与雷蒙德氏棉基因中结构域信息包含PF00931或NB-ARC的基因，即为NBS基因家族成员，获取的陆地棉、亚洲棉与雷蒙德氏棉的NBS基因家族分别包含682、588和365个NBS基因。

第四，将陆地棉(参考物种)的NBS基因与雷蒙德氏棉(第一考察物种)的NBS基因使用ClustalX软件进行两两序列比对，生成的nj文件包含两两序列间的匹配序列长度和匹配序列距离值，匹配序列相似性值即为1减去匹配序列距离值的得数，再通过Perl程序编程分别获取大于第三步所述的7个阈值条件的雷蒙德氏棉的非冗余基因个数及其占雷蒙德氏棉NBS基因家族基因总数的比例，并将此比例叫做陆地棉的雷蒙德氏棉非冗余基因比例，由于第三步所述阈值条件为7个，因此获得相应的7个陆地棉的雷蒙德氏棉非冗余基因个数及其比例，结果如表2所示，从相似性90％至30％，陆地棉的雷蒙德氏棉非冗余基因个数分别为189、259、280、284、293、296和315，占雷蒙德氏棉NBS基因总数的比例分别为52％、71％、77％、78％、80％、81％和86％。

第五，将陆地棉(参考物种)的NBS基因与亚洲棉(第二考察物种)的NBS基因进行两两序列比对，并获取陆地棉的亚洲棉非冗余基因及其占亚洲棉NBS基因家族基因总数的比例，将此比例叫做陆地棉的亚洲棉非冗余基因比例，比对方法和获取方法与第四步所述的方法相同，获得相应7个陆地棉的亚洲棉非冗余基因个数及比例，结果如表2所示，从相似性90％至30％，陆地棉的亚洲棉非冗余基因个数分别为140、202、211、214、219、221和225，占亚洲棉NBS基因总数的比例分别为57％、82％、86％、87％、89％、90％和91％。

第六，计算相同阈值条件时陆地棉的亚洲棉非冗余基因比例与雷蒙德氏棉非冗余基因比例的差值，如表2所示，差值绝对值最大的是相似性大于80％条件下陆地棉的亚洲棉非冗余基因比例比雷蒙德氏棉非冗余基因比例大11％，并且，在各阈值条件下陆地棉的亚洲棉非冗余基因比例均比雷蒙德氏棉非冗余基因比例大，表明亚洲棉的NBS基因家族与陆地棉的NBS基因家族的关系更密切。

表2不同阈值条件陆地棉的的雷蒙德氏棉非冗余基因与亚洲棉非冗余基因数量

由于NBS基因为抗病相关基因，事实上，研究表明，亚洲棉和陆地棉均易受黄萎病菌的侵染，而雷蒙德氏棉却能够免疫，本发明方法的分析结果显示陆地棉的NBS基因家族与亚洲棉的NBS基因家族的关系更密切，表明陆地棉与亚洲棉具有更相近的NBS基因家族对应的形态发育特征，即抗病水平，与棉花抗病现实情况一致。另外，由于陆地棉均为雷蒙德氏棉和亚洲棉间异源杂交形成的新物种，根据本发明的结果表明，陆地棉从亚洲棉中遗传到了更多的NBS基因家族基因，推测陆地棉的抗病水平与亚洲棉更相近，与棉花抗病的现实情况一致。同样表明，不同物种间同一基因家族的关系密切程度跟基因家族功能表现水平和物种相应形态发育特征水平具有紧密的关系，表明本发明方法具有较高的应用可信度。

实验例3：可可、雷蒙德氏棉与亚洲棉间NBS基因家族的亲缘关系计算

1材料来源：雷蒙德氏棉和亚洲棉的基因组数据来自于中国农业科学院棉花研究所(http://cgp.genomics.org.cn/)，分别包含40,976和40,134个基因；可可的基因组数据来自于可可基因组中心(http://cocoagendb.cirad.fr./)，包含46,143个基因。

2方法与步骤

第一，将可可、雷蒙德氏棉与亚洲棉基因组数据中的全部基因的蛋白序列(也是氨基酸序列)提交EMBL-EBI的Pfam搜索工具(http://pfam.xfam.org/search)预测每一条基因所包含的结构域信息。

第二，由于NBS基因家族为包含NB-ARC结构域(结构域编号或Pfam编号为PF00931)的基因，于是通过Perl程序编程分别从第一步中获取可可、雷蒙德氏棉与亚洲棉基因中结构域信息包含PF00931或NB-ARC的基因，即为NBS基因家族成员，获取的海岛棉、陆地棉与亚洲棉的NBS基因家族分别包含298、365和246个NBS基因。

第四，将可可(参考物种)的NBS基因与雷蒙德氏棉(第一考察物种)的NBS基因使用ClustalX软件进行两两序列比对，生成的nj文件包含两两序列间的匹配序列长度和匹配序列距离值，匹配序列相似性值即为1减去匹配序列距离值的得数，再通过Perl程序编程分别获取大于第三步所述的7个阈值条件的雷蒙德氏棉的非冗余基因个数及其占雷蒙德氏棉NBS基因家族基因总数的比例，并将此比例叫做可可的雷蒙德氏棉非冗余基因比例，由于第三步所述阈值条件为7个，因此获得相应的7个可可的雷蒙德氏棉非冗余基因个数及其比例，结果如表3所示，从相似性90％至30％，可可的雷蒙德氏棉非冗余基因个数分别为0、14、43、91、262、338和346，占雷蒙德氏棉NBS基因总数的比例分别为0％、4％、12％、25％、72％、93％和95％。

第五，将可可(参考物种)的NBS基因与亚洲棉(第二考察物种)的NBS基因进行两两序列比对，并获取可可的亚洲棉非冗余基因个数及其占亚洲棉NBS基因家族基因总数的比例，比对方法和获取方法与第四步所述的方法相同，获得相应7个可可的亚洲棉非冗余基因个数及其比例，结果如表3所示，从相似性90％至30％，可可的亚洲棉非冗余基因个数分别为0、10、32、99、171、217和222，占亚洲棉NBS基因总数的比例分别为0％、4％、13％、40％、69％、88％和90％。

第六，计算相同阈值条件时可可的雷蒙德氏棉非冗余基因比例与亚洲棉非冗余基因比例的差值，如表3所示，差值最大的是相似性大于60％条件下可可的亚洲棉非冗余基因比例比雷蒙德氏棉非冗余基因比例大15％，虽然在相似性大于50、40、30％的阈值条件下可可的雷蒙德氏棉非冗余基因比例均比亚洲棉非冗余基因比例大，但差值相当小，没有显著性，因此确定以差值绝对值最大且非冗余基因比例较大的亚洲棉为与可可的NBS基因家族的关系更密切。

表3不同阈值条件下雷蒙德氏棉与亚洲棉的可可非冗余基因的数量

由于NBS基因为抗病相关基因，事实上，研究表明，亚洲棉和可可均易受黄萎病菌的侵染，而雷蒙德氏棉却能够免疫，本发明方法的分析结果显示可可的NBS基因家族与亚洲棉的NBS基因家族的关系更密切，表明可可与亚洲棉具有更相近的NBS基因家族对应的形态发育特征，即抗病水平，与棉花抗病现实情况一致，表明本发明方法具有较高的应用可信度。

Claims

1.基于相似基因比例的基因家族关系判定方法，其特征在于：

确定一个待分析的基因家族作为指定基因家族；

选择参考物种，获取其所述指定基因家族所包含的基因，作为参考基因；

确定第一考察物种，获取其所述指定基因家族所包含的基因作为第一考察基因；

确定第二考察物种，获取其所述指定基因家族所包含的基因作为第二考察基因；

将参考基因与第一考察基因进行两两序列比对，从该序列比对结果中获取大于阈值条件的第一考察物种的非冗余基因个数及该个数占第一考察基因总数的比例，将该比例作为参考物种的第一考察物种非冗余基因比例，所述阈值条件为匹配序列长度和匹配序列相似性值的组合，所述大于阈值条件为大于匹配序列长度且大于匹配序列相似性值，所述阈值条件的数量大于等于2个；

将参考基因与第二考察基因进行两两序列比对，从该序列比对结果中获取大于所述阈值条件的第二考察物种的非冗余基因个数及该个数占第二考察基因总数的比例，将该比例作为参考物种的第二考察物种非冗余基因比例；

2.如权利要求1所述的基于相似基因比例的基因家族关系判定方法，其特征在于，参考物种和考察物种均已被全基因组测序。

3.如权利要求1所述的基于相似基因比例的基因家族关系判定方法，其特征在于，所述序列为氨基酸序列。

4.如权利要求1所述的基于相似基因比例的基因家族关系判定方法，其特征在于，所述阈值条件中的匹配序列长度为所述指定基因家族的结构域序列长度的三分之二。

5.如权利要求1所述的基于相似基因比例的基因家族关系判定方法，其特征在于，所述阈值条件中的匹配序列相似性值为至少30％以上。

6.如权利要求5所述的基于相似基因比例的基因家族关系判定方法，其特征在于，所述匹配序列相似性值为至少30％以上包括30％，40％，50％，60％，70％，80％和90％。

7.如权利要求1所述的基于相似基因比例的基因家族关系判定方法，其特征在于，所述阈值条件的数量大于等于2个为7个阈值条件，所述7个阈值条件的匹配序列相似性值分别为30％，40％，50％，60％，70％，80％和90％，匹配序列长度均为所述指定基因家族的结构域序列长度的三分之二。