專利名稱：一種基因序列數據的篩選方法
一種基因序列數據的篩選方法技術領域
本發明屬于應用生物信息學技術領域，尤其涉及一種基因序列數據篩選方法，主要應用于生物系統發育、生物條形碼、生物物種鑒定等相關領域的基因數據篩選和質量控制。
背景技術：
生物分子系統發育研究在不同水平和層次上依賴于對數據的使用從簡單的檢索到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日
發明者周園春, 孟珍, 黎建輝 申請人:中國科學院計算機網絡信息中心