本发明属于物种鉴定领域，涉及一种基于新一代测序的古生物物种鉴定方法。
背景技术：
：古生物样本对现代生物种群的进化史研究至关重要，古人类基因组的研究成果使人们重新认识到现代人的遗传组成并非只有非洲祖先成分，而是在走出非洲之后又与古尼安德特人和古丹尼索尔人发生过基因交流，颠覆了以往人们对现代人进化史的认识。同时，古生物基因组的研究对现代生物种群，尤其是人类的自然选择和疾病的研究也有着无法替代的重要作用，藏族人的高原适应性基因被证明是来自于古丹尼索尔人的基因组之间的渗透作用。古生物样本作为一种无法复制的遗传资源，对现代生物群体的进化、选择和疾病等研究具有巨大的促进作用而且无法替代。我国是一个古生物资源大国，有着极其丰富的化石和亚化石资源，如热河古生物群、湖南长沙马王堆古墓、以及近期发现并发表在国际顶级期刊《自然》上的古人类骨骼和牙齿样本等。对这些古生物资源的遗传信息的解读，对于汉族人的起源进化和复杂疾病的研究以及动植物的分类和进化研究都将起到革命性的推动作用。但是，在古dna领域的研究中，由于所有样本的采集几乎都不会直接来自于生物体本身，而是千百年来幸存下来的动物残骸，包括骨骼、毛发、牙齿等，而且长期暴露在恶劣的环境中，很多古生物样本的形态学特征早已不复存在。如果无法准确辨别古生物残骸所属物种，将无法进行后续生物学信息的解读和挖掘。古dna提取和二代测序技术在古生物样本中的测序应用已经趋于成熟，而目前并没有基于二代测序技术的古生物物种鉴定的方法，并且利用现代生物的物种鉴定的方法并不能做到对古生物样本的准确鉴定，因此，古生物物种鉴定成为古生物基因组研究的重要瓶颈，如果能突破这一技术瓶颈，势必会给整个行业带来更迅速的发展。目前考古学进行物种鉴定还是主要通过形态学方法。一般都要经过下述步骤：①熟悉标本外部形态和内部构造，对大化石的细微构造或微体化石，一般需要借助显微镜或电子显微镜进行观察。有时要将化石做连续切片，以便于了解其内部构造特征；②利用所具有的知识并查阅有关文献，确定较大的分类阶元，一般确定到科；③利用检索表、图版等文献资料，将标本进一步检索到属、种；④选择有代表性的种群标本或典型的单个标本进行特征描记，度量各种性状要素及照相。除此之外，基于pcr的物种鉴定技术可以对古生物样本的线粒体dna进行扩增，然后在数据库中进行比对。现有的古生物物种鉴定方法存在以下缺点：第一、基于古生物保存的形态进行物种鉴定，只能鉴定形态保存比较完好的古生物样品，对于一些历史比较久远的，没有完整形态的古生物残骸，用现有的方法无法实现古生物物种鉴定。第二、每个物种存在连续的个体变异。形态学方法要参考模式标本，而这些标本往往只是某个物种首次发现的标本，不一定全面地反映该种的特征。由于古生物样本获得不易，在鉴定物种时，不能尽可能多地、全面地采集标本，从而使得这些标本不能全面反映真正种群的总面貌，最后无法用各种统计方法来区分种间变异和种内变异，导致鉴定物种时出现一定的误差。第三、如果进行大批量古生物样品鉴定，依赖于形态特征的鉴定方法需要大量的人力物力，人力成本过高。第四、基于pcr的物种鉴定方法对古生物样本的保存条件要求很高，大多数古生物样本无法达到。此外，古生物样本有其特有的dna分子特征，通过pcr的方法会增加犯错误的概率。技术实现要素：本发明的目的是提供一种基于新一代测序的古生物物种鉴定方法。本发明所提供的鉴定古生物物种的方法，具体可包括如下步骤：(a)从来源于古生物(如古动物或古植物)的待测样本中提取dna，进行高通量测序，得到原始测序数据。(b)对所述原始测序数据进行过滤处理，去除长度大于阈值a的读长；所述阈值a为70～150bp，具体如99bp。因为古dna序列高度片段化，平均长度一般在50-70bp，如果读长过长(大于所述阈值a)，则很有可能是来自于现代dna的污染，因此，为了最大限度保留古dna，这些读长将被删除。这一步极其重要，如果不删除大于所述阈值a的读长，将会影响后续物种鉴定的准确性，这也是跟现代生物样本物种鉴定的一个很大的不同点。(c)构建动物或植物线粒体核酸数据库，作为对比数据库。即，若所述待测样本为动物样本，则构建动物线粒体核酸数据库作为对比数据库；若所述待测样本为植物样本，则构建植物线粒体核酸数据库作为对比数据库。本发明之所以选择线粒体核酸库而没有选择基因组数据库进行比对，是为了防止比对过程中出现随机偏差，这也是本发明的重要发明点。因为不同物种的基因组数据量在数据库中相差巨大，有些物种在国际上进行过大量的基因组测序计划(如国际千人基因组计划)，这些物种在数据库中具有海量的dna序列，而很多物种并没有进行过全基因组测序，这种情况下，数据库中不同物种的序列储量就会相差极大。而古dna的序列高度片段化，平均长度50-70bp，如此小的片段很容易发生随机错误比对，而这种错误比对会随着基因组数据量的增加而增加。也就是说，某个物种在数据库中的数据越多，由于随机错误导致的比对到这个物种基因组上的概率就大，自然比对上的读长也就越多，这会对物种鉴定结果产生致命的错误。(d)将步骤(b)经过过滤处理后的测序数据与步骤(c)所建立的比对数据库进行比对，从比对结果中去除所有的具有100％相似度读长的部分，然后根据剩余比对结果确定所述待测样本所属的物种。在本发明中，步骤(a)中，所述高通量测序可为二代测序，具体如illuminahiseq测序。在步骤(a)中，进行所述高通量测序时构建dna文库的过程中，由于古dna的片段化程度很高，因此不需进行dna片段化处理，dna提取完成后直接进行dna文库构建。在步骤(a)中，进行所述高通量测序时，可将读长控制在100bp以内，不宜选择长片段测序，因为古dna平均长度在50-70bp左右，如果测序时读长的长度超过100bp，一方面会引入大量的接头污染，另一方面会造成大量的数据浪费。另外，测序数据量每个样本10gb足以支撑本发明进行古生物物种鉴定。在步骤(b)中，对所述原始测序数据进行过滤处理时，除了去除长度大于所述阈值a的读长外，还包括如下处理中的至少一种；(1)若发现读长中包含有接头序列，则将接头序列部分切除，读长的其余部分保留。(2)若质量值q≤10的碱基数占整条读长总碱基数的50％以上，则删除整条读长；若质量值q≤10的碱基在读长的端头，且数量不超过整条读长总碱基数的50％，则仅切除质量值q≤10的碱基。(3)去除n区位于中间且含n比例大于10％的读长；如果n区仅存在于读长两端，无论读长多长仅切除读长两端的n区，其余碱基进行保留；其中，n的意思是表示任何一个碱基，表明测序时软件无法分辨是哪个碱基。(4)去除长度小于阈值b的读长；所述阈值b为10～50bp，具体如30bp。因为如果读长小于所述阈值b，在后续的比对过程中会造成较多的错误比对。在步骤(c)中，构建所述动物或植物线粒体核酸数据库时，从ncbi网站下载所有动物或植物线粒体dna的测序数据，并以100000条读长为一个单位进行切分，得到若干个子库，所有子库构成所述动物或植物线粒体核酸数据库，即所述比对数据库。步骤(d)中所进行的比对通常为blast比对，由于blast比对执行运算速度很慢，因此在使用过程中建议将对比数据库进行切分运算，以100000条读长为一个单位进行切分，然后将每一个切分文件进行独立的blast比对，在计算资源充足的前提下，可以在几个小时内完成比对。在步骤(d)中，从比对结果中去除所有的具有100％相似度读长的部分后，是将比对的相似度控制在大于或者等于95％，从而得到所述剩余比对结果的。具体的，从比对结果中去除所有的具有100％相似度读长的部分后，若所述经过过滤处理后的测序数据中存在读长a与所述比对数据库中物种a的线粒体dna序列具有95％或95％以上的相似度，则视为所述读长a来自于所述物种a；据此，将所述经过过滤处理后的测序数据中所有读长能够确定物种来源的均做相应确认，所得结果即为所述剩余比对结果。在步骤(d)中，在所述剩余比对结果中，被确认频次最高的的物种即视为所述待测样本所属的物种。古生物样本在长期的保存过程中，双链dna会受到两种最主要的化学损伤，即胞嘧啶脱氨基作用和脱嘌呤作用。脱氨基作用主要发生在dna片段的端头部位，这种脱氨基作用会使胞嘧啶转换成尿嘧啶，因此在文库构建和测序的时候会引入c->t的突变。可以通过判断dna两端是否存在c->t的突变来对古dna进行认定。但由于blast比对结果没有显示错配发生的位置信息，因此无法通过查找dna在5’和3’端的c->t错配信息来区分古dna和现代dna污染。因此，为了在最大限度保留古dna的同时，保证比对结果的准确性，本发明在对结果的分析过程中去除掉了所有的具有100％相似度的reads，保证在比对过程中是有错配的，这样可以最大限度的去掉现代dna污染和保留古dna。由于过滤后的读长长度约为50bp，我们在允许有2个错配的前提下，将比对的相似度控制在大于或者等于95％。在本发明中，步骤(a)中，所述原始测序数据为fastq格式；步骤(c)中，构建所述动物线粒体核酸数据库时，从ncbi网站下载的所有动物线粒体dna的测序数据保存成fasta格式；步骤(d)中，将所述经过过滤处理后的测序数据与所述比对数据库进行比对时，先将所述经过过滤处理后的测序数据的格式由fastq格式转换为fasta格式。本发明的关键技术发明点在于选择动物线粒体核酸库作为比对数据库，并针对古dna的特征，对大于阈值(所述阈值为70～150bp，具体如99bp)的读长进行过滤，同时去除了100％匹配的读长等，本发明大大降低了现代dna污染率，保证了物种鉴定的可靠性和准确性，任何通过本发明的关键技术发明点进行古dna物种鉴定的技术都可以达到类似技术效果。此外，将本发明应用于具有微量dna特征的法医鉴定领域，特别是经过鞣制后的毛皮等的物种鉴定也视为本发明的保护范围。本发明基于高通量测序数据和生物信息学方法进行古生物物种鉴定，可以解决现有技术的不足。第一、对于一些无法通过形态学特征鉴定物种的古生物残骸(包括骨骼、毛发和组织等样本)，只要能够提取到古dna并最终获得高通量测序数据，根据古生物dna特征，可以用一系列生物信息分析方法进行这类古生物的物种鉴定。第二、虽然同一物种形态学上会有较大差异，但是基因组水平上的差异微乎其微，而且不同物种之间的基因组数据差异又较大，因此，通过基因组测序数据可以巧妙的避免因个体差异造成的误判，从而有效的进行物种鉴定。第三、高通量测序技术已经实现了自动化，生物信息学方法通过编写好的程序也可以一次性处理大批量的数据，均可以同时对大量古生物样品进行物种鉴定，从而节省了时间和成本，提高了效率。第四、高通量测序技术可以对残留在古生物样本中极微量的dna序列进行捕获测序，对于样本保存条件依赖性低，而且可以根据古dna的数据特征来挑选有效序列进行物种鉴定。有效地避免了基于pcr技术物种鉴定的不足。附图说明图1为基于高通量测序的古生物物种鉴定方法的总体流程图。具体实施方式下述实施例中所使用的实验方法如无特殊说明，均为常规方法。下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。实施例1、基于新一代测序的古生物物种鉴定方法的建立本发明基于illumina二代测序数据和生物信息学分析方法，研发出了一种鉴定古生物物种的方法。主要包括以下几大步骤：获取illumina二代测序数据；原始下机数据的质控；比对数据库的构建；blast比对分析；结果的统计(图1)。具体技术细节如下：一、illumina二代测序数据的获取本发明基于illumina二代测序数据，因此，该发明所涉及的第一步就是获取数据，获取数据包括dna提取、dna文库构建和上机测序。其中dna提取、建库和上机测序不是本发明的主要内容，在此不进行详细描述，可以参考已经发表的文献中描述的试验方法完成。但是有一些细节需要在实验过程中注意：由于古dna的片段化程度很高，因此在构建dna文库的过程中不需对dna进行片段化处理，dna提取完成后可直接进行文库构建。此外，在上机测序时，不宜选择长片段测序，读长长度控制在100bp以内，因为古dna平均长度在50-70bp左右，如果测序时读长的长度超过100bp，一方面会引入大量的接头污染，另一方面会造成大量的数据浪费。测序数据量每个样本10gb足以支撑本发明进行古生物物种鉴定。二、原始下机fastq数据的质控原始fastq数据下机后的最重要的一步就是根据illumina数据特征以及古dna的序列特征对数据进行过滤，目的在于最大限度的去除低质量的序列以及被外源污染的dna序列。数据过滤主要包括4方面：第一，对接头进行过滤。在这一步中，如果发现读长中包含有接头序列，就将接头序列部分切除，读长的其余部分保留。第二，对低质量碱基进行过滤。质量值q≤10的碱基数占整条读长总碱基数量的50％以上时，删除整条读长。如果质量值q≤10的碱基在读长的端头，且数量不超过整条读长的50％，则仅切除质量值q≤10部分的碱基。第三，对n区进行过滤。去除n区位于中间且含n比例大于10％的读长；如果n区仅存在于读长两端，无论读长多长仅切除读长两端的n区，其余碱基进行保留。其中，n的意思是表示任何一个碱基，表明测序时软件无法分辨是哪个碱基。第四，去除长度小于30bp以及长度大于99bp的读长。如果读长小于30bp，在后续的比对过程中会造成较多的错误比对。因为古dna序列高度片段化，平均长度一般在50-70bp，如果读长过长(大于99bp)，则很有可能是来自于现代dna的污染，因此，为了最大限度保留古dna，这些读长将被删除。这一步及其重要，如果不删除大于99bp的读长，将会影响后续物种鉴定的准确性，这也是跟现代生物样本物种鉴定的一个很大的不同点。三、blast比对分析1、动物线粒体dna核酸数据库构建数据完成质控后，要将保留下来的读长与已知物种的核酸数据库进行比对，因此，需要构建一个尽可能包含所有物种序列信息的核酸库。本发明选择构建线粒体dna核酸库。构建数据库的步骤如下：①从ncbi网站下载所有动物线粒体dna的测序数据，并保存成fasta格式；②下载并安装blastall软件，此软件可从ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.3.0/免费下载安装，其中包含用于数据库构建的formatdb。2、blast比对本技术使用动物线粒体核酸库进行比对分析。比对前要先将fastq文件进行格式转换成fasta格式。由于blast比对执行运算速度很慢，因此在使用过程中建议将待比对fasta数据进行切分运算。本发明建议将fasta文件以100，000条读长为一个单位进行切分，然后将每一个切分文件进行独立的blast比对，在计算资源充足的前提下，可以在几个小时内完成比对。本技术之所以选择线粒体核酸库而没有选择基因组数据库进行比对，是为了防止比对过程中出现随机偏差。因为不同物种的基因组数据量在数据库中相差巨大，有些物种在国际上进行过大量的基因组测序计划(如国际千人基因组计划)，这些物种在数据库中具有海量的dna序列，而很多物种并没有进行过全基因组测序，这种情况下，数据库中不同物种的序列储量就会相差极大。而古dna的序列高度片段化，平均长度50-70bp，如此小的片段很容易发生随机错误比对，而这种错误比对会随着基因组数据量的增加而增加。也就是说，某个物种在数据库中的数据越多，由于随机错误导致的比对到这个物种基因组上的概率就大，自然比对上的读长也就越多，这会对物种鉴定结果产生致命的错误。本发明曾尝试使用全基因组数据作为数据库进行比对，结果发现，比对结果错误率极高。如在本案例中曾对古牛的测序数据进行全基因组数据库比对，最后发现大部分测序读长比对到了北美安乐蜥，经查证，北美安乐蜥曾发起过一个庞大的基因组测序计划，在数据库中拥有海量的物种测序数据(表1)。表1使用全基因组数据库和使用线粒体基因组数据库比对结果比较注：此表中结果为实施例2中9-1样本的比较结果相比全基因组数据，研究人员对线粒体基因组的研究在不同物种之间的差异就小多了。因此，本发明采用线粒体核酸库进行比对的策略，将测序读长与动物线粒体核酸数据库进行比对来降低由物种之间数据量大小对研究造成的影响。实际测试结果表明，用线粒体核酸数据库进行比对可以极大提高比对的准确性和可靠性。四、结果的统计古生物样本在长期的保存过程中，双链dna会受到两种最主要的化学损伤，即胞嘧啶脱氨基作用和脱嘌呤作用。脱氨基作用主要发生在dna片段的端头部位，这种脱氨基作用会使胞嘧啶转换成尿嘧啶，因此在文库构建和测序的时候会引入c->t的突变。可以通过判断dna两端是否存在c->t的突变来对古dna进行认定。但由于blast比对结果没有显示错配发生的位置信息，因此无法通过查找dna在5’和3’端的c->t错配信息来区分古dna和现代dna污染。因此，为了在最大限度保留古dna的同时，保证比对结果的准确性，我们在对结果的分析过程中去除掉了所有的具有100％相似度的reads，保证在比对过程中是有错配的，这样可以最大限度的去掉现代dna污染和保留古dna。由于过滤后的读长长度约为50bp，我们在允许有2个错配的前提下，将比对的相似度控制在大于或者等于95％。若某一读长与所述比对数据库中某一物种的线粒体dna序列具有95％或95％以上的相似度，则视为该读长来自于该物种；据此，将所有读长能够确定物种来源的均做相应确认，最终被确认频次最高的的物种即视为所述待测样本所属的物种。实施例2、基于新一代测序的古生物物种鉴定方法的验证为证明本发明的可行性和准确性，本发明对4例未知古生物样本进行了物种鉴定，样本编号分别为：9-1；9-2；10-1；11-1。同时，为了证明本发明所得结果的可靠性，本发明对5例已知物种信息的古人类样本对本发明进行了验证。一、4例古生物样本的物种鉴定本发明中涉及到的4例未知古生物样本中，包括3例古生物骨骼样本(分别记为9-2、10-1和11-1)和1例古生物组织样本(记为9-1)，所有4例古生物样本距今约2000多年。所有样本均由中国科学院古脊椎动物与古人类研究所提供。根据实施例1中所描述的方法对4例古生物进行物种鉴定：(1)基于illumina二代测序数据，4例古生物样本的dna提取和建库方法详见以下两篇参考文献：[1]n.rohland，m.hofreiter.ancientdnaextractionfrombonesandteeth[j].natureprotocols，2007，2(7):1756-1762.[2]m.t.gansauge，m.meyer.single-strandeddnalibrarypreparationforthesequencingofancientordamageddna[j].natureprotocols，2013，8(3):737-748.其中9-2采用的是illuminahiseq2000pe50测序策略，其余样本的测序均采用illuminahiseq2000pe101，每个样本的原始下机fastq格式的测序数据量详见表2。(2)按照实施例1中描述的数据过滤方法，对4例古生物样本的原始下机fasta数据进行了严格的过滤。具体执行标准如下：1)如果发现读长中包含有接头序列，切除接头序列部分；2)如果质量值q≤10的碱基数占整条读长总碱基数量的50％以上时，删除整条读长，如果质量值q≤10碱基在读长的端头，且数量不超过整条读长的50％，则仅切除质量值q≤10部分的碱基；3)去除n区位于中间且含n比例大于10％的读长；如果n区仅存在于读长两端，无论读长多长仅切除读长两端的n区，其余碱基进行保留；其中n的意思是任意一个碱基，表明测序时软件没有识别出是哪个碱基；4)去除长度小于30bp以及长度大于99bp的读长。过滤后的数据量详见表2。表24例未知物种样本在数据质控前后数据量(3)按照本发明实施例1中的技术流程，首先构建动物线粒体基因组数据库，然后将过滤后的数据比对到动物线粒体基因组核酸库。在比对过程中，为了在最大限度的排除现代dna的污染的同时，保证比对结果的准确性，本发明去除掉了所有的具有100％匹配的比对，保证在比对过程中是有错配的，这样可以最大限度的去掉现代dna污染并保留古dna。在本发明中已经做过尝试，由于本实施例中的骨骼样本很显然来自大型哺乳动物，但是，当加入100％匹配的比对后，发现大量的非目标物种序列会被比对出来，出现在最终结果，导致目标物种在结果中所占比例大大降低(表3)，这些都很可能来自于现代dna的污染，而去掉100％匹配的读长后，这些背景噪音明显降低。由于过滤后的读长长度约为50bp，我们在允许有2个错配的前提下，将比对的相似度控制在了≥95％。表3去除100％比对结果之后对物种鉴定结果的影响注：此表中结果为10-1样本的比较结果最后，根据结果统计出比对次数排在前10名的物种，具体的比对结果详见表4。表44例未知物种样本的物种鉴定结果根据最后的统计结果，我们得到的结论是：编号为9-2和11-1的古生物样本归属于马(equuscaballus)，编号为10-1和9-1的古生物样本归属于牛(bostaurus)。在所有的4个样本的物种鉴定结果中，排名第一的物种比对次数均远远大于排名第二的物种，其中二者比值有3例在10倍以上，其中一例为较小，但也大于4倍，均未出现难以分辨的结果。二、对本发明的验证为了验证本发明实施例1给出的技术方案对古生物物种鉴定的可靠性，本发明应用此技术方案对5例已知物种信息的古人类样本进行了物种鉴定，进而验证此方法是否准确可靠。该5例古生物样本均由中国科学院古脊椎动物与古人类研究所提供，样本年代大约在3000-8000年前，其中4例为古人类骨骼样本(human_1～human_4)，1例为古人类毛发样本(human_5)。在验证过程中，每个样本分别从质控后的数据中提取5，000，000条读长进行物种鉴定。所使用的实验和分析方法与步骤一中4例未知样本的物种鉴定方法完全一致，最后统计结果详见表5。表55例古人类样本物种鉴定结果根据最后的统计结果，即表5可知，我们得到的结论是：所有5例样本均归属于人类(homosapiens)。在所有的5个样本的物种鉴定结果中，均未出现难以分辨的结果，比对次数排名第一的占比在97％以上，均远远超过后面的结果，证明了此种方法进行物种鉴定的准确性和可靠性。当前第1页12当前第1页12