基于混合統(tǒng)計模型的DNA序列壓縮算法*

孫季豐 仝雪珂 譚麗

(華南理工大學電子與信息學院，廣東廣州510640)

DNA序列數(shù)據在諸多領域具有重要的研究和應用價值，如醫(yī)學、遺傳科學等.以獲取DNA序列數(shù)據為目標的生物測序工程近年來得到廣泛重視.隨著各種測序項目的展開，DNA序列數(shù)據以指數(shù)規(guī)模急劇增長，在有限的存儲資源內有效存儲和應用成本日益增大，因此需對DNA序列進行有效的壓縮.由于DNA序列的特殊性，普通壓縮算法效果并不理想，需要使用專門的算法進行壓縮［1］.

目前主要有兩類專門針對DNA序列數(shù)據的壓縮算法.第1類基于替代的壓縮算法在過去是主流算法，自1993年Grumbach等［2］學者首次提出Bio-Compress算法后，相繼出現(xiàn) GenCompress算法［3］、DNACompression算法［4］、DNAPack算法［5］、GeNML算法［6］、在線DNA序列的水平與垂直壓縮［7］、使用基于壓縮算法的近似重復和提供隨機存取能力的系統(tǒng)［8］、LSBD算法［9］、BioLZMA算法［10］以及改進的游程編碼算法［11］等基于替代的壓縮算法.第2類基于統(tǒng)計信息的壓縮算法只是其得不到滿意結果的補充.Allison等［12］在1998年提出ARM算法，Loewenstern等［13］在1999年提出CDNA算法，兩種算法在當時都產生了比替代算法更好的壓縮率.Korodi等［14］在2007年提出基于NML的算法，Cao等［15］在2007年提出專家模型算法(XM算法).NML算法和XM算法可單獨處理DNA序列，使統(tǒng)計信息壓縮算法的應用效果有了顯著提高.

文中主要研究第2類基于統(tǒng)計信息的數(shù)據壓縮算法.基于XM算法理論，Pinho等［16-18］使用有限上下文模型進行序列的統(tǒng)計壓縮，實驗結果表明，改進后的算法能產生更好的壓縮率.不同有限上下文模型的區(qū)別在于模型階數(shù)不同，鑒于其特點，文中將每一階模型改為含有某特定規(guī)律的多個模型的混合，可稱為混合統(tǒng)計模型.另外，文獻［16-18］中所述算法并未應用于標準DNA序列集，因此不具備普遍意義上的一般性.因此，基于XM算法理論，文中使用有限上下文混合統(tǒng)計模型，對標準DNA序列集進行算術編碼統(tǒng)計壓縮;并將文中算法與基于統(tǒng)計信息的當前較先進的XM算法以及其他壓縮算法進行了比較.

1 DNA序列的壓縮

DNA序列數(shù)據的壓縮通過兩個步驟完成，首先使用描述的模型完成序列的統(tǒng)計概率計算，其次使用所得概率計算結果進行算術編碼壓縮.

1.1 DNA序列的概率統(tǒng)計模型

文中使用的概率統(tǒng)計模型基于XM算法，提出的混合統(tǒng)計模型代替XM算法中使用的模型，應用于XM算法，進行序列符號的統(tǒng)計概率計算.

1.1.1 基本XM算法描述

首先基于已知符號的統(tǒng)計情況預測每一個符號位的概率分布，而后進行編碼來壓縮符號.預測概率需要設置專門的模型，XM算法結合使用少量復制和反相模型，提供了一個將各種模型混合使用預測符號出現(xiàn)概率的基本原理.

統(tǒng)計n位后，模型θk(k∈E)預測第n+1位符號的概率為［15］p(xn+1|θk，x1，2，…，n).其中θk代表使用的不同模型，k為模型的階數(shù)(θ2即2階該模型)，E為算法使用的模型階數(shù)集合.令 x1，2，…，n=x1x2…xn.基于貝葉斯平均，混合使用不同階數(shù)模型進行預測［15］:

權重ωθk，n為統(tǒng)計n位后k階模型θk估計xn+1出現(xiàn)概率的參與度，依照分析可由p(θk|x1，2，…，n)表示.評估模型的使用價值，需要使用距離當前位置最近的歷史序列.當前位置前的20個符號值能夠得到最好實驗結果的窗口值，因此，XM算法取窗口值w=20.

1.1.2 基本有限上下文模型

有限上下文模型描述如下(以4階為例)［16］.

圖1中輸出結果概率一般由下式所得［17］:

式(2)分母中∑na為計算到第n+1位時，前n位序列中上下文序列中(如圖1中ACAG)出現(xiàn)的次數(shù)，分子中ns為上下文序列后出現(xiàn)某符號s(如G)的次數(shù)，δ控制分配給未知事件的可能性，在階數(shù)小時默認為1.此式含義為計算某上下文序列(如ACAG)后一位出現(xiàn)某符號(如G)的概率，此式亦為公式(1)中當θk使用有限上下文模型時概率p(xn+1|θk，x1，2，…，n)的計算結果.

圖1 4階上下文模型統(tǒng)計概率計算Fig 1 Statistical probability calculation of 4-order context model

表1 上下文模型統(tǒng)計結果Table 1 Statistics results of context models

文中總結此模型的概率求解特點如下:①圖1所示的統(tǒng)計表是實時更新的;②當前位的統(tǒng)計概率基于更新到當前位置的統(tǒng)計表;③可設置不同的階數(shù)，某一位置符號針對不同階數(shù)存在不同的統(tǒng)計表和相應的統(tǒng)計概率.

將有限上下文模型做為XM算法中模型θk，即k(k∈E)階模型θk使用k階有限上下文模型，進行統(tǒng)計概率計算.

1.1.3 各階上下文模型的貢獻權重

(1)權重求解

使用貝葉斯理論估計式(1)中權重ωθk，n的值，計算基于下式［15］:

式(3)代表正相關比例關系.k階模型θk用k階有限上下文表示，令ωk，n為k階有限上下文模型應用于k階模型θk所對應的ωθk，n，重新表達式(3):

式(4)右邊是編碼xi(i=1，2，…，n)共消耗的二進制比特數(shù)，其中概率部分求解同式(2).式(4)右邊正比于權重值ωk，n，根據其所求各階模型的權重比例以及∑ωk，n=1，即可求出各階模型的權重值.由于上下文模型所得概率是實時更新的，因此所得權重也是實時更新的，即每運算到一位，都要計算當前位各階模型的貢獻權重.

至此求出上下文模型用于模型θk所需的統(tǒng)計量:各階模型的符號實時統(tǒng)計概率和貢獻權重，以此綜合計算符號的實時統(tǒng)計概率.其中，貢獻權重的求解基于相應階數(shù)模型的符號實時統(tǒng)計概率.

(2)k階上下文模型權重

有限上下文模型并不適用于數(shù)據序列的任意區(qū)域，而是依賴距離當前處理符號位最近的區(qū)域.借鑒學者Armando J.Pinho的基于遺忘策略的遞歸關系，可得較好的權重預測結果.

為式(4)引入參量γ［18］:

式(5)中參量γ能夠使符號的概率估計值受最接近當前位置的符號影響更大，較符合實際情況，因此能夠產生更準確的預測.

定義［18］:重新表達式(5)［18］:

最終可得［18］

上式僅為所求得的權重比例，再根據∑ωk，n=1，即可求得改進的權重值.參量γ的選取須進行實驗，根據最優(yōu)結果值(能夠使標準DNA序列產生最小的信息熵)進行選取.

1.1.4 混合統(tǒng)計模型

針對具體某一階上下文模型提出改進模型.以4階模型為例，原有限上下文模型的上下文Ci是4個連續(xù)符號序列(如ACAG).現(xiàn)將4個連續(xù)符號序列改為每2個連續(xù)符號取一個符號(隔一位取位模型)和每3個連續(xù)符號取一個符號(隔兩位取位模型)來進行上下文模型的選擇.如圖2所示.

圖2 原模型以及兩種改進模型Fig 2 Original model and two improved models

隔位取位模型處理更多的歷史數(shù)據，有助于取得更準確的預測概率.在原模型中，預測當前位符號的概率要處理4位歷史數(shù)據(見圖2(a)所示);而在隔一位取位模型中，需要處理7位歷史數(shù)據(如圖2(b)所示)，將7位數(shù)據隔一位取位，得到加粗表示的四位符號，作為隔一位取位模型的4階上下文，非加粗表示的數(shù)據略去不用;同理，圖2(c)所表示的隔兩位取位模型需要處理10位歷史數(shù)據，將加粗表示的4位符號作為此改進模型的4階上下文，非加粗表示的數(shù)據略去不用.

為使用更多歷史數(shù)據，可繼續(xù)生成隔三位取位模型，隔四位取位模型，…;文中只采用了原模型、隔一位取位模型與隔兩位取位模型.后續(xù)生成的一系列改進模型，可用做后續(xù)研究.

使用連續(xù)取位或隔一位取位、隔兩位取位，并沒有相關的先驗概率.故使用此3種模型，則分別取它們的概率各為1/3，將同一階數(shù)的3種取位模型混合使用，作為文中所提出的混合統(tǒng)計模型.將混合統(tǒng)計模型應用于XM算法.

使用4階原模型預測序列第n位出現(xiàn)符號A的概率用P4a1(n)表示，使用隔一位、隔兩位取位模型預測的概率分別為P4a2(n)、P4a3(n)，由于無先驗概率，得到

P4a(n)為4階混合統(tǒng)計模型在第n位上出現(xiàn)符號A的概率，同理可計算出P2a(n)、P6a(n).依XM算法原理，應用模型都要確定模型階數(shù)集合(包括原模型、混合統(tǒng)計模型)，在文中選擇k=2、4、6，即E={2，4，6}.可根據1.1.3部分實時計算出使用2、4、6階模型的權重ω2(n)、ω4(n)、ω6(n).據此可由XM算法部分公式(1)得

上式為第n位出現(xiàn)A的概率.同理可得此位出現(xiàn)C、G、T的概率分別Pc(n)、Pg(n)、Pt(n)，并存在

這即為進行算術編碼計算實時編碼比例所需的4種符號概率.當編碼序列長度為N時，需求解n在1～N中所有符號位上的4種符號概率.

1.2 DNA序列編碼

混合統(tǒng)計模型應用于XM算法，即可求出算術編碼的實時編碼比例，這既不同于傳統(tǒng)算術編碼中已知固定編碼比例，也不同于自適應算術編碼中編碼比例的實時求出，而是已求得的不同符號位置的不同編碼比例，以此為基礎進行算術編碼.

計算到第n位，對ACAG進行算術編碼:

(1)1.1.4部分實時計算出第n位上4種符號概率Pa(n)、Pc(n)、Pg(n)、Pt(n)，進一步轉化如圖3所示.

圖3 第n位實時編碼比例Fig 3 Real-time encoding proportion in the n bits

其中整個線段長度代表全部概率1，計算PA(n)、PC(n)、PG(n)和PT(n)，使得:

事實上，由于Pa(n)+Pc(n)+Pg(n)+Pt(n)=1，故PT(n)=1.第n位為符號A，則編碼取概率區(qū)間0～PA(n).

(2)再根據實時計算的第n+1位出現(xiàn)的4種符號概率Pa(n+1)、Pc(n+1)、Pg(n+1)、Pt(n+1)進行在區(qū)間0～PA(n)上的轉化，如圖4所示.

圖4 第n+1位實時編碼比例Fig 4 Real-time encoding proportion in the n+1 bits

計算PA(n+1)、PC(n+1)、PG(n+1)和PT(n+1)，使得:

又由于Pa(n+1)+Pc(n+1)+Pg(n+1)+ Pt(n+1)=1，故PT(n+1)=PA(n).第n+1位為符號C，則編碼取概率區(qū)間PA(n+1)～PC(n+1).

(3)下面繼續(xù)進行的編碼步驟同上，直到編碼完成得到一個區(qū)間，取區(qū)間中小數(shù)表示此段序列的編碼值.以此值編碼此一系列符號.所得小數(shù)數(shù)值由數(shù)字0～9表示，以下分析將其轉換為二進制表示.

在沒有先驗概率的前提下，首先假設編碼后10個數(shù)字0～9出現(xiàn)的概率是相同的.由于23＜10＜24，所以表示每一個數(shù)字的二進制符號位數(shù)在3～4之間.三位二進制數(shù)表示值有000、001、010、011、100、101、110、111，共8個.為使增添的四位二進制數(shù)在解碼時不與三位二進制數(shù)相混淆.去掉000，增添0000和 0001;去掉 111，增添 1110和1111.則10個數(shù)字可以由001、010、011、100、101、110、0000、0001、1110、1111表示.在解碼時，先觀察前3個符號位，當都相同(為000或111)時，則使用前4個符號位對數(shù)字進行解碼.當不全相同時，則使用3個符號位對數(shù)字解碼.則可一一對應解碼數(shù)字0～9.編碼每一個數(shù)字所需的平均二進制符號位為

(3×6+4×4)/10=3.4.

(4)隨著編碼位數(shù)的增加，所需的小數(shù)位數(shù)也在逐步增加，因此需要設定一個固定的長度，在計算機軟件的計算精度范圍以內，經試驗取25最為合適.即每次取25位進行算術編碼，二進制轉換.

混合統(tǒng)計模型應用于XM算法的流程圖如圖5所示.

圖5 完整算法流程圖Fig 5 Flowchart of full algorithm

2 實驗數(shù)據與結果分析

2.1 參數(shù)設置

1.1.2 部分上下文模型由序列第1位開始進行實時統(tǒng)計.在初始位，令表1中數(shù)據全部為零;各符號(A、C、G、T)出現(xiàn)的初始概率均為1/4;模型階數(shù)集選取E={2，4，6};各階模型初始貢獻權重相同(文中令ω2(1)、ω4(1)、ω6(1)均為1/3);式(2)中δ=1.

2.2 結果與分析

實驗選取美國國家生物技術信息中心(NCBI)收集的DNA序列數(shù)據.為進行對照，將算法應用于DNA壓縮領域文獻中最常使用的標準DNA序列集.

壓縮之前，分析1.1.3部分的權重改進引入的參數(shù)γ.經驗獲得的γ在0.970～0.998之間.混合統(tǒng)計模型所得標準DNA序列符號平均信息熵見表2.

表2中數(shù)據為表示一個堿基符號所需的二進制比特數(shù)(單位為信息熵單位Bits Per Base，縮寫為BPB)，加粗數(shù)據為使用不同參量值所得最優(yōu)值(最小信息熵).當γ取0.995時，能夠取得在各標準序列上的最小信息熵，故γ=0.995.

表3示出了應用于標準DNA序列集時，使用原模型和文中提出的混合統(tǒng)計模型進行壓縮的結果比較(原模型與混合統(tǒng)計模型均需混合使用3個階數(shù)有限上下文模型，階數(shù)k分別為2、4、6).顯然有:混合統(tǒng)計模型比原模型需要更小的二進制比特表示(9段序列使用混合統(tǒng)計模型顯示更好的壓縮效果)，因此混合統(tǒng)計模型能夠顯示更好的壓縮效果.

表2 使用不同參量值在標準序列上的信息熵Table 2 Information entropies of the standard sequences by using different parameter values BPB

表3 不同模型使用方式在標準序列上的壓縮率Table 3 Compression ratios in the standard sequences by using different models BPB

表4示出了文中提出的混合統(tǒng)計模型壓縮算法同其他應用于DNA序列壓縮算法的比較.文中提出的多階混合統(tǒng)計模型和一系列用于對照的壓縮算法相比更有優(yōu)勢.在序列 CHMPXX、CHNTXX、HUMGHCSA上，XM算法產生最優(yōu)壓縮結果.在序列HUMDYSTROP、HUMHBB、HUMHDABCD、HUMHPRTB、MPOMTCG和MTPACGA上，混合統(tǒng)計模型產生最優(yōu)壓縮結果.在序列 VACCG上 DNACompress算法有最好的壓縮比，序列HEHCMV上是DNAPack算法.

表4 不同壓縮算法在標準序列上的壓縮率Table 4 Compression ratios of compression algorithms in the standard sequences BPB

根據數(shù)據分析，和基于統(tǒng)計信息的當前較先進的XM算法相比，文中所使用的混合統(tǒng)計模型能夠產生更好的壓縮效果.對于使用XM算法壓縮效果較差的標準DNA序列，混合統(tǒng)計模型算法表現(xiàn)出更好的壓縮效果，能夠彌補XM算法應用于一些序列的不足，完成進一步壓縮。

將文中算法壓縮人類染色體的結果與XM算法進行對比(見表5).實驗結果表明，文中算法壓縮效果比XM算法差，未彰顯算法優(yōu)勢，因此僅在此列出兩條染色體(FASTA格式)的對比結果(序列下載鏈接為 http:∥ www.ncbi.nlm.nih.gov/nuccore/ 1577282 28?report=fasta和http:∥www.ncbi.nlm. nih.gov/nuccore/157729478?report=fasta).

表5 對兩條人類染色體的壓縮率Table 5 Compression ratios in two Homo sapiens chromosomes BPB

由表5可知:文中提出的算法對人類染色體這類長度超過千萬的高通量DNA數(shù)據［19-20］的壓縮效果有待提高.

3 結語

文中在多個階數(shù)有限上下文模型綜合應用的基礎上，對每一階模型進行多種取位的混合，再使用混合模型壓縮數(shù)據，同原模型進行對比，并同其他一系列經典DNA序列壓縮算法(包含當前基于統(tǒng)計理論壓縮效果最優(yōu)的XM算法)進行對照.實驗結果表明:改進算法優(yōu)于其他壓縮算法，能夠彌補XM算法應用于序列標準DNA集部分數(shù)據的不足.但對高通量DNA數(shù)據(如人類染色體)的壓縮有待提高.

由于混合模型只是原模型與隔一位、隔兩位取位模型的未知先驗概率的混合使用，因此還可以進行模型的后續(xù)改進:①增加新模型的種類，可依據文中的兩種取位模型的趨勢進一步改進;②統(tǒng)計估計每種模型的應用效果，進一步為其設置使用權重值，此步驟的進行被限制在某一階模型中.可依據以上提出的兩個方向進行進一步研究，分析壓縮率.

［1］ 紀震，周家銳，姜來，等.DNA序列數(shù)據壓縮技術綜述［J］.電子學報，2010，38(5):1113-1121. Ji Zhen，Zhou Jia-rui，Jiang Lai，et al.Overview of DNA sequence data compression techniques［J］.Acta Electronica Sinica，2010，38(5):1113-1121.

［2］ Grumbach S，Tahi F.Compression of DNA sequences［C］∥Proceedings of 1993 Data Compression Conference.Snowbird:IEEE，1993:340-350.

［3］ Chen X，Kwong S，Li M.A compression algorithm for DNA sequences and its applications in genome comparison［J］.Genome Informatics Series，1999(10):51-61.

［4］ Chen X，Li M，Ma B，et al.DNACompress:fast and effective DNA sequence compression［J］.Bioinformatics，2002，18(12):1696-1698.

［5］ Behzadi B，Le Fessant F.DNA compression challenge revisited:a dynamic programming approach［C］∥Proceedings of Combinatorial Pattern Matching.Berlin-Heidelberg:Springer-Verlag，2005:190-200.

［6］ Korodi G，Tabus I.An efficient normalized maximum likelihood algorithm for DNA sequence compression［J］. ACM Transactions on Information Systems，2005，23(1): 3-34.

［7］ Mishra K N，Aaggarwal D A，Abdelhadi D E，et al.An efficient horizontal and vertical method for online DNA sequence compression［J］.International Journal of Computer Applications，2010，3(1):39-46.

［8］ Kaipa K K，Lee K，Ahn T，et al.System for random access DNA sequence compression［C］∥Proceedings of 2010 IEEE International Conference on Bioinformatics and Biomedicine Workshops.Hong Kong:IEEE，2010:853-854.

［9］ Mridula T V，Samuel P.Lossless segment based DNA compression［C］∥Proceedings of 2011 3rd International Conference on Electronics Computer Technology.Kanyakumari:IEEE，2011:298-302.

［10］ 紀震，周家銳，朱澤軒，等.基于生物信息學特征的DNA序列數(shù)據壓縮算法［J］.電子學報，2011，39 (5):991-995. Ji Zhen，Zhou Jia-rui，Zhu Ze-xuan，et al.Bioimformatics features based DNA sequence data compression algorithm［J］.Acta Electronica Sinica，2011，39(5):991-995.

［11］ Ouyang J，F(xiàn)eng P，Kang J.Fast compression of huge DNA sequence data［C］∥Proceedings of 2012 5th IEEE International Conference on Biomedical Engineering and Informatics.Chongqing:IEEE，2012:885-888.

［12］ Allison L，Edgoose T，Dix T I.Compression of strings with approximate repeats［C］∥Proceedings of the 6th International Conference on Intelligent Systems for Molecular Biology.Montréal:AAAI，1998:8-16.

［13］ Loewenstern D，Yianilos P N.Significantly lower entropy estimates for natural DNA sequences［J］.Journal of Computational Biology，1999，6(1):125-142.

［14］ Korodi G，Tabus I.Normalized maximum likelihood model of order-1 for the compression of DNA sequences［C］∥Proceedings of 2007 Data Compression Conference. Snowbird:IEEE，2007:33-42.

［15］ Cao M D，Dix T I，Allison L，et al.A simple statistical algorithm for biological sequence compression［C］∥Proceedings of 2007 Data Compression Conference. Snowbird:IEEE，2007:43-52.

［16］ Pinho A J，Neves A J R，F(xiàn)erreira P J S G.Inverted-repeatsaware finite-context models for DNA coding［C］∥Proceedings of the 16th European Signal Processing Conference.Lausanne:Technical Program Committee，2008.

［17］ Pinho A J，Neves A J R，Bastos C A C，et al.DNA coding using finite-context models and arithmetic coding［C］∥Proceedings of 2009 IEEE International Conference on Acoustics，Speech and Signal Processing.Taipei:IEEE，2009:1693-1696.

［18］ Pinho A J，Pratas D，F(xiàn)erreira P J S G.Bacteria DNA sequence compression using a mixture of finite-context models［C］∥Proceedings of 2011 IEEE Statistical Signal Processing Workshop(SSP).Aveiro:IEEE，2011: 125-128.

［19］ Wang C，Zhang D.A novel compression tool for efficient storage of genome resequencing data［J］.Nucleic Acids Research，2011，39(7):e45/1-6.

［20］ Pinho A J，Pratas D，Garcia S P.GReEn:A tool for efficient compression of genome resequencing data［J］.Nucleic Acids Research，2012，40(4):e27/1-8.