序列搜索算法的進(jìn)化2

Next?generation sequencing era

第二代測序技術(shù)（NGS）出現(xiàn)后發(fā)生了重大的技術(shù)轉(zhuǎn)變谭贪。測序儀產(chǎn)生非常短的DNA片段，稱為reads（35-400nt）植酥，數(shù)量巨大（每次運(yùn)行數(shù)十GB的測序數(shù)據(jù)），成本低弦牡。其應(yīng)用包括：研究全基因組（genome-wide）友驮、RNA轉(zhuǎn)錄本（RNA-seq）、DNA甲基化（Methyl-seq）驾锰、蛋白質(zhì)-DNA相互作用（ChIP-seq）和蛋白質(zhì)組學(xué)圖譜（proteomic profiles）（外顯子測序卸留，exome sequencing）。

處理NGS測序數(shù)據(jù)最常見的任務(wù)是讀段映射（read mapping）：在參考基因組序列中定位reads并計(jì)算相應(yīng)的比對（locating reads within a reference genomic sequence and computing corresponding alignments）椭豫。數(shù)以百萬計(jì)的reads必須與一個(gè)包含數(shù)十億個(gè)字母的序列進(jìn)行比對耻瑟，類似BLAST的工具無法在合理的時(shí)間內(nèi)處理這樣的數(shù)據(jù)旨指。

常用的read mapping軟件包括：BWA-MEM、Bowtie2喳整、GEM谆构，它們比類似BLAST的工具快幾個(gè)數(shù)量級。這些算法的一個(gè)主要效率（efficiency）來源是被稱為BWT-或FM-index的索引結(jié)構(gòu)（indexing structure）算柳。

BWT-index基于Burrows-Wheeler transform低淡，從組合學(xué)（combinatorics）的角度BWT與Gessel-Reutenauer bijection（雙射姓言，一一對應(yīng)關(guān)系）有著密切的聯(lián)系瞬项。雖然BWT的主要應(yīng)用是文本壓縮（text compression），但它也適用于構(gòu)造一個(gè)簡潔的文本索引（compact text index）何荚。與以BLAST為代表的基于哈希（hash）的索引基礎(chǔ)結(jié)構(gòu)支持種子和擴(kuò)展策略（seed-and-extend strategies）不同囱淋，BWT索引是一個(gè)全文索引（full-text index），因?yàn)樗С炙阉魅我獯笮〉淖址吞痢WT-index以位（bits）為單位的大小接近于序列的無損表示所需的信息理論最小值妥衣。其他流行的全文索引（如后綴樹（suffix tree）、后綴數(shù)組（suffix array）及其變體）也被用于生物信息學(xué)算法戒傻，但需要更多的空間税手。

BWT算法是一種數(shù)據(jù)轉(zhuǎn)換算法，它將一個(gè)字符串中的相似字符放在相鄰的位置需纳，以便于后續(xù)的壓縮芦倒。BWT算法可以分為編碼和解碼兩部分。編碼后不翩，原始字符串中的相似字符會(huì)處在比較相鄰的位置兵扬；解碼就是將編碼后的字符串重新恢復(fù)成原始字符串的過程。BWT的一個(gè)特點(diǎn)就是經(jīng)過編碼后的字符串可以完全恢復(fù)成原始字符串（無損壓縮）口蝠。

Burrows和Wheeler在1994年提出了字符串匹配和壓縮索引的BWT方法器钟，應(yīng)用于無損壓縮。BWT通過對字符串循環(huán)移位得到一個(gè)字符矩陣妙蔗，然后通過排序和變換得到一個(gè)新的字符串傲霸，只改變了字符串中字符的順序，未改變字符本身眉反。轉(zhuǎn)換后的字符串中許多連續(xù)相同的字符會(huì)被重組在一起昙啄，更容易被壓縮。解碼時(shí)有序地按字典序排列BWT轉(zhuǎn)換后的字符串禁漓，即可還原未經(jīng)轉(zhuǎn)換前的字符串跟衅。在BWT結(jié)構(gòu)和技術(shù)的基礎(chǔ)上，F(xiàn)erragina和Manzini在2000年提出了一種后綴查找算法FM-index播歼，其搜索查找的時(shí)間復(fù)雜度為O(m)（m為查找串的長度）伶跷，并且與參考庫中序列長度無關(guān)掰读。

編碼：輸入字符串→加標(biāo)記$（$比字符串里的所有字符都要小）→循環(huán)移位（重復(fù)將最后一個(gè)字符移到開頭）→得到的新字符串從小到大排序→排序后的矩陣中最后一個(gè)字符構(gòu)成L列叭莫。

解碼：輸入L列蹈集；對L列進(jìn)行從小到大排序得到F列；根據(jù)L列的字符Li找到F列中的相同字符Fj雇初，然后得到Fj所在行的最后一個(gè)字符Lj拢肆。將Lj記錄下來；重復(fù)上面一步靖诗，直到Fj等于標(biāo)記字符為止郭怪；按照上述步驟找到的各個(gè)Lj進(jìn)行反向排列，得到字符串r刊橘。一般地鄙才，如果Li (i=1,2,…,n)這個(gè)字符在F列中出現(xiàn)過多次，分別是Fj, Fj+1, …促绵；并且假設(shè)L1到Li里這i個(gè)字符中一共有k個(gè)(k<=i)字符等于Li這個(gè)字符攒庵，那么我們要為Li在F列中找的對應(yīng)的字符就是Fj+k-1。

read mapping的結(jié)果：得到SAM（sequence alignment/map）格式败晴，或它的二進(jìn)制（binary）BAM格式浓冒。可以上傳到UCSC Genome Browser尖坤。使用SAMtools和BEDTools可以識別變異（call variants）稳懒、堆疊讀段（pile reads）、將SAM/BAM轉(zhuǎn)換為BED格式（瀏覽器擴(kuò)展數(shù)據(jù)糖驴，Browser Extensible Data）僚祷。BEDTools能夠發(fā)現(xiàn)重疊（finding overlaps）、計(jì)算覆蓋度（computing coverage）贮缕。SAMTools能夠?qū)AM文件進(jìn)行排序辙谜、合并和索引（sorting, merging, and indexing），在每個(gè)位置生成比對感昼，能夠發(fā)現(xiàn)SNPs（single-nucleotide polymorphisms）装哆。

無比對方法（alignment-free methods）

宏基因組學(xué)（metagenomics）任務(wù)可能是將數(shù)百萬reads比對到數(shù)千計(jì)的微生物基因組（microbial genomes），以闡明正在研究的環(huán)境樣本的組成定嗓。即使是專門的比對軟件也可能不滿足這項(xiàng)任務(wù)的實(shí)際時(shí)間要求蜕琴。另外，宏基因組學(xué)可能不需要計(jì)算比對宵溅，通常只需要識別給定read可能起源于哪個(gè)基因組凌简。一種方法是放棄比對，降低精度恃逻，以加快速度和節(jié)約內(nèi)存雏搂。

一個(gè)常見的想法是將序列看作是出現(xiàn)在其中的一（多）組模式（a (multi-)set of patterns）藕施。在最簡單的情況下，這些模式是固定長度k的子字符串（substrings）（k-mers）凸郑。然后通過相應(yīng)集合之間的相似性來表示序列之間的相似性裳食。這種方法被稱為無比對（alignment-free）或基于組合的比較（composition-based comparison）。大多數(shù)用于全基因組測序數(shù)據(jù)的現(xiàn)代宏基因組分類器（classifiers）都是基于k-mer分析的芙沥。并且可以通過將k-mers替換為spaced k-mers來提高精度（accuracy）诲祸。

k-mers是包含在生物序列中的長度為k的子序列。mer即為monomeric unit（單體單元）而昨。長度為L的序列對于一個(gè)給定的K可以得到L-k+1個(gè)k-mers救氯。每個(gè)位點(diǎn)的堿基可以為（A、T配紫、C径密、G）中的任意一個(gè)，因此k-mer理論上有4^k個(gè)不同的序列躺孝。由k-mer構(gòu)成的特征向量（x1, x2, …, xn），其中理論的n = 4^k底桂。例如植袍，如果k=4，x1=AAAA籽懦，x2=AAAT于个，...，x256=GGGG暮顺。

基因組組裝中厅篓，二代測序技術(shù)reads可能有100多個(gè)堿基，如果直接用reads進(jìn)行組裝捶码，原始reads中可能出現(xiàn)的錯(cuò)誤就會(huì)累加到最終的序列中羽氮，導(dǎo)致結(jié)果的不準(zhǔn)確性。通過將reads切割成以更短的k為單位的k-mer惫恼，由于測序錯(cuò)誤具有隨機(jī)性档押，這些包含測序錯(cuò)誤的k-mer絕大多數(shù)都是原測序物種中不存在的k-mer，只出現(xiàn)了1次祈纯，將這些k-mer去掉可以使組裝的結(jié)果更加可靠令宿。切割成相同長度的較小k-mer也有助于緩解不同起始read長度的問題。另外腕窥，能夠利用k-mer對物種基因組特征進(jìn)行評估粒没。例如：評估物種基因組大小、雜合度和重復(fù)序列比例簇爆、物種樣品污染等癞松。

De Bruijn graph是組裝基因組最常用的算法之一（其他還有Overlap-Layout-Consensus倾贰，String Graph等）。De Bruijn graph就是通過將reads打斷成k-mer后拦惋，利用k-mer之間的重復(fù)部分構(gòu)建圖匆浙，得到最優(yōu)化路徑從而拼接contig。尋找k-mer之間的重疊關(guān)系厕妖，建立De Bruijn圖首尼，即對于任意兩個(gè)k-mer，如果u的后k-1個(gè)堿基序列與w的前k-1個(gè)堿基序列相同言秸，則建立一條由u指向w的有向邊（edge）软能。尋找歐拉路徑（Eulerian path）來獲得結(jié)果序列Contigs。實(shí)際生成的De Bruijn graph非常復(fù)雜举畸，例如由于測序錯(cuò)誤查排、雜合或者高重復(fù)序列產(chǎn)生的tips翼尖（a）和bubbles氣泡（c），由于低覆蓋率造成的鏈接low coverage links（b）和（d）抄沮，微小重復(fù)造成的壓縮（e）跋核。因此需要對DBG圖進(jìn)行簡化，移除錯(cuò)誤鏈接叛买、刪除低覆蓋鏈接砂代、解開短重復(fù)序列等，最終輸出Contigs序列率挣。

根據(jù)overlap構(gòu)圖并搜尋contigs僅是第一步刻伊，完整的基因組組裝流程還包含Scaffold搭建、gap修補(bǔ)等過程椒功。雙末端測序中捶箱，相同序列名字的read1和read2來自同一條序列（中間不一定有overlap），可以根據(jù)paired-end信息將不同的Contigs搭建成Scaffold动漾。得到的Scaffold中間會(huì)有較多的gap丁屎，為了使組裝的序列更完整，需再次利用測序的雙末端數(shù)據(jù)之間的配對關(guān)系連接Contigs谦炬，并利用測序數(shù)據(jù)與已經(jīng)組裝的Contig之間的覆蓋關(guān)系對Contig之間空隙進(jìn)行補(bǔ)洞悦屏，延長Contigs。

k-mer大小的選擇對基因組組裝有多種影響键思。較小的k-mer將減少圖中存儲(chǔ)的edges數(shù)础爬，節(jié)約內(nèi)存；增加所有k-mer重疊的機(jī)會(huì)吼鳞；但是會(huì)面臨多頂點(diǎn)通向單個(gè)k-mer的風(fēng)險(xiǎn)看蚜，信息也會(huì)丟失；無法解決DNA中出現(xiàn)小微衛(wèi)星或重復(fù)序列問題赔桌。較大的k-mer會(huì)增加存儲(chǔ)DNA序列所需的內(nèi)存供炎；頂點(diǎn)的數(shù)目減少有助于基因組的構(gòu)建渴逻，因?yàn)槁窂阶兩倭耍惠^大的k-mer也會(huì)有較高的風(fēng)險(xiǎn)音诫，即沒有從每個(gè)k-mer出發(fā)的向外頂點(diǎn)惨奕；導(dǎo)致大量較小的contigs。

一般來說竭钝，選用17-mer來估算基因組大小梨撞，因?yàn)锳TCG四種不同的堿基組成長度為17的核苷酸有4^17=17G，足以覆蓋一般大小的基因組香罐；如果選擇15則只有1G卧波，對正常的基因組來說可能覆蓋度不夠，導(dǎo)致估計(jì)不準(zhǔn)庇茫。越長的k-mer片段越具很強(qiáng)的物種特異性港粱。當(dāng)基因組中有較多的重復(fù)序列時(shí)，可以使用較大的k-mer來跨過高重復(fù)的區(qū)域旦签，從而獲得更加準(zhǔn)確及完整的基因組草圖查坪。由于reads上的堿基錯(cuò)誤率的存在，選擇較長的k-mer會(huì)帶來較高的錯(cuò)誤率顷霹，但這也可以加大測序深度來彌補(bǔ)咪惠。用奇數(shù)的k-mer進(jìn)行分析，是為了防止組裝時(shí)正反鏈混淆淋淀，奇數(shù)的k-mer已經(jīng)被證明是不能夠匹配其反義互補(bǔ)鏈的。

k-mer常用的工具：JELLYFISH覆醇、KAT(The K-mer Analysis Toolkit)朵纷、KmerGenie、KMC永脓、Gerbil袍辞、Velvet、ABySS常摧、AllPath-LG搅吁、SOAPdenovo、EulerUSR落午、IDBA谎懦、SPAdes、Ray溃斋。

除了BWT-index外界拦，還有幾種數(shù)據(jù)結(jié)構(gòu)建立在Bloom filters的基礎(chǔ)上。Bloom filters已被應(yīng)用于表示de Bruijn graphs梗劫。Bloom filters的擴(kuò)展包括Cascading Bloom Filter享甸、Bloom Filter Trie截碴、Sequence Bloom tree及其變體。Bloom filters的替代方案有節(jié)省空間（space-efficient）和緩存友好（cache-friendly）的無損哈希方案（lossless hashing schemes）蛉威，例如counting quotient filters日丹。

Bloom過濾器是一種空間高效的（space-efficient）數(shù)據(jù)結(jié)構(gòu)，具有快速的查找時(shí)間和產(chǎn)生假陽性的可控風(fēng)險(xiǎn)蚯嫌。它最初是由Burton Bloom在20世紀(jì)70年代開發(fā)的哲虾，目的是減少容納哈希編碼信息所需的空間。Bloom filter可以用于檢索一個(gè)元素是否在一個(gè)集合中齐帚，它的優(yōu)點(diǎn)是空間效率和查詢時(shí)間都遠(yuǎn)遠(yuǎn)超過一般的算法妒牙，缺點(diǎn)是有一定的誤識別率和刪除困難（有假陽性，無假陰性）对妄。

Bloom Filter 原理：當(dāng)一個(gè)元素被加入集合時(shí)湘今，通過K個(gè)散列函數(shù)將這個(gè)元素映射成一個(gè)位數(shù)組（Bit Array）中的K個(gè)點(diǎn)，把它們置為1（它實(shí)際上是一個(gè)很長的二進(jìn)制向量和一系列隨機(jī)映射函數(shù)）剪菱。檢索時(shí)摩瞎，我們只要看看這些點(diǎn)是不是都是1就（大約）知道集合中有沒有它了：如果這些點(diǎn)有任何一個(gè)0，則被檢元素一定不在孝常；如果都是1旗们，則被檢元素很可能在。

Hash（散列构灸、哈希）函數(shù)上渴，就是把任意長度的輸入通過算法變成固定長度的輸出。直觀解釋起來喜颁，就是對一串?dāng)?shù)據(jù)m進(jìn)行雜糅稠氮，輸出另一段固定長度的數(shù)據(jù)h，作為這段數(shù)據(jù)的特征（指紋）半开。例如隔披，可以簡單地為字符串中的每個(gè)字符添加代碼值，并將結(jié)果mod（取余數(shù)）返回給定值寂拆。散列函數(shù)總是從相同的數(shù)據(jù)生成相同的散列值奢米，但實(shí)際上可以為兩個(gè)不同的數(shù)據(jù)值生成相同的散列值。也就是說纠永，哈希值對于給定的數(shù)據(jù)項(xiàng)不是唯一的鬓长，并且無法反轉(zhuǎn)哈希函數(shù)來獲取數(shù)據(jù)值。Bloom Filter與單哈希函數(shù)Bit-Map不同之處在于：Bloom Filter使用了k個(gè)哈希函數(shù)渺蒿，每個(gè)字符串跟k個(gè)bit對應(yīng)痢士。從而降低了沖突的概率。

所謂的BitMap就是用一個(gè)bit位來標(biāo)記某個(gè)元素所對應(yīng)的value，而key即是該元素怠蹂，由于BitMap使用了bit位來存儲(chǔ)數(shù)據(jù)善延，因此可以大大節(jié)省存儲(chǔ)空間。一個(gè)byte是占8個(gè)bit城侧，每一個(gè)bit的值是有或者沒有易遣，也就是二進(jìn)制的1或者0。