CN102521528A

CN102521528A - 一种基因序列数据的筛选方法

Info

Publication number: CN102521528A
Application number: CN2011104001226A
Authority: CN
Inventors: 黎建辉; 孟珍; 周园春
Original assignee: Computer Network Information Center of CAS
Current assignee: Computer Network Information Center of CAS
Priority date: 2011-12-05
Filing date: 2011-12-05
Publication date: 2012-06-27

Abstract

本发明公开了一种基因序列数据的筛选方法。首先利用基因的注释信息抽提初始数据集，然后逐条对基因序列进行N/R/K/M/S/Y/W/H/D/V/B的含量计算、终止密码子和自定义序列串(如污染序列片段)的检测、与模板蛋白的相似性计算，根据预设条件决定是否选取。该方法克服了现有基因序列数据筛选时存在的注释错误或模糊、测序精度参差不齐等质量问题继而导致无法正确构建系统发育树的问题，可以用于生物系统发育、生物条形码、生物物种鉴定等相关领域的基因数据筛选。

Description

一种基因序列数据的筛选方法

技术领域

本发明属于应用生物信息学技术领域，尤其涉及一种基因序列数据筛选方法，主要应用于生物系统发育、生物条形码、生物物种鉴定等相关领域的基因数据筛选和质量控制。

背景技术

生物分子系统发育研究在不同水平和层次上依赖于对数据的使用：从简单的检索到PCR污染物的检查，到寻找一个给定序列的类群同源性序列，到更全面的基于大量数据进行的类群和位点的数据挖掘(McMahon，M.M.，and M.J.Sanderson.2006.”Phylogeneticsupermatrix analysis of GenBank sequences from 2228papilionoid legumes”.Syst.Biol.55：818-836；Ciccarelli，F.D.，T.Doerks，C.von Mering，C.J.Creevey，B.Snel，and P.Bork.2006.”Toward automatic reconstruction of a highly resolved tree of life.Science 311：1283-1287；Bininda-Emonds，O.R.P.，M.Cardillo，K.E.Jones，R.D.E.MacPhee，R.M.D.Beck，R.Grenyer， S.A.Price，R.A.Vos，J.L.Gittleman，and A.Purvis.2007.”The delayed rise ofpresent-day mammals”.Nature 446：507-512；Li，C.H.，G.Orti，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 BOSS，et al.2008.“ThePhyLoTA Browser：Processing GenBank for Molecular Phylogenetics Research”，Syst.Biol.57(3)：335-346.)。

分子生物学的早期研究积累了大量的基因序列数据。以国际核算序列数据库联盟(International Nucleotide Sequence Database Collaboration，INSDC)成员之一的GenBank为例(Michael Y.Galperin.2011.“The Molecular Biology Database Collection：2011updae”.Nucl.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，M.M.，and M.J.Sanderson.2006.“Phylogenetic 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中任一种字符的含量：

Pi = \frac{Ni}{Nall}

其中每种字符(N/R/K/M/S/Y/W/H/D/V/B)的含量为Pi，每条序列的字符总数为Nall，字符i的个数为Ni(i＝N，R，K，M，S，Y，W，H，D，V或B)。字符N，R，K，M，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％)，在步骤5)淘汰Pi高于其阈值的序列；对于步骤4)得到的比对结果，设定一致性值和期望值的阈值(如Identities＞93％，evalue＜1.0×10^-10)，在步骤5)淘汰一致性和期望值不符合条件的序列。

在步骤5)，针对步骤2)、3)和4)的结果，根据预设的条件对序列进行取舍，通常所选取的序列应该同时满足下述预设条件：步骤2)的计算结果Pi符合所设定的阈值要求；在步骤3)检测不到终止密码子和不应具有的自定义序列串；步骤4)的计算结果符合设定的同源性要求。满足预设条件的序列被选取，否则被淘汰。

为了进一步确保没有出现漏筛或错筛的情况，在步骤5)之后增加下列步骤：

6)将各序列与已确定的同源基因序列进行多重序列比对，验证序列同源与否，并将此同源性结果与步骤4)的结果进行比较。

为了提高筛选效率，本发明可以对步骤1)获得的数据集中的序列进行分组并行处理，也就是说，在步骤1)之后首先依据文件块大小对序列进行分组，然后检测每组文件的头和尾并做相应调整，以保证每一条fasta格式的基因序列被分到一个文件块中；随后对各文件块并行进行步骤2)、3)、4)和5)的处理；最后将各文件块的处理结果汇总。

上述并行处理可选用基于Map/Reduce的并行运算方式。将已经调整为.fasta格式的数据集输入文件按大小切分成文件块；然后检测文件块的头尾并进行调整，使得每个文件块中所包含的文件均为完整的fasta格式的文件；对每个文件块再进行切分产生键值对，将键值对发送到Map计算节点；每个Map节点接收计算信息，进行步骤2)～5)的运算，并产生结果发送到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.0×10^-10，其中Identities指一致性，evalue指期望值，接下来执行步骤E；

E、根据预设条件(阀值)，综合B、C、D的结果决定是否选取，并生成每条序列的结果报告，接下来执行步骤F；

F、将每条序列与已知的同源基因序列进行多重序列比对，验证序列同源与否。

本具体实施方式中的初始数据检索通过调用NCBI的API得到，我们的检索词是：(matk[Gene Name]AND″Embryophyta″[Organism])AND″ddbj embl genbank″[Filter]，得到相关的数据集，并保存为fasta的格式。

本具体实施方式中所述步骤B的具体过程为：针对每条基因序列中的N/R/K/M/S/Y/W/H/D/V/B字符数进行计数，然后与序列长度即此条序列的字符总数相除，得到各字符的含量Pi；

本具体实施方式中所述步骤C的具体过程为：每条基因序列分为6个阅读框(正向3个、反向3个)分别匹配检索字符串TAG、TAA、TGA、ACCCAGTCCATCT、GGAAATCTTGGTCC。

本具体实施方式中所述步骤D的具体过程为：首先每条待筛选基因序列分为6个阅读框先翻译成蛋白质序列：应用EMBOSS平台(http://emboss.sourceforge.net/)的transeq实现6个阅读框的翻译，其密码子编码使用标准密码子(Standard)；然后分别选取6个阅读框中蛋白质序列延伸最长的序列作为BLAST算法的输入序列；应用NCBIncbi-blast-2.2.22版本(ftp://ftp.ncbi.nlm.nih.gov/blast/)运行blastp，进行待筛选基因序列的翻译蛋白和已测蛋白模板LIB之间的序列比对，设置阀值为Identities＞93％、evalue(Expectation value)＜1.0×10^-10。

本具体实施方式中的蛋白模板选取来自Swiss-Prot的此种基因对应的已测序蛋白序列(DBSOURCE UniProtKB：locus MATK_ANTLI，accession Q7YJG1；)，然后格式转化为BLAST算法需要的LIB。

本具体实施方式中所述步骤E的具体过程为：根据预设阀值，综合B∪C∪D的结果决定是否选取，其条件为：(1)目标序列Pi＜0.01(i＝N，R，K，M，S，Y，W，H，D，V，B)；并且(2)待筛选基因序列6个阅读框检测不含有“TAG”、“TAA”、“TGA”、“ACCCAGTCCATCT”、“GGAAATCTTGGTCC”任一字符串；并且(3)blastp设定的阀值Identities＞93％；evalue＜1.0×10^-10。同时生成针对每条待筛选基因序列的报告，其报告内容定义如下表所示。

表1

报告项	说明	阀值
			name	每条基因序列的名字，如来自公共序列筛选为序列的“accession”	-
Length	基因序列长度	-
			score	序列和模板蛋白的比对打分	-
identities	序列和模板蛋白的一致性(百分数表示)	＞93％
			evalue	序列和模板蛋白期望值	＜1.0×10^-10
P_N	“N”字符在序列中的百分含量，N代表A/G/C/T任一碱基；	＜1％
			P_R	“R”字符在序列中的百分含量，R代表A/G任一碱基；	＜1％
P_K	“K”字符在序列中的百分含量，K代表G/T任一碱基；	＜1％
			P_M	“M”字符在序列中的百分含量，M代表A/C任一碱基；	＜1％
P_S	“S”字符在序列中的百分含量，S代表G/C任一碱基；	＜1％
			P_Y	“Y”字符在序列中的百分含量，Y代表C/T任一碱基；	＜1％
P_W	“W”字符在序列中的百分含量，W代表A/T任一碱基；	＜1％
			P_H	“H”字符在序列中的百分含量，H代表A/C/T任一碱基；	＜1％
P_D	“D”字符在序列中的百分含量，D代表A/G/T任一碱基；	＜1％
			P_V	“V”字符在序列中的百分含量，V代表A/G/C任一碱基；	＜1％
P_B	“B”字符在序列中的百分含量，B代表G/C/T任一碱基；	＜1％
			TAG	终止密码子“TAG”个数，为整数；	＜1
TAA	终止密码子“TAA”个数，为整数；	＜1
			TGA	终止密码于“TGA”个数，为整数；	＜1
SC1	自定义检测串1“ACCCAGTCCATCT”个数，为整数；	＜1
			SC2	自定检测串2“GGAAATCTTGGTCC”个数，为整数；	＜1
Result	标记T(选取)或F(淘汰)；	-

本具体实施方式中所述步骤F的具体过程为：针对每条基因序列和已验证的同源序列簇，应用clustalx1.83(J.D.Thompson，D.G.Higgins，T.J.Gibson.1994.“CLUSTAL W：improving the sensitivity of progressive multiple sequence alignment through sequence weighting，position-specific gap penalties and weight matrix choice”.Nucleic Acids Research.22：4673.)，选取Blosum80矩阵，进行多重序列比对；然后应用Jalview Java Alignment Editor(M.Clamp，J.Cuff，S.M.Searle，G.J.Barton.2004.“The Jalview Java alignment editor”.Bioinformatics.20：426-427.)进行序列排列状况的查看，应用Archaeopteryx 0.954beta(Han MV and Zmasek CM.2009.″phyloXML：XML for evolutionary biology and comparative genomics″.BMCBioinformatics.10：356.)进行序列聚类情况的查看，对判断结果进行序列同源与否验证。

本实施方式中每条基因序列(T或F)均可以和已验证的同源基因序列簇(即已验证的多条同源基因序列)进行多重序列比对。

为提高计算效率，本实施方式引入了并行计算的概念，将含有待筛选序列数据集的文件进行分组(数据切分)，并行处理待筛选序列的筛选计算。

本实施方式进一步引入了基于Map/Reduce的并行运算方式。Map/Reduce(Jeffrey Deanand Sanjay Ghemawat.2004.“Map/Reduce：Simplified Data Processing on Large Clusters”.Proceedings of USENIX OSDI’04.137-150.)是Google提出的一个软件架构，用于大规模数据集的并行运算。

基于Map/Reduce计算模型、采用Hadoop(APACHE[HADOOP].http://hadoop.apache.org/)框架、通过切分待筛选序列文件的方法来实现基因序列数据集的快速筛选。

本实施方式对待筛选序列数据集文件(.fasta)进行分组，然后进行步骤B、C、D、E的筛选计算，整个的流程如图2所示。对基因序列数据集进行筛选时，首先把输入文件按大小切分成文件块，检测文件块的头尾并进行调整，使得每个文件块中所包含的文件均为完整的fasta格式的文件；对每个文件块再进行切分产生键值对，将键值对发送到Map计算节点；每个Map节点接收计算信息，进行(B、C、D、E)步骤的运算，并产生结果发送到Reduce节点；Reduce节点接收Map输出的结果信息并组织报告输出。

本实施方式的进行筛选计算的具体流程如图2所示，

1.任务及数据分割：在每个分组当中的JobTracker上，读取输入文件由FileInputFormat按一定份数切分成大小相同的文件块(FileSplits)，生成同样数目的子任务，每个子任务对应处理一个文件块；然后，将该子任务的配置信息(如处理的文件块信息)发送到不同的节点去运算。每个FileSplits均要检测头尾是否是fasta的格式，如图2一个fasta格式的序列数据被分成了两部分，分别位于part-000的尾部和part-001的头部，因此在生成键值对<Key，Value>之前要先作调整，保证一条序列数据在一个part中，然后再对每个文件块切分产生键值对。

2.读取输入文件，执行Map运算：计算节点接受由JobTracker发送的子任务，读取相应的文件块进行处理；读取文件块的操作由RecordReader来完成，它将文件块组织成<Key，Value>的形式提交给Map任务执行，Map任务将这些键值对与该节点的(B、C、D、E)步骤的运算。

3.Reduce：当有部分节点的Map运算执行完毕之后，即由部分计算节点对结果进行Reduce运算，将所有计算结果汇总到输出文件返回。

Claims

1.一种基因序列数据的筛选方法，包括如下步骤：

2)针对数据集中的每条序列进行N/R/K/M/S/Y/W/H/D/V/B各字符含量的计算；

3)对数据集中的每条序列进行终止密码子或自定义序列串的检测；

2.如权利要求1所述的筛选方法，其特征在于，步骤2)用Nall代表序列的字符总数，Ni代表该序列中字符i的个数，其中i指N，R，K，M，S，Y，W，H，D，V或B，则字符i的含量Pi为：步骤5)淘汰Pi高于预设阈值的序列。

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)淘汰含有终止密码子和不希望出现的自定义序列串的序列。

5.如权利要求1所述的筛选方法，其特征在于，步骤4)采用正向3个、反向3个共6个阅读框的方式将待筛选基因序列翻译为蛋白序列，再应用BLAST算法将其与对应的模板蛋白序列进行序列相似性比对，获得序列的一致性值和期望值；步骤5)淘汰一致性值和期望值不符合预设条件的序列。

6.如权利要求5所述的筛选方法，其特征在于，步骤5)所选取的序列一致性＞93％，期望值＜1.0×10^-10。

7.如权利要求1所述的筛选方法，其特征在于，在步骤5)之后增加下列步骤：

8.如权利要求1所述的筛选方法，其特征在于，对步骤1)获得的数据集中的序列文件进行分组，然后各组并列进行步骤2)～5)的处理，最后将处理结果汇总。

9.如权利要求8所述的筛选方法，其特征在于，采用基于Map/Reduce的并行运算方式对数据集中的序列进行并行处理。

10.如权利要求9所述的筛选方法，其特征在于，将已经调整为.fasta格式的数据集输入文件按大小切分成文件块；然后检测文件块的头尾并进行调整，使得每个文件块中所包含的文件均为完整的fasta格式的文件；对每个文件块再进行切分产生键值对，将键值对发送到Map计算节点；每个Map节点接收计算信息，进行步骤2)～5)的运算，并产生结果发送到Reduce节点；Reduce节点接收Map输出的结果信息并组织报告输出。