1 .TH samtools 1 "27 October 2010" "samtools-0.1.9" "Bioinformatics tools"
4 samtools - Utilities for the Sequence Alignment/Map (SAM) format
7 samtools view -bt ref_list.txt -o aln.bam aln.sam.gz
9 samtools sort aln.bam aln.sorted
11 samtools index aln.sorted.bam
13 samtools idxstats aln.sorted.bam
15 samtools view aln.sorted.bam chr2:20,100,000-20,200,000
17 samtools merge out.bam in1.bam in2.bam in3.bam
19 samtools faidx ref.fasta
21 samtools pileup -vcf ref.fasta aln.sorted.bam
23 samtools mpileup -C50 -agf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam
25 samtools tview aln.sorted.bam ref.fasta
29 Samtools is a set of utilities that manipulate alignments in the BAM
30 format. It imports from and exports to the SAM (Sequence Alignment/Map)
31 format, does sorting, merging and indexing, and allows to retrieve reads
32 in any regions swiftly.
34 Samtools is designed to work on a stream. It regards an input file `-'
35 as the standard input (stdin) and an output file `-' as the standard
36 output (stdout). Several commands can thus be combined with Unix
37 pipes. Samtools always output warning and error messages to the standard
38 error output (stderr).
40 Samtools is also able to open a BAM (not SAM) file on a remote FTP or
41 HTTP server if the BAM file name starts with `ftp://' or `http://'.
42 Samtools checks the current working directory for the index file and
43 will download the index upon absence. Samtools does not retrieve the
44 entire alignment file unless it is asked to do so.
46 .SH COMMANDS AND OPTIONS
50 samtools view [-bhuHS] [-t in.refList] [-o output] [-f reqFlag] [-F
51 skipFlag] [-q minMapQ] [-l library] [-r readGroup] [-R rgFile] <in.bam>|<in.sam> [region1 [...]]
53 Extract/print all or sub alignments in SAM or BAM format. If no region
54 is specified, all the alignments will be printed; otherwise only
55 alignments overlapping the specified regions will be output. An
56 alignment may be given multiple times if it is overlapping several
57 regions. A region can be presented, for example, in the following
58 format: `chr2' (the whole chr2), `chr2:1000000' (region starting from
59 1,000,000bp) or `chr2:1,000,000-2,000,000' (region between 1,000,000 and
60 2,000,000bp including the end points). The coordinate is 1-based.
66 Output in the BAM format.
69 Only output alignments with all bits in INT present in the FLAG
70 field. INT can be in hex in the format of /^0x[0-9A-F]+/ [0]
73 Skip alignments with bits present in INT [0]
76 Include the header in the output.
79 Output the header only.
82 Only output reads in library STR [null]
88 Skip alignments with MAPQ smaller than INT [0]
91 Only output reads in read group STR [null]
94 Output reads in read groups listed in
99 Input is in SAM. If @SQ header lines are absent, the
104 This file is TAB-delimited. Each line must contain the reference name
105 and the length of the reference, one line for each distinct reference;
106 additional fields are ignored. This file also defines the order of the
107 reference sequences in sorting. If you run `samtools faidx <ref.fa>',
108 the resultant index file
115 Output uncompressed BAM. This option saves time spent on
116 compression/decomprssion and is thus preferred when the output is piped
117 to another samtools command.
122 samtools tview <in.sorted.bam> [ref.fasta]
124 Text alignment viewer (based on the ncurses library). In the viewer,
125 press `?' for help and press `g' to check the alignment start from a
126 region in the format like `chr10:10,000,000' or `=10,000,000' when
127 viewing the same reference sequence.
131 samtools pileup [-2sSBicv] [-f in.ref.fasta] [-t in.ref_list] [-l
132 in.site_list] [-C capMapQ] [-M maxMapQ] [-T theta] [-N nHap] [-r
133 pairDiffRate] [-m mask] [-d maxIndelDepth] [-G indelPrior]
136 Print the alignment in the pileup format. In the pileup format, each
137 line represents a genomic position, consisting of chromosome name,
138 coordinate, reference base, read bases, read qualities and alignment
139 mapping qualities. Information on match, mismatch, indel, strand,
140 mapping quality and start and end of a read are all encoded at the read
141 base column. At this column, a dot stands for a match to the reference
142 base on the forward strand, a comma for a match on the reverse strand,
143 a '>' or '<' for a reference skip, `ACGTN' for a mismatch on the forward
144 strand and `acgtn' for a mismatch on the reverse strand. A pattern
145 `\\+[0-9]+[ACGTNacgtn]+' indicates there is an insertion between this
146 reference position and the next reference position. The length of the
147 insertion is given by the integer in the pattern, followed by the
148 inserted sequence. Similarly, a pattern `-[0-9]+[ACGTNacgtn]+'
149 represents a deletion from the reference. The deleted bases will be
150 presented as `*' in the following lines. Also at the read base column, a
151 symbol `^' marks the start of a read. The ASCII of the character
152 following `^' minus 33 gives the mapping quality. A symbol `$' marks the
153 end of a read segment.
157 is applied, the consensus base, Phred-scaled consensus quality, SNP
158 quality (i.e. the Phred-scaled probability of the consensus being
159 identical to the reference) and root mean square (RMS) mapping quality
160 of the reads covering the site will be inserted between the `reference
161 base' and the `read bases' columns. An indel occupies an additional
162 line. Each indel line consists of chromosome name, coordinate, a star,
163 the genotype, consensus quality, SNP quality, RMS mapping quality, #
164 covering reads, the first alllele, the second allele, # reads supporting
165 the first allele, # reads supporting the second allele and # reads
166 containing indels different from the top two alleles.
168 The position of indels is offset by -1.
174 Disable the BAQ computation. See the
179 Call the consensus sequence using SOAPsnp consensus model. Options
180 .BR -T ", " -N ", " -I " and " -r
181 are only effective when
186 Coefficient for downgrading the mapping quality of poorly mapped
189 command for details. [0]
194 reads in the pileup for indel calling for speed up. Zero for unlimited. [1024]
197 The reference sequence in the FASTA format. Index file
203 Generate genotype likelihood in the binary GLFv3 format. This option
204 suppresses -c, -i and -s. This option is deprecated by the
209 Only output pileup lines containing indels.
212 Phred probability of an indel in sequencing/prep. [40]
215 List of sites at which pileup is output. This file is space
216 delimited. The first two columns are required to be chromosome and
217 1-based coordinate. Additional columns are ignored. It is
218 recommended to use option
221 Filter reads with flag containing bits in
226 Cap mapping quality at INT [60]
229 Number of haplotypes in the sample (>=2) [2]
232 Expected fraction of differences between a pair of haplotypes [0.001]
235 Print the mapping quality as the last column. This option makes the
236 output easier to parse, although this format is not space efficient.
239 The input file is in SAM.
242 List of reference names ane sequence lengths, in the format described
245 command. If this option is present, samtools assumes the input
247 is in SAM format; otherwise it assumes in BAM format.
251 as in the default format we may not know the mapping quality.
254 The theta parameter (error dependency coefficient) in the maq consensus
260 samtools mpileup [-Bug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]
262 Generate BCF or pileup for one or multiple BAM files. Alignment records
263 are grouped by sample identifiers in @RG header lines. If sample
264 identifiers are absent, each input file is regarded as one sample.
270 Disable probabilistic realignment for the computation of base alignment
271 quality (BAQ). BAQ is the Phred-scaled probability of a read base being
272 misaligned. Applying this option greatly helps to reduce false SNPs
273 caused by misalignments.
276 Coefficient for downgrading mapping quality for reads containing
277 excessive mismatches. Given a read with a phred-scaled probability q of
278 being generated from the mapped position, the new mapping quality is
279 about sqrt((INT-q)/INT)*INT. A zero value disables this
280 functionality; if enabled, the recommended value is 50. [0]
283 The reference file [null]
286 Compute genotype likelihoods and output them in the binary call format (BCF).
291 except that the output is uncompressed BCF, which is preferred for pipeing.
294 File containing a list of sites where pileup or BCF is outputted [null]
297 Minimum mapping quality for an alignment to be used [0]
300 Minimum base quality for a base to be considered [13]
303 Only generate pileup in region
310 samtools reheader <in.header.sam> <in.bam>
312 Replace the header in
316 This command is much faster than replacing the header with a
317 BAM->SAM->BAM conversion.
321 samtools sort [-no] [-m maxMem] <in.bam> <out.prefix>
323 Sort alignments by leftmost coordinates. File
325 will be created. This command may also create temporary files
326 .I <out.prefix>.%d.bam
327 when the whole alignment cannot be fitted into memory (controlled by
334 Output the final alignment to the standard output.
337 Sort by read names rather than by chromosomal coordinates
340 Approximately the maximum required memory. [500000000]
345 samtools merge [-nur] [-h inh.sam] [-R reg] <out.bam> <in1.bam> <in2.bam> [...]
347 Merge multiple sorted alignments.
348 The header reference lists of all the input BAM files, and the @SQ headers of
350 if any, must all refer to the same set of reference sequences.
351 The header reference list and (unless overridden by
357 and the headers of other files will be ignored.
365 as `@' headers to be copied to
367 replacing any header lines that would otherwise be copied from
370 is actually in SAM format, though any alignment records it may contain
374 Merge files in the specified region indicated by
378 Attach an RG tag to each alignment. The tag value is inferred from file names.
381 The input alignments are sorted by read names rather than by chromosomal
385 Uncompressed BAM output
390 samtools index <aln.bam>
392 Index sorted alignment for fast random access. Index file
398 samtools idxstats <aln.bam>
400 Retrieve and print stats in the index file. The output is TAB delimited
401 with each line consisting of reference sequence name, sequence length, #
402 mapped reads and # unmapped reads.
406 samtools faidx <ref.fasta> [region1 [...]]
408 Index reference sequence in the FASTA format or extract subsequence from
409 indexed reference sequence. If no region is specified,
411 will index the file and create
413 on the disk. If regions are speficified, the subsequences will be
414 retrieved and printed to stdout in the FASTA format. The input file can
421 samtools fixmate <in.nameSrt.bam> <out.bam>
423 Fill in mate coordinates, ISIZE and mate related flags from a
424 name-sorted alignment.
428 samtools rmdup [-sS] <input.srt.bam> <out.bam>
430 Remove potential PCR duplicates: if multiple read pairs have identical
431 external coordinates, only retain the pair with highest mapping quality.
432 In the paired-end mode, this command
434 works with FR orientation and requires ISIZE is correctly set. It does
435 not work for unpaired reads (e.g. two ends mapped to different
436 chromosomes or orphan reads).
442 Remove duplicate for single-end reads. By default, the command works for
443 paired-end reads only.
446 Treat paired-end reads and single-end reads.
451 samtools calmd [-eubSr] [-C capQcoef] <aln.bam> <ref.fasta>
453 Generate the MD tag. If the MD tag is already present, this command will
454 give a warning if the MD tag generated is different from the existing
455 tag. Output SAM by default.
461 Convert a the read base to = if it is identical to the aligned reference
462 base. Indel caller does not support the = bases at the moment.
465 Output uncompressed BAM
468 Output compressed BAM
471 The input is SAM with header lines
474 Coefficient to cap mapping quality of poorly mapped reads. See the
476 command for details. [0]
479 Perform probabilistic realignment to compute BAQ, which will be used to
485 SAM is TAB-delimited. Apart from the header lines, which are started
486 with the `@' symbol, each alignment line consists of:
492 Col Field Description
494 1 QNAME Query (pair) NAME
496 3 RNAME Reference sequence NAME
497 4 POS 1-based leftmost POSition/coordinate of clipped sequence
498 5 MAPQ MAPping Quality (Phred-scaled)
499 6 CIAGR extended CIGAR string
500 7 MRNM Mate Reference sequence NaMe (`=' if same as RNAME)
501 8 MPOS 1-based Mate POSistion
502 9 ISIZE Inferred insert SIZE
503 10 SEQ query SEQuence on the same strand as the reference
504 11 QUAL query QUALity (ASCII-33 gives the Phred base quality)
505 12 OPT variable OPTional fields in the format TAG:VTYPE:VALUE
509 Each bit in the FLAG field is defined as:
517 0x0001 p the read is paired in sequencing
518 0x0002 P the read is mapped in a proper pair
519 0x0004 u the query sequence itself is unmapped
520 0x0008 U the mate is unmapped
521 0x0010 r strand of the query (1 for reverse)
522 0x0020 R strand of the mate
523 0x0040 1 the read is the first read in a pair
524 0x0080 2 the read is the second read in a pair
525 0x0100 s the alignment is not primary
526 0x0200 f the read fails platform/vendor quality checks
527 0x0400 d the read is either a PCR or an optical duplicate
532 Import SAM to BAM when
534 lines are present in the header:
536 samtools view -bS aln.sam > aln.bam
542 samtools faidx ref.fa
543 samtools view -bt ref.fa.fai aln.sam > aln.bam
547 is generated automatically by the
554 tag while merging sorted alignments:
556 perl -e 'print "@RG\\tID:ga\\tSM:hs\\tLB:ga\\tPL:Illumina\\n@RG\\tID:454\\tSM:hs\\tLB:454\\tPL:454\\n"' > rg.txt
557 samtools merge -rh rg.txt merged.bam ga.bam 454.bam
561 tag is determined by the file name the read is coming from. In this
574 Call SNPs and short indels for one diploid individual:
576 samtools pileup -vcf ref.fa aln.bam > var.raw.plp
577 samtools.pl varFilter -D 100 var.raw.plp > var.flt.plp
578 awk '($3=="*"&&$6>=50)||($3!="*"&&$6>=20)' var.flt.plp > var.final.plp
582 option of varFilter controls the maximum read depth, which should be
583 adjusted to about twice the average read depth. One may consider to add
587 if mapping quality is overestimated for reads containing excessive
588 mismatches. Applying this option usually helps
590 but may not other mappers. It also potentially increases reference
594 Call SNPs (not short indels) for multiple diploid individuals:
596 samtools mpileup -augf ref.fa *.bam | bcftools view -bcv - > snp.raw.bcf
597 bcftools view snp.raw.bcf | vcfutils.pl filter4vcf -D 2000 | bgzip > snp.flt.vcf.gz
599 Individuals are identified from the
603 header lines. Individuals can be pooled in one alignment file; one
604 individual can also be separated into multiple files. Similarly, one may
609 SNP calling in this way also works for single sample and has the
610 advantage of enabling more powerful filtering. The drawback is the lack
611 of short indel calling, which may be implemented in future.
614 Derive the allele frequency spectrum (AFS) on a list of sites from multiple individuals:
616 samtools mpileup -gf ref.fa *.bam > all.bcf
617 bcftools view -bl sites.list all.bcf > sites.bcf
618 bcftools view -cGP cond2 sites.bcf > /dev/null 2> sites.1.afs
619 bcftools view -cGP sites.1.afs sites.bcf > /dev/null 2> sites.2.afs
620 bcftools view -cGP sites.2.afs sites.bcf > /dev/null 2> sites.3.afs
625 contains the list of sites with each line consisting of the reference
626 sequence name and position. The following
628 commands estimate AFS by EM.
631 Dump BAQ applied alignment for other SNP callers:
633 samtools calmd -br aln.bam > aln.baq.bam
635 It adds and corrects the
639 tags at the same time. The
641 command also comes with the
643 option, the same as the on in
652 Unaligned words used in bam_import.c, bam_endian.h, bam.c and bam_aux.c.
654 In merging, the input files are required to have the same number of
655 reference sequences. The requirement can be relaxed. In addition,
656 merging does not reconstruct the header dictionaries
657 automatically. Endusers have to provide the correct header. Picard is
660 Samtools paired-end rmdup does not work for unpaired reads (e.g. orphan
661 reads or ends mapped to different chromosomes). If this is a concern,
662 please use Picard's MarkDuplicate which correctly handles these cases,
663 although a little slower.
667 Heng Li from the Sanger Institute wrote the C version of samtools. Bob
668 Handsaker from the Broad Institute implemented the BGZF library and Jue
669 Ruan from Beijing Genomics Institute wrote the RAZF library. Various
670 people in the 1000 Genomes Project contributed to the SAM format
675 Samtools website: <http://samtools.sourceforge.net>
projects / samtools.git / blob