基于多參考序列的基因序列分級壓縮方法

基于多參考序列的基因序列分級壓縮方法
【專利摘要】本發明提出了一種基于多參考序列的基因序列分級壓縮方法，該方法首先將BAM格式文件轉化成SAM格式的文件，SAM格式的基因序列由11個強制域和多個可選域構成，將原文件按域提取成12個獨立文件，然后對12個文件進行并行壓縮：（1）對‘Sequence’域采用基于多個參考序列逐步減半序列長度的分級壓縮方法；（2）對于‘Quality?Value’域采用k均值聚類結合上下文建模PPMVC壓縮的方法；該方案相對于現有的同格式的壓縮方案既提高或保證了壓縮效率，又提供了壓縮等級的多選擇性，使得其更有適應性與擴展性。
【專利說明】基于多參考序列的基因序列分級壓縮方法
【技術領域】
[0001]本發明涉及一種面向超大規模SAM/BAM (兩種序列比對格式標準)格式基因序列的信息壓縮方法，具體是一種利用目標序列與相同物種的不同參考序列間的相似重復性，多次改變參考序列并減少短序列長度的分級基因壓縮方法。
【背景技術】
[0002]DNA是生物生存、延續和發展的重要物質基礎，具有重大的科學價值和社會價值。目前，DNA的研究廣泛應用于生物學、醫學、遺傳科學等諸多重要領域，如通過收集和保存DNA信息以保護瀕臨滅絕的生物物種、基于人類基因序列的信息預測以及找到基因變異規律以治療癌癥腫瘤等。為這些學科研究提供基礎實驗數據的各種DNA序列測定工程已成為各國重點發展的研究項目。隨著這些測序項目的展開，每天都有海量的DNA序列數據產生，相關數據量呈指數方式增長，生物信息數據這種急速的積累增長在人類的科學研究歷史中是空前的。例如，公共測序序列的寄存器“序列檔案”(SRA, the Sequence Read Archive)的儲存量截至2013年年底將預計達到1000TB。存儲和使用這些數據的成本已越來越面臨著無法承擔的規模，如何在有限的存儲資源內有效儲存急劇膨脹的DNA序列數據成為了計算機專家和生物學家面臨的新課題，也是國內外諸多重大計劃所面臨的前進障礙。因此，采用更有效的壓縮編碼方式，用較小的存儲空間存放較大的基因信息序列是必然的選擇。
[0003]迄今為止，所有的基因壓縮研究主要圍繞三種序列格式展開:(I)最初的DNA基集合形式的FASTA格式(一種基于文本用于表不核苷酸序列或氨基酸序列的格式)；(2)未比對短序列形式的FASTQ格式(一種基于文本用于表示核苷酸序列或氨基酸序列和相應質量信息的格式)；(3)比對后短序列形式的SAM/BAM格式。由于序列比對是進行序列分析和處理的首要步驟，且SAM/BAM格式包含了豐富完整的基因比對信息，因此近幾年的基因分析與壓縮研究均聚焦于SAM (由于BAM是SAM的二進制表示，因此壓縮BAM格式時可先將其解壓，然后還原為壓縮SAM格式的序列)格式。2011年Muhammad等人在PLOS One期刊的 “Improving transmission efficiency of large sequence alignment/map (SAM)files”上提出了專門針對于SAM格式和特征的壓縮方法SAMZIP，利用基本壓縮技巧，如游程編碼(Run-Length Encoding, RLE),差分(Delta)編碼,哈弗曼(Huffman)編碼和字典編碼，對SAM格式中的每一列獨立處理。同年,Kozanitis等人在Journal of ComputationalBiology 期干丨J的“Compressing genomic sequence fragments using SlimGene，，上提出了壓縮方法SLIMGENE，其中首次嘗試了對SAM格式中的Quality Value (質量值)項進行保證后續處理的有損壓縮。考慮到相同物種核苷酸差異性的微小不同(人類核苷酸差異性僅為0.1%左右)，研究學者們開始將已知的參考序列引入SAM格式的基因壓縮中。2011年，Fritz 等人在 Genome research 期干丨J 的 “Efficient storage of high throughputDNA sequencing datausing reference-based compression，，上發表了一種基于參考序列的壓縮方法CRAM，將SAM格式序列中的每一子序列(Read)與參考基因做比對，然后壓縮其比對結果，用較少比特表示相應的差別。2012年Hach等人在Bioinformatics期刊的“SCALCE:boosting sequence compression algorithms using locally consistentencoding”上提出了另一種基于參考序列的壓縮方法SCALCE，該方法基于短序列的重組以適應參考序列的局部特性，便于進一步比對壓縮。以上所述的這些方法將SAM格式分離成多個獨立項并行處理。在壓縮DNA序列時將其視作由特殊字符(即‘A’，‘C’，‘T’，‘G’，‘N’)構成的長字符串，從數據的構成特點和自身冗余性出發進行整體處理，有效的提高了壓縮效率和壓縮時間。但總體而言SAM格式的大基因壓縮技術仍處于起步階段，已知公開的參考基因的信息并沒有被充分發揮利用，目標信息中的準確比對率也仍待提高。

【發明內容】

[0004]針對現有技術中的缺陷，本發明的目的是提供一種更加有效的基于多參考序列的SAM/BAM格式基因序列分級壓縮方法。
[0005]本發明是通過以下技術方案實現的:該方法通過利用多個公開的參考基因序列，并將短序列長度逐步減半，多次比對目標序列以提聞被壓縮序列的比對準確率，進而提聞壓縮效率。另外，對于SAM格式中的Quality Value項,本發明采用了用戶可指定壓縮等級的k均值(k-means)聚類后再壓縮的策略，提高了方法的廣泛實用性。由于BAM格式是SAM格式的簡單轉化，因此下面僅討論針對SAM格式的壓縮技術，對于BAM格式，只需將其轉換成SAM格式然后再按照SAM格式的壓縮方法處理即可。
[0006]本發明所述的基因序列分級壓縮方法，該方法首先將BAM格式文件轉化成SAM格式的文件，SAM格式的基因序列是由比對工具輸出的文本形式文件，它由11個強制域(Field)和一系列可選域(看作第12個域)構成，因此本發明首先將原文件按域提取成12個獨立文件，然后對12個文件進行并行壓縮:
[0007](I)對‘Sequence’域采用基于多個參考序列逐步減半序列長度的分級壓縮方法；
[0008](2)對于‘Quality Value’域采用k均值聚類結合上下文建模PPMVC (部分串匹配的預測(Prediction with Partial string Matching)法的一種擴展)壓縮的方法；
[0009](3)對于剩下的十個域采用基于域內特征和域間相關性的壓縮方法。
[0010]進一步的，所述對‘Sequence’域采用基于多個參考序列逐步減半序列長度的分級壓縮方法，具體為:利用快速比對工具S0AP3將SAM/BAM格式基因序列文件的‘Sequence’域中的短序列分線程地與參考序列作比對，對于準確匹配序列高效壓縮，對于非準確匹配和未匹配的短序列，將其序列長度減半，即一個序列分成長度相同的兩個序列，并改變參考序列，再進行第二次比對，得到比對結果，如此重復三至四次結束，剩余的非準確匹配和未匹配的短序列進行PPMVC編碼即可。這樣多次比對后需要處理壓縮的非準確匹配和未匹配序列變少。
[0011]優選的，所述對于準確匹配序列高效壓縮，該壓縮方法基于比對結果的特點，富有針對性。具體為:對于準確比對的子序列(Read)，使用〈Read編號，參考序列上重復發生的染色體號，參考序列上重復發生的偏移位置，重復類型 > 這四個量來替代目標序列上重復的子序列(Read),分別使用差分編碼+Huffman編碼、游程編碼、差分編碼+Huffman編碼和游程編碼來壓縮這四個分量。
[0012]進一步的，所述的對于‘Quality Value’域采用k均值聚類結合上下文建模PPMVC壓縮的方法，具體為:采用k均值聚類法將η個Qashi值(代表對應基的比對質量分值的ASCII碼值)聚成k類，使得聚類后每類內所有Quality value的值與聚類前的值差值平方最小，然后采用基于上下文建模和統計的自適應文本壓縮方法PPMVC壓縮聚類后的‘Quality Value’文件。這樣同等壓縮條件下失真較小。
[0013]進一步的，所述的對于剩下的十個域采用基于域內特征和域間相關性的壓縮方法，具體為:
[0014]對于‘QNAME’域，用‘0’表示之前未出現過的QNAME，用逐漸遞增的數字編號與當前位置只差表示之前已經出現的QNAME，然后采用Huffman編碼壓縮這些非均勻分布的小型數值；
[0015]對于‘FLAG’域，用一個字節表示I?255之間的數值，用三個字節(即，0，x/256andx%256)表示其它數值，然后采用Huffman編碼壓縮變換后的數值；
[0016]對于‘RNAME’域，用相同的數字標記整個SAM文件中的相同的參考序列名字，記錄下來所有參考序列，然后用游程編碼進行壓縮；
[0017]對于‘P0S’域，采用差分編碼+Huffman編碼；
[0018]對于‘MAPQ’域,采用游程編碼；
[0019]對于‘CIGAR’ 域，采用 Lempel-Ziv-Welch LZW 字典壓縮方法；
[0020]對于‘MRNM’域，采用游程編碼；
[0021]對于‘MP0S’域，結合‘MRNM’域的字符串，采用差分編碼+Huffman編碼；
[0022]對于‘TLEN’域，‘TLEN’域的值與‘MP0S’域減去‘P0S’域的值的差(即，TLEN-(MPOS-POS))的絕對值服從于一個有限的字符集，因此，對于該域的壓縮，本發明結合‘POS’，‘MP0S’和‘MRNM’三個域的信息采用Huffman編碼壓縮變換后的值；
[0023]對于‘OPTIONAL’域，使用bzip2壓縮工具。
[0024]與現有技術相比，本發明的有益效果是:
[0025]本發明所提出的基于多參考序列的SAM/BAM格式基因序列分級壓縮方法，提高了基因壓縮的效率和完整性。本發明將多個公開的基因序列作為參考結合使用，充分利用了相同物種間的基因相似性；將非準確匹配的子序列進行截斷式再次比對，克服了以往方法無針對性地對準確匹配和非準確匹配序列統一處理的缺點，提高了準確比對率，進而節省了壓縮比特；對‘Quality Value’域采用用戶可指定壓縮級的k均值聚類法，在提高壓縮效率的同時保證了質量值的準確度；對剩余域既考慮了域間獨立性也考慮了域內特征和部分域間的分布相關性，針對性的進行變換然后再編碼，挖掘出了 SAM序列格式中的潛在信息，進一步提高了壓縮效率。
【專利附圖】

【附圖說明】
[0026]通過閱讀參照以下附圖對非限制性實施例所作的詳細描述，本發明的其它特征、目的和優點將會變得更明顯:
[0027]圖1是本發明實施例的壓縮流程圖；
[0028]圖2是本發明實施例的‘Sequence’域經過不同參考序列不同子序列長度多次比對的不意圖；
[0029]圖3是本發明實施例的‘Sequence’域在某一實例中的壓縮效果圖。【具體實施方式】
[0030]下面結合具體實施例對本發明進行詳細說明。以下實施例將有助于本領域的技術人員進一步理解本發明，但不以任何形式限制本發明。應當指出的是，對本領域的普通技術人員來說，在不脫離本發明構思的前提下，還可以做出若干變形和改進。這些都屬于本發明的保護范圍。
[0031]本發明首先將BAM格式文件轉化成SAM格式的文件，由于SAM格式的基因序列是由比對工具輸出的文本形式文件，它由11個強制域(Field)和一系列可選域(看作第12個域)構成，因此本發明首先將原文件按域提取成12個獨立文件(這12個獨立文件對應的是上述11個強制域和第12個域，本實施例中即以下所描述的各個域)，然后對12個文件進行并行壓縮。對‘Sequence’域的處理為本發明的第一部分；對‘Quality Value’域的處理為本發明的第二部分；對剩余10個域的壓縮為本發明的第三部分。其中，‘Sequence’域和‘Quality Value’域占據了整個基因序列的50%左右內容且不易壓縮，是本發明的處理重點。每個部分的編碼經過如下:
[0032]1、‘Sequence’域的多參考分級式壓縮方案
[0033]對由‘A’，‘C’，‘T’，‘G’，‘N，五個基構成的‘Sequence’域，將經過以下編碼過程:
[0034](I)通過快速比對工具S0AP3將每一 Read分別與指定的參考序列進行比對；
[0035](2)對于準確比對的子序列(Read)，使用〈Read編號，參考序列上重復發生的染色體號，參考序列上重復發生的偏移位置，重復類型〉這四個量來替代目標序列上重復的子序列(Read)；
[0036](3)分別使用差分編碼+哈弗曼編碼、游程編碼、差分編碼+哈弗曼編碼和游程編碼壓縮〈Read編號，參考序列上重復發生的染色體號，參考序列上重復發生的偏移位置，重復類型 > 中的四個分量；
[0037](4)對于未準確比對或未匹配的子序列(Read)，將其長度減半，并選擇其他參考序列(可選)作為參考，然后重復(I) - (4)，共重復三次。
[0038]2、‘Quality Value’域的用戶可指定編碼級的k均值聚類法壓縮方案
[0039]對由多個(一般情況下有51個)不同字符構成的‘Quality Value’域，每個字符的ASCII值與當前位置所對應的基的誤差概率的關系為:



Qasc 丨丨一33
[0040]Qascii = Qphrcd + 33, Pe = 10 10 > Qphrcd = -1Olog10Pc
[0041]其中Pe是當前位置基的測序錯誤率,Qphred是其轉化為整數后的值，Qascti是可以在文本中顯示出來的質量值。
[0042]采用k均值聚類法將η個Qascq值聚成k類,使得聚類后每類內所有Quality value的值與聚類前的值差值平方最小，即
iq 廠 33Uj-33


y 10 ?ο" , nurrij

i = l qjESj
[0044]以保證聚類后的誤差概率Pe的變化最小。其中Ui是類別Si的均值，numj是整個序列中q」出現的次數。
[0045]最后，采用基于上下文建模和統計的自適應文本壓縮方法PPMVC壓縮聚類后的‘Quality Value，文件。
[0046]3、剩余域的基于域內特征和域間聯系的壓縮方案
[0047]根據剩余域域內和域間的數值關聯及特征性，分別對每一域采用如下編碼方法:
[0048]QNAME:用‘0’表示之前未出現過的QNAME，用逐漸遞增的數字編號與當前位置只差表示之前已經出現的QNAME，然后采用哈弗曼編碼壓縮這些非均勻分布的小型數值。
[0049]FLAG:用一個字節表示I?255之間的數值，用三個字節(即，0，x/256and x%256)表示其它數值，然后采用哈弗曼編碼壓縮變換后的數值。
[0050]RNAME:用相同的數字標記整個SAM文件中的相同的參考序列名字，記錄下來所有參考序列，然后用游程編碼進行壓縮。
[0051]P0S:差分編碼+哈弗曼編碼。
[0052]MAPQ:游程編碼。
[0053]CIGAR:Lempel-Ziv-Welch(LZW)字典壓縮方法。
[0054]MRNM:游程編碼。
[0055]MPOS:結合‘MRNM’域的字符串，采用差分編碼+哈弗曼編碼。
[0056]TLEN:因為(TLEN-(MPOS-POS))的絕對值服從于一個有限的字符集，因此，對于該域的壓縮，結合‘POS’，‘MPOS’and ‘MRNM’三個域的信息采用哈弗曼編碼壓縮變換后的值。
[0057]0PT10NAL:bzip2 壓縮工具。
[0058]解壓縮時，首先按照上述同樣的步驟并行恢復出原基因序列中每一域的部分，然后將所有域合并得到完整無損的SAM或BAM基因文件。
[0059]以下提供本發明的一個具體應用實例:
[0060]如圖1所示，本實施例的壓縮過程包括如下步驟:
[0061]步驟一，利用samtools工具將SAM/BAM文件提取成12個獨立文件，其中每個域為一單獨文件，所有的可選域(即Optional Fields)提取成一個文件,然后送入每個域單獨的編碼器；
[0062]步驟二，對有規律性的十個域采用上述‘剩余域的基于域內特征和域間聯系的壓縮方案’中描述的方法進行壓縮；
[0063]步驟三，對重點的‘Sequence’域采用圖1所示的分級壓縮結構，其目的是為了提高短序列的準確匹配率。如圖2所示，以單一序列為參考進行一次比對時的匹配率遠遠低于多次比對時的匹配率，更低于以多個序列為參考時多次比對的匹配率。圖中的hgl9和HuRef 是兩個具有代表性的參考序列,EMRs (Exact Mapped Reads), IMR (Inexact MappedReads) , UMR (Unmapped Reads)分別指準確匹配序列、非準確匹配序列、未匹配序列；
[0064]步驟四，對重點的‘Quality Value’域先采用k均值聚類法,然后采用基于上下文建模和統計的自適應文本壓縮方法PPMVC (命令行參數為‘e-04-rl’)壓縮聚類后的‘Quality Value，文件；
[0065]步驟五，將12個域的壓縮文件打包成一個完整的壓縮文件，作為最終的壓縮結果輸出。
[0066]如圖3所示，本實施例壓縮過程中的步驟三具體實施包括如下細節:
[0067]1、首先利用快速比對工具S0AP3將SAM/BAM文件的‘Sequence’域中的上百萬個短序列分線程地與參考序列作比對，得到以〈Read編號，參考序列上重復發生的染色體號，參考序列上重復發生的偏移位置，重復類型 > 標注的比對結果；
[0068]2、對于非準確匹配和未匹配的短序列，將其序列長度減半，即一個序列分成長度相當的兩個序列，并改變參考序列，再進行第二次比對，得到比對結果；
[0069]3、返回第一步，如此重復三到四次直到此時的序列長度低于15，剩下的非準確匹配和未匹配的短序列數已非常少，如圖3所示，經過分別以hgl9.fa，HuRef.fa, hgl9.fa和hgl9.fa為參考序列的四次比對后，剩余的非準確匹配和未匹配的短序列比例僅為0.01%,對它們再進行PPMVC (命令行參數為‘e-o8’)編碼即可。
[0070]實施效果:
[0071]依據上述步驟，選取了多組不同種類的下一代測序數據(Next GenerationSequencing, NGS):7 組來自千人基因組項目(IOOOGenomes Project), 一組來自于ChIP-Seq的mouse數據，一組來自于RNA-Seq的E.coli數據，還有一組來自于癌癥基因集(theCancer Genome Atlas, TCGA)的基因數據。本實施例比較了采用本發明所述的方法,MarkusHs1-Yang Fritz等人提出的基于目標序列與參考序列間的區別編碼的CRAM方法和JamesK.Bonfield等人提出的基于序列中每個base的位置坐標建模的Samcomp方法的性能:
[0072]對于來自不同項目組織的NGS基因數據，本發明所提出的方法均取得了明顯優于CRAM和與Samcomp相當的壓縮率(壓縮后的文件大小/壓縮前文件大小)。相對于CRAM的無損壓縮方法，本發明的無損模式在BAM文件上所產生的壓縮率降低了 6%-20%，達到了
0.5-0.65的壓縮率，節省了 35%-50%的存儲空間。在壓縮時間相當的情況下，其解壓時間小于CRAM。相對于Samcomp,本發明在保證壓縮效率相當的情況下,提供了 ‘Quality Value’壓縮等級的多選擇性，使得其更有適應性與擴展性。
[0073]以上對本發明的具體實施例進行了描述。需要理解的是，本發明并不局限于上述特定實施方式，本領域技術人員可以在權利要求的范圍內做出各種變形或修改，這并不影響本發明的實質內容。
【權利要求】
1.一種基于多參考序列的基因序列分級壓縮方法，其特征是，首先將BAM格式文件轉化成SAM格式的文件，SAM格式的基因序列由11個強制域和多個可選域構成，將可選域作為第12個域，原文件按域提取成12個獨立文件，然后對12個文件進行并行壓縮: (1)對‘Sequence’域采用基于多個參考序列逐步減半序列長度的分級壓縮方法； (2)對于‘QualityValue’域采用k均值聚類結合上下文建模PPMVC壓縮的方法； (3)對于剩下的十個域采用基于域內特征和域間相關性的壓縮方法。
2.根據權利要求1所述的基于多參考序列的基因序列分級壓縮方法，其特征是，所述對‘Sequence’域采用基于多個參考序列逐步減半序列長度的分級壓縮方法，具體為:利用快速比對工具S0AP3將SAM/BAM文件的‘Sequence’域中的短序列分線程地與參考序列作比對，對于準確匹配序列高效壓縮，對于非準確匹配和未匹配的短序列，將其序列長度減半，即一個序列分成長度相同的兩個序列，并改變參考序列，再進行第二次比對，得到比對結果,如此重復三至四次結束,剩余的非準確匹配和未匹配的短序列進行PPMVC編碼。
3.根據權利要求2所述的基于多參考序列的基因序列分級壓縮方法，其特征是，所述的對于準確匹配序列高效壓縮，具體為:對于準確比對的子序列Read，使用〈Read編號，參考序列上重復發生的染色體號，參考序列上重復發生的偏移位置，重復類型 > 這四個量來替代目標序列上重復的子序列，分別使用差分編碼+哈弗曼編碼、游程編碼、差分編碼+哈弗曼編碼和游程編碼來壓縮這四個分量。
4.根據權利要求1-3任一項所述的基于多參考序列的基因序列分級壓縮方法，其特征是，所述的對于‘Quality Value’域采用k均值聚類結合上下文建模PPMVC壓縮的方法，具體為:采用k均值 聚類法將η個Qasqi值聚成k類,使得聚類后每類內所有Quality value的值與聚類前的值差值平方最小，然后采用基于上下文建模和統計的自適應文本壓縮方法PPMVC壓縮聚類后的‘Quality Value’文件。
5.根據權利要求1-3任一項所述的基于多參考序列的基因序列分級壓縮方法，其特征是，所述的對于剩下的十個域采用基于域內特征和域間相關性的壓縮方法，具體為: 對于‘QNAME’域，用‘0’表示之前未出現過的QNAME，用逐漸遞增的數字編號與當前位置只差表示之前已經出現的QNAME，然后采用哈弗曼編碼壓縮這些非均勻分布的小型數值； 對于‘FLAG’域，用一個字節表示I~255之間的數值，用三個字節即0，x/256和x%256表示其它數值，然后采用哈弗曼編碼壓縮變換后的數值； 對于‘RNAME’域，用相同的數字標記整個SAM文件中的相同的參考序列名字，記錄下來所有參考序列，然后用游程編碼進行壓縮； 對于‘POS’域，采用差分編碼+哈弗曼編碼； 對于‘MAPQ’域，采用游程編碼； 對于‘CIGAR’域，采用LZW字典壓縮方法； 對于‘MRNM’域，采用游程編碼； 對于‘MPOS’域，結合‘MRNM’域的字符串，采用差分編碼+哈弗曼編碼； 對于‘TLEN’域，‘TLEN’域的值與‘MPOS’域減去‘POS’域的值的差即TLEN-(MPOS-POS))的絕對值服從于一個有限的字符集，對于該域的壓縮，結合‘POS’，‘MPOS’and ‘MRNM’三個域的信息采用Huffman編碼壓縮變換后的值；對于‘OP TIONAL’域，使用bzip2壓縮工具。
【文檔編號】H03M7/30GK103546160SQ201310433248
【公開日】2014年1月29日 申請日期:2013年9月22日 優先權日:2013年9月22日
【發明者】熊紅凱, 李平好 申請人:上海交通大學