被忽视的Samtools参数

来源:tiramisutes评论7,307

Samtools是一个用于操作序列比对结果sam和bam文件的工具合集。

sam文件格式

SAM格式由两部分组成:头部区和比对区,都以tab分列。

头部区:以’@’开始,体现了比对的一些总体信息。比对的SAM格式版本,比对的参考序列,比对使用的软件等。

比对区: 比对结果,每一个比对结果是一行,有11个主列和1个可选列

@HD VN:1.0 SO:unsorted  
头部区第一行:VN是格式版本;SO表示比对排序的类型,有unkown(default),unsorted,
queryname和coordinate几种。samtools软件在进行排序后不能自动更新bam文件的SO值。
picard却可以。
@SQ SN:A.auricula_all_contig_1 LN:9401
参考序列名。这些参考序列决定了比对结果sort的顺序。SN是参考序列名;LN是参考序列
长度;
@RG ID:sample01
Read Group. 1个sample的测序结果为1个Read Group;该sample可以有多个library
的测序结果。
@PG ID:bowtie2 PN:bowtie2 VN:2.0.0-beta7
比对所使用的软件。

比对区11个列和可选列的解释
1  QNAME  比对的序列名
2  FLAG   Bwise FLAG(表明比对类型:pairing,strand,mate strand等)
3  RNAME  比对上的参考序列名
4  POS    1-Based的比对上的最左边的定位
5  MAPQ   比对质量,255表示没有map
6  CIGAR  Extended CIGAR string (操作符:MIDNSHP) 比对结果信息:匹配碱基数,可变剪接等,*表示不可用。
7  MRNM   相匹配的另外一条序列,比对上的参考序列名
8  MPOS   1-Based leftmost Mate POsition
9  ISIZE  插入片段长度
10 SEQ    和参考序列在同一个琏上的比对序列(若比对结果在负意链上,则序列是其反向重复序列)
11 QUAL   比对序列的质量(ASCII-33=Phred base quality)
12 可选的行,以TAG:TYPE:VALUE的形式提供额外的信息

 

比对区解释

sam/bam比对区包含有此次比对的结果信息,其中主要信息解释如下:

 

  • FLAG部分

被忽视的Samtools参数-图片1

0x800 表明相应位置的比对属于嵌合体比对;

0x4 没有map上的reads;

  • CIGAR部分

被忽视的Samtools参数-图片2

对于mRNA到基因组的比对,N表示内含子。

More: http://samtools.github.io/hts-specs/SAMv1.pdf

 

view

-c 计数

-f 返回指定区间/flags比对结果

-q 返回比对质量大于等于指定值的比对数目

-F 4:统计map 上的 reads总数;

-f 4:统计没有map 上的 reads总数;

To get the unmapped reads from a bam file use :

samtools view -f 4 file.bam > unmapped.sam

the output will be in sam

to get the output in bam use :

samtools view -b -f 4 file.bam > unmapped.bam

To get only the mapped reads use the parameter ‘F’, which works like -v of grep and skips the alignments for a specific flag.

samtools view -b -F 4 file.bam > mapped.bam
samtools view -b -F 4 -f 8 file.bam > onlyThisEndMapped.bam
samtools view -b -F 8 -f 4 file.bam > onlyThatEndMapped.bam
samtools view -b -F12 file.bam > bothEndsMapped.bam
samtools merge merged.bam onlyThisEndMapped.bam onlyThatEndMapped.bam bothEndsMapped.bam

 

对于tophat比对结果:

samtools view -b -f 2 accepted_hits.bam > mappedPairs.bam

Better with:

samtools view -b -f 0x2 accepted_hits.bam > mappedPairs.bam

sort

-m 指定运算内存,支持K,M,G等缩写
-@ 并行运算核数

index

必须对bam文件进行默认情况下的排序后,才能进行index。否则会报错。

建立索引后将产生后缀为.bai的文件,用于快速的随机处理。很多情况下需要有bai文件的存在,特别是显示序列比对情况下。比如samtool的tview命令就需要;gbrowse2显示reads的比对图形的时候也需要。

faidx

对fasta文件建立索引,生成的索引文件以.fai后缀结尾。该命令也能依据索引文件快速提取fasta文件中的某一条(子)序列。

See more: bedtools 使用小结

flagstat

给出BAM文件的比对结果

$ samtools flagstat example.bam
11945742 + 0 in total (QC-passed reads + QC-failed reads)
#总共的reads数
0 + 0 duplicates
7536364 + 0 mapped (63.09%:-nan%)
#总体上reads的匹配率
11945742 + 0 paired in sequencing
#有多少reads是属于paired reads
5972871 + 0 read1
#reads1中的reads数
5972871 + 0 read2
#reads2中的reads数
6412042 + 0 properly paired (53.68%:-nan%)
#完美匹配的reads数:比对到同一条参考序列,并且两条reads之间的距离符合设置的阈值
6899708 + 0 with itself and mate mapped
#paired reads中两条都比对到参考序列上的reads数
636656 + 0 singletons (5.33%:-nan%)
#单独一条匹配到参考序列上的reads数,和上一个相加,则是总的匹配上的reads数。
469868 + 0 with mate mapped to a different chr
#paired reads中两条分别比对到两条不同的参考序列的reads数
243047 + 0 with mate mapped to a different chr (mapQ>=5)

mpileup

samtools还有个非常重要的命令mpileup,以前为pileup。该命令用于生成bcf文件,再使用bcftools进行SNP和Indel的分析。bcftools是samtool中附带的软件,在samtools的安装文件夹中可以找到。

-f 来输入有索引文件的fasta参考序列;

-g 输出到bcf格式。

贡献来源

http://www.plob.org/2014/01/26/7112.html

http://blog.sina.com.cn/s/blog_670445240101l30k.html

https://www.biostars.org/p/56246/

https://www.biostars.org/p/110039/

https://www.biostars.org/p/95929/

https://www.biostars.org/p/110157/

https://groups.google.com/forum/#!forum/bedtools-discuss

https://code.google.com/p/hydra-sv/

发表评论

匿名网友