CN105740650B - 一种快速准确鉴定高通量基因组数据污染源的方法 - Google Patents

Abstract

本发明公开了一种快速准确鉴定高通量基因组数据污染源的方法，该方法首先组装denovo测序的原始基因组测序数据，得到组装结果，将组装结果进行基因预测，并翻译得到基因所对应蛋白的氨基酸序列，将组装的基因组序列和氨基酸序列分别与NCBI的NT数据库和NR数据库进行blast比对，得到有同源性的序列，作为原始比对数据库；从原始比对数据库中，提取序列对应的物种信息并排序，将序列对应的物种从多到少进行排序，结合基因数据的结果和氨基酸数据的结果，综合判断是否存在外源污染。本发明的方法可大限度的降低基因组denovo项目中，外来污染源对高通量基因组测序数据的污染和对后续生物信息学分析的影响，并提高污染源鉴定的速度和效率。

Description

一种快速准确鉴定高通量基因组数据污染源的方法

技术领域

本发明属于分子生物学技术领域，涉及一种快速准确鉴定高通量基因组数据污染源的方法。

背景技术

高通量测序技术(High-throughput sequencing)又称“下一代”测序技术，可以一次对几十万到几百万条DNA分子进行序列测定。

近年来，随着高通量测序技术的测序通不断提高、运行时间不断缩短、测序片段不断增长、成本不断降低，使高通量测序技术的应用范围越来越广，越来越多的团队选择通过高通量测序方法开展科学研究、辅助育种等，随着海量基因数据被测序，越来越多的物种公布了全基因组数据，为更进一步的研究打下了很好的基础。

但是测序样品的外源污染问题一直是不容忽视的问题，给后续数据分析造成很大的影响和障碍。

造成污染的原因有很多，主要的原因有以下几点：一是由于负责样品制备的试验人员的不规范操作，引入外源污染；二是很多样品的收集并不是也无法在无菌的环境下进行，空气中本身就存在着很多微生物；三是在文库制备及测序的过程中，也很可能因为实验人员的不规范操作会造成不同样品间的交叉污染。

在全基因组denovo从头测序中，由于测序通量通常比较高，因此小部分的污染也可能使最终产出的raw data中存在相当可观的污染数据量，对后续的基因组组装结果造成很大影响，进而影响后续的基因结构注释、基因功能注释和后续的生物信息学分析，直接关系到项目的成败。因此，除了在样品准备、文库制备及测序阶段要尽量规范操作避免外源污染外，在拿到测序平台的下机raw data数据开展后续的生物信息学分析之前，还需要采取高效准确的方法从GB甚至TB级的海量测序原始数据中快速对污染情况进行定性，并快速锁定外源污染源，进而对外源污染数据进行去除，得到去污染后的clean data后再开展后续的生物信息学分析。

由于生物信息学是新兴学科，很多数据分析人员由于刚接触高通量基因组测序数据的生物信息学分析工作，缺乏足够的知识储备，对外源污染源污染问题缺乏认识，或对污染问题对后续生物信息学分析工作的影响问题认识不足，所以很多时候对测序的原始数据不采取任何质控或污染鉴定和去除的工作，或仅仅进行简单的质控分析去除低质量的reads和接头后，就直接进行后续的生物信息学分析。

一般情况下，对高通量测序数据的去污染原理是利用blast程序，设置特定的参数和比对阈值(主要包括identity和e value两个参数，分别设置阈值为80％和1e-5)，与NCBI的NT数据库(非冗余核酸库)进行比对，通过比对上的序列对应的物种信息，来确定测序数据是否有污染，以及具体是何种物种导致的污染。由于高通量测序平台下机的raw data数据量非常大(达到GB甚至TB级别)，而NT数据库所包含的数据信息更是巨大的。因此，对原始数据逐条进行比对就意味着原始数据的每一条都要与每条NT数据库中的序列进行一一比对，从而找到最佳匹配结果，然后才能进行下一条测序数据的比对。采用该算法所花费CPU资源巨大，时间消耗巨大，并随着测序数据量的增加呈正比例增加，少则数周，多则数月都无法完成。

为了控制与NT数据库比对的时间成本，目前较普遍的方法是对每份数据进行随机抽样，再将抽样的数据进行与NT数据库的比对工作(对于双端测序的数据，一般是随机抽取10000对数据)。该方法能够较好的降低污染鉴定工作的时间成本，但是也存在明显的问题。就是由于抽样是随机的，导致基于抽样数据进行的污染鉴定分析很难准确反映测序数据整体的污染情况。特别是对于测序深度非常深、测序数据量非常大的项目，抽样数据在全部数据中所占的比例非常有限，几乎不可避免出现污染鉴定结论与实际真实污染的偏差，甚至出现结论完全错误的情况，比如实际上是存在某种污染源物种造成的污染，但是由于测序数据本身比较大，导致抽样数据中没有涵盖该污染数据，导致无法正确地识别该污染源物种造成的污染。

不管是抽样还是不抽样，两个方法都普遍存在的问题是，采用二代测序获得的reads都非常短，一般在100-250bp之间，为了污染鉴定的准确性，比对参数中设置的阈值一般都比较高(主要包括identity和e value两个参数，分别设置阈值为90％和1e-05)，比对结果中低于这个阈值的序列就会被认为不是该污染源。而对于突变频繁的位点，本身的遗传多样性就比较高，这样就会导致很多时候，污染情况被低估。

总之，目前存在对污染问题对分析工作的影响认识不足的情况；而目前所普遍采用的污染鉴定及排除方法中，全部测序数据与NT数据库比对存在对CPU资源占用巨大，时间耗费长的缺点；抽样的方法存在因抽样导致对污染的真实情况评估不准甚至评估错误的风险；两种方法都存在因阈值过高导致的污染情况被低估的风险，进而影响后续的污染去除和后续的生物信息学分析工作。

发明内容

为了解决上述技术问题，本发明的目的在于提供一种快速准确鉴定高通量基因组数据污染源的方法。本发明的方法可大限度的降低基因组denovo项目中，外来污染源对高通量基因组测序数据的污染和对后续生物信息学分析的影响，并提高污染源鉴定的速度和效率。

为实现上述目的，本发明采取的技术方案是一种快速准确鉴定高通量基因组数据污染源的方法，包括以下步骤：

(3)组装denovo测序的原始基因组测序数据，得到组装结果；

(4)将组装结果与NCBI的NT数据库进行blast比对，得到同源性的序列，作为原始比对数据库；

(3)从原始比对数据库中，提取序列对应的物种信息并排序，将序列对应的物种从多到少进行排序，判断是否存在外源污染；

(4)将组装结果进行基因结构的注释并获得对应的基因所翻译蛋白的氨基酸序列，将氨基酸序列与NCBI的NR数据库进行blast比对；

(5)得到步骤(4)的比对结果后，再依照步骤(3)的方法，提取物种信息和排序，将氨基酸序列对应的物种从多到少进行排序，判断是否存在外源污染；

(6)结合步骤(3)和步骤(5)的结果，根据两步分析统计中均存在的结果，最终确定污染情况及具体污染源。

其中，步骤(1)中，优选采用基于k-mer的开源基因组组装软件进行组装；更优选地，所述基于k-mer的开源基因组组装软件为ABYSS、SOAPdenovo、Velvet或ALLPATHS-LG软件；

其中，步骤(2)中，进行blast比对时使用的具体命令为：blastall-p blastn-dBlastDB-i in_file.fasta-m8–b 20-e 1e-5>blast_output；

其中，步骤(3)中，提取物种信息和排序的功能是通过基于Java程序的算法来实现的，具体包含以下两个步骤：

①java

PollutedIdentificationgetAnnotationFromNCBIDBByBlastOutblastOutResultncbiDB>result1，本步骤从NCBI的数据库中获取原始比对数据库中每条序列对应的注释信息，注释信息中就包含对应的物种信息；②java PollutedIdentificationgetPollutedSpeciesInfoByNCBIDBBlastOut result1ThreadHold>result2，本步骤设定e value≤1e-5、identity≥80％，将第①步的结果中的物种信息进行统计，凡是大于设定阈值的物种都会被统计，从而得到最终的物种统计结果，最终统计结果的格式为：每一个物种信息占一行，每行两列，第一列为物种信息，第二列为对应的比对到此物种的query序列总数，统计结果以文本文档的格式用记事本或Notepad++软件打开，然后将内容拷贝到excel表格中，按第二列统计数字按从大到小排序，排序后，就可以清晰直观的观察到最终的组装序列对应的物种从多到少的排序，进而根据物种的信息，判断是否有污染，以及主要污染源及其污染情况。

其中，在完成步骤(3)后，根据排序结果判断，如排在最前面的不是测序物种及近源物种，而是细菌或其他明显为污染的物种(比如测的是植物，而排在前面的是细菌、动物等)，且e value≤1e-5、identity≥80％，就说明污染较严重；或排在前面的是近源物种，但其他非近源物种的e value≤1e-5、identity≥80％，且每个占比或总量≥1％，也说明有污染。

其中，在步骤(4)中，使用AUGUSTUS软件将组装结果进行基因结构的注释并获得对应的基因所翻译蛋白的氨基酸序列。

其中，在步骤(4)中，进行blast比对时使用的具体命令为：blastall-p blastp-dBlastDB-i in_file.fasta-m8–b 20-e 1e-5>blast_output2。

通过以上技术方案，本实用新型的有益效果如下：

(1)本发明基于基因组组装结果和基因预测结果进行分析，因为组装后的序列和预测出的基因结果要比原始测序数据的条数少了非常多，首先克服了使用原始测序数据进行分析所耗时间太长的缺点；

(2)本发明方案的基因组组装结果和基因预测结果是采用基因组denovo组装软件和基因预测软件基于所有的原始测序数据得出，因此能够全面的反映测序数据，克服了基于对原始测序数据进行抽样从而导致无法真实反映污染的实际情况的缺点。

(3)常规技术中因原始测序序列长度过短，因此会出现变异大的片段因达不到阈值条件而被过滤掉，从而造成污染情况估计不准确的缺点；而由于组装结果和预测的基因序列较长，可以在更广的范围内与目标序列进行匹配，从而规避了这一缺点。

附图说明

图1为本发明的方法流程图。

具体实施方式

下面结合实施例对本发明的具体实施方式作进一步描述，本发明的优点和特点将会随着描述而更为清楚。但这些实施例仅是范例性的，并不对本发明的范围构成任何限制。本领域技术人员应该理解的是，在不偏离本发明的精神和范围下可以对本发明技术方案的细节和形式进行修改或替换，但这些修改和替换均落入本发明的保护范围内。

下述实施例中所使用的实验方法如无特殊说明，均为常规方法。

下述实施例中所用的材料、试剂等，如无特殊说明，均可从商业途径得到。

实施例1

某致病真菌(向日葵霜霉Plasmopara halstedii)基因组denovo测序，，二代illumina平台共2个文库180bp和500bp，测序深度分别是35X和34X，每条reads长度为100bp，每个文库的reads总数分别为46308070和43435185，共89743255条,总数据量为8.36G，使用以下方法鉴定污染源：

(1)使用ABYSS软件组装(k-mer参数设置为k＝50，其他参数为软件默认参数)，组装结果中的scaffold个数共30428，N50为10506，最长为479848，大小为80M；可以很容易的看出：①组装后的序列条数共30428，仅仅为原始的序列总条数89743255的0.03％；②总数据量118M，仅为原8.36G总数据量的1.38％。③序列长度有100bp提升到N50为10506，为100bp的105倍，且最长可达479848。

上述结果①和②可大大降低比对的工作量，而结果③可保证在更长的序列范围内进行比对分析。

(2)将组装结果与NCBI的NT数据库进行blast比对，进行blast比对时使用的具体命令为：blastall-p blastn-d BlastDB-i in_file.fasta-m8–b 20-e1e-5>blast_output，得到有同源性的序列，作为原始比对数据库。

(3)从原始比对数据库中，提取序列对应的物种信息并排序，将序列对应的物种从多到少进行排序，判断是否存在外源污染；提取物种信息和排序的功能是通过基于Java程序的算法来实现的，具体包含以下两个步骤：

①java

PollutedIdentificationgetAnnotationFromNCBIDBByBlastOutblastOutResultncbiDB>result1，本步骤从NCBI的数据库中获取原始比对数据库中每条序列对应的注释信息，注释信息中就包含对应的物种信息；②java PollutedIdentificationgetPollutedSpeciesInfoByNCBIDBBlastOutresult1ThreadHold>result2，本步骤设定e value≤1e-5、identity≥80％，将第①步的结果中的物种信息进行统计，凡是大于设定阈值的物种都会被统计，从而得到最终的物种统计结果，最终统计结果的格式为：每一个物种信息占一行，每行两列，第一列为物种信息，第二列为对应的比对到此物种的query序列总数，统计结果以文本文档的格式用记事本或Notepad++软件打开，然后将内容拷贝到excel表格中，按第二列统计数字按从大到小排序，排序后，就可以清晰直观的观察到最终的组装序列对应的物种从多到少的排序，进而根据物种的信息，判断是否有污染，以及主要污染源及其污染情况。

(4)使用AUGUSTUS软件将组装结果进行基因结构的注释并获得对应的基因所翻译蛋白的氨基酸序列，共预测出22028个基因，将氨基酸序列与NCBI的NR数据库进行blast比对；进行blast比对时使用的具体命令为：blastall-p blastp-d BlastDB-i in_file.fasta-m8–b 20-e 1e-5>blast_output2；

(5)得到步骤(4)的比对结果后，再依照步骤(3)的方法，提取物种信息和排序，将氨基酸序列对应的物种从多到少进行排序，判断是否存在外源污染；

(6)结合步骤(3)和步骤(5)的结果发现，物种信息占比最多的为该真菌及其近源物种。另外，从比对结果发现，紧接着排在近源物种后面的分别是Pseudomonasalcaligenes(16个scaffold/1324个基因)和Chryseobacterium gleum(10个scaffold/593个基因)，且组装结果和预测基因结果的比对结果是一致的。从而可以快速鉴定出该数据中存在较多的两种细菌污染，该结论可为后续的建库测序提供参考，提示后续样品准备、文库制备及测序过程中应特别留意细菌的污染问题。

上述步骤的参考时间(仅供参考，还与具体数据量、软件、服务器性能、同一时间并行运行的程序等因素有关)：

ABYSS组装：约1天；组装结果的数据库比对：约3天；AUGUSTUS基因预测：约1天；基因结果的数据库比对：约2天；

共计约7天

该方案可在一周左右快速全面的对污染源进行鉴定分析，比对工作量和时间成本缩减至少100倍以上，比抽样数据更准确全面，可有效避免误判漏判。

一些概念：

Contig:基于reads之间的overlap关系，拼接获得的中间没有gap的序列称为Contig。

Scaffold:通过reads拼接获得Contigs后，根据paired-end或mate-pair信息，确定一些Contig之间的顺序关系，这些先后顺序已知的Contigs组成Scaffold。

Contig/Scaffold N50：假设组装得到总长为1M的基因组，把contig和scaffold从长到短进行排列，然后相加，当恰好加到1M的50％，也就是500k的时候，那一条contig或者scaffold的长度就叫做Contig N50和Scaffold N50。很明显这个数值越大说明组装的质量越好。

Claims (6)

1.一种快速准确鉴定高通量基因组数据污染源的方法，其特征在于，包括以下步骤：

(1)组装denovo测序的原始基因组测序数据，得到组装结果；

(2)将组装结果与NCBI的NT数据库进行blast比对，得到有同源性的序列，作为原始比对数据库；

(3)从原始比对数据库中，提取序列对应的物种信息并排序，将序列对应的物种从多到少进行排序，判断是否存在外源污染；

(4)将组装结果进行基因结构的注释并获得对应的基因所翻译蛋白的氨基酸序列，将氨基酸序列与NCBI的NR数据库进行blast比对；

(5)得到步骤(4)的比对结果后，再依照步骤(3)的方法，提取物种信息和排序，将氨基酸序列对应的物种从多到少进行排序，判断是否存在外源污染；

(6)结合步骤(3)和步骤(5)的结果，根据两步分析统计中均存在的结果，最终确定污染情况及具体污染源。

2.按照权利要求1所述的方法，其特征在于：步骤(1)中，采用基于k-mer的开源基因组组装软件进行组装。

3.按照权利要求2所述的方法，其特征在于：步骤(1)中，所述基于k-mer的开源基因组组装软件为ABYSS、SOAPdenovo、Velvet或ALLPATHS-LG软件。

4.按照权利要求1所述的方法，其特征在于：步骤(2)中，进行blast比对时使用的具体命令为：blastall-p blastn-d BlastDB-i in_file.fasta-m8–b 20-e 1e-5>blast_output。

5.按照权利要求1所述的方法，其特征在于：在步骤(4)中，使用AUGUSTUS软件将组装结果进行基因结构的注释并获得对应的基因所翻译蛋白的氨基酸序列。

6.按照权利要求1所述的方法，其特征在于：在步骤(4)中，进行blast比对时使用的具体命令为：blastall-p blastp-d BlastDB-i in_file.fasta-m8–b 20-e 1e-5>blast_output2。