专利名称：一种基因序列数据的筛选方法
一种基因序列数据的筛选方法技术领域
本发明属于应用生物信息学技术领域，尤其涉及一种基因序列数据筛选方法，主要应用于生物系统发育、生物条形码、生物物种鉴定等相关领域的基因数据筛选和质量控制。
背景技术：
生物分子系统发育研究在不同水平和层次上依赖于对数据的使用从简单的检索到PCR污染物的检查，到寻找一个给定序列的类群同源性序列，到更全面的基于大量数据进行的类群和位点的数据挖掘(McMahon，Μ. Μ.，and Μ. J. Sanderson. 2006. ” Phylogenetic supermatrix analysis of GenBank sequences from 2228papilionoid legumes". Syst. Biol. 55 :818-836 ；Ciccarelli, F. D. , Τ. Doerks, C. von Mering, C. J. Creevey, B. Sne 1, and P. Bork. 2006. "Toward automatic reconstruction of a highly resolved tree of life. Science 311 1283-1287 ；Bininda-Emonds, 0. R. P. , Μ. Cardillo, K.Ε. Jones, R. D. Ε. MacPhee, R. Μ. D. Beck, R. Grenyer, S. Α. Price, R. Α. Vos, J. L. Gittleman, and Α. Purvis. 2007. "The delayed rise ofpresent-day mammals". Nature 446 :507-512 ；Li, C. H. ,G. 0rti,G. Zhang,and G. Q. Lu. 2007.，，Apractical approach to phylogenomics :The phylogeny of ray-finned fish (Actinopterygii) as a casestudy" . BMC Evol.Biol.7 44 ；MICHAEL J.SANDERSON, 1 DARREN B0SS，et al.2008."ThePhyLoTA Browser processing GenBank for Molecular Phylogenetics Research", Syst.Biol. 57(3) :335-346.)。
分子生物学的早期研究积累了大量的基因序列数据。以国际核算序列数据库 (International Nucleotide Sequence Database Collaboration, INSDC) t一的 GenBank 为例(Michael Y. Galperin.2011. "The Molecular Biology Database Collection 201 lupdaeNuc 1. Acids Res. 35 :D3_D4)，截至 2010 年 9 月统计的数字，传统的GenBank版本中在720，000, 000条序列纪录中有75，000, 000, 000碱基对数据；在WGS版本中有92，369，977，826碱基对的海量数据。
与生物分子系统发育学相关的最重要的注释是类群的名称和基因或序列区域的名称的注释，但在其发布的数据中呈现明显的问题，同时，其中还存在注释错误或模糊、一条数据重复提交的问题(Vilgalys，R. 2003. “Taxonomic misidentification in public DNA databases". New Phytol. 160 :4_5 ；McMahon,Μ. Μ. ,and Μ. J. Sanderson. 2006,Phylog enetic supermatrixanalysis of GenBank sequences from 2228papilionoid legumes". Syst. Biol. 55 :818-836.)。
即使从INSDC拿到的序列，不存在注释错误的问题，但是其测序的质量却不一定符合相关系统发育学研究的需要。如在BARCODE Data Standards v. 2. 3 (26March 2009) 中就建议做为潜在物种条形码的序列是在测序中双向覆盖无N碱基且序列谱图文件的 PHRED scores 不能低于 40%。
所以，需要提供一种方法对现有基因序列数据进行筛选，摈弃注释错误或模糊、测序精度参差的不符合后续挖掘要求的数据。随后，当在已测公开数据中没有找到符合条件的基因序列数据时进行补充测序。发明内容
从上面的分析可以看出，由于历史数据积累的原因，基因序列数据存在注释错误或模糊、测序精度参差不齐等质量问题，继而导致无法正确构建系统发育树的问题。本发明的目的在于提供一种基于注释信息和同源性比对以及特定筛选序列片段相结合的基因序列数据筛选方法。
另外，由于基因序列数据筛选计算属于数据密集型计算，筛选效率问题也是一个需要重点考虑的对象。
因此，本发明的基因序列数据筛选方法首先要解决基因序列数据的质量问题，进一步要提高目标基因序列数据集的筛选效率。本发明的基因序列数据筛选方法也可以作为自测数据的质量控制筛选方法。
本发明的基因序列数据筛选方法，其步骤包括
1)基于基因注释信息的初始数据检索得到数据集，并将其调整为.fasta的格式；
2)针对数据集中的每条序列进行N/R/K/M/S/Y/W/H/D/V/B含量的计算；
3)对数据集中的每条序列进行终止密码子(TAG、TAA、TGA)或其它自定义序列串的检测；
4)将数据集中的每条序列翻译为蛋白序列，将它们与基因对应的模板蛋白序列进行相似性比对计算；
5)根据预设条件，综合步骤2)、3)和4)的结果对每条序列进行评判，决定是否选取。
上述步骤中，步骤2)、3)、4)的执行顺序可以互换或并列进行。
在步骤2)，本发明通过对目标基因序列中N/R/K/M/S/Y/W/H/D/V/B含量的计算， 保证筛选者选取合适测序质量的序列。
本发明按照下面公式计算N/R/K/M/S/Y/W/H/D/V/B中任一种字符的含量γ π …Ni
Pi =-Nall
其中每种字符(N/R/K/M/S/Y/W/H/D/V/B)的含量为Pi，每条序列的字符总数为 Nall，字符 i 的个数为 Ni(i = N，R，K，Μ, S，Y，W, H，D，V 或 B)。字符 N，R，K，Μ, S，Y，W, H， D，V，B代表序列表中不确定的核苷酸残基，其具体含义参见表1。
在步骤3)，本发明通过终止密码子(TAG、TAA、TGA)的检测以排除目标基因序列是假基因序列的可能；通过自定义序列串的检测以排除在各自研究领域内常见的污染序列串或是引物去除不净的序列等不希望出现的序列串。
本发明优选采用正反共6个阅读框(正向3个、反向3个)的方式检测基因序列中是否含有以上终止密码子和自定义序列串。
本发明在步骤4)将待筛选序列与对应的模板蛋白序列进行相似性比对计算，得到一致性值(identity)和期望值(evalue,即Expectation value)，以此作为去除注释错误的非同源序列的依据。
本发明优选采用正反6个阅读框(正向3个、反向3个)的方式将待筛选基因序列翻译为蛋白序列，应用BLAST (Basic Local Alignment Search Tool)算法进行序列相似性比对。
考虑到各具体数据挖掘领域对数据精度要求的不同，在进行具体的基因序列筛选时，对步骤2)和步骤4)的计算结果，可以通过预设阀值来确定筛选范围。例如，对于步骤 2)得到的各字符的含量Pi，设定Pi阈值(例如Pi < 1 % )，在步骤幻淘汰Pi高于其阈值的序列；对于步骤4)得到的比对结果，设定一致性值和期望值的阈值(如Identities > 93%, evalue < 1. 0X 10_1(l)，在步骤5)淘汰一致性和期望值不符合条件的序列。
在步骤幻，针对步骤2、,3)和4)的结果，根据预设的条件对序列进行取舍，通常所选取的序列应该同时满足下述预设条件步骤幻的计算结果Pi符合所设定的阈值要求； 在步骤；3)检测不到终止密码子和不应具有的自定义序列串；步骤4)的计算结果符合设定的同源性要求。满足预设条件的序列被选取，否则被淘汰。
为了进一步确保没有出现漏筛或错筛的情况，在步骤幻之后增加下列步骤
6)将各序列与已确定的同源基因序列进行多重序列比对，验证序列同源与否，并将此同源性结果与步骤4)的结果进行比较。
为了提高筛选效率，本发明可以对步骤1)获得的数据集中的序列进行分组并行处理，也就是说，在步骤1)之后首先依据文件块大小对序列进行分组，然后检测每组文件的头和尾并做相应调整，以保证每一条fasta格式的基因序列被分到一个文件块中；随后对各文件块并行进行步骤幻、3)、4)和幻的处理；最后将各文件块的处理结果汇总。
上述并行处理可选用基于Map/Reduce的并行运算方式。将已经调整为.fasta格式的数据集输入文件按大小切分成文件块；然后检测文件块的头尾并进行调整，使得每个文件块中所包含的文件均为完整的fasta格式的文件；对每个文件块再进行切分产生键值对，将键值对发送到Map计算节点；每个Map节点接收计算信息，进行步骤幻 幻的运算， 并产生结果发送到Reduce节点；Reduce节点接收Map输出的结果信息并组织报告输出。
本发明的基因序列数据的筛选方法克服了现有基因序列数据筛选时存在的注释错误或模糊、测序精度参差不齐等质量问题。本发明的方法首先利用基因的注释信息抽提初始数据集，然后通过逐条对基因序列进行N/R/K/M/S/Y/W/H/D/V/B含量计算、终止密码子(TAG、TAA、TGA)以及污染序列片段等(自定义序列串)的检测、与模板蛋白的相似性计算，最后根据预设条件决定是否选取。进一步的，可将待筛选序列与同源基因序列进行相似性计算验证所筛选出的序列的同源性。本方法可以用于生物系统发育、生物条形码、生物物种鉴定等相关领域的基因数据筛选。
图1是本发明具体实施方式
中基因序列数据筛选方法的工作流程图2是具体实施方式
中针对陆地植物系统发育分析所需MatK基因序列数据进行并行筛选的流程图。
具体实施方式
参见图1，本具体实施方式
所述的基因序列数据筛选方法的具体过程为
A、基于基因注释信息的初始数据检索得到数据集并调整为.fasta的格式，接下来执行步骤B;
B、针对每条序列进行N/R/K/M/S/Y/W/H/D/V/B含量的计算，每种字符的含量为 Pi，设定阀值Pi < 1% (i = N，R，K，M，S，Y，W，H，D，V或B)，接下来执行步骤C;
C、针对每条序列进行终止密码子(TAG、TAA、TGA)和其它自定义序列串 (ACCCAGTCCATCT和GGAAATCTTGGTCC)的检测，接下来执行步骤D ；
D、将每条序列与基因对应的模板蛋白序列进行相似性比对计算，设定阀值 Identities >93%，设定阀值 evalue (Expectation value) < 1. OXKT1 ，其中 Identities 指一致性，evalue指期望值，接下来执行步骤E ；
E、根据预设条件(阀值)，综合B、C、D的结果决定是否选取，并生成每条序列的结果报告，接下来执行步骤F;
F、将每条序列与已知的同源基因序列进行多重序列比对，验证序列同源与否。
本具体实施方式
中的初始数据检索通过调用NCBI的API得到，我们的检索词是(matk[Gene NameJAND “ Embryophyta “ [Organism])AND “ ddbj embl genbank" [Filter]，得到相关的数据集，并保存为fasta的格式。
本具体实施方式
中所述步骤B的具体过程为针对每条基因序列中的N/R/K/M/S/ Y/ff/H/D/V/Β字符数进行计数，然后与序列长度即此条序列的字符总数相除，得到各字符的含量Pi ；
本具体实施方式
中所述步骤C的具体过程为每条基因序列分为6个阅读框(正向 3 个、反向 3 个)分别匹配检索字符串 TAG、TAA、TGA, ACCCAGTCCATCT、GGAAATCTTGGTCC。
本具体实施方式
中所述步骤D的具体过程为首先每条待筛选基因序列分为6 个阅读框先翻译成蛋白质序列应用EMBOSS平台(http://emboss, sourceforge.net/) 的transeq实现6个阅读框的翻译，其密码子编码使用标准密码子(Mandard)；然后分别选取6个阅读框中蛋白质序列延伸最长的序列作为BLAST算法的输入序列；应用 NCBIncbi-blast-2. 2. 22 版本(ftp://ftp. ncbi. nlm. nih. gov/blast/)运行 blastp,进行待筛选基因序列的翻译蛋白和已测蛋白模板LIB之间的序列比对，设置阀值为Identities > 93% > evalue (Expectation value) < L OX ICT100
本具体实施方式
中的蛋白模板选取来自Swiss-Prot的此种基因对应的已测序蛋白序歹丨J (DBS0URCE UniProtKB =Iocus MATK_ANTLI, accession Q7YJG1 ；)，然后格式转化为 BLAST算法需要的LIB。
本具体实施方式
中所述步骤E的具体过程为根据预设阀值，综合B U C U D的结果决定是否选取，其条件为(1)目标序列Pi < 0.01(i = N，R，K，M，S，Y，W，H，D，V，B)； 并且⑵待筛选基因序列6个阅读框检测不含有“TAG”、“TAA”、“TGA”、“ACCCAGTCCATCT”、 “GGAAATCTTGGTCC” 任一字符串；并且(3) blastp 设定的阀值 Identities >93% ；evalue < 1.0X10,。同时生成针对每条待筛选基因序列的报告，其报告内容定义如下表所示。
表 1
name每条基因序列的名字，如来自公共序列筛选为序列的“accession”Length基因序列长度score序列和模板蛋白的比对打分identities序列和模板蛋白的一致性(百分数表示)> 93%evalue序列和模板蛋白期望值< 1. OX1CT10Pn“N”字符在序列中的百分含量，N代表A/G/C/T任一碱基■’< 1%Pr"R"字符在序列中的百分含量，R代表A/G任一碱基；< 1%Pk“K”字符在序列中的百分含量，K代表G/T任一碱基；< 1%Pm“M”字符在序列中的百分含量，M代表A/C任一碱基；< 1%Ps“S”字符在序列中的百分含量，S代表G/C任一碱基；< 1%Py“Y”字符在序列中的百分含量，Y代表c/τ任一碱基；< 1%Pw“W”字符在序列中的百分含量，W代表A/T任一碱基；< 1%Ph“H”字符在序列中的百分含量，H代表A/C/T任一碱基；< 1%Pd“D”字符在序列中的百分含量，D代表A/G/T任一碱基；< 1%
权利要求
1.一种基因序列数据的筛选方法，包括如下步骤1)基于基因注释信息的初始数据检索得到数据集，并将其调整为.fasta的格式；2)针对数据集中的每条序列进行N/R/K/M/S/Y/W/H/D/V/B各字符含量的计算；3)对数据集中的每条序列进行终止密码子或自定义序列串的检测；4)将数据集中的每条序列翻译为蛋白序列，将它们与基因对应的模板蛋白序列进行相似性比对计算；5)根据预设条件，综合步骤2)、3)和4)的结果对每条序列进行评判，决定是否选取。
2.如权利要求1所述的筛选方法，其特征在于，步骤2)用Nall代表序列的字符总数， Ni代表该序列中字符i的个数，其中i指N，R，K，M，S，Y，W，H，D，V或B，则字符i的含量PiNi为= 步骤5)淘汰Pi高于预设阈值的序列。 Nall
3.如权利要求2所述的筛选方法，其特征在于，步骤5)所选取的序列各字符含量Pi < 1%，其中 i 指 N，R，K，M，S，Y，W，H，D，V 或 B。
4.如权利要求1所述的筛选方法，其特征在于，步骤3)所述终止密码子是指TAG、TAA 和TGA，采用正向3个、反向3个共6个阅读框的方式检测序列中是否含有终止密码子和自定义序列串；步骤幻淘汰含有终止密码子和不希望出现的自定义序列串的序列。
5.如权利要求1所述的筛选方法，其特征在于，步骤4)采用正向3个、反向3个共6个阅读框的方式将待筛选基因序列翻译为蛋白序列，再应用BLAST算法将其与对应的模板蛋白序列进行序列相似性比对，获得序列的一致性值和期望值；步骤幻淘汰一致性值和期望值不符合预设条件的序列。
6.如权利要求5所述的筛选方法，其特征在于，步骤幻所选取的序列一致性>93%， 期望值< 1.0X10,。
7.如权利要求1所述的筛选方法，其特征在于，在步骤5)之后增加下列步骤6)将各序列与已确定的同源基因序列进行多重序列比对，验证序列同源与否，并将此同源性结果与步骤4)的结果进行比较。
8.如权利要求1所述的筛选方法，其特征在于，对步骤1)获得的数据集中的序列文件进行分组，然后各组并列进行步骤幻 5)的处理，最后将处理结果汇总。
9.如权利要求8所述的筛选方法，其特征在于，采用基于Map/Reduce的并行运算方式对数据集中的序列进行并行处理。
10.如权利要求9所述的筛选方法，其特征在于，将已经调整为.fasta格式的数据集输入文件按大小切分成文件块；然后检测文件块的头尾并进行调整，使得每个文件块中所包含的文件均为完整的fasta格式的文件；对每个文件块再进行切分产生键值对，将键值对发送到Map计算节点；每个Map节点接收计算信息，进行步骤2) 5)的运算，并产生结果发送到Reduce节点；Reduce节点接收Map输出的结果信息并组织报告输出。
全文摘要
本发明公开了一种基因序列数据的筛选方法。首先利用基因的注释信息抽提初始数据集，然后逐条对基因序列进行N/R/K/M/S/Y/W/H/D/V/B的含量计算、终止密码子和自定义序列串(如污染序列片段)的检测、与模板蛋白的相似性计算，根据预设条件决定是否选取。该方法克服了现有基因序列数据筛选时存在的注释错误或模糊、测序精度参差不齐等质量问题继而导致无法正确构建系统发育树的问题，可以用于生物系统发育、生物条形码、生物物种鉴定等相关领域的基因数据筛选。
文档编号G06F19/22GK102521528SQ20111040012
公开日2012年6月27日 申请日期2011年12月5日 优先权日2011年12月5日
发明者周园春, 孟珍, 黎建辉 申请人:中国科学院计算机网络信息中心