Dean's blog
0%

samtools-view-简单介绍

Posted on
Symbols count in article: 2k Reading time ≈ 2 mins.

整理高通量测序数据分析中使用的软件。

samtools view

samtools是常用的对sam/bam文件操作的工具，其中samtools view命令可以实现查看序列、sam-bam文件转换、过滤序列等功能。下面用实际的例子简单介绍个人使用samtools view过程中的一些经验。

查看序列

如果只是查看sam文件的序列比对结果的前几行，可以用该命令简单的查看.

1
2
3
4
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

默认情况，samtools view输出序列到屏幕，可以重导向到别的文件。

1
samtools view test.sam > test2.sam

sam-bam文件格式转换

1
samtools view -@ 16 -bh test.sam -o test.bam

-b：声明输出为bam格式的文件。
-h：保留sam文件的header（如果有的话），header信息常包括比对的参考基因组的染色体信息和比对的命令
-@ ：指明使用的线程数

过滤序列

以下简单介绍samtools view 过滤序列的参数

根据MAPQ过滤

1
samtools view -@ 8 -bh -q 30 test.bam -o test.q30.bam

-q参数只输出MAPQ（比对质量）大于等于该值的序列，上述命令过滤了MAPQ小于30的序列

根据FLAG过滤

另外，通过-f和-F参数，我们可以根据FLAG的值过滤序列。两者的区别在于：
-f：保留该flag值的序列（相当于grep ）
-F：保留除了该flag值以外的所有序列（相当于grep -v）

对于单端测序的比对结果

因此，假设我们需要取比对上的reads，可以使用-F 4

1
samtools view -@ 8 -bh -F 4 test.bam -o test.mapped.bam

或者取出比对不上的reads（前提是比对软件也输出了比对不上的reads）

1
samtools view -@ 8 -bh -f 4 test.bam -o test.unmapped.bam

对于双端测序的比对结果

1
2
samtools view -@ 8  -bh -f 2 test.bam > test.mapped.bam
samtools view -@ 8  -bh -F 2 test.bam > test.unmapped.bam

这里提供一个小工具可以快速查询flag值所对应的含义 (https://broadinstitute.github.io/picard/explain-flags.html)

在左下侧的框中勾选会给出相应的flag值
在上方搜索栏中输入数值会在右侧给出flag值对应的含义。

根据染色体坐标进行过滤

如果bam文件已经使用samtools index建好index的话，可以输出特定染色体坐标内的reads

1
samtools view -@ 8 -b test.bam chr1:10420000-10421000 > subset.bam

如果想取出多个染色体区域的reads的话，就不再建议使用上述的方法了，可以使用bedtools之类的工具根据bed文件进行提取。对samtools view命令的简单介绍就到此结束，以后使用有心得再作更新。

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/

完。