CN101957892B - 一种全基因组复制事件的检测方法和系统 - Google Patents
一种全基因组复制事件的检测方法和系统 Download PDFInfo
- Publication number
- CN101957892B CN101957892B CN2010102849664A CN201010284966A CN101957892B CN 101957892 B CN101957892 B CN 101957892B CN 2010102849664 A CN2010102849664 A CN 2010102849664A CN 201010284966 A CN201010284966 A CN 201010284966A CN 101957892 B CN101957892 B CN 101957892B
- Authority
- CN
- China
- Prior art keywords
- gene
- homologous
- protein sequence
- gene set
- genome
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000001514 detection method Methods 0.000 title claims abstract description 43
- 230000010076 replication Effects 0.000 title abstract description 12
- 108090000623 proteins and genes Proteins 0.000 claims abstract description 480
- 102000004169 proteins and genes Human genes 0.000 claims abstract description 131
- 230000002068 genetic effect Effects 0.000 claims abstract description 38
- 239000002773 nucleotide Substances 0.000 claims description 72
- 125000003729 nucleotide group Chemical group 0.000 claims description 72
- 210000000349 chromosome Anatomy 0.000 claims description 11
- 238000006243 chemical reaction Methods 0.000 claims description 8
- 230000004069 differentiation Effects 0.000 claims description 8
- 238000005259 measurement Methods 0.000 claims description 8
- 238000012216 screening Methods 0.000 claims description 8
- 108020004705 Codon Proteins 0.000 claims description 6
- 108020004707 nucleic acids Proteins 0.000 claims description 4
- 102000039446 nucleic acids Human genes 0.000 claims description 4
- 150000007523 nucleic acids Chemical class 0.000 claims description 4
- 238000000034 method Methods 0.000 abstract description 33
- 230000008901 benefit Effects 0.000 abstract description 2
- 108091028043 Nucleic acid sequence Proteins 0.000 abstract 1
- 241000894007 species Species 0.000 description 24
- 244000061456 Solanum tuberosum Species 0.000 description 14
- 235000002595 Solanum tuberosum Nutrition 0.000 description 14
- 238000010586 diagram Methods 0.000 description 10
- 230000035772 mutation Effects 0.000 description 10
- 230000008569 process Effects 0.000 description 10
- 230000006870 function Effects 0.000 description 7
- 241000196324 Embryophyta Species 0.000 description 5
- 238000013519 translation Methods 0.000 description 5
- 101150112970 up gene Proteins 0.000 description 5
- 241000219194 Arabidopsis Species 0.000 description 4
- 241000208292 Solanaceae Species 0.000 description 4
- 230000008711 chromosomal rearrangement Effects 0.000 description 4
- 239000011159 matrix material Substances 0.000 description 4
- 108091081024 Start codon Proteins 0.000 description 3
- 230000008859 change Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000015654 memory Effects 0.000 description 3
- 210000003205 muscle Anatomy 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 101150040974 Set gene Proteins 0.000 description 2
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 230000010429 evolutionary process Effects 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000003780 insertion Methods 0.000 description 2
- 230000037431 insertion Effects 0.000 description 2
- 230000035800 maturation Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000007935 neutral effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000012549 training Methods 0.000 description 2
- 101150096316 5 gene Proteins 0.000 description 1
- 244000189799 Asimina triloba Species 0.000 description 1
- 235000006264 Asimina triloba Nutrition 0.000 description 1
- 241000219193 Brassicaceae Species 0.000 description 1
- 235000009467 Carica papaya Nutrition 0.000 description 1
- 240000008067 Cucumis sativus Species 0.000 description 1
- 235000010799 Cucumis sativus var sativus Nutrition 0.000 description 1
- 244000068988 Glycine max Species 0.000 description 1
- 235000010469 Glycine max Nutrition 0.000 description 1
- 235000007688 Lycopersicon esculentum Nutrition 0.000 description 1
- 241000218922 Magnoliophyta Species 0.000 description 1
- 240000003768 Solanum lycopersicum Species 0.000 description 1
- 235000009754 Vitis X bourquina Nutrition 0.000 description 1
- 235000012333 Vitis X labruscana Nutrition 0.000 description 1
- 240000006365 Vitis vinifera Species 0.000 description 1
- 235000014787 Vitis vinifera Nutrition 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000008034 disappearance Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 244000144992 flock Species 0.000 description 1
- 230000004545 gene duplication Effects 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 238000009396 hybridization Methods 0.000 description 1
- 238000009342 intercropping Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000008521 reorganization Effects 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
Images
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)
Abstract
本发明公开一种全基因组复制的检测方法和系统。该方法包括：获取与目标物种基因组上基因集的蛋白质序列以及基因在基因组上的位置信息；对基因集的蛋白质序列进行蛋白比蛋白两两比对获得基因集的同源基因对；利用动态规划获得基因组上具有共线性的同源区块；获得同源区块中同源基因对的核酸序列的同源位点对位排列；计算同源区块中同源基因对的遗传距离，获得年龄分布；根据年龄分布确定全基因组复制事件。本发明的方法和系统，全面考虑了全基因组上的所有基因的序列信息，排除了基因组上非保守和选择压力受到释放的区域的干扰，并利用动态规划算法实现最佳复制区域的寻找，确定的同源基因对和同源区块全面、准确，从而能够更准确地实现全基因组复制事件的检测。
Description
技术领域
本发明涉及生物信息学技术领域，尤其涉及一种全基因组复制事件的检测方法和系统。
背景技术
全基因组复制事件(Whole Genome Duplication，WGD)，也称古多倍化事件，在植物界是一个频繁发生和普遍存在的现象，这些全基因组范围内的复制对植物的多样化和适应性产生了广泛而深远的影响。随着大量的植物物种全基因组被测序，以及对其基因组内功能元件的研究，推断基因组的演变历史及其对物种形成和多样化的影响已越来越受到人们的关注。基因组上的序列信息蕴含了一个物种所有进化历史所遗留下来的痕迹，因此，鉴定基因组上序列的进化过程和变化模型为人们了解该物种提供了很好的信息。在基因组测序技术成熟之前，一般的方法是通过估计相关物种基因组大小，染色体数目的变化，或通过染色体荧光定位杂交的方法来推断某植物物种是否经历过全基因组复制。但是，由于全基因组复制事件是一个特别复杂的动态过程，并且随后会伴随大量基因的丢失和染色体重排等复杂变化，该方法检测出来的结果只能提供一小部分信息或仅针对少数物种，或只能针对较近期发生的比较明显的全基因组复制事件。随着测序技术的逐渐成熟，出现基于EST(序列表达标签)数据鉴定全基因组复制事件的方法，但该方法也只是基于能够表达的在基因组上比较保守的基因区域，由于表达数据的获取不全面或者说表达数据只是基因组序列水平上的一个间接产物，该方法鉴定出来的植物基因组内全基因组复制事件的次数和时间往往模糊而不能确定。
发明内容
本发明要解决的一个技术问题是提供一种全基因组复制事件的检测方法，使得确定的全基因组复制事件更准确。
本发明提供一种全基因组复制事件的检测方法，包括：
获取与目标物种基因组上基因集的核酸序列和对应的基因集的蛋白质序列以及基因在基因组上的位置信息；
对基因集的蛋白质序列进行蛋白比蛋白两两比对获得基因集的同源基因对；
根据基因集中基因在基因组上的位置信息以及基因集的同源基因对，利用动态规划获得基因组上具有共线性的同源区块；
获得同源区块中同源基因对的核酸序列的同源位点对位排列；
计算同源区块中同源基因对的遗传距离，获得年龄分布；
根据年龄分布确定全基因组复制事件。
根据本发明方法的一个实施例，根据基因集中基因在基因组上的位置信息以及基因集的同源基因对利用动态规划获得基因组的具有共线性的同源区块的步骤包括：根据基因集的同源基因对，形成所有基因间的同源关系矩阵；利用基因集中每个基因在染色体上的位置信息，利用聚类算法，将相互连锁的共线性的同源基因块聚集在一起，形成具有共线性的同源区块。
根据本发明方法的一个实施例，获取与目标物种基因组上基因集的核酸序列对应的基因集的蛋白质序列的步骤包括：获取目标物种基因组上基因集的核酸序列；根据密码子表将基因集的核酸序列转化为基因集的蛋白质序列。
根据本发明方法的一个实施例，获取与目标物种基因组上基因集的核酸序列对应的基因集的蛋白质序列的步骤还包括：当基因集中的基因有多个转录本时，选择最长的转录本对应的蛋白质序列作为基因的核酸序列对应蛋白质序列。
根据本发明方法的一个实施例，将基因集的蛋白质序列进行蛋白比蛋白两两比对获得基因集的同源基因对的步骤包括：通过blastp对基因集的蛋白质序列进行蛋白比蛋白两两比对，设置p-value＜＝1e-7；对比对输出的同源基因对通过e-value＜＝1e-10进行筛选；当同一基因对有多个比对结果仅以不同的e-value值给出时，选择e-value值最小的比对结果作为同源基因对。
根据本发明方法的一个实施例，获得同源区块中同源基因对的核酸序列的同源位点对位排列的步骤包括：对同源区块上的同源基因对的蛋白质序列进行全局序列比对获得同源基因对的蛋白质序列的同源位点对位排列；根据同源基因对的蛋白质序列的同源位点对位排列构建基因集的核酸序列的同源位点对位排列。
根据本发明方法的一个实施例，根据同源基因对的蛋白质序列的同源位点对位排列构建基因集的核酸序列的同源位点对位排列的步骤包括：根据同源基因对的蛋白质序列的同源位点对位排列构建基因集的核酸序列的同源位点对位排列；其中，利用三联密码子将蛋白质“反翻译”(CDS-back-translation)成相应核酸序列后，关注的同源位点中不包含起始密码子和终子密码子。
根据本发明方法的一个实施例，计算同源区块中同源基因对的遗传距离建立年龄分布图的步骤包括：每一个同源比对块上的每对同源基因对，计算遗传距离，对同一个同源区块上的所有基因对的遗传距离计算一个平均值；将遗传距离作为其分化年龄的一种衡量，建立年龄分布图。
通过本发明方法的实施例，全面考虑了全基因组上的所有基因的序列信息，排除了基因组上非保守和选择压力受到释放的区域的干扰，并利用动态规划算法实现最佳复制区域的寻找，确定的同源基因对和同源区块全面、准确，从而能够更准确地实现全基因组复制事件的检测。
本发明要解决的一个技术问题是提供一种全基因组复制事件的检测系统，使得确定的全基因组复制事件更准确。
本发明提供一种全基因组复制事件的检测系统，包括：
蛋白质序列获取模块，用于获取与目标物种基因组上基因集的核酸序列对应的基因集的蛋白质序列以及基因在基因组上的位置信息；
蛋白质序列比对模块，用于将基因集的蛋白质序列进行蛋白比蛋白两两比对，获得基因集的同源基因对；
同源区块确定模块，用于根据基因集中基因在基因组上的位置信息以及基因集的同源基因对利用动态规划获得基因组的具有共线性的同源区块；
对位排列确定模块，用于获得同源区块中同源基因对的核酸序列的同源位点对位排列；
年龄分布获取模块，用于计算同源区块中同源基因对的遗传距离，获得年龄分布；
全基因组复制确定模块，用于根据年龄分布确定全基因组复制事件。
根据本发明系统的一个实施例，蛋白质序列获取模块包括：核酸序列获取单元，用于获得目标物种基因序列的基因集的核酸序列以及基因在基因组上的位置信息，发送基因集的核酸序列；蛋白质序列转化单元，用于接收基因集的核酸序列，根据密码子表将基因集中各个基因的核酸序列转化为蛋白质序列，获得基因集的蛋白质序列。
根据本发明系统的一个实施例，蛋白质序列比对模块包括：
同源基因对获取单元，用于通过BLASTP蛋白质组比对对基因集的蛋白质序列进行蛋白比蛋白两两比对，获得同源基因对，其中blastp蛋白质组比对的参数设置为p-value＜＝1e-7；
同源基因对过滤单元，用于同源基因对获取单元获得的同源基因对进行筛选，筛选的参数为设置为e-value＜＝1e-10。
根据本发明系统的一个实施例，对位排列确定模块包括：
蛋白质对位排列获取单元，用于对同源区块上的同源基因对的蛋白质序列进行全局序列比对获得同源基因对的蛋白质序列的同源位点对位排列；
核酸对位排列获取单元，用于根据同源基因对的蛋白质序列的同源位点对位排列构建基因集的核酸序列的同源位点对位排列。
通过本发明系统的实施例，蛋白质序列获取模块获取与目标物种基因组上基因集的核酸序列对应的基因集的蛋白质，全面考虑了全基因组上的所有基因的序列信息，排除了基因组上非保守和选择压力受到释放的区域的干扰，同源区块确定模块利用动态规划算法和有效的聚类算法实现最佳复制区域的寻找，确定的同源基因对和同源区块全面、准确，能够更准确地实现全基因组复制事件的检测。
附图说明
图1示出本发明全基因组复制事件的检测方法的一个实施例的流程图；
图2示出本发明全基因组复制事件的检测方法的另一个实施例的流程图；
图3示出本发明全基因组复制事件的检测系统的一个实施例的框图；
图4示出本发明全基因组复制事件的检测系统的另一个实施例的框图。
图5示出本发明全基因组复制事件的检测方法的一个应用例的流程图；
图6示出根据图5的应用例获得的土豆全基因组复制时间的4DTV图。
具体实施方式
下面参照附图对本发明进行更全面的描述，其中说明本发明的示例性实施例。
图1示出本发明全基因组复制事件的检测方法的一个实施例的流程图。
如图1所示，在步骤102，获取与目标物种基因组上基因集的核酸序列对应的基因集的蛋白质序列以及基因在基因组上的位置信息。获取目标物种基因组序列的基因注释信息，例如，基因集中各个基因的核酸序列和各个基因在基因组上的位置信息；根据密码子表将基因集中各个基因的核酸序列转化为该基因的蛋白质序列。基因集指物种基因组包括的所有基因的集合。
步骤104，对基因集的蛋白质序列进行蛋白比蛋白两两比对获得基因集的同源基因对。将基因集中每个基因的蛋白质序列和其他每个基因的蛋白质序列进行比对，获得基因集的同源基因对。
步骤106，根据基因集中基因的位置信息以及基因集的同源基因对，利用动态规划获得基因组上具有共线性的同源区块。通过步骤104的基因的同源比对，获得基因集的同源基因对，形成了一个所有基因间的同源关系矩阵。然后，利用每个基因在染色体上的位置信息，利用聚类算法，将相互连锁的共线性的同源基因块聚集在一起，即将单对的同源基因对聚集成多对的共线性的连锁基因块，即具有共线性的同源区块。由于基因组的进化过程中会涉及到大量的染色体重组，缺失等因素，所以要保证所鉴定出来的共线性区块确实是由于全基因组复制所遗留下来的信号，在这里，相邻基因对之间的距离不能太远(一般平均设为10k)，整个共线性基因区块的长度要足够大(例如，至少＞＝100K)来确定一个基因区块是否为同源区块。
步骤108，获得同源区块中同源基因对的核酸序列的同源位点对位排列。对于所有同源区块上的同源基因对，用其蛋白质序列进行全局序列比对(MUSCLE)，获得同源基因对的蛋白质序列的同源位点对位排列方式，然后将蛋白质序列比对反溯至核酸序列水平上去，构建对应的同源基因对的核酸序列的同源位点对位排列方式。
步骤110，计算同源区块中同源基因对的遗传距离，获得年龄分布。例如，对于同源区块上的每对同源基因对，计算其遗传距离(也即分歧度，ks或4dtv值)，对同一个同源区块上的所有同源基因对的遗传距离进行平均，以其平均值作为该同源区块的遗传距离。也可以通过例如加权平均的方法根据每个同源基因对的遗传距离获得同源区块的遗传距离。
步骤112，根据年龄分布确定全基因组复制事件。由于每一次全基因组复制(随后可能会发生基因丢失、染色体重排等事件)都是同源基因对的一次大量涌现，在分布统计上会对应地形成一个显著的峰，据此可对全基因组复制事件发生的次数和时间形成判断。
在上述实施例中，全面考虑了全基因组上的所有基因的序列信息，排除了基因组上非保守和选择压力受到释放的区域的干扰，并利用动态规划算法和有效的聚类算法实现最佳复制区域的寻找，确定的同源基因对和同源区块全面、准确，建立基因序列最有可能的进化模型，从而能够更准确地实现全基因组复制事件的检测。
图2示出本发明全基因组复制事件的检测方法的另一个实施例的流程图。
如图2所示，在步骤202，获取目标物种基因组序列的基因注释信息，包括基因集核酸序列；根据密码子表将基因集的核酸序列转化为基因集的蛋白质序列。当基因集中的基因有多个转录本时，选择最长的转录本对应的蛋白质序列作为该基因的核酸序列对应蛋白质序列。
在步骤204，将基因集的蛋白质序列进行蛋白比蛋白两两比对，对比对结果进行过滤后获得基因集的同源基因对。该过滤包括去除基因的蛋白质序列自身的比对结果，去除同一对基因中多重比对的结果，只保留一个基因对的多个比对结果中的最佳同源基因对。在blastp比对结果中，可以根据e-value值参数进行过滤，有时也需要结合相似度和覆盖度。通过对比对结果进行过滤，使得最终结果可靠性更高。
步骤206，获取每个基因在基因组上的位置信息，若基因在负链上，按5’-＞3’的方向给出位置。以基因为锚点利用动态规划算法找出所有基因组内部具有共线性的同源区块。对步骤204中所有的同源基因对利用基于图论的聚类算法，根据打分矩阵，将局部匹配的比对区域聚集在一起，形成最高分值的大比对区块。在确定同源区块的过程中考虑相邻基因之间的距离，基因缺失和插入所带来的影响，以及最小区块的设置。可以根据植物基因组的进化特点，进行了有效的模型设计。
步骤208，对于所有同源区块上的同源基因对，用其蛋白质序列例如通过MUSCLE进行全局序列比对，获得最佳同源位点对位排列方式。同源位点对位排列方式是指将可能是由同一个祖先传下来的同源区域或同源碱基位点对齐。在MUSLCE算法，可以根据打分获得一个最高分的同源位点对位排列结果，这个最高分值的同源位点对位排列结果作为最佳同源位点对位排列方式。
步骤210，将蛋白质序列比对反溯至核酸序列水平上去，构建对应的核酸序列对位排列方式。根据密码子将蛋白质“反翻译”回核酸序列(CDS-back-translation)，充分利用核酸序列的蛋白序列之间的对应密码子，去除起始密码子和终子密码子。对于残缺基因(有移码或中间有终子密码子的基因)在序列准备阶段给以抛弃。
步骤212，对于每个同源区块上的每对基因对，计算其遗传距离(也即分歧度，ks或4dtv值)，对同一同源区块上的所有基因对的遗传距离计算一个平均值，作为给同源区块的遗传距离。遗传距离的衡量，一个是通过ks(同义突变位点上的平均同义突变率)，它主要是基于分子钟假设，在完全中性的条件下基因的突变速率的衡量；而由于平等突变、回复突变等的影响，ks的衡量标准往往会受到影响而失真，所以需要对其进行矫正，第三位密码子四重简并位点的计算(4dtv，并考虑转换、颠换等参数)会更加接近真实的进化速率。对遗传距离的计算时，可以针对物种基因的特点进行参数训练，建立基因序列最有可能的进化模型，并对距离计算中可能会带来噪音的一些假设进行矫正。
步骤214，将同源区块的遗传距离作为其分化年龄的一种衡量，建立年龄分布图。
步骤216，每一次全基因组复制(伴随着随后的基因丢失，染色体重排等事件)都是同源基因对的一次大量涌现，在分布图上会对应地形成一个显著的峰，据此可对全基因组复制事件发生的次数和时间形成一个直观的认识和判断。
随后，可根据需要，逐步分析每次全基因组复制事件中产生的同源基因的进化趋势或命运，分析它们的功能分化或表达差异，分析它们分析所受到的选择压力和丢失、死亡等事件。
在上述实施例中，在鉴定同源区块时，首先以基因为锚点利用了高效的动态规划算法找出所有可能的最佳同源基因对，并利用基于图论的马尔可夫聚类算法对各对同源基因进行聚类，以得到所有可能的同源区块。处理过程中利用了一系列加权采样值来表示每对基因可能所在位置的后验概率分布，算法每一步都包括基因位置预测和位置更新两个阶段。位置预测阶段是利用m个加权采样值对后验概率分布进行描述的过程，位置更新阶段则是通过重要性采样操作对其进行及时不断更新，最终达到最佳同源区块的鉴定。此外，对遗传距离的计算时，针对物种基因的特点进行参数训练，建立基因序列最有可能的进化模型，并对距离计算中可能会带来噪音的一些假设进行矫正(如充分排除平行突变、回复突变，考虑转换巅换比及考虑基因所受到的选择压力等)，进一步保证了全基因组复制事件检测的准确性。
图3示出本发明全基因组复制事件的检测系统的一个实施例的框图。如图3所示，该实施例的检测系统包括蛋白质序列获取模块31、蛋白质序列比对模块32、同源区块确定模块33、对位排列确定模块34、年龄分布获取模块35和全基因组复制确定模块36。
在图3的实施例中，蛋白质序列获取模块31获取与目标物种基因组上基因集的核酸序列对应的基因集的蛋白质序列以及基因在基因组上的位置信息。例如，蛋白质序列获取模块31获得目标物种基因序列的注释信息，其中包括基因集的核酸序列，然后根据密码子表将基因集中各个基因的核酸序列转化为该基因的蛋白质序列。当基因集中的基因有多个转录本时，选择最长的转录本对应的蛋白质序列作为基因的核酸序列对应蛋白质序列。蛋白质序列比对模块32将基因集的蛋白质序列进行蛋白比蛋白两两比对，获得基因集的同源基因对。蛋白质序列比对模块32将基因集中每个基因的蛋白质序列和分别和其他基因的蛋白质序列进行比对，从而获得基因集的同源基因对。同源区块确定模块33根据基因集基因在基因组上的位置信息以及基因集的同源基因对利用动态规划获得基因组的具有共线性的同源区块。例如，同源区块确定模块33根据蛋白质序列比对模块32获得的基因集的同源基因对，形成一个所有基因间的同源关系矩阵；然后利用每个基因在染色体上的位置信息，利用聚类算法，将相互连锁的共线性的同源基因块聚集在一起，即将单对的同源基因对聚集成多对的共线性的连锁基因块，即具有共线性的同源区块。对位排列确定模块34获得同源区块中同源基因对的核酸序列的同源位点对位排列。对位排列确定模块34对所有同源区块上的同源基因对，利用其蛋白质序列进行全局序列比对，获得同源基因对的蛋白质序列的同源位点对位排列方式，然后将蛋白质序列比对反溯至核酸序列水平上去，构建对应的彤云基因对的核酸序列的同源位点对位排列方式。年龄分布获取模块35计算同源区块中同源基因对的遗传距离，获得年龄分布。例如，年龄分布获取模块35对于同源区块上的每对同源基因对，计算其遗传距离，对同一个同源区块上的所有同源基因对的遗传距离进行平均，以其平均值作为该同源区块的遗传距离。全基因组复制确定模块36根据年龄分布确定全基因组复制事件。
在上述实施例中，蛋白质序列获取模块获取与目标物种基因组上基因集的核酸序列对应的基因集的蛋白质，全面考虑了全基因组上的所有基因的序列信息，排除了基因组上非保守和选择压力受到释放的区域的干扰，同源区块确定模块利用动态规划算法和有效的聚类算法实现最佳复制区域的寻找，确定的同源基因对和同源区块全面、准确，建立基因序列最有可能的进化模型，从而能够更准确地实现全基因组复制事件的检测。
图4示出本发明全基因组复制事件的检测系统的另一个实施例的框图。如图4所示，该实施例的检测系统包括蛋白质序列获取模块41、蛋白质序列比对模块42、同源区块确定模块33、对位排列确定模块44、年龄分布获取模块35和全基因组复制确定模块36。同源区块确定模块33、年龄分布获取模块35和全基因组复制确定模块36的功能和实现可以参见图3中对应模块的描述，为简洁起见在此不再详细描述。下面具体描述蛋白质序列获取模块41、蛋白质序列比对模块42、和对位排列确定模块44的实现。
在图4的实施例中，蛋白质序列获取模块41包括核酸序列获取单元411和蛋白质序列转化单元412。核酸序列获取单元411获得目标物种基因序列的基因集的核酸序列以及基因在基因组上的位置信息，将基因集的核酸序列发送给蛋白质序列转化单元412。蛋白质序列转化单元412接收基因集的核酸序列，根据密码子表将基因集中各个基因的核酸序列转化为蛋白质序列，从而获得基因集的蛋白质序列。蛋白质序列比对模块42包括同源基因对获取单元421和同源基因对过滤单元422。同源基因对获取单元421用于通过BLASTP蛋白质组比对对基因集的蛋白质序列进行蛋白比蛋白两两比对，获得同源基因对，其中blastp蛋白质组比对的参数设置为p-value＜＝1e-7。同源基因对过滤单元422对同源基因对获取单元421获得的同源基因对进行筛选，筛选的参数为设置为e-value＜＝1e-10。对位排列确定模块44包括蛋白质对位排列获取单元441和核酸对位排列获取单元442。蛋白质对位排列获取单元441对同源区块上的同源基因对的蛋白质序列进行全局序列比对获得同源基因对的蛋白质序列的同源位点对位排列。核酸对位排列获取单元442根据蛋白质对位排列获取单元411输出的同源基因对的蛋白质序列的同源位点对位排列来构建基因集的核酸序列的同源位点对位排列。
对于图3至图4中各个装置或单元的功能，可以参考上文中关于本发明方法的实施例中对应部分的说明，为简洁起见，在此不再详述。
本领域的技术人员应当理解，对于图3至图4中的各个装置，可以通过单独的计算处理设备实现，或者将其集成为一个独立的设备实现。在图3至图4中用框示出以说明它们的功能。这些功能块可以用硬件、软件、固件、中间件、微代码、硬件描述语音或者它们的任意组合来实现。举例来说，一个或者两个功能块都可以利用运行在微处理器、数字信号处理器(DSP)或任何其他适当计算设备上的代码实现。代码可以表示过程、功能、子程序、程序、例行程序、子例行程序、模块或者指令、数据结构或程序语句的任意组合。代码可以位于计算机可读介质中。计算机可读介质可以包括一个或者多个存储设备，例如，包括RAM存储器、闪存存储器、ROM存储器、EPROM存储器、EEPROM存储器、寄存器、硬盘、移动硬盘、CD-ROM或本领域公知的其他任何形式的存储介质。计算机可读介质还可以包括编码数据信号的载波。
本领域技术人员将意识到硬件、固件和软件配置在这些情况下的可替换性，以及如何最好地实现每个特定应用地所述功能。
图5示出本发明全基因组复制事件的检测方法的一个应用例的流程图。该应用例通过本发明的检测方法鉴定土豆全基因组复制事件。土豆是茄科植物中非常重要的一个物种，其基因组大小约760M，功能基因有约4320个。每一次全基因组复制事件都会是同源基因数目的大爆发，其后伴随的大量基因的丢失和染色体重排虽然会抹去许多全基因组复制的痕迹，但功能基因在整体上会保留这些曾经发生过大爆发的证据。当把所有的今天依然存留下来的同源区块都鉴定出来后，计算每对同源区块之间分化的可能时间(遗传距离)，根据遗传距离累积所反映出来的“年龄”分布，就可判断时全基因复制事件是否发生和可能发生的时间。具体步骤如下：
步骤502，获得土豆基因集的核酸序列和对应蛋白质序列。土豆基因集包括4320个基因。如果一个基因有多个转录本，则选择其最长的转录本对应的蛋白质序列。
步骤504，blastp蛋白质组比对，比对时一般可设置p-value＜＝1e-7；比对输出结果中可以e-value＜＝1e-10来筛选，进一步保证比对结果的可信度。如果包括自身比对，则要过滤自身比对结果；若同一对基因对有多个比对结果只是以不同的e-value值给出，选择其最优比对结果(即e-value值最小)。土豆的蛋白质组比对结果过滤后有1,716,931条纪录。
步骤506，获取每个基因在染色体上的位置信息，若基因在负链上，按5’-＞3’的方向给出位置。一般可从基因注释结果中提取基因在染色体上的位置信息。
步骤508，以基因为锚点，利用以上基因对比上的结果和位置信息，实现动态规划算法获取最优局部匹配。对步骤504得到的比对结果利用基于图论的聚类算法，根据打分矩阵，将局部匹配的比对区域聚集在一起，形成两边延伸的大比对区块，即同源区块。将重点考虑相邻基因之间的距离，基因缺失和插入所带来的影响，以及最小区块的设置。土豆分析中相邻基因之间的最小距离是10k，每个区块中至少要有5个基因对。结果有284个大区块被鉴定出来，平均每个区块有9个基因，共2556个基因对涉及到3200个土豆基因。
步骤510，对2556个同源基因对的蛋白质序列进行全局比对，形成蛋白质的同源位点对位排列，随之转换为基于核酸序列的全局同源位点对位排列方式。蛋白质比对时选用MUSCLE全局比对软件包处理；CDS-back-translation中，要充分利用核酸序列的蛋白序列之间的对应密码子，并去除起始密码子和终子密码子。对于残缺基因(有移码或中间有终子密码子的基因)在序列准备阶段给以抛弃。
步骤512，基于核酸序列的比对结果，计算每个区块中每对同源基因的遗传距离。遗传距离的衡量，一个是通过ks(同义突变位点上的平均同义突变率)，它主要是基于分子钟假设，在完全中性的条件下基因的突变速率的衡量；而由于平等突变、回复突变等的影响，ks的衡量标准往往会受到影响而失真，所以需要对其进行矫正，第三位密码子四重简并位点的计算(4dtv，并考虑转换、颠换等参数)会更加接近真实的进化速率。
步骤514，对每个区块中的多对遗传距离值计算一个平均值，作为该对区块分化时间的一个衡量。对所有区块分化时间作一个累积分布图。
步骤516，每一次全基因组复制事件都是重复区块的一个大涌现，根据分布图上的峰值可以判断全基因组复制事件发生的时间和发生的次数。
通过图5所示的应用例，土豆全基因组复制事件至少有两次被鉴定出来(如图6所示)。一次是在4dtv＝～0.33处，这是最近的一次全基因组复制事件；一次是在4dtv＝～0.75处，这是古老的一次全基因组复制事件。从图6中也可以看出，土豆古老的那次全基因组复制事件发生在与拟南芥分化之前，可能是与拟南芥甚至十字花科共有的一次全基因组复制事件；而土豆最近的那次全基因组复制事件发生在与拟南芥、葡萄分化之后，可能是土豆或说茄科所特有的。
上述应用例中，通过本发明的方法鉴定出来的两次全基因组复制事件非常清楚，比现有技术中只根据染色体数目进行的估计有足够的数据说服力，对土豆/茄科中确实有一次特有的全基因组复制事件给予了充分的回答。过去几年曾有不少研究根据EST数据对土豆、番茄等茄科植物的全基因组复制事件进行过推断，但都由于数据不全或方法不够有效而归于失败，或得到模凌两可的处于猜测阶段的结果。本发明实施例的方法得力于土全基因组数据及其注释基因的获得，更得力于高效的鉴定方法。本发明的方法可应用于其它所有的被子植物，已在木瓜、大豆、拟南芥、黄瓜等物种中做过测试，都取得显著的效果。
本发明的描述是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显然的。选择和描述实施例是为了更好说明本发明的原理和实际应用，并且使本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。
Claims (14)
1.一种全基因组复制事件的检测方法，其特征在于，包括：
获取与目标物种基因组上基因集的核酸序列对应的所述基因集的蛋白质序列以及基因在基因组上的位置信息；
对所述基因集的蛋白质序列进行蛋白比蛋白两两比对获得所述基因集的同源基因对；
根据所述基因集中基因在基因组上的位置信息以及所述基因集的同源基因对，利用动态规划算法获得所述基因组上具有共线性的同源区块；
获得所述同源区块中同源基因对的核酸序列的同源位点对位排列；
计算所述同源区块中同源基因对的遗传距离，获得年龄分布；
根据所述年龄分布确定全基因组复制事件。
2.根据权利要求1所述的检测方法，其特征在于，所述根据所述基因集中基因在基因组上的位置信息以及所述基因集的同源基因对利用动态规划获得所述基因组的具有共线性的同源区块的步骤包括：
根据所述基因集的同源基因对，形成所有基因间的同源关系矩阵；
利用基因集中每个基因在染色体上的位置信息，利用聚类算法，将相互连锁的共线性的同源基因块聚集在一起，形成具有共线性的同源区块。
3.根据权利要求1所述的检测方法，其特征在于，所述获取与目标物种基因组上基因集的核酸序列对应的所述基因集的蛋白质序列的步骤包括：
获取所述目标物种基因组上基因集的核酸序列；
根据密码子表将所述基因集的核酸序列转化为所述基因集的蛋白质序列。
4.根据权利要求3所述的检测方法，其特征在于，所述获取与目标物种基因组上基因集的核酸序列对应的所述基因集的蛋白质序列的步骤还包括：
当所述基因集中的基因有多个转录本时，选择最长的转录本对应的蛋白质序列作为所述基因的核酸序列对应蛋白质序列。
5.根据权利要求1所述的检测方法，其特征在于，所述将所述基因集的蛋白质序列进行蛋白比蛋白两两比对获得所述基因集的同源基因对的步骤包括：
通过blastp对所述基因集的蛋白质序列进行蛋白比蛋白两两比对，其中blastp蛋白质组比对的参数设置为p-value＜＝1e-7；
对比对输出的同源基因对通过将筛选的参数为设置为e-value＜＝1e-10进行筛选；
当同一基因对有多个比对结果仅以不同的e-value值给出时，选择e-value值最小的比对结果作为同源基因对。
6.根据权利要求1所述的检测方法，其特征在于，所述获得所述同源区块中同源基因对的核酸序列的同源位点对位排列的步骤包括：
对所述同源区块上的同源基因对的蛋白质序列进行全局序列比对获得所述同源基因对的蛋白质序列的同源位点对位排列；
根据所述同源基因对的蛋白质序列的同源位点对位排列构建所述基因集的核酸序列的同源位点对位排列。
7.根据权利要求6所述的检测方法，其特征在于，所述根据所述同源基因对的蛋白质序列的同源位点对位排列构建所述基因集的核酸序列的同源位点对位排列的步骤包括：
利用核酸序列的蛋白序列之间的对应密码子根据所述同源基因对的蛋白质序列的同源位点对位排列构建所述基因集的核酸序列的同源位点对位排列。
8.根据权利要求1所述的检测方法，其特征在于，所述计算所述同源区块中同源基因对的遗传距离建立年龄分布图的步骤包括：
每一个同源比对块上的每对同源基因对，计算遗传距离，对同一个同源区块上的所有基因对的遗传距离计算一个平均值；
将遗传距离作为其分化年龄的一种衡量，建立年龄分布图。
9.一种全基因组复制的检测系统，其特征在于，包括：
蛋白质序列获取模块，用于获取与目标物种基因组上基因集的核酸序列对应的所述基因集的蛋白质序列以及基因在基因组上的位置信息；
蛋白质序列比对模块，用于将所述基因集的蛋白质序列进行蛋白比蛋白两两比对，获得所述基因集的同源基因对；
同源区块确定模块，用于根据所述基因在基因组上的位置信息以及所述基因集的同源基因对利用动态规划获得所述基因组的具有共线性的同源区块；
对位排列确定模块，用于获得所述同源区块中同源基因对的核酸序列的同源位点对位排列；
年龄分布获取模块，用于计算所述同源区块中同源基因对的遗传距离，获得年龄分布；
全基因组复制确定模块，用于根据所述年龄分布确定全基因组复制事件。
10.根据权利要求9所述的全基因组复制的检测系统，其特征在于，所述蛋白质序列获取模块包括：
核酸序列获取单元，用于获得目标物种基因序列的基因集的核酸序列以及基因在基因组上的位置信息，发送所述基因集的核酸序列；
蛋白质序列转化单元，用于接收所述基因集的核酸序列，根据密码子表将所述基因集中各个基因的核酸序列转化为蛋白质序列，获得所述基因集的蛋白质序列。
11.根据权利要求9所述的全基因组复制的检测系统，其特征在于，所述蛋白质序列比对模块包括：
同源基因对获取单元，用于通过BLASTP蛋白质组比对对所述基因集的蛋白质序列进行蛋白比蛋白两两比对，获得同源基因对，其中所述blastp蛋白质组比对的参数设置为p-value＜＝1e-7；
同源基因对过滤单元，用于同源基因对获取单元获得的同源基因对进行筛选，筛选的参数为设置为e-value＜＝1e-10。
12.根据权利要求9所述的全基因组复制的检测系统，其特征在于，所述对位排列确定模块包括：
蛋白质对位排列获取单元，用于对所述同源区块上的同源基因对的蛋白质序列进行全局序列比对获得所述同源基因对的蛋白质序列的同源位点对位排列；
核酸对位排列获取单元，用于根据所述同源基因对的蛋白质序列的同源位点对位排列构建所述基因集的核酸序列的同源位点对位排列。
13.根据权利要求9所述的全基因组复制的检测系统，其特征在于，所述同源区块确定模块根据所述基因集的同源基因对，形成所有基因间的同源关系矩阵；利用基因集中每个基因在染色体上的位置信息，利用聚类算法，将相互连锁的共线性的同源基因块聚集在一起，形成具有共线性的同源区块。
14.根据权利要求9所述的全基因组复制的检测系统，其特征在于，所述年龄分布获取模块对于同源区块上的每对同源基因对，计算其遗传距离，对同一个同源区块上的所有同源基因对的遗传距离进行平均，以其平均值作为所述同源区块的遗传距离。
Priority Applications (2)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2010102849664A CN101957892B (zh) | 2010-09-17 | 2010-09-17 | 一种全基因组复制事件的检测方法和系统 |
| HK11102472.0A HK1148373A1 (en) | 2010-09-17 | 2011-03-11 | A method and a system for detecting genome-wide duplication |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN2010102849664A CN101957892B (zh) | 2010-09-17 | 2010-09-17 | 一种全基因组复制事件的检测方法和系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN101957892A CN101957892A (zh) | 2011-01-26 |
| CN101957892B true CN101957892B (zh) | 2012-09-05 |
Family
ID=43485218
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN2010102849664A Active CN101957892B (zh) | 2010-09-17 | 2010-09-17 | 一种全基因组复制事件的检测方法和系统 |
Country Status (2)
| Country | Link |
|---|---|
| CN (1) | CN101957892B (zh) |
| HK (1) | HK1148373A1 (zh) |
Families Citing this family (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104134016B (zh) * | 2014-07-30 | 2017-12-15 | 北京诺禾致源科技股份有限公司 | 分子水平上的系谱重建的装置和方法 |
| CN104134018B (zh) * | 2014-07-30 | 2017-09-26 | 北京诺禾致源科技股份有限公司 | 系谱中染色体区段的来源推断的装置和方法 |
| CN104298892B (zh) * | 2014-09-18 | 2017-05-10 | 天津诺禾致源生物信息科技有限公司 | 基因融合的检测装置和方法 |
| CN106951729A (zh) * | 2017-03-19 | 2017-07-14 | 中国海洋大学 | 一种利用细胞器基因组的同源模块进行系统发生分析的方法 |
| CN107330303B (zh) * | 2017-06-12 | 2020-06-30 | 浙江工业大学 | 一种多域蛋白模板无缝比对方法 |
| CN107832583B (zh) * | 2017-11-08 | 2021-04-16 | 武汉大学 | 一种基于图匹配的跨物种生物通路发现方法 |
| CN108733977A (zh) * | 2018-05-31 | 2018-11-02 | 中国人民解放军军事科学院军事医学研究院 | 真核生物保守转录因子结合位点聚集区tfcr的识别方法与应用 |
| CN108920901B (zh) * | 2018-07-24 | 2019-10-01 | 中国医学科学院北京协和医院 | 一种测序数据突变分析系统 |
| CN109949868B (zh) * | 2019-03-01 | 2020-10-16 | 深圳乐土生物科技有限公司 | 基于耐受性分析的基因等级排序方法和装置 |
| CN111508561B (zh) * | 2019-07-04 | 2024-02-06 | 北京希望组生物科技有限公司 | 同源序列和同源序列中串联重复序列的检测方法、计算机可读介质和应用 |
| CN111445953B (zh) * | 2020-03-27 | 2022-04-26 | 武汉古奥基因科技有限公司 | 一种利用全基因组比对拆分四倍体鱼类亚基因组的方法 |
| CN111564180A (zh) * | 2020-05-12 | 2020-08-21 | 西藏自治区农牧科学院水产科学研究所 | 一种鮡科鱼类古染色体进化比较分析的方法 |
| CN112270954B (zh) * | 2020-11-17 | 2023-09-26 | 中山大学 | 一种大规模连锁不平衡遗传位点快速分区方法 |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101089184A (zh) * | 2006-06-12 | 2007-12-19 | 国家海洋局第三海洋研究所 | 高温双链dna噬菌体gbsv1的全基因组序列 |
-
2010
- 2010-09-17 CN CN2010102849664A patent/CN101957892B/zh active Active
-
2011
- 2011-03-11 HK HK11102472.0A patent/HK1148373A1/xx not_active Application Discontinuation
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101089184A (zh) * | 2006-06-12 | 2007-12-19 | 国家海洋局第三海洋研究所 | 高温双链dna噬菌体gbsv1的全基因组序列 |
Also Published As
| Publication number | Publication date |
|---|---|
| HK1148373A1 (en) | 2011-09-02 |
| CN101957892A (zh) | 2011-01-26 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN101957892B (zh) | 一种全基因组复制事件的检测方法和系统 | |
| CN106202991B (zh) | 一种基因组多重扩增测序产物中突变信息的检测方法 | |
| Parks et al. | Phylogenomics reveals an extensive history of genome duplication in diatoms (Bacillariophyta) | |
| CN104164479B (zh) | 杂合基因组处理方法 | |
| CN101930502B (zh) | 表型基因的检测及生物信息分析的方法及系统 | |
| CN104302781B (zh) | 一种检测染色体结构异常的方法及装置 | |
| CN104762402B (zh) | 超快速检测人类基因组单碱基突变和微插入缺失的方法 | |
| CN110648721B (zh) | 针对外显子捕获技术检测拷贝数变异的方法及装置 | |
| CN103080333A (zh) | 一种基因组结构性变异检测方法和系统 | |
| CN108256292A (zh) | 一种拷贝数变异检测装置 | |
| CN106480221B (zh) | 基于基因拷贝数变异位点对林木群体基因型分型的方法 | |
| CN110111843A (zh) | 对核酸序列进行聚类的方法、设备及存储介质 | |
| CN110055306A (zh) | 一种基于转录组测序挖掘玉米耐低氮基因的方法 | |
| CN108004302A (zh) | 一种转录组参考的关联分析方法及其应用 | |
| CN111161797B (zh) | 一种基于三代测序检测多样本量比较转录组分析方法 | |
| Linde et al. | Rates and patterns of molecular evolution in bryophyte genomes, with focus on complex thalloid liverworts, Marchantiopsida | |
| Chira et al. | A cluster merging method for time series microarray with production values | |
| Lv et al. | Diverse phylogenomic datasets uncover a concordant scenario of laurasiatherian interordinal relationships | |
| CN110956155B (zh) | 基于co数据的综采工作面作业工序模糊聚类识别方法 | |
| CA3140963A1 (en) | Deep learning based system and method for prediction of alternative polyadenylation site | |
| CN111192636A (zh) | 一种适用于oligodT富集的mRNA二代测序结果分析方法 | |
| CN103184275A (zh) | 一种水稻基因组基因标识的新方法 | |
| Izuno et al. | Demography and selection analysis of the incipient adaptive radiation of a Hawaiian woody species | |
| CN112086128B (zh) | 一种适用于Sequel测序的三代全长转录组测序结果分析方法 | |
| Yu et al. | Genomic insights into biased allele loss and increased gene numbers after genome duplication in autotetraploid Cyclocarya paliurus |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| REG | Reference to a national code |
Ref country code: HK Ref legal event code: DE Ref document number: 1148373 Country of ref document: HK |
|
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| REG | Reference to a national code |
Ref country code: HK Ref legal event code: GR Ref document number: 1148373 Country of ref document: HK |
|
| EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20110126 Assignee: BGI TECH SOLUTIONS Co.,Ltd. Assignor: BGI SHENZHEN Co.,Ltd. Contract record no.: 2012440020389 Denomination of invention: Whole-genome replication event detection method and system Granted publication date: 20120905 License type: Exclusive License Record date: 20121219 |
|
| LICC | Enforcement, change and cancellation of record of contracts on the licence for exploitation of a patent or utility model | ||
| ASS | Succession or assignment of patent right |
Owner name: BGI TECHNOLOGY SOLUTIONS CO., LTD. Free format text: FORMER OWNER: BGI-SHENZHEN CO., LTD. Effective date: 20130422 |
|
| C41 | Transfer of patent application or patent right or utility model | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20130422 Address after: 518083 science and Technology Pioneer Park, comprehensive building, Beishan Industrial Zone, Yantian District, Guangdong, Shenzhen 201 Patentee after: BGI TECH SOLUTIONS Co.,Ltd. Address before: Beishan Industrial Zone Building in Yantian District of Shenzhen city of Guangdong Province in 518083 Patentee before: BGI SHENZHEN Co.,Ltd. |
|
| TR01 | Transfer of patent right |
Effective date of registration: 20221214 Address after: No. 128, Hengtong Road, huankeyuan, Yixing, Wuxi, Jiangsu, 214205 Patentee after: Huada Qinglan Biotechnology (Wuxi) Co.,Ltd. Patentee after: BGI TECH SOLUTIONS Co.,Ltd. Address before: 518083 science and Technology Pioneer Park 201, Beishan Industrial Park, Yantian District, Shenzhen City, Guangdong Province Patentee before: BGI TECH SOLUTIONS Co.,Ltd. |
|
| TR01 | Transfer of patent right |