本发明属于生物信息学领域，具体涉及一种基于蛋白质质谱数据注释真核生物基因组的生物信息学方法。
背景技术：
基因组测序只能测出整个dna的碱基对排列顺序，不能直接测出dna上的基因及其功能，必须通过生物信息学方法，结合蛋白组学、转录组学，对测出来的序列进行分析，将基因及其功能加以挖掘、注释，这称作基因注释。
基因组注释(genomeannotation)是利用生物信息学方法和工具，对基因组所有基因的生物学功能进行高通量注释，是当前功能基因组学研究的一个热点。
基因结构的预测对于发现新基因、了解基因组结构规律具有重要意义，是各类基因组计划的重要内容。目前，无论是对于原核生物还是真核生物，对基因组的预测注释方法主要有三种：第一种是利用已有的转录数据来识别基因组中基因编码区的转录作图法；第二种方法是利用计算算法从统计学的角度来预测基因和调节元件的基因结构统计建模法；第三种是比较其它物种基因组序列中的保守区域进行基因预测的同源建模法。
1、转录作图法(transcriptionmapping)
转录作图法是利用已有的转录数据来识别基因组中的基因编码区。将反转录的cdna进行测序分析，以mrna转录本作为基因结构注释的依据，被认为是基因序列结构注释的“金标准”。但该方法存在一些不足：如由于mrna具有不稳定性，部分基因很难获得全长mrna序列；某些物种中转录本编码的是操纵子而不是基因；由于翻译框的不确定性以及翻译起始位点易受多种因素的影响，依靠mrna来决定起始密码子和终止密码子比较困难。也是由于上述原因，转录信息在大多数基于基因序列特征的统计建模的基因预测软件中仅作为补充。
2、基因结构统计建模法
基因结构统计建模法是利用计算算法从统计学的角度来预测基因和调节元件，是以从头计算(abinitio)为基础的软件直接从序列本身提取有关基因的信息，基于广义隐马尔可夫模型(generalizedhiddenmarkovmodel,ghmm)的数据结构，并通过使用机器学习所构建的方法，例如glimmer、gfenesh、genemarks、mgene、evigan和genemark等方法。
3、同源建模法
同源建模法是基于基因组中的功能元件在进化过程中保守的特点，通过比较其它物种基因组序列中的保守区域进行基因预测，是整合了数据库中已经存在的外部信息来进行预测注释，例如augustus、jigsaw、eugene、geneid和critica等方法。
现今注释流水线(annotationpipeline)已经被普遍采用，例如ensembl的流水线，ncbi的gnomon和ucsc的knowngenes。但是当同时使用多个基因组注释软件时就会出现注释的差异性。排除基因组测序质量因素，此类差异首先与给定的注释的证据的支持强度有关，如某一段序列的长度和上下游信息不充足就会使得小基因、短外显子和读码框移码的预测变得很困难；某些特殊密码子的使用也可能造成基因预测的错误；非标准的起始和终止密码子可以用来标定orfs，然而在某些情况下，标准终止密码在真核生物中编码硒代半胱氨酸或古细菌中编码吡咯赖氨酸。某些时候，即使预测是正确的，有些类型的基因仅仅编码rna而不会翻译，其最终产物可能是非编码的rna(ncrna)或者可以部分激活的假基因，其rnas没有编码蛋白质的功能并会很快被降解。最后，由于不同的物种具有不同的基因组，导致使用从头计算工具需要为每一个新的基因组进行机器学习训练。
而传统的基因组注释方法(例如生物信息学预测和同源信息比对等)在鉴定精度和深度上受到了较大的限制。相对于基因组学和转录组学，基于质谱技术的蛋白质组学直接研究编码基因翻译出的蛋白质产物，比转录组学注释基因组获得的结果更直接，而且可以发现由于知识不足导致的基因从头预测算法遗漏的基因和基因结构注释的错误。因此，为了完善基因组注释，使用蛋白质组学数据注释基因组，即蛋白基因组学(proteogenomics)，是继基因组学和转录组学之后发展起来的又一项新兴高通量基因组注释技术。”蛋白质基因组学“一词由jaffe等于2004年首次提出，作者采用串联质谱数据匹配dna翻译得到氨基酸序列的研究方法成功应用于运动型支原体(mycoplasmamobile)基因组的注释过程中。后来这种分析方法被应用到注释病毒基因组，原核生物基因组以及真核生物基因组。
蛋白基因组学直接对编码基因的表达产物--蛋白质进行研究，不仅能在蛋白质水平上验证基因表达和基因精细模式，还能提供蛋白质组层面特有的信息，如翻译后修饰、信号肽等。该方法有别于基因组学、转录组学和功能基因组学，对解释基因组注释结果和理解生命现象来说是非常重要的。在无法获得转录组数据的研究内容中，比如动物的血液和其他组织液中，蛋白质基因组学更是起到了不可替代的作用。另外，基于蛋白质基因组学的策略还可以被用在通过分析蛋白质水平来鉴定一些基于个体的序列变异导致的癌症，并确定其中最相关的变异上。2011年国际人类蛋白质研究组织(hupo)启动人类蛋白质组计划，通过整合多组学数据对可变剪接、snp以及三类主要翻译后修饰(磷酸化、乙酰化、糖基化)进行注释以加深人们对蛋白质(基因)功能的理解，从而指导疾病的研究，这是蛋白质基因组学一个重要的尝试和应用。而2014年在nature杂志上发布的人类蛋白质组草图的工作标志着蛋白基因组学已成为目前基因组注释和功能研究的有力工具。越来越多的研究人员提倡在基因组注释工程中加入标准蛋白质组分析作为互补，甚至直接采用蛋白质基因组学的方法对完成测序的基因组进行注释，体现了蛋白质基因组学在蛋白质组层面上注释基因组特有的优势。
自2004年蛋白基因组学的概念被首次提出，10多年的时间已经过去了。这些年间，得益于质谱技术的重大突破及蛋白质分离和富集方法的进步和rna-seq等转录组学研究的飞速发展，使得蛋白质组学研究中的关键问题灵敏性和准确性得到很大的解决。这也推动蛋白基因组学分析方法、流程设计的发展，以及各种流程化和系统化分析软件的产生。这些方法流程已被应用于一些物种基因组的重注释。尽管能从这些方法流程的建立中获益，但目前蛋白基因组学分析方法的研究尚处于起步阶段，其分析方法在计算上仍具有挑战，不可避免的存在一些问题：a.数据搜索的局限性:直接使用基因组构建真核生物蛋白基因组数据库，加上其他搜索条件带来的搜索空间的迅速膨胀，不可避免的会拉长搜索时间并降低质谱匹配的敏感性(鉴定数目)和特异性(正确鉴定数目)。面对真核生物庞大的基因组数据，通过何种方法构建蛋白基因组搜索数据库，才能应对数据膨胀带来的负面效应？b.精准鉴定:数据分析后端缺乏可靠而深入的信息学分析和评估策略，错误注释没有被有效地过滤，直接影响了基因组注释质量。尤其对于真核蛋白基因组分析，关注重点在新肽段的鉴定，搜索结果混合有已注释肽段和新肽段，传统的全局发现错误率评估未考虑新肽段类别错误率问题，致使研究结果带有较高的假阳性。虽然考虑了新肽段类别错误率问题，但数据分析端仍缺少准确评估新肽段类别错误率的策略。因此，如何正确的评估新肽段与已注释肽段的发现错误率、提高鉴定结果的精准度仍然是一个具有挑战性的问题。c.深度分析：目前的蛋白基因组学研究主要集中在核酸层(即编码基因注释和基因结构注释)，而更深入的、对蛋白功能研究更重要的蛋白组层面特有的现象却很少触及。比如翻译后修饰-对基因组功能注释具有重要的作用。翻译后修饰的分析能开启蛋白基因组学研究新的方向吗？如何实现翻译后修饰的大规模鉴定？d.标准化的应用:除了过高的新肽段发现错误率，数据分析端鉴定结果筛选标准的宽松性，也会使研究结果中混入过多的错误注释。同时，由于筛选标准的多样化，使得同一种分析手段确得到不同的研究结果。因此，建立蛋白质组基因学研究共同的数据分析的标准是极其有必要的。标准的蛋白基因组分析流程能提升研究结果的准确性和一致性，然而研究人员常常需要手动综合多种工具来完成整个流程的数据分析，极大加重了研究人员的工作负担。自动化的流程处理能为新物种的蛋白基因组分析带来便利，如何将理论和实际分析的经验想结合，开发一套完整而标准的适用于真核生物的蛋白基因组分析软件来支持分析工作，是目前从事相关分析的研究人员急切希望解决的。此外，该领域内还存在着一些亟待解决的问题，比如，哪些新的剪切变异体片段可以被翻译成稳定的有功能的蛋白质？dna变异的数据如何在蛋白质组水平调控发挥提供有价值的信息？
技术实现要素：
本发明的目的是为了克服上述现有技术的缺陷和问题而提供一种基于蛋白质质谱数据注释真核生物基因组的生物信息学方法。具体方法包括(1)构建高覆盖度的真核生物多组学序列数据库；(2)去除真核蛋白序列数据库冗余；(3)质谱原始数据格式转换；(4)采用多种不同算法的数据库搜索引擎，分别对质谱数据进行检索；(5)对检索及处理后的结果分别进行肽段图谱匹配打分；(6)对经类别fdr体系评估后的结果数据进行筛选；(7)验证已注释编码基因；(8)鉴定未注释新基因；(9)可变剪接的鉴定；(10)功能性点突变的鉴定；(11)针对蛋白质翻译后修饰进行大规模鉴定；(12)新基因和翻译后修饰的功能性注释。
1.构建高覆盖度的真核生物多组学序列数据库
(a)建立六阅读框翻译数据库
使用真核生物的基因组序列进行六阅读框翻译。建库时充分考虑真核生物的翻译特性，例如采用真核生物的标准密码子表和物种特异性的密码子表进行数据库序列翻译。按照中心法则和终止子-终止子(stop-to-stop)翻译方式，在真核基因组上三个连续碱基翻译成一种氨基酸，以三位为周期，加上dna双链因素，翻译出六种不同的、长度大于38的全酶切肽段，构建六阅读框翻译序列数据库。
(b)建立n端肽段数据库
已有文献报道真核生物中存在除了以atg翻译起始密码子外的其他翻译起始密码子(比如gtg和ttg)。三连密码子作为翻译起始密码子，只会翻译成甲硫氨酸。所以同一种密码子在一条蛋白序列中可能对于两种不同的氨基酸，即起始的氨基酸和中间序列的氨基酸。且翻译过程中无法确定蛋白的翻译起始位点。为了校正已注释基因的翻译起始位点注释错误或鉴定新的翻译起始位点，我们在六阅读框翻译数据库建立的基础上，考虑多种起始密码子，枚举所有长度大于38且以甲硫氨酸开头的全酶切肽段。
(c)建立从头预测蛋白序列数据库
为了提高鉴定翻译起始位点的鉴定，同时为新基因的鉴定提供证据，我们使用开源的genemarks软件包的隐马尔可夫模型算法结合蛋白编码模型、非编码区域集和基因调控位点集，产生编码蛋白基因的从头预测序列，来构建从头预测蛋白序列数据库。
(d)建立转录组翻译序列数据库
考虑到不同的数据来源，拟通过两种方式构建转录本翻译序列数据库。第一种方式，ncbi、ensemble、ucsc等公共数据库上有该物种组装好的转录本数据，直接下载到本地。逐一对其转录本序列，按照中心法则和终止子-终止子翻译方式，考虑真核生物的特异性的密码子，仅考虑单链因素，翻译出三种不同的、长度大于38的全酶切肽段，构建三阅读框翻译的转录组序列数据库。第二种方式，公共数据库上没有该物种组装好的转录本数据，则可在ncbi上下载该物种的转录组测序数据，使用开源软件trinity将转录组测序片段自组装成转录组序列。根据(a)步骤的建库策略，继续构建三阅读框翻译的转录组序列数据库。
(e)整合多组学数据库
建立多组学数据库会出现不同数据库之间序列的重复，则需将其重复序列去除，降低数据库的冗余性。氨基酸序列相同而命名编号不同的序列，则保留一条氨基酸序列，命名编号需按单斜杠重新拼接成。同时由于基因组测序和组装中存在准确性和完整性的局限，基因组序列中不可避免的会出现“n”，建库过程中被翻译成“x”。为了提高后续分析的精准性，出现三个以上“x”数据库序列被过滤。
2.去除真核蛋白序列数据库冗余
(a)从对应的公共数据库上下载该物种的非编码rna、假基因、未编码基因序列和est序列数据，分别按照中心法则和终止子-终止子翻译方式，考虑双链因素，翻译成六种不同的、长度大于38的全酶切肽段序列。
(b)按照第一点(e)步骤的整合策略，将上一步骤中创建的四类数据库，整合成去冗余数据库。
(c)对构建的真核蛋白序列数据库进行过滤。假如真核蛋白序列数据库中出现与去冗余数据库相同的序列，这些序列将从真核蛋白序列数据库中移除，缩减其数据库容量。
(d)对去冗余的真核蛋白序列数据库，根据序列长度和氨基酸字母排序，创建数据字典、建立序列索引，重建并优化其数据库的存储结构，为后续进一步的高效率的数据库搜索和数据分析打下基础。
3.质谱原始数据格式转换
由于不同品牌的质谱仪产生的下机原始数据格式不同，为了提升后续数据库搜索和数据分析的统一性与效用性，我们拟使用开源的proteowizard质谱数据转换软件，将每一个质谱原始数据转换成标准的mgf数据格式。proteowizard软件可以处理来自不同格式的质谱原始数据，也方便集成至现有的流程软件中(比如gape软件)。转换后的mgf格式数据，使用来自openms的peakpickerhires工具进行质心校正，然后被openms中的filemerger工具统一进行数据融合。
4.复合式搜索引擎策略
为了提升蛋白基因组学肽段鉴定的灵敏性，拟采用多种不同算法的数据库搜索引擎，分别对质谱数据进行检索。比如，x！tandem(搜索算法基于超几何分布模型)，msamanda(基于累积二项分布函数来确定显著性的谱图肽段匹配)，ms-gf+(使用标量积打分功能去评估计算谱图肽段匹配)，omssa(使用泊松分布算法来确定显著性的谱图肽段匹配)，comet(则使用基于不同属性的参数的经验打分算法)等等。数据库检索完成之后，分别对五种搜索引擎的结果进行格式归一化，处理并整合成统一的结果数据格式。
5.建立新肽段类别发现错误率评估方法
对检索及处理后的结果分别进行肽段图谱匹配打分。将打分后的结果肽段使用内建的blastp算法回帖到物种蛋白库和基因组上，并将结果肽段分类为已注释肽段谱图匹配(能够比对到蛋白库上)和广义新肽段谱图匹配(未能比对到蛋白库上)。针对已注释肽段谱图匹配结果，根据打分结果利用靶标-诱饵数据库搜索来计算全局fdr，对鉴定结果进一步的评价，以获得可靠的肽段谱图匹配结果；针对鉴定到的广义新肽段谱图匹配结果(突变肽段、可变剪接肽段、未注释肽段等)建立假发现率评估技术体系，主要结合全局fdr去推算各类别fdr，建立更加可信的适合蛋白质基因组学研究的肽谱匹配打分和评价方法，提高鉴定结果精准度。
6.建立严苛过滤标准
我们拟采用一套严苛的筛选标准体系，对经类别fdr体系评估后的结果数据进行筛选。首先，拟采用以下标准对其结果数据进行全局性过滤：(a)过滤同一张谱图匹配到被不同搜索引擎鉴定到的不同肽段序列的肽段谱图匹配结果；(b)一个肽段比对到不同的基因组位置的肽段谱图匹配结果(共享肽段)，也被过滤掉；(c)过滤掉长度小于7的肽段结果；(d)过滤污染库序列以及反库序列。其次，针对鉴定到的广义新肽段结果(突变肽段、可变剪接肽段、未注释肽段等)拟使用以下过滤标准：(e)过滤掉长度大于29aa的肽段结果；(f)仅保留全酶切的肽段结果；(g)含有脱酰胺化修饰或者脲甲基化修饰的肽段谱图匹配被过滤掉；(h)过滤单肽段序列(一个肽段序列只对应一个蛋白)。
7.验证已注释编码基因
已注释肽段的结果经过严苛过滤标准过滤后，被回溯到基因组上用于验证已注释编码基因的表达和结构。
8.鉴定未注释新基因
建立转录本数据与基因组数据之间的索引文件，将经过第六步骤严苛筛选标准过滤后的广义新肽段的结果比对到索引文件上，比对上的结果被用来筛选未注释新肽段。该类肽段不包含在原注释蛋白数据库中，但可与核酸序列信息匹配(genomesearchspecificpeptide,gssp)，它们被回溯至基因组进行重注释并分类，同时通过自主研发的基因预测算法鉴定新基因和校正已注释基因的结构。最后，综合多组学的信息对重注释的新基因结果进行验证(如利用mrna表达数据辅助蛋白确认)，提高注释结果的可信度。
9.可变剪接的鉴定
广义新肽段的结果经上一步筛选后，另一部分的新肽段集合首先被用来在肽段水平上鉴定跨越剪接位点的肽段。拟对该部分新肽段按位置从左至右依次循环分割成两段氨基酸片段，通过直接枚举的方法比对到六阅读框翻译数据库上，参考标准剪接位点的组合形式(gt/ag,gc/ag或at/ac)在规定长度内枚举可能的外显子剪接位点，并确认跨越剪接位点的肽段。该种算法鉴定到的跨越剪接位点的肽段，可能会因字符排列上的随机性和偶然性导致错配。因此我们通过延长新肽段分割的两段氨基酸片段分别匹配到的六阅读框翻译序列，建立一种评估方法来从鉴定的跨越剪接位点肽段中，筛选高可信度的可变剪接肽段。之后，基于自主研发的基因预测算法，针对筛选的高可信度可变剪接肽段集合进行新基因鉴定，以发现因可变剪接而未能预测到的新基因以及校正已注释基因的结构。同时也会综合多组学的信息(如mrna表达数据)对重注释的新基因结果进行验证。
10.功能性点突变的鉴定
拟对广义新肽段集合过滤未注释新肽段集合与跨越剪接位点肽段集合，来获得可能存在功能性点突变的新肽段集合。点突变新肽段集合中长度小于10的肽段被过滤掉。之后，点突变的新肽段被比对到六阅读框翻译数据库上，通过自主研发的功能性点突变搜索算法，搜索最多任意两个点突变的新肽段。该点突变搜索算法在比对中可能会因氨基酸字符排列上的随机性和偶然性导致错配。所以我们通过上一步所述的评估方法从鉴定到的功能性点突变肽段中，筛选高可信度的点突变肽段。
11.针对蛋白质翻译后修饰进行大规模鉴定
拟建立海量质谱数据中蛋白质非限制性翻译后修饰分析的新算法，结合unimod蛋白质翻译后修饰类型数据库，通过开放式蛋白序列数据库检索，鉴定不同类型的已知和未知的蛋白质翻译后修饰，并对翻译后修饰位点实施精准定位评估，提升质谱数据鉴定的准确度，实现真核生物翻译后修饰的快速大规模发现和精准解析。同时，建立质谱数据中蛋白质非特异酶切肽段鉴定的新算法，结合非特异酶切肽段检索鉴定结果和预测工具(如signalp、predisi等)，对蛋白质n端的信号肽进行注释与验证。
12.新基因和翻译后修饰的功能性注释
新基因(包括未被注释肽段和可变剪接肽段注释的新基因)的注释通过执行go(geneontology)功能注释来完成。鉴定到的新基因通过使用blastp算法比对到来自于ncbi的uniref数据库上，并通过得到的每个新基因的e-value值来筛选。然后，再将其结果比对到go功能数据库上，挑选最合适的功能标签对新基因进行注释。翻译后修饰结果的注释使用cog(clusteroforthologousgroupsofproteins)功能注释来完成。鉴定到的翻译后修饰蛋白被比对到cog数据库上，同时进行数据转换、处理和整合来实现翻译后修饰的注释。
附图说明
为了更清楚地说明本发明专利的技术方案，下面将对本方法的技术路线以附图的形式作简单地介绍。
图1为一种基于蛋白质质谱数据注释真核生物基因组的生物信息学方法的路线图。
具体实施方式
展示一下实例来具体说明发明的某些实施例，且不应解释为限制本发明的范围。对本发明公开的内容可以同时从材料、方法和反应条件进行改进，所有这些改进，均应落入本发明的的精神和范围之内。
实施案例
黄曲霉菌质谱数据10g,建立六阅读框翻译数据库、n端肽段数据库、从头预测蛋白序列数据库、转录组翻译序列数据库、整合多组学数据库,构建高覆盖度的真核生物多组学序列数据库。
从对应的公共数据库上下载该物种的非编码rna、假基因、未编码基因序列和est序列数据，分别按照中心法则和终止子-终止子翻译方式，考虑双链因素，翻译成六种不同的、长度大于38的全酶切肽段序列。
(b)按照第一点(e)步骤的整合策略，将上一步骤中创建的四类数据库，整合成去冗余数据库。
(c)对构建的真核蛋白序列数据库进行过滤。假如真核蛋白序列数据库中出现与去冗余数据库相同的序列，这些序列将从真核蛋白序列数据库中移除，缩减其数据库容量。
(d)对去冗余的真核蛋白序列数据库，根据序列长度和氨基酸字母排序，创建数据字典、建立序列索引，重建并优化其数据库的存储结构。
用x！tandem(搜索算法基于超几何分布模型)，msamanda(基于累积二项分布函数来确定显著性的谱图肽段匹配)，ms-gf+(使用标量积打分功能去评估计算谱图肽段匹配)，omssa(使用泊松分布算法来确定显著性的谱图肽段匹配)，comet(则使用基于不同属性的参数的经验打分算法)进行数据库检索。数据库检索完成之后，分别对五种搜索引擎的结果进行格式归一化，处理并整合成统一的结果数据格式。
对检索及处理后的结果分别进行肽段图谱匹配打分。将打分后的结果肽段使用内建的blastp算法回帖到物种蛋白库和基因组上，并将结果肽段分类为已注释肽段谱图匹配(能够比对到蛋白库上)和广义新肽段谱图匹配(未能比对到蛋白库上)。针对已注释肽段谱图匹配结果，根据打分结果利用靶标-诱饵数据库搜索来计算全局fdr，对鉴定结果进一步的评价，以获得可靠的肽段谱图匹配结果；
(a)过滤同一张谱图匹配到被不同搜索引擎鉴定到的不同肽段序列的肽段谱图匹配结果；(b)一个肽段比对到不同的基因组位置的肽段谱图匹配结果(共享肽段)，也被过滤掉；(c)过滤掉长度小于7的肽段结果；(d)过滤污染库序列以及反库序列。其次，针对鉴定到的广义新肽段结果(突变肽段、可变剪接肽段、未注释肽段等)拟使用以下过滤标准：(e)过滤掉长度大于29aa的肽段结果；(f)仅保留全酶切的肽段结果；(g)含有脱酰胺化修饰或者脲甲基化修饰的肽段谱图匹配被过滤掉；(h)过滤单肽段序列(一个肽段序列只对应一个蛋白)。
最后，通过该方法鉴定到黄曲霉菌93.8％的全部预测的编码基因，鉴定到的肽段序列覆盖了81.2％的编码基因序列，发现了39种蛋白质的翻译后修饰，发现了219个新基因，校正了47个预测的基因。
综上所述，本发明解决了常规方法注释真核生物基因组效率低，准确度低、覆盖度低等一系列问题，全面提升蛋白质质谱数据分析的准确度和灵敏度，实现了对真核生物基因组的深度解析和注释，具体有高效、准确、全面的特点。
虽然已经针对具体特征对本发明作了详细描述，然而本领域技术人员明显可知，该描述仅是优选的实施方式，并不限制本发明的范围，因此，本发明的实质范围将通过所附权利要求及其等同体来限定。