CN114882949A - 基于氨基酸序列比对的蛋白家族系统发育分析方法 - Google Patents

基于氨基酸序列比对的蛋白家族系统发育分析方法 Download PDF

Info

Publication number
CN114882949A
CN114882949A CN202210749152.6A CN202210749152A CN114882949A CN 114882949 A CN114882949 A CN 114882949A CN 202210749152 A CN202210749152 A CN 202210749152A CN 114882949 A CN114882949 A CN 114882949A
Authority
CN
China
Prior art keywords
sequence
amino acid
analysis
matrix
site
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.)
Pending
Application number
CN202210749152.6A
Other languages
English (en)
Inventor
郑波
张哲�
施雪萍
朱苗苗
谢琪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Huazhong Agricultural University
Original Assignee
Huazhong Agricultural University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Huazhong Agricultural University filed Critical Huazhong Agricultural University
Priority to CN202210749152.6A priority Critical patent/CN114882949A/zh
Publication of CN114882949A publication Critical patent/CN114882949A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Health & Medical Sciences (AREA)
  • Evolutionary Biology (AREA)
  • Biotechnology (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Medical Informatics (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biophysics (AREA)
  • General Health & Medical Sciences (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • General Physics & Mathematics (AREA)
  • Genetics & Genomics (AREA)
  • Molecular Biology (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于氨基酸序列比对的蛋白家族系统发育分析方法，包括以下步骤：基于氨基酸序列比对融合方法，获得合并的多序列比对结果；对合并的多序列比对结果进行数字化，构建分数矩阵；对分数矩阵进行降维、聚类处理，获得输入序列；鉴定输入序列的特异位点及保守位点；对输入序列进行拟时间分析，获得输入序列的轨迹排序；基于轨迹排序，获得输入序列的发育轨迹。本发明通过序列位点特征进行分数矩阵构建及其降维分析，从而推断基因家族间的聚类和进化关系，在保证序列聚类稳定性的情况下，有效提高了序列基因聚类速度，为基因系统发育分析、发育轨迹分析提供了新的工具和方法。

Description

基于氨基酸序列比对的蛋白家族系统发育分析方法
技术领域
本发明属于聚类分析及生物学技术领域，特别是涉及基于氨基酸序列比对的蛋白家族系统发育分析方法。
背景技术
依据多序列比对融合方法对一组同源蛋白序列进行系统发育分析，推断这些同源蛋白序列之间的进化关系，是蛋白功能分析的第一步。在获得同源蛋白序列的多序列比对结果之后，通常有两类方法对其进行系统发育分析，分别是基于序列位点特征的系统发育分析方法，包括极大似然法，最大简约法和贝叶斯推断法；以及基于序列之间的进化距离的系统发育分析方法，包括邻接法，最小进化法和非加权组平均法。非加权组平均法假设碱基或氨基酸的替换速率是均等且恒等的，不利于模拟氨基酸序列之间的进化关系。最小进化法在邻接法的基础上搜索分支长度最小的树，本质与邻接法相同。相比于基于序列之间的进化距离的分析方法，极大似然法、贝叶斯推断法和最大简约法能保留更多的序列信息，因此能够得到更为准确的结果，但由于极大似然法、贝叶斯推断法计算量过大，而最大简约法只适用于近缘序列，其适用性不如邻接法广。
发明内容
为了解决上述问题，我们开发了一种基于氨基酸序列比对的蛋白家族系统发育分析方法，通过氨基酸序列比对的融合和数字化对同源蛋白序列进行系统发育分析。本发明通过序列位点特征进行分数矩阵构建及其降维分析，从而推断基因家族间的聚类和进化关系。
为实现上述目的，本发明提供了基于氨基酸序列比对的蛋白家族系统发育分析方法，包括以下步骤：
基于氨基酸序列比对融合方法，获得合并的多序列比对结果；
对所述合并的多序列比对结果进行数字化，构建分数矩阵；
对所述分数矩阵进行降维、聚类处理，获得输入序列；
鉴定所述输入序列的特异位点、保守位点；
对所述输入序列进行拟时间分析，获得所述输入序列的轨迹排序；
基于所述轨迹排序，获得所述输入序列的发育轨迹。
优选地，基于氨基酸序列比对融合方法，获得合并的多序列比对结果的过程包括：
过滤掉待合并的第一多序列比对结果、第二多序列比对结果中的噪音位点，获得第一剩余位点、第二剩余位点；
获取所述第一剩余位点、第二剩余位点的氨基酸组成；
获取所述第一多序列比对结果的所有位点、第二多序列比对结果的所有位点的匹配分数；
基于所述匹配分数，通过动态规划算法合并所述第一多序列比对结果、第二多序列比对结果，并计算所述第一多序列比对结果、第二多序列比对结果的合并质量分数矩阵；
对所述合并质量分数矩阵进行层次聚类，确定待合并多序列比对结果的合并顺序，对所述待合并多序列比对结果依序合并，获得合并的多序列比对结果。
优选地，获取所述第一多序列比对结果的所有位点、第二多序列比对结果的所有位点的匹配分数的过程包括：
通过比较所述第一剩余位点、第二剩余位点的氨基酸组成，获得位点匹配比例、位点不匹配比例、氨基酸匹配比例、氨基酸不匹配比例；基于所述位点匹配比例、位点不匹配比例、氨基酸匹配比例、氨基酸不匹配比例，通过整合氨基酸替代矩阵，获得所述匹配分数。
优选地，对所述合并的多序列比对结果进行数字化，构建分数矩阵的过程包括：
获取所述合并的多序列比对结果的位点氨基酸组成，构建氨基酸频率矩阵；基于所述氨基酸频率矩阵构建所述分数矩阵。
优选地，对所述分数矩阵进行降维、聚类处理的过程包括：
基于R语言单细胞测序分析软件包的主成分分析函数和降维函数，对所述分数矩阵进行主成分分析和统一流形逼近和投影分析，获得输入序列；
基于共享最近邻算法，对所述输入序列进行聚类；
基于R语言层次聚类分析软件包的随机引用函数计算调整兰德系数，评估聚类效果。
优选地，鉴定所述输入序列的特异位点及保守位点的过程包括：
基于威尔科克森符号秩检验方法鉴定所述输入序列的特异位点；
基于氨基酸种类、位点信息熵、校正数、缺口比例，计算所述输入序列的每个位点的保守分数。
优选地，对所述输入序列进行拟时间分析的过程包括：
基于R语言单细胞拟时间分析软件包中的降维函数对所述分数矩阵进行树形判别降维分析和聚类；
基于细胞排序函数对所述输入序列进行轨迹推断，并沿着轨迹进行排序；
基于R语言层次聚类分析软件包随机引用函数计算调整兰德系数，评估聚类效果。
优选地，获得所述输入序列的发育轨迹的过程包括：
基于矩阵实验室软件单细胞能量路径软件包中的能量预测函数和家族推断函数，计算获得所述输入序列的能量景观及转移概率，并构建最小有向生成树；
基于每一类基因的物种出现顺序及家族信息对所述输入序列进行排序；
如果所述转移概率大于0.2，则使用双向箭头连接对应基因簇，获得发育路径；
如果所述发育路径出现在所述最小有向生成树中，但所述转移概率小于0.2，则使用单向虚线箭头连接对应基因簇，最终获得所述输入序列的发育轨迹。
本发明的技术效果为：
本发明通过序列位点特征进行分数矩阵构建及其降维分析，从而推断基因家族间的聚类和进化关系，在保证序列聚类稳定性的情况下，有效提高了序列基因聚类速度，为基因系统发育分析、发育轨迹分析提供了新的工具和方法。
本发明从蛋白氨基酸序列入手，研究同源蛋白之间的进化关系，进而分别从物种和家族进化角度探究蛋白氨基酸序列与其功能之间的关系，为后续的蛋白质功能鉴定提供理论依据，是一种兼具速度快、适用范围广、结果准确、对用户友好等优势的氨基酸序列系统发育分析方法，
附图说明
构成本申请的一部分的附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：
图1为本发明实施例中的方法流程图；
图2为本发明实施例中的多序列比对结果合并模块工作流程图；
图3为本发明实施例中的应用于植物HB家族的聚类效果图；
图4为本发明实施例中的推断植物HB超家族中每个基因簇的发育轨迹图；
图5为本发明实施例中的推断植物HB超家族中的WOX家族每个基因簇的发育轨迹图；
图6为本发明实施例中的构建植物HB超家族的聚类树图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
如图1-2所示，本发明提供了一种基于氨基酸序列比对的蛋白家族系统发育分析方法，包括如下步骤：
1)合并多个家族的多序列比对结果MSA；
2)对于已合并的蛋白家族多序列比对结果MSA进行数字化，构建分数矩阵；
3)对分数矩阵进行降维和聚类；
4)鉴定各家族的保守位点和特异位点；
5)对输入序列进行拟时间分析；
6)推断各家族输入序列的发育轨迹。
进一步地，在步骤1)中，使用氨基酸序列比对融合方法M-MSA合并多个家族的多序列比对结果MSA，步骤如下：
(1)设待合并的两个多序列比对结果MSA分别为X和Y，首先过滤掉两个多序列比对结果MSA中的噪音位点。本研究中，噪音位点即多序列比对结果MSA中缺口比例超过90％的位点。
(2)计算X和Y中每个剩余位点的氨基酸组成(图3，b)。
(3)计算X中所有位点与Y中所有位点之间的匹配分数(图3)。以X中的位点x1和Y中的位点y1为例。首先比较x1和y1的氨基酸组成，计算x1和y1的匹配比例pmatch(x1,y1)和不匹配比例pmismatch(x1,y1)(图3)。接着，计算x1和y1中每种氨基酸的匹配比例pmx1,y1和不匹配比例pnx1,y1。然后整合氨基酸替代矩阵M，计算x1和y1的匹配分数，计算公式如下：
Figure BDA0003720589020000071
其中，Mi,i代表氨基酸i在氨基酸替代矩阵M中的替代分数，Mi,j代表氨基酸i和氨基酸j在氨基酸替代矩阵M中的替代分数，pmx1,y1(i)代表氨基酸i在x1和y1中的匹配比例，pnx1,y1(i)和pnx1,y1(j)分别代表氨基酸i和氨基酸j在x1和y1中的不匹配比例，ss(x1,y1)即位点x1和位点y1的匹配分数。
(4)使用动态规划算法Needleman–Wunsch，根据X和Y中的所有位点分数合并X和Y，同时计算X和Y之间的合并质量分数Q(X,Y)。
(5)计算所有待合并多序列比对结果MSA两两之间合并质量分数矩阵Q，并对其进行层次聚类，确定所有待合并多序列比对结果MSA的合并顺序，依序合并多个多序列比对结果MSA。
进一步地，在步骤2)中使用多序列比对结果MSA构建分数矩阵，步骤如下：
(1)计算每个位点的氨基酸组成C，构建氨基酸频率矩阵Q：
Figure BDA0003720589020000072
其中，Cai(n)代表基本氨基酸ai在位点n出现的频率；
(2)根据对应位点的氨基酸频率矩阵构建分数矩阵S：
Figure BDA0003720589020000073
其中，am，n代表序列m在位点n的对应氨基酸残基，
Figure BDA0003720589020000081
即序列m在位点n的对应氨基酸残基在位点n出现的频率，S即使用MSA构建的分数矩阵。
进一步地，在步骤3)中，对分数矩阵S进行降维及聚类分析，步骤如下：
(1)调用R语言单细胞测序分析软件包Seurat的主成分分析函数RunPCA和降维函数RunUMAP分别对分数矩阵进行主成分分析(Principal Component Analysis，PCA)和统一流形逼近和投影(Uniform Manifold Approximation and Projection，UMAP)分析，将输入序列投影到二维平面上。
(2)使用共享最近邻算法(Shared Nearest Neighbor，SNN)对输入序列进行聚类；
(3)使用R语言层次聚类分析软件包flexclust的随机引用函数randIndex计算调整兰德系数(Adjusted Rand Index，ARI)，评估聚类效果。
进一步地，在步骤4)中，鉴定每一类输入序列的特异位点及保守位点，步骤如下：
(1)使用威尔科克森符号秩检验方法检验鉴定每一类输入序列的特异位点；
(2)计算每一类HB基因超家族氨基酸序列每个位点的保守分数，计算公式如下：
保守分数＝Bits值×(1-gap比例) (4)
Bits值计算公式如下：
Bits＝Iog2(m)-H+em (5)
其中m为氨基酸种类，H为该位点信息熵，em为校正数；
进一步地，在步骤5)中，对输入氨基酸序列进行拟时间分析，步骤如下：
(1)调用R语言单细胞拟时间分析软件包monocle中的降维函数reduceDimension对分数矩阵进行树形判别降维(Discriminative Dimensionality Reduction withTrees，DDRTree)分析和聚类；
(2)使用细胞排序函数orderCells对输入氨基酸序列进行轨迹推断，并沿着轨迹进行排序。使用R语言聚类分析软件包flexclust的随机引用函数randIndex计算调整兰德系数(Adjusted Rand Index,ARI)，评估聚类效果。
进一步地，在步骤6)中，对输入氨基酸序列进行发育轨迹分析，步骤如下：
(1)分别使用矩阵实验室软件MATLAB单细胞能量路径软件包scEpath中的能量预测函数estimatingscEnergy和家族推断函数inferingLineage计算每一类集群(Cluster)输入序列的能量景观及转移概率(Transition Probabilities，TP)，并构建最小有向生成树(Minimum Directed Spanning Tree，MDST)；
(2)依据每一类基因的物种出现顺序及家族信息对其进行排序；
(3)如果两类基因之间的转移概率(TP)大于0.2，则使用双向蓝色箭头连接这两个基因簇；
(4)如果两类基因之间的转移概率(TP)大于0.2，则使用双向蓝色箭头连接这两个基因簇；
(5)如果连接两类基因之间的路径出现在最小有向生成树MDST中，但两类基因之间的转移概率TP小于0.2，则使用单向蓝色虚线箭头连接这两类基因，最终获得输入氨基酸序列的可能发育轨迹。
实施例一
如图1-2所示，本实施例中提供一种基于氨基酸序列比对的蛋白家族系统发育分析方法，包括：
1.植物HB蛋白家族氨基酸序列提取：我们在PlantTFDB数据库中调取植物HB家族的全长蛋白(共18147条)和同源异形盒结构域(Homeodomain，HD)(共15184条)氨基酸序列，以及其家族注释信息。
2.获取植物HB蛋白超家族基因物种注释信息：针对物种发育程度，依照物种进化程度的从高到低将其依次分为被子植物(Angiospermae)，针叶植物(Coniferophyta)，石松植物(Lycopodiophyta)，苔藓植物(Bryophyta)，苔类植物(Marchantiophyta)，轮藻(Charophyta)和绿藻(Chlorophytae)。针对双子叶植物，将物种信息分为紫菀科(Asterids)、基底木兰科(Basal Magnoliophyta)、豆科(Fabids)、锦葵科(Malvids)、其他真双子叶植物(Other Eudicots)、单子叶植物(Monocots)和其他植物(Other plants)。
3.获取植物HB蛋白超家族基因亚家族注释信息：通过查找文献，获得拟南芥HD-ZIP家族各亚家族(I、II、III和IV)、TALE家族各亚家族(BEL、KNOX_I和KNOX_II)和WOX家族各亚家族(Ancient、Intermediate和WUS)的注释信息，进一步使用局部序列比对基本检索工具(Basic Local Alignment Search Tool，BLAST)(BLAST-v2.7.1)找出每个HD-ZIP家族、TALE家族和WOX家族基因最同源的2个拟南芥基因(max_target_seqs＝2，其余参数为默认参数)。如果最同源的2个拟南芥基因属于同一个亚家族，将该基因注释为对应家族，如果最同源的2个拟南芥基因分别属于不同亚家族，将该基因注释为最同源基因对应的亚家族，同时对其进行标注。
4.获取植物Homeobox蛋白超家族及其各基因亚家族的多序列比对结果MSA：使用MAFFT软件渐进式方法FFT-NS-1默认参数对HB超家族全长蛋白序列及HD蛋白序列进行序列比对，同时对HB超家族中各个家族全长蛋白序列及HD蛋白序列进行序列比对。
5.合并HB蛋白超家族各基因亚家族的多序列比对结果：
使用氨基酸序列比对融合方法M-MSA合并HB蛋白超家族各亚家族的多序列比对结果，步骤如下：
(1)设待合并的两个多序列比对结果分别为X和Y，首先过滤掉两个多序列比对结果中的噪音位点。本研究中，噪音位点即多序列比对结果MSA中缺口比例gap超过90％的位点。
(2)计算X和Y中每个剩余位点的氨基酸组成。
(3)计算X中所有位点与Y中所有位点之间的匹配分数(图3)。以X中的位点x1和Y中的位点y1为例。首先比较x1和y1的氨基酸组成，计算x1和y1的匹配比例pmatch(x1,y1)和不匹配比例pmismatch(x1,y1)(图3)。接着，计算x1和y1中每种氨基酸的匹配比例pmx1,y1和不匹配比例pnx1,y1。然后整合氨基酸替代矩阵M，计算x1和y1的匹配分数，计算公式如下：
Figure BDA0003720589020000121
其中，Mi,i代表氨基酸i在氨基酸替代矩阵M中的替代分数，Mi,j代表氨基酸i和氨基酸j在氨基酸替代矩阵M中的替代分数，pmx1,y1(i)代表氨基酸i在x1和y1中的匹配比例，pnx1,y1(i)和pnx1,y1(j)分别代表氨基酸i和氨基酸j在x1和y1中的不匹配比例，ss(x1,y1)即位点x1和位点y1的匹配分数。
(4)使用动态规划算法Needleman–Wunsch，根据X和Y中的所有位点分数合并X和Y，同时计算X和Y之间的合并质量分数Q(X,Y)。
(5)计算所有待合并多序列比对结果两两之间合并质量分数矩阵Q，并对其进行层次聚类，确定所有待合并多序列比对结果的合并顺序，依序合并多个多序列比对结果。
6.HB基因超家族分数矩阵的构建：
对合并的氨基酸序列进行数字化，构建分数矩阵，步骤如下：
(1)计算每个位点的氨基酸组成C，构建氨基酸频率矩阵Q：
Figure BDA0003720589020000122
其中，Cai(n)代表基本氨基酸ai在位点n出现的频率；
(2)根据对应位点的氨基酸频率矩阵构建分数矩阵S：
Figure BDA0003720589020000131
其中，am，n代表序列m在位点n的对应氨基酸残基，
Figure BDA0003720589020000132
即序列m在位点n的对应氨基酸残基在位点n出现的频率，S即使用氨基酸序列构建的分数矩阵。
7.HB基因超家族分数矩阵降维及聚类分析：
对HB基因超家族的分数矩阵进行降维及聚类分析，步骤如下：
(1)调用R语言单细胞测序分析软件包Seurat的主成分分析函数RunPCA和降维函数RunUMAP分别对分数矩阵进行主成分分析(Principal Component Analysis，PCA)和统一流形逼近和投影(Uniform Manifold Approximation and Projection，UMAP)分析，将HB基因超家族序列投影到二维平面上。
(2)使用共享最近邻算法(Shared Nearest Neighbor，SNN)对HB基因超家族序列进行聚类；
(3)使用R语言层次聚类分析软件包flexclust的随机引用函数randIndex计算调整兰德系数(Adjusted Rand Index，ARI)，评估聚类效果。
8.HB基因超家族每一类氨基酸序列特异位点及保守位点鉴定
鉴定HB基因超家族每一类氨基酸序列的特异位点及保守位点，步骤如下：
(1)使用威尔科克森符号秩检验方法检验鉴定每一类HB基因超家族氨基酸序列的特异位点；
(2)计算每一类HB基因超家族氨基酸序列每个位点的保守分数，计算公式如下：
保守分数＝Bits值×(1-gap比例) (4)
Bits值计算公式如下：
Bits＝log2(m)-H+em (5)
其中m为氨基酸种类，H为该位点信息熵，em为校正数。
9.HB基因超家族的拟时间分析：
对HB基因超家族氨基酸序列进行拟时间分析，步骤如下：
(1)调用R语言单细胞拟时间分析软件包monocle中的降维函数reduceDimension对分数矩阵进行树形判别降维(Discriminative Dimensionality Reduction withTrees，DDRTree)分析和聚类；
(2)使用细胞排序函数orderCells对HB基因超家族氨基酸序列进行轨迹推断，并沿着轨迹进行排序。使用R语言层次聚类分析软件包flexclust的随机引用函数randIndex计算调整兰德系数(Adjusted Rand Index,ARI)，评估聚类效果。
10.分别对HB基因超家族氨基酸序列和HB基因超家族WOX家族氨基酸序列进行发育轨迹分析，步骤如下：
(1)分别使用矩阵实验室软件MATLAB单细胞能量路径软件包scEpath中的能量预测函数estimatingscEnergy和家族推断函数inferingLineage计算每一类集群(Cluster)HB基因超家族氨基酸序列和HB基因超家族WOX家族氨基酸序列的能量景观及转移概率(Transition Probabilities，TP)，并构建最小有向生成树(Minimum Directed SpanningTree，MDST)(如图6所示)；
(2)依据HB基因超家族氨基酸序列和HB基因超家族WOX家族氨基酸序列中的每一类基因的物种出现顺序及家族信息对其进行排序；
(3)如果两类基因之间的转移概率TP大于0.2，则使用双向蓝色箭头连接这两个基因簇；
(4)如果连接两类基因之间的路径出现在最小有向生成树MDST中，但两类基因之间的转移概率TP小于0.2，则使用单向蓝色虚线箭头连接这两类基因，最终获得输入氨基酸序列的可能发育轨迹(如图4-图5所示)。
本发明公开了一种基于氨基酸序列比对的蛋白家族系统发育分析方法。包括如下步骤：(1)合并多个蛋白家族或亚家族的多序列比对结果(Multiple Sequence Alignment，MSA)；(2)使用合并后的MSA构建分数矩阵；(3)使用氨基酸序列比对融合方法UMAP和树形判别降维算法DDRTree对分数矩阵进行降维，之后使用共享最近邻算法SNN对输入序列进行聚类；(4)鉴定每一类输入序列保守位点及特异位点；(5)对分数矩阵进行拟时间分析；(6)推断各家族输入序列的发育轨迹。本发明有效提高了序列基因聚类速度，为基因系统发育分析、发育轨迹分析提供了新的工具和方法。
以上所述，仅为本申请较佳的具体实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应该以权利要求的保护范围为准。

Claims (8)

1.基于氨基酸序列比对的蛋白家族系统发育分析方法，其特征在于，包括以下步骤：
基于氨基酸序列比对融合方法，获得合并的多序列比对结果；
对所述合并的多序列比对结果进行数字化，构建分数矩阵；
对所述分数矩阵进行降维、聚类处理，获得输入序列；
鉴定所述输入序列的特异位点、保守位点；
对所述输入序列进行拟时间分析，获得所述输入序列的轨迹排序；
基于所述轨迹排序，获得所述输入序列的发育轨迹。
2.根据权利要求1所述的基于氨基酸序列比对的蛋白家族系统发育分析方法，其特征在于，基于氨基酸序列比对融合方法，获得合并的多序列比对结果的过程包括：
过滤掉待合并的第一多序列比对结果、第二多序列比对结果中的噪音位点，获得第一剩余位点、第二剩余位点；
获取所述第一剩余位点、第二剩余位点的氨基酸组成；
获取所述第一多序列比对结果的所有位点、第二多序列比对结果的所有位点的匹配分数；
基于所述匹配分数，通过动态规划算法合并所述第一多序列比对结果、第二多序列比对结果，并计算所述第一多序列比对结果、第二多序列比对结果的合并质量分数矩阵；
对所述合并质量分数矩阵进行层次聚类，确定待合并多序列比对结果的合并顺序，对所述待合并多序列比对结果依序合并，获得合并的多序列比对结果。
3.根据权利要求2所述的基于氨基酸序列比对的蛋白家族系统发育分析方法，其特征在于，获取所述第一多序列比对结果的所有位点、第二多序列比对结果的所有位点的匹配分数的过程包括：
通过比较所述第一剩余位点、第二剩余位点的氨基酸组成，获得位点匹配比例、位点不匹配比例、氨基酸匹配比例、氨基酸不匹配比例；基于所述位点匹配比例、位点不匹配比例、氨基酸匹配比例、氨基酸不匹配比例，通过整合氨基酸替代矩阵，获得所述匹配分数。
4.根据权利要求1所述的基于氨基酸序列比对的蛋白家族系统发育分析方法，其特征在于，对所述合并的多序列比对结果进行数字化，构建分数矩阵的过程包括：
获取所述合并的多序列比对结果的位点氨基酸组成，构建氨基酸频率矩阵；基于所述氨基酸频率矩阵构建所述分数矩阵。
5.根据权利要求1所述的基于氨基酸序列比对的蛋白家族系统发育分析方法，其特征在于，对所述分数矩阵进行降维、聚类处理的过程包括：
基于R语言单细胞测序分析软件包的主成分分析函数和降维函数，对所述分数矩阵进行主成分分析和统一流形逼近和投影分析，获得输入序列；
基于共享最近邻算法，对所述输入序列进行聚类；
基于R语言层次聚类分析软件包的随机引用函数计算调整兰德系数，评估聚类效果。
6.根据权利要求1所述的基于氨基酸序列比对的蛋白家族系统发育分析方法，其特征在于，鉴定所述输入序列的特异位点及保守位点的过程包括：
基于威尔科克森符号秩检验方法鉴定所述输入序列的特异位点；
基于氨基酸种类、位点信息熵、校正数、缺口比例，获取所述输入序列的每个位点的保守分数。
7.根据权利要求1所述的基于氨基酸序列比对的蛋白家族系统发育分析方法，其特征在于，对所述输入序列进行拟时间分析的过程包括：
基于R语言单细胞拟时间分析软件包中的降维函数对所述分数矩阵进行树形判别降维分析和聚类；
基于细胞排序函数对所述输入序列进行轨迹推断，并沿着轨迹进行排序；
基于R语言层次聚类分析软件包随机引用函数计算调整兰德系数，评估聚类效果。
8.根据权利要求1所述的基于氨基酸序列比对的蛋白家族系统发育分析方法，其特征在于，获得所述输入序列的发育轨迹的过程包括：
基于矩阵实验室软件单细胞能量路径软件包中的能量预测函数和家族推断函数，获取所述输入序列的能量景观及转移概率，并构建最小有向生成树；
基于每一类基因的物种出现顺序及家族信息对所述输入序列进行排序；
如果所述转移概率大于0.2，则使用双向箭头连接对应基因簇，获得发育路径；
如果所述发育路径出现在所述最小有向生成树中，但所述转移概率小于0.2，则使用单向虚线箭头连接对应基因簇，最终获得所述输入序列的发育轨迹。
CN202210749152.6A 2022-06-29 2022-06-29 基于氨基酸序列比对的蛋白家族系统发育分析方法 Pending CN114882949A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210749152.6A CN114882949A (zh) 2022-06-29 2022-06-29 基于氨基酸序列比对的蛋白家族系统发育分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210749152.6A CN114882949A (zh) 2022-06-29 2022-06-29 基于氨基酸序列比对的蛋白家族系统发育分析方法

Publications (1)

Publication Number Publication Date
CN114882949A true CN114882949A (zh) 2022-08-09

Family

ID=82682835

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210749152.6A Pending CN114882949A (zh) 2022-06-29 2022-06-29 基于氨基酸序列比对的蛋白家族系统发育分析方法

Country Status (1)

Country Link
CN (1) CN114882949A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115312122A (zh) * 2022-10-12 2022-11-08 之江实验室 一种CRISPR-Cas酶可突变位点推荐方法和装置
CN116206675A (zh) * 2022-09-05 2023-06-02 北京分子之心科技有限公司 用于预测蛋白质复合物结构的方法、设备、介质及程序产品

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116206675A (zh) * 2022-09-05 2023-06-02 北京分子之心科技有限公司 用于预测蛋白质复合物结构的方法、设备、介质及程序产品
CN116206675B (zh) * 2022-09-05 2023-09-15 北京分子之心科技有限公司 用于预测蛋白质复合物结构的方法、设备、介质及程序产品
CN115312122A (zh) * 2022-10-12 2022-11-08 之江实验室 一种CRISPR-Cas酶可突变位点推荐方法和装置
CN115312122B (zh) * 2022-10-12 2022-12-16 之江实验室 一种CRISPR-Cas酶可突变位点推荐方法和装置

Similar Documents

Publication Publication Date Title
CN114882949A (zh) 基于氨基酸序列比对的蛋白家族系统发育分析方法
Kim et al. Introducing EzAAI: a pipeline for high throughput calculations of prokaryotic average amino acid identity
Allio et al. Whole genome shotgun phylogenomics resolves the pattern and timing of swallowtail butterfly evolution
Zhang et al. Assembly of allele-aware, chromosomal-scale autopolyploid genomes based on Hi-C data
Kellis et al. Methods in comparative genomics: genome correspondence, gene identification and regulatory motif discovery
Patwardhan et al. Molecular markers in phylogenetic studies-a review
CN105740650B (zh) 一种快速准确鉴定高通量基因组数据污染源的方法
Rijzaani et al. The pangenome of banana highlights differences between genera and genomes
Li et al. Characterization and comparison of the mitochondrial genomes from two Lyophyllum fungal species and insights into phylogeny of Agaricomycetes
Wu et al. Evolution at the subgene level: domain rearrangements in the Drosophila phylogeny
Barley et al. Sun skink landscape genomics: assessing the roles of micro‐evolutionary processes in shaping genetic and phenotypic diversity across a heterogeneous and fragmented landscape
CA2696843A1 (en) Method, system and software arrangement for comparative analysis and phylogeny with whole-genome optical maps
Rivera-Colón et al. Population genomics analysis with RAD, reprised: Stacks 2
Garcia-Erill et al. Warthog genomes resolve an evolutionary conundrum and reveal introgression of disease resistance genes
Feng et al. Evolution and multiple origins of zona pellucida genes in vertebrates
Meredith et al. Relationships and divergence times among the orders and families of Marsupialia
Pollo et al. Benchmarking hybrid assemblies of Giardia and prediction of widespread intra-isolate structural variation
Qiu et al. Defining honeybee subspecies in an evolutionary context warrants strategized conservation
Xia et al. Lycophyte transcriptomes reveal two whole-genome duplications in Lycopodiaceae: Insights into the polyploidization of Phlegmariurus
Balboa et al. African bushpigs exhibit porous species boundaries and appeared in Madagascar concurrently with human arrival
Liu et al. Extensive mitogenomic heteroplasmy and its implications in the phylogeny of the fish genus Megalobrama
Earley et al. Genomic regions associate with major axes of variation driven by gas exchange and leaf construction traits in cultivated sunflower (Helianthus annuus L.)
Jiang et al. Interspecific delimitation and relationships among four Ostrya species based on plastomes
Lendemer et al. Bacidia depriestiana (Ramalinaceae), a new species from the southern Appalachian Mountains of eastern North America
Huang et al. The revised complete mitogenome sequence of the tree frog Polypedatesmegacephalus (Anura, Rhacophoridae) by next-generation sequencing and phylogenetic analysis

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination