整理高通量測(cè)序中使用的軟件嗦明。
samtools是常用的對(duì)sam/bam文件操作的工具库北，其中samtools view命令可以實(shí)現(xiàn)查看序列朴恳、sam-bam文件轉(zhuǎn)換欺栗、過(guò)濾序列等功能毫痕。下面用實(shí)際的例子簡(jiǎn)單介紹個(gè)人使用samtools view過(guò)程中的一些經(jīng)驗(yàn)。
查看序列
如果只是查看sam文件的序列比對(duì)結(jié)果的前幾行迟几，可以用該命令簡(jiǎn)單的查看.
$samtools view test.sam | head -n2
SRR6924174.1 163 chr11 106049782 255 80M = 106049950 248 GCCAAACACAGAACAAACTTTTTAGAATCTTGGTAATAATGACAACTGACATCAGGGTGGAAGATAAGAATATCTCACTA @BBBCGGGGGEFGGG;EFGGGGGFGGGCGGCDGG>FFEGGGEGGGGBF>1EFB:@GGGGGGGDEGGG@GGDFGGGF1@FC NH:i:1 HI:i:1 AS:i:157 nM:i:0
SRR6924174.1 83 chr11 106049950 255 80M = 106049782 -248 TTGTCTACAGATCAAGTTCCAGGACAGCCAGGAACAGAGAAATGCTGTCTTGAAAACTAAAAAAACACCTTTAATCCCAN GGGGGGGGGGGEGCEGEBFGGBDGFEFGGGGGGGGGEGGGGGFGGGGGGGGGGGEGGFGGGGGGGGFGGGGGGGGBA<<# NH:i:1 HI:i:1 AS:i:157 nM:i:0
默認(rèn)情況消请，samtools view輸出序列到屏幕，可以重導(dǎo)向到別的文件类腮。
$samtools view test.sam > test2.sam
sam-bam文件格式轉(zhuǎn)換
$samtools view -@ 16 -bh test.sam -o test.bam
-b：聲明輸出為bam格式的文件臊泰。
-h：保留sam文件的header（如果有的話），header信息常包括比對(duì)的參考基因組的染色體信息和比對(duì)的命令
-@ ：指明使用的線程數(shù)
過(guò)濾序列
以下簡(jiǎn)單介紹samtools view 過(guò)濾序列的參數(shù)
根據(jù)MAPQ過(guò)濾
$samtools view -@ 8 -bh -q 30 test.bam -o test.q30.bam
-q參數(shù)只輸出MAPQ（比對(duì)質(zhì)量）大于等于該值的序列蚜枢，上述命令過(guò)濾了MAPQ小于30的序列
根據(jù)FLAG過(guò)濾
另外因宇，通過(guò)-f和-F參數(shù)七婴，我們可以根據(jù)FLAG的值過(guò)濾序列。兩者的區(qū)別在于：
-f：保留該flag值的序列（相當(dāng)于grep）
-F：保留除了該flag值以外的所有序列（相當(dāng)于grep -v）
對(duì)于單端測(cè)序的比對(duì)結(jié)果
因此察滑，假設(shè)我們需要取比對(duì)上的reads打厘，可以使用-F 4
$samtools view -@ 8 -bh -F 4 test.bam -o test.mapped.bam
或者取出比對(duì)不上的reads（前提是比對(duì)軟件也輸出了比對(duì)不上的reads）
$samtools view -@ 8 -bh -f 4 test.bam -o test.unmapped.bam
對(duì)于雙端測(cè)序的比對(duì)結(jié)果
$samtools view -@ 8 -bh -f 2 test.bam > test.mapped.bam
$samtools view -@ 8 -bh -F 2 test.bam > test.unmapped.bam
這里提供一個(gè)小工具可以快速查詢flag值所對(duì)應(yīng)的含義，https://broadinstitute.github.io/picard/explain-flags.html
在左下側(cè)的框中勾選會(huì)給出相應(yīng)的flag值
在上方搜索欄中輸入數(shù)值會(huì)在右側(cè)給出flag值對(duì)應(yīng)的含義贺辰。
根據(jù)染色體坐標(biāo)進(jìn)行過(guò)濾
如果bam文件已經(jīng)使用samtools index建好index的話户盯，可以輸出特定染色體坐標(biāo)內(nèi)的reads
samtools view -@ 8 -b test.bam chr1:10420000-10421000 > subset.bam
如果想取出多個(gè)染色體區(qū)域的reads的話，就不再建議使用上述的方法了饲化，可以使用bedtools之類的工具根據(jù)bed文件進(jìn)行提取莽鸭。對(duì)samtools view命令的簡(jiǎn)單介紹就到此結(jié)束，以后使用有心得再作更新吃靠。
ref
samtools view manual: http://www.htslib.org/doc/samtools-view.html
How To Filter Mapped Reads With Samtools: https://www.biostars.org/p/56246/
Decoding SAM flags: https://broadinstitute.github.io/picard/explain-flags.html
Extract Alignment by genomic coordinates: https://www.biostars.org/p/45936/
完硫眨。
