BWT輪換矩陣之完全桶排序演算法

BWT輪換矩陣之完全桶排序演算法Luo Weifeng 2011-7-3聲明： 演算法不專業，高手請繞行。        首先，慶祝一下順利到達北京，入住到798藝術區，這裡的人兒還真是一個個藝術到家了。住處挺不錯，處處顯示出文藝細胞。一路上的行程都很順利，幾乎沒有什麼停頓。又一次來到帝都，有種很特別的感覺。還有一個灰常值得慶祝的事情便是csdn的新版blog了，因為之前我的很多同學都因為csdn的不好用而投身其他地方了，只有我還在這裡掙紮著，新版blog很給力，就這個線上編輯器有了很大的改觀，為了慶祝新版blog的上線，準備這幾天每天堅持寫blog。        BWT(Burrows-Wheeler Transformation)演算法在人類基因組測序方面有很重要的應用，開放源碼的bzip就是bwt壓縮演算法的成功案例。關於這裡的其他知識可以去看維基百科。這篇文章主要介紹我設計的一個產生BWT演算法所需的 L（也就是輪換矩陣最後一列）的桶排序的方法。       問題輸入：一個2  M的人類基因組一條基因的序列chr1.fa。檔案內容是NACGT字母的亂序。整個檔案是一個完整的序列。       要求輸出： BWT演算法的L串。演算法分析： 資料模型：

  演算法所需要的資料儲存在堆中的一個地區，這個地區的頭為指標P，地區長度為N。那麼第一個字串便是從P到P+（N-1）。第a個字串的第b個元素為 P+（a + b）%N。

  F a(b)= P+( a + b)%N.

核心思想：

演算法的核心思想就是利用人類基因組串內有用字元只有ACGT的特點，對數百萬個基因串進行桶排序，並且對每一個子桶也進行相同的桶排序，直到桶裡邊的元素只剩下1個或0個，觸法桶的收集。桶的收集在檢查到父桶可以合并時觸法合并，直到合并完成。

詳細：

桶存放的是串標識。桶裡邊的元素4，就表示第5個串，按照上頁得公式便表示從P+（5+0）%N到P+（5+N-1）%N的串。

預先處理：首先由預先處理程式從檔案中讀入未經處理資料檔案，對所有的N進行過濾，並將其他內容載入到內容堆中的一個字元數組中。然後初始化一個0桶，這個桶包含了所有的字串（桶裡邊只存的是字串所在的位移，比方說，我存了4，那麼就表示自 P+4 到P + （4 + N-1）%N的那個字元序列）。

分發：由控製程序將0桶傳遞給分發程式進行分發處理。分發程式按照字元序列的首個字母是A還是C、G、T，把這些字元標示分發到0A、0C、0G、0T桶裡邊。並且一邊分發一邊統計分發到這些桶裡的entry數量，若數量大於等於2則繼續下一層分發，否則對桶做標記，並觸發收集。（比如：如果檢測到桶0AA的資料條目少於兩條，將桶0AA重新命名為桶$AA）。

收集：收集時先檢測同目錄的桶是不是都被收集，若有桶沒有收集，則放棄執行，否則產生父收集桶，然後將四個子桶的資料收集到父桶之中。如果父桶不是$（最終桶），則觸發父桶的收集。（同樣接著前面的例子：現在觸發了對桶$AA的收集，則首先檢測$AC、$AG、$AT是否存在，如果有一個不存在則放棄執行並退出，如果都存在則建立父收集桶$A,並把$AA、$AC、$AG、$AT裡邊的資料收集到父桶$A中)

收尾：程式到這裡，頂層收集桶$已經產生。依次提取桶中的串，輸出此串的最後一個字元。即可得到bwt輪換矩陣的最後一列L。比如在裡邊拿到一個int資料a，則需要輸出的是P+(a+N-1)%N。

問題： 這裡最直接的一個問題就是，桶作為一個檔案時，檔案的名稱的長度會隨著未經處理資料的增長而同步增長。解決的辦法是，在邏輯上還是使用剛才提到的方法，而在實際使用桶時，使用此桶名的32位md5值作為實際的檔案名稱。

           還有一個問題，就是如果使用上面介紹的方法時在操作上會出現記憶體不夠、堆疊溢位的情況，解決的方法也很直接，通過第三方控制函數對這兩個原子操作進行控制，徹底拋棄遞迴方法。這個在以後慢慢改過來，現在貼一個遞迴的實現，注意，這個程式有bug，主程式我會在有時間了把控製程序加進去，就可以避免記憶體問題了。

想了想還是直接上傳工程吧，工程是code：：blocks的工程，使用標準C，現在的工程還很爛，等過幾天改改。

連結地址：http://download.csdn.net/source/3414631

演算法的優缺點： 缺點很明顯了，就是需要大量的檔案操作，需要動態建立和刪除檔案，但是優點也是很明顯的，整個產生過程不需要比較操作，操作被原子的分為了分發和收集，醬紫的模型與分布式hadoop幾乎完全一樣，也就是說，這樣的模型很容易被擴充到分散式處理應用當中，所以，儘管在單機看來，不如其他一些成型的演算法，也不如桶排加快排的做法，但是可以很容易的擴充到分散式處理當中，那時候效能的提高就非常明顯了。