經(jīng)過了第四節(jié)的長文，我想大家基本上已經(jīng)知道了一個(gè)WGS流程該如何構(gòu)建起來了吧同窘。但在那一節(jié)中限于篇幅有兩個(gè)很重要的文件我沒能展開來講喷鸽，分別是：BAM和VCF文件。這篇我們先說BAM文件稚晚。
什么是BAM
BAM是目前基因數(shù)據(jù)分析中最通用的比對(duì)數(shù)據(jù)存儲(chǔ)格式，它既適合于短read也適合于長read型诚，最長可以支持128Mbp的超大read客燕！除了后綴是.bam之外，有些同學(xué)可能還會(huì)看到.cram狰贯，甚至.sam后綴的文件也搓，其實(shí)它們一個(gè)是BAM的高壓縮格式(.cram)——IO效率比原來的BAM要略差；另一個(gè)是BAM的純文本格式（.sam）涵紊。當(dāng)然格式都是一樣的傍妒，因此為了描述上的清晰，我下面都統(tǒng)一用BAM摸柄。
BAM文件格式
其實(shí)一開始它的名字是SAM（The Sequencing Alignment/Map Format的英文簡稱）颤练，第一次出現(xiàn)的時(shí)候，它是bwa比對(duì)軟件的標(biāo)準(zhǔn)輸出文件驱负，但原生的SAM是純文本文件嗦玖，十分巨大（比如：一個(gè)人30x全基因組測序的sam大小超過600G），非常容易導(dǎo)致存儲(chǔ)空間爆滿跃脊！為了解決這個(gè)問題宇挫，bwa的開發(fā)者李恒（lh3）設(shè)計(jì)了一種比gz更加高效的壓縮算法，對(duì)其進(jìn)行壓縮匾乓，這就是我們所說的BAM捞稿，它的文件大小差不多只有原來的1/6。
在2007年拼缝，NGS技術(shù)剛剛興起之時(shí)娱局，各類短序列比對(duì)軟件層出不窮，輸出格式也是各有特點(diǎn)咧七，各家各有一套衰齐，并沒有什么真正的標(biāo)準(zhǔn)可言，可以說那是一個(gè)誰都說我最好的時(shí)期继阻。
但逐漸的耻涛，研究者們發(fā)現(xiàn)BAM格式對(duì)Mapping信息的記錄是最全面的，用起來也是最靈活的瘟檩。bwa的作者還為BAM文件開發(fā)了一個(gè)非常好用的工具包——Samtools抹缕，使得人們對(duì)BAM文件的處理變得十分便利，拓展性也變得非常強(qiáng)墨辛，后來還有類似于IGV等專門支持BAM的工具也越來越多卓研，因此它就逐漸成為了主流。
現(xiàn)在基本上所有的比對(duì)數(shù)據(jù)都是用BAM格式存儲(chǔ)的，儼然已經(jīng)成為了業(yè)內(nèi)的默認(rèn)標(biāo)準(zhǔn)奏赘。
在2013年寥闪，研究者們還專門將Samtools的處理核心剝離出來，并將其打包成為一個(gè)專門用于處理高通量數(shù)據(jù)的API——htslib磨淌，除了C語言版本之外還有Java和Python版本疲憋，這些在github上都能直接找到。后續(xù)許多與NGS數(shù)據(jù)處理有關(guān)的工具基本都會(huì)使用這個(gè)API進(jìn)行相關(guān)功能的開發(fā)梁只，可見其影響力缚柳。
ok，背景的介紹就先到此為止了敛纲，我們回歸主題喂击。下面這個(gè)圖是我從一份剛剛完成比對(duì)的bam文件中截取出來的內(nèi)容：
由于屏幕所限，無法把全部的內(nèi)容都包含進(jìn)來淤翔，特別是header信息翰绊，貼在這里僅是為了讓還沒見過BAM文件的同學(xué)們能夠?qū)λ幸粋€(gè)總體的感覺。
如果是SAM文件旁壮，同時(shí)你也熟悉linux操作的話监嗜，直接在linux終端用less打開即可（注意：不要試圖在本地使用文本編輯器，如vim等直接打開文件抡谐，會(huì)撐死機(jī)子的）裁奇，但如果我們要查看的是BAM，那么必須通過Samtools（可以到samtools的網(wǎng)站下載并安裝）麦撵。
$ less -SN in.sam # 打開sam文件
$ samtools view -h in.bam # 打開bam文件
BAM文件分為兩個(gè)部分：header和record刽肠。這里額外說一句，許多NGS組學(xué)數(shù)據(jù)的存儲(chǔ)格式都是由header和record兩部分組成的免胃。
以上例子音五，在samtools view中加上-h參數(shù)目的是為了同時(shí)把它的header輸出出來，如果沒有這個(gè)參數(shù)羔沙，那么header默認(rèn)是不顯示的躺涝。
header內(nèi)容不多，也不會(huì)太復(fù)雜扼雏，每一行都用‘@’ 符號(hào)開頭坚嗜，里面主要包含了版本信息，序列比對(duì)的參考序列信息诗充，如果是標(biāo)準(zhǔn)工具（bwa苍蔬，bowtie，picard）生成的BAM蝴蜓，一般還會(huì)包含生成該份文件的參數(shù)信息（如上圖）碟绑，@PG標(biāo)簽開頭的那些這里需要重點(diǎn)提一下的是header中的@RG也就是Read group信息，這是在做后續(xù)數(shù)據(jù)分析時(shí)專門用于區(qū)分不同樣本的重要信息。它的重要性還體現(xiàn)在蜈敢，如果原來樣本的測序深度比較深，一般會(huì)按照不同的lane分開比對(duì)汽抚，最后再合并在一起抓狭，那么這個(gè)時(shí)候你會(huì)在這個(gè)BAM文件中看到有多個(gè)RG，里面記錄了不同的lane造烁，甚至測序文庫的信息否过，唯一不變的一定是SM的sample信息，這樣合并后才能正確處理惭蟋。
其實(shí)苗桂，關(guān)于這一點(diǎn)我在上一篇文章中講序列比對(duì)時(shí)的也特意強(qiáng)調(diào)了這些方面，不記得的同學(xué)們也可以翻看上一篇的相關(guān)內(nèi)容告组。
接下來重點(diǎn)要說的是BAM的核心：record（有時(shí)候也叫alignment section煤伟，即，比對(duì)信息）木缝。這是我們通常所說的序列比對(duì)內(nèi)容便锨，每一行都是一條read比對(duì)信息，它的記錄看起來是這樣的：
我這里借用了網(wǎng)上的一張圖片來輔助說明我碟，recoed中的每一個(gè)信息都是用制表符tab分開的放案。
下面我們就來仔細(xì)瞧瞧這里的每一個(gè)信息分別都是什么。
以上矫俺，前11列是所有BAM文件中都必須要有的信息吱殉，而且從描述中我們也能夠比較清楚地知道其所代表的含義。但其中厘托，有幾個(gè)信息實(shí)在太重要了友雳，以至于我認(rèn)為有必要對(duì)其進(jìn)行詳細(xì)說明。
第一催烘，F(xiàn)lag信息
這是一個(gè)非常特別并且重要的數(shù)字沥阱，也是一個(gè)容易被忽視的數(shù)字，這可能和許多生信工程師也并不完全理解這個(gè)值有關(guān)伊群。許多同學(xué)在第一次看到其官方文檔中的描述之后依然會(huì)覺得十分困惑考杉，但它里面實(shí)際上記錄了許多有關(guān)read比對(duì)情況的信息。想要讀懂它的一個(gè)關(guān)鍵點(diǎn)是我們不能夠?qū)⑵湟暈橐粋€(gè)數(shù)字舰始，而是必須將其轉(zhuǎn)換為一串由0和1組成的二進(jìn)制碼崇棠，這一串二進(jìn)制數(shù)中的每一個(gè)位(注意是“位”，bit的意思)都代表了一個(gè)特定信息丸卷，它一共有12位（以前只有8位）枕稀，所以一般會(huì)用一個(gè)16位的整數(shù)來代表，這個(gè)整數(shù)的值就是12個(gè)0和1的組合計(jì)算得來的，因此它的數(shù)值范圍是0~2048（2的12次方萎坷，計(jì)算機(jī)科學(xué)的同學(xué)對(duì)這種計(jì)算應(yīng)該不陌生）凹联。
那么下面我就結(jié)合其文檔和自己的實(shí)踐經(jīng)驗(yàn)對(duì)這12個(gè)位的含義用更加通俗易懂的語言來重新描述，如下表：
所以哆档，通過上面這個(gè)表的信息蔽挠，我們就可以清楚地知道每一個(gè)FLAG中都包含了什么信息。比如看到FLAG = 77時(shí)瓜浸，我們第一步要做的就是將其分解為二進(jìn)制序列（也可以理解為分解成若干個(gè)2的n次方之和）：
77 = 000001001101 = 1 + 4 + 8 +64澳淑，這樣就得到了這個(gè)FLAG包含的意思：PE read，read比對(duì)不上參考序列插佛，它的配對(duì)read也同樣比不上參考序列杠巡，它是read1。
當(dāng)然雇寇，如果你希望自己在程序中寫一段處理FLAG的代碼氢拥，那么顯然是不會(huì)像我們這個(gè)例子那樣去分解這個(gè)整數(shù)的，多麻煩跋呛睢兄一！那么該如何做呢？其實(shí)也很簡單识腿，比如我們要獲得其中某個(gè)位（假設(shè)第N位）的值——只需要將這個(gè)FLAG值和2的N次方做與的運(yùn)算即可出革。在與運(yùn)算時(shí)，F(xiàn)LAG值首先會(huì)被轉(zhuǎn)換成一串二進(jìn)制序列（如77=000001001101）渡讼，而2的N次方除了第N位是1之外骂束，其它的都是0，“與”了之后其它信息就會(huì)被屏蔽掉成箫。比如展箱，我們想知道該read是否比對(duì)上了參考序列，那么只需要計(jì)算FLAG & 4 的值就行了蹬昌，如果結(jié)果是1那么就是比對(duì)上了混驰，如果是0則代表沒有比上。
不過皂贩，在實(shí)際工作中栖榨，除非遇到特殊的情況，否則我一般更推薦調(diào)用官方的htslib這個(gè)包來協(xié)助處理明刷，它是一個(gè)C語言庫婴栽，如果你用Python，則是pysam——htslib的python包（Java則是htsjdk）辈末，包中已經(jīng)幫我們做了這些處理愚争，可以直接得到結(jié)果映皆，下一篇文章里我會(huì)用pysam舉例說明如何用它來操作bam文件。
另外轰枝，下面這一段代碼是htslib（samtools的核心庫）中定義的12個(gè)與flag值進(jìn)行與操作獲取對(duì)應(yīng)位信息的變量捅彻，感興趣的同學(xué)可以再htslib里面的sam.h文件中找到，在做一些需要觸達(dá)基礎(chǔ)性原理的開發(fā)時(shí)或許你會(huì)用到鞍陨。
/*! @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */
#define BAM_FPAIRED 1
/*! @abstract the read is mapped in a proper pair */
#define BAM_FPROPER_PAIR 2
/*! @abstract the read itself is unmapped; conflictive with BAM_FPROPER_PAIR */
#define BAM_FUNMAP 4
/*! @abstract the mate is unmapped */
#define BAM_FMUNMAP 8
/*! @abstract the read is mapped to the reverse strand */
#define BAM_FREVERSE 16
/*! @abstract the mate is mapped to the reverse strand */
#define BAM_FMREVERSE 32
/*! @abstract this is read1 */
#define BAM_FREAD1 64
/*! @abstract this is read2 */
#define BAM_FREAD2 128
/*! @abstract not primary alignment */
#define BAM_FSECONDARY 256
/*! @abstract QC failure */
#define BAM_FQCFAIL 512
/*! @abstract optical or PCR duplicate */
#define BAM_FDUP 1024
/*! @abstract supplementary alignment */
#define BAM_FSUPPLEMENTARY 2048
第二沟饥，CIGAR
CIGAR是CompactIdiosyncraticGappedAlignmentReport的首字母縮寫，稱為“雪茄”字符串湾戳。
作為一個(gè)字符串，它用數(shù)字和幾個(gè)字符的組合形象記錄了read比對(duì)到參考序列上的細(xì)節(jié)情況广料，讀起來要比FLAG直觀友好許多砾脑，只是記錄的是不同的信息。比如艾杏，一條150bp長的read比對(duì)到基因組之后韧衣，假如看到它的CIGAR字符串為：33S117M，其意思是說在比對(duì)的時(shí)候這條read開頭的33bp在被跳過了（S）购桑，緊接其后的117bp則比對(duì)上了參考序列（M）畅铭。這里的S代表軟跳過（Soft clip），M代表匹配（Match）勃蜘。CIGAR的標(biāo)記字符有“MIDNSHP=XB”這10個(gè)硕噩，分別代表read比對(duì)時(shí)的不同情況：
除了最后‘=XB’非常少見之外，其它的標(biāo)記符通常都會(huì)在實(shí)際的BAM文件中碰到缭贡。另外炉擅，對(duì)于M還是再強(qiáng)調(diào)一次，CIGAR中的M阳惹，不能覺得它代表的是匹配就以為是百分百?zèng)]有任何miss-match谍失，這是不對(duì)的，多態(tài)性堿基或者單堿基錯(cuò)配也是用M標(biāo)記莹汤！
第三快鱼，MAPQ，比對(duì)質(zhì)量值
這個(gè)值同樣非常重要纲岭，它告訴我們的是這個(gè)read比對(duì)到參考序列上這個(gè)位置的可靠程度抹竹，用錯(cuò)誤比對(duì)到該位置的概率值（轉(zhuǎn)化為Phred scale）來描述：-10logP{錯(cuò)比概率}。因此MAPQ（mapping quality）值大于30就意味著錯(cuò)比概率低于0.001（千分之一）止潮，這個(gè)值也是我們衡量read比對(duì)質(zhì)量的一個(gè)重要因子柒莉。
剩下的幾列在上面的格式表中描述的也比較清楚，基本沒有過于隱藏的信息沽翔，因此我就不打算再一一細(xì)說了兢孝，如果大家依然有困惑可以到后臺(tái)留言窿凤。
此外，細(xì)心的同學(xué)可能也已經(jīng)發(fā)現(xiàn)了：fastq的所有信息都被涵蓋到了BAM文件中了跨蟹，包括比對(duì)不上的read也在雳殊，因此獲得了BAM其實(shí)也等于獲得了所有的read。而且窗轩，fastq有時(shí)也會(huì)被轉(zhuǎn)換成一種uBam文件夯秃，指的就是un-mapping BAM——沒有做過比對(duì)的BAM文件。它相比于Fastq可以用metadata存儲(chǔ)更多有用的信息痢艺，不過這不是我們這篇文章想說的內(nèi)容仓洼。
最后，還是再說明一次：BAM文件中除了必須的前11列信息之外堤舒，不同的BAM文件中后面記錄metadata的列是不固定的色建，在不同的處理軟件中輸出時(shí)也會(huì)有所不同，我們也可以依據(jù)實(shí)際的情況增刪不同的metadata信息舌缤。
使用samtools view查看BAM文件
BAM文件由于是特殊的二進(jìn)制格式箕戳，因此沒辦法通過文本的形式直接打開，要用samtools的view功能在終端上進(jìn)行查看（上文也已經(jīng)說到這里在進(jìn)行系統(tǒng)補(bǔ)充）国撵，如：
$ samtools view in.bam
如果不想從頭開始看陵吸，希望快速地跳轉(zhuǎn)到基因組的其它位置上，比如chr22染色體介牙，那么可以先用samtools index生成BAM文件的索引(如果已經(jīng)有索引文件則不需該步驟)壮虫，然后這樣操作：
$ samtools index in.bam # 生成in.bam的索引文件in.bam.bai
$ samtools view in.bam chr22 # 跳轉(zhuǎn)到chr22染色體
$ samtools view in.bam chr22:16050103 # 跳轉(zhuǎn)到chr22:16050103位置
$ samtools view in.bam chr22:16050103-16050103 # 只查看該位置
IGV或者samtools tview查看比對(duì)情況
以上，我基本上列舉了我們會(huì)在終端上如何查看BAM文件的幾個(gè)最常用操作环础。但如果你想更直觀查看的BAM文件旨指，IGV是目前最好的一個(gè)選擇，但僅適合于文件還比較小的情況喳整，效果如下：
如果你的BAM文件很大谆构，都超過了你的本地電腦磁盤了，你還是想看該怎么辦框都？你有兩個(gè)選擇：
第一搬素，把你想查看的那部分區(qū)域用samtools view提取出來，生成一份小一些的BAM魏保，然后下載下來熬尺，在導(dǎo)入到IGV中。
$ samtools view -h in.bam chr22:16050103-16050203 | samtools view -Sb - > small.bam
第二谓罗，不下載粱哼，直接在終端用samtools tview進(jìn)行查看。samtools tview有類似于IGV的功能檩咱，雖然體驗(yàn)會(huì)稍差一些揭措。
$ samtools tview --reference hg38.fa in.bam
在該模式下胯舷，按下鍵盤‘g’后，會(huì)跳出一個(gè)Goto框绊含，在里面輸入想要調(diào)整過去的位置桑嘶，就行了，比如：
按下esc鍵則可以取消躬充。另外逃顶，為了節(jié)省空間，加快查詢效率充甚，read中與參考序列相同的部分被用一串串不同顏色的點(diǎn)表示以政，只留下miss-match的堿基和發(fā)生indel變異的區(qū)域。其中圓點(diǎn)表示正鏈比對(duì)伴找，逗號(hào)表示負(fù)鏈比對(duì)盈蛮。不同的顏色代表不同的比對(duì)質(zhì)量值：白色>=30，黃色20-29疆瑰，綠色10-19，藍(lán)色0-9昙啄。如果你還想知道的其他的功能穆役，可以在tview模式里按下“?”問號(hào)，就會(huì)彈出類似下面這樣的幫助窗口梳凛，然后按照指引做就行了耿币。
雖然看起來不如IGV體驗(yàn)?zāi)菢雍茫δ芤脖容^單一（僅可以查看比對(duì)情況）韧拒，但可貴之處在于可以在終端里面直接操作淹接，當(dāng)需要快速查看某個(gè)位置的比對(duì)情況時(shí)，操作效率非常高叛溢。而如果要退出該模式塑悼，也非常簡單，按下q鍵就可以了楷掉。
小結(jié)
那么厢蒜，有關(guān)BAM格式的內(nèi)容我們就暫且先到這里吧，大家如果有疑惑或者感興趣的內(nèi)容都可以到后臺(tái)留言烹植，我都會(huì)定時(shí)進(jìn)行回復(fù)斑鸦。在下一篇文章中，我們將重點(diǎn)介紹如何使用pysam來操作bam文件了草雕。
