samtools基本命令实例

Converting a SAM file to a BAM file

First, if you use the Unix command

head test.sam

The first 10 lines on your terminal after typing "head test.sam", should be lines starting with an "@" sign, which is an indicator for a header line. If you don't see lines starting with the "@" sign, the header information is most likely missing.

If the header information is absent from the SAM file use the command below, where reference.fa is the reference fasta file used to map the reads:

samtools view -bT reference.fa test.sam > test.bam

If the header information is available:

samtools view -bS test.sam > test.bam

Sorting a BAM file

samtools sort test.bam test_sorted

Creating a BAM index file

samtools index test_sorted.bam test_sorted.bai

Converting a BAM file to a SAM file

Note: remember to use -h to ensure the converted SAM file contains the header information. Generally, I suggest storing only sorted BAM files as they use much less disk space and are faster to process.

samtools view -h NA06984.chrom16.ILLUMINA.bwa.CEU.low_coverage.20100517.bam > NA06984.chrom16.ILLUMINA.bwa.CEU.low_coverage.20100517.sam

Simple stats

samtools flagstat NA06984.chrom16.ILLUMINA.bwa.CEU.low_coverage.20100517.bam

10182494 in total
0 QC failure
223627 duplicates
9861117 mapped (96.84%)
10095646 paired in sequencing
5049066 read1
5046580 read2
8174084 properly paired (80.97%)
9452892 with itself and mate mapped
321377 singletons (3.18%)
215316 with mate mapped to a different chr
126768 with mate mapped to a different chr (mapQ>=5)
For more statistics of SAM or BAM files have a look at the SAMStat program.

Interpreting the BAM flags

(也可以利用FLAG值含义解释工具:https://www.plob.org/2012/02/04/1697.html

Here are some common BAM flags:

163: 10100011 in binary

147: 10010011 in binary

99: 1100011 in binary

83: 1010011 in binary

Interpretation of 10100011 (reading the binary from left to right):

1the read is paired in sequencing, no matter whether it is mapped in a pair
1the read is mapped in a proper pair (depends on the protocol, normally inferred during alignment)
0the query sequence itself is unmapped
0the mate is unmapped
0strand of the query (0 for forward; 1 for reverse strand)
1strand of the mate
0the read is the first read in a pair
1the read is the second read in a pair

163 second read of a pair on the positive strand with negative strand mate

147 second read of a pair on the negative strand with positive strand mate

99 first read of a pair on the forward strand with negative strand mate

83 first read of a pair on the reverse strand with positive strand mate

Extracting only the first read from paired end BAM files

samtools view -h -f 0x0040 test.bam > test_first_pair.sam

0x0040 is hexadecimal for 64 (i.e. 16 * 4), which is binary for 1000000, corresponding to the read in the first read pair.

Filtering out unmapped reads in BAM files

samtools view -h -F 4 blah.bam > blah_only_mapped.sam

Creating FASTQ files from a BAM file

I found this great tool at http://www.hudsonalpha.org/gsl/software/bam2fastq.php

For example to extract ONLY unaligned from a bam file:

bam2fastq -o blah_unaligned.fastq --no-aligned blah.bam

To extract ONLY aligned reads from a bam file:

bam2fastq -o blah_aligned.fastq --no-unaligned blah.bam

More information

http://samtools.sourceforge.net/

http://sourceforge.net/apps/mediawiki/samtools/index.php?title=SAM_FAQ

  • 本文由 整理发表
  • 网站部分文章源自互联网,若未正确标注来源,请联系管理员更新。文章转载,请务必保留本文链接

您必须才能发表评论!

评论:1   其中:访客  1   博主  0
    • Levin

      【疑问】:看了提取only the first read from paired end BAM files之后我想问一下0x0040是什么?如果要提取second reads呢?
      还有个问题:如果只有一条read map上了(不论是read1还是read2),我想把它对应的没有map上的read(1或者2)提取出来,该用什么参数?