1 .TH samtools 1 "21 April 2011" "samtools-0.1.16" "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 -gf 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 [-bchuHS] [-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 Instead of printing the alignments, only count them and print the
105 total number. All filter options, such as
110 , are taken into account.
113 This file is TAB-delimited. Each line must contain the reference name
114 and the length of the reference, one line for each distinct reference;
115 additional fields are ignored. This file also defines the order of the
116 reference sequences in sorting. If you run `samtools faidx <ref.fa>',
117 the resultant index file
124 Output uncompressed BAM. This option saves time spent on
125 compression/decomprssion and is thus preferred when the output is piped
126 to another samtools command.
131 samtools tview <in.sorted.bam> [ref.fasta]
133 Text alignment viewer (based on the ncurses library). In the viewer,
134 press `?' for help and press `g' to check the alignment start from a
135 region in the format like `chr10:10,000,000' or `=10,000,000' when
136 viewing the same reference sequence.
140 samtools mpileup [-EBug] [-C capQcoef] [-r reg] [-f in.fa] [-l list] [-M capMapQ] [-Q minBaseQ] [-q minMapQ] in.bam [in2.bam [...]]
142 Generate BCF or pileup for one or multiple BAM files. Alignment records
143 are grouped by sample identifiers in @RG header lines. If sample
144 identifiers are absent, each input file is regarded as one sample.
150 Do not skip anomalous read pairs in variant calling.
153 Disable probabilistic realignment for the computation of base alignment
154 quality (BAQ). BAQ is the Phred-scaled probability of a read base being
155 misaligned. Applying this option greatly helps to reduce false SNPs
156 caused by misalignments.
159 Coefficient for downgrading mapping quality for reads containing
160 excessive mismatches. Given a read with a phred-scaled probability q of
161 being generated from the mapped position, the new mapping quality is
162 about sqrt((INT-q)/INT)*INT. A zero value disables this
163 functionality; if enabled, the recommended value for BWA is 50. [0]
166 At a position, read maximally
168 reads per input BAM. [250]
171 Output per-sample read depth
174 Phred-scaled gap extension sequencing error probability. Reducing
176 leads to longer indels. [20]
179 Extended BAQ computation. This option helps sensitivity especially for MNPs, but may hurt
180 specificity a little bit.
183 The reference file [null]
186 Compute genotype likelihoods and output them in the binary call format (BCF).
189 Coefficient for modeling homopolymer errors. Given an
192 run, the sequencing error of an indel of size
199 Do not perform INDEL calling
202 File containing a list of sites where pileup or BCF is outputted [null]
205 Skip INDEL calling if the average per-sample depth is above
210 Phred-scaled gap open sequencing error probability. Reducing
212 leads to more indel calls. [40]
215 Comma dilimited list of platforms (determined by
217 from which indel candidates are obtained. It is recommended to collect
218 indel candidates from sequencing technologies that have low indel error
219 rate such as ILLUMINA. [all]
222 Minimum mapping quality for an alignment to be used [0]
225 Minimum base quality for a base to be considered [13]
228 Only generate pileup in region
233 Output per-sample Phred-scaled strand bias P-value
238 except that the output is uncompressed BCF, which is preferred for piping.
243 samtools reheader <in.header.sam> <in.bam>
245 Replace the header in
249 This command is much faster than replacing the header with a
250 BAM->SAM->BAM conversion.
254 samtools cat [-h header.sam] [-o out.bam] <in1.bam> <in2.bam> [ ... ]
256 Concatenate BAMs. The sequence dictionary of each input BAM must be identical,
257 although this command does not check this. This command uses a similar trick
260 which enables fast BAM concatenation.
264 samtools sort [-no] [-m maxMem] <in.bam> <out.prefix>
266 Sort alignments by leftmost coordinates. File
268 will be created. This command may also create temporary files
269 .I <out.prefix>.%d.bam
270 when the whole alignment cannot be fitted into memory (controlled by
277 Output the final alignment to the standard output.
280 Sort by read names rather than by chromosomal coordinates
283 Approximately the maximum required memory. [500000000]
288 samtools merge [-nur1f] [-h inh.sam] [-R reg] <out.bam> <in1.bam> <in2.bam> [...]
290 Merge multiple sorted alignments.
291 The header reference lists of all the input BAM files, and the @SQ headers of
293 if any, must all refer to the same set of reference sequences.
294 The header reference list and (unless overridden by
300 and the headers of other files will be ignored.
306 Use zlib compression level 1 to comrpess the output
309 Force to overwrite the output file if present.
314 as `@' headers to be copied to
316 replacing any header lines that would otherwise be copied from
319 is actually in SAM format, though any alignment records it may contain
323 The input alignments are sorted by read names rather than by chromosomal
327 Merge files in the specified region indicated by
332 Attach an RG tag to each alignment. The tag value is inferred from file names.
335 Uncompressed BAM output
340 samtools index <aln.bam>
342 Index sorted alignment for fast random access. Index file
348 samtools idxstats <aln.bam>
350 Retrieve and print stats in the index file. The output is TAB delimited
351 with each line consisting of reference sequence name, sequence length, #
352 mapped reads and # unmapped reads.
356 samtools faidx <ref.fasta> [region1 [...]]
358 Index reference sequence in the FASTA format or extract subsequence from
359 indexed reference sequence. If no region is specified,
361 will index the file and create
363 on the disk. If regions are speficified, the subsequences will be
364 retrieved and printed to stdout in the FASTA format. The input file can
371 samtools fixmate <in.nameSrt.bam> <out.bam>
373 Fill in mate coordinates, ISIZE and mate related flags from a
374 name-sorted alignment.
378 samtools rmdup [-sS] <input.srt.bam> <out.bam>
380 Remove potential PCR duplicates: if multiple read pairs have identical
381 external coordinates, only retain the pair with highest mapping quality.
382 In the paired-end mode, this command
384 works with FR orientation and requires ISIZE is correctly set. It does
385 not work for unpaired reads (e.g. two ends mapped to different
386 chromosomes or orphan reads).
392 Remove duplicate for single-end reads. By default, the command works for
393 paired-end reads only.
396 Treat paired-end reads and single-end reads.
401 samtools calmd [-EeubSr] [-C capQcoef] <aln.bam> <ref.fasta>
403 Generate the MD tag. If the MD tag is already present, this command will
404 give a warning if the MD tag generated is different from the existing
405 tag. Output SAM by default.
411 When used jointly with
413 this option overwrites the original base quality.
416 Convert a the read base to = if it is identical to the aligned reference
417 base. Indel caller does not support the = bases at the moment.
420 Output uncompressed BAM
423 Output compressed BAM
426 The input is SAM with header lines
429 Coefficient to cap mapping quality of poorly mapped reads. See the
431 command for details. [0]
434 Compute the BQ tag (without -A) or cap base quality by BAQ (with -A).
437 Extended BAQ calculation. This option trades specificity for sensitivity, though the
443 samtools targetcut [-Q minBaseQ] [-i inPenalty] [-0 em0] [-1 em1] [-2 em2] [-f ref] <in.bam>
445 This command identifies target regions by examining the continuity of read depth, computes
446 haploid consensus sequences of targets and outputs a SAM with each sequence corresponding
447 to a target. When option
449 is in use, BAQ will be applied. This command is
451 designed for cutting fosmid clones from fosmid pool sequencing [Ref. Kitzman et al. (2010)].
456 samtools phase [-AF] [-k len] [-b prefix] [-q minLOD] [-Q minBaseQ] <in.bam>
458 Call and phase heterozygous SNPs.
463 Drop reads with ambiguous phase.
466 Prefix of BAM output. When this option is in use, phase-0 reads will be saved in file
470 Phase unknown reads will be randomly allocated to one of the two files. Chimeric reads
471 with switch errors will be saved in
472 .BR STR .chimeric.bam.
476 Do not attempt to fix chimeric reads.
479 Maximum length for local phasing. [13]
482 Minimum Phred-scaled LOD to call a heterozygote. [40]
485 Minimum base quality to be used in het calling. [13]
490 samtools pileup [-2sSBicv] [-f in.ref.fasta] [-t in.ref_list] [-l
491 in.site_list] [-C capMapQ] [-M maxMapQ] [-T theta] [-N nHap] [-r
492 pairDiffRate] [-m mask] [-d maxIndelDepth] [-G indelPrior]
495 Print the alignment in the pileup format. In the pileup format, each
496 line represents a genomic position, consisting of chromosome name,
497 coordinate, reference base, read bases, read qualities and alignment
498 mapping qualities. Information on match, mismatch, indel, strand,
499 mapping quality and start and end of a read are all encoded at the read
500 base column. At this column, a dot stands for a match to the reference
501 base on the forward strand, a comma for a match on the reverse strand,
502 a '>' or '<' for a reference skip, `ACGTN' for a mismatch on the forward
503 strand and `acgtn' for a mismatch on the reverse strand. A pattern
504 `\\+[0-9]+[ACGTNacgtn]+' indicates there is an insertion between this
505 reference position and the next reference position. The length of the
506 insertion is given by the integer in the pattern, followed by the
507 inserted sequence. Similarly, a pattern `-[0-9]+[ACGTNacgtn]+'
508 represents a deletion from the reference. The deleted bases will be
509 presented as `*' in the following lines. Also at the read base column, a
510 symbol `^' marks the start of a read. The ASCII of the character
511 following `^' minus 33 gives the mapping quality. A symbol `$' marks the
512 end of a read segment.
516 is applied, the consensus base, Phred-scaled consensus quality, SNP
517 quality (i.e. the Phred-scaled probability of the consensus being
518 identical to the reference) and root mean square (RMS) mapping quality
519 of the reads covering the site will be inserted between the `reference
520 base' and the `read bases' columns. An indel occupies an additional
521 line. Each indel line consists of chromosome name, coordinate, a star,
522 the genotype, consensus quality, SNP quality, RMS mapping quality, #
523 covering reads, the first alllele, the second allele, # reads supporting
524 the first allele, # reads supporting the second allele and # reads
525 containing indels different from the top two alleles.
528 Since 0.1.10, the `pileup' command is deprecated by `mpileup'.
534 Disable the BAQ computation. See the
539 Call the consensus sequence. Options
540 .BR -T ", " -N ", " -I " and " -r
541 are only effective when
546 Coefficient for downgrading the mapping quality of poorly mapped
549 command for details. [0]
554 reads in the pileup for indel calling for speed up. Zero for unlimited. [1024]
557 The reference sequence in the FASTA format. Index file
563 Generate genotype likelihood in the binary GLFv3 format. This option
564 suppresses -c, -i and -s. This option is deprecated by the
569 Only output pileup lines containing indels.
572 Phred probability of an indel in sequencing/prep. [40]
575 List of sites at which pileup is output. This file is space
576 delimited. The first two columns are required to be chromosome and
577 1-based coordinate. Additional columns are ignored. It is
578 recommended to use option
581 Filter reads with flag containing bits in
586 Cap mapping quality at INT [60]
589 Number of haplotypes in the sample (>=2) [2]
592 Expected fraction of differences between a pair of haplotypes [0.001]
595 Print the mapping quality as the last column. This option makes the
596 output easier to parse, although this format is not space efficient.
599 The input file is in SAM.
602 List of reference names ane sequence lengths, in the format described
605 command. If this option is present, samtools assumes the input
607 is in SAM format; otherwise it assumes in BAM format.
611 as in the default format we may not know the mapping quality.
614 The theta parameter (error dependency coefficient) in the maq consensus
620 SAM is TAB-delimited. Apart from the header lines, which are started
621 with the `@' symbol, each alignment line consists of:
627 Col Field Description
629 1 QNAME Query (pair) NAME
631 3 RNAME Reference sequence NAME
632 4 POS 1-based leftmost POSition/coordinate of clipped sequence
633 5 MAPQ MAPping Quality (Phred-scaled)
634 6 CIAGR extended CIGAR string
635 7 MRNM Mate Reference sequence NaMe (`=' if same as RNAME)
636 8 MPOS 1-based Mate POSistion
637 9 ISIZE Inferred insert SIZE
638 10 SEQ query SEQuence on the same strand as the reference
639 11 QUAL query QUALity (ASCII-33 gives the Phred base quality)
640 12 OPT variable OPTional fields in the format TAG:VTYPE:VALUE
644 Each bit in the FLAG field is defined as:
652 0x0001 p the read is paired in sequencing
653 0x0002 P the read is mapped in a proper pair
654 0x0004 u the query sequence itself is unmapped
655 0x0008 U the mate is unmapped
656 0x0010 r strand of the query (1 for reverse)
657 0x0020 R strand of the mate
658 0x0040 1 the read is the first read in a pair
659 0x0080 2 the read is the second read in a pair
660 0x0100 s the alignment is not primary
661 0x0200 f the read fails platform/vendor quality checks
662 0x0400 d the read is either a PCR or an optical duplicate
667 Import SAM to BAM when
669 lines are present in the header:
671 samtools view -bS aln.sam > aln.bam
677 samtools faidx ref.fa
678 samtools view -bt ref.fa.fai aln.sam > aln.bam
682 is generated automatically by the
689 tag while merging sorted alignments:
691 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
692 samtools merge -rh rg.txt merged.bam ga.bam 454.bam
696 tag is determined by the file name the read is coming from. In this
709 Call SNPs and short indels for one diploid individual:
711 samtools mpileup -ugf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf
712 bcftools view var.raw.bcf | vcfutils.pl varFilter -D 100 > var.flt.vcf
716 option of varFilter controls the maximum read depth, which should be
717 adjusted to about twice the average read depth. One may consider to add
721 if mapping quality is overestimated for reads containing excessive
722 mismatches. Applying this option usually helps
724 but may not other mappers.
727 Generate the consensus sequence for one diploid individual:
729 samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq
732 Phase one individual:
734 samtools calmd -AEur aln.bam ref.fa | samtools phase -b prefix - > phase.out
738 command is used to reduce false heterozygotes around INDELs.
741 Call SNPs and short indels for multiple diploid individuals:
743 samtools mpileup -P ILLUMINA -ugf ref.fa *.bam | bcftools view -bcvg - > var.raw.bcf
744 bcftools view var.raw.bcf | vcfutils.pl varFilter -D 2000 > var.flt.vcf
746 Individuals are identified from the
750 header lines. Individuals can be pooled in one alignment file; one
751 individual can also be separated into multiple files. The
753 option specifies that indel candidates should be collected only from
758 Collecting indel candidates from reads sequenced by an indel-prone
759 technology may affect the performance of indel calling.
762 Derive the allele frequency spectrum (AFS) on a list of sites from multiple individuals:
764 samtools mpileup -Igf ref.fa *.bam > all.bcf
765 bcftools view -bl sites.list all.bcf > sites.bcf
766 bcftools view -cGP cond2 sites.bcf > /dev/null 2> sites.1.afs
767 bcftools view -cGP sites.1.afs sites.bcf > /dev/null 2> sites.2.afs
768 bcftools view -cGP sites.2.afs sites.bcf > /dev/null 2> sites.3.afs
773 contains the list of sites with each line consisting of the reference
774 sequence name and position. The following
776 commands estimate AFS by EM.
779 Dump BAQ applied alignment for other SNP callers:
781 samtools calmd -bAr aln.bam > aln.baq.bam
783 It adds and corrects the
787 tags at the same time. The
789 command also comes with the
791 option, the same as the one in
800 Unaligned words used in bam_import.c, bam_endian.h, bam.c and bam_aux.c.
802 In merging, the input files are required to have the same number of
803 reference sequences. The requirement can be relaxed. In addition,
804 merging does not reconstruct the header dictionaries
805 automatically. Endusers have to provide the correct header. Picard is
808 Samtools paired-end rmdup does not work for unpaired reads (e.g. orphan
809 reads or ends mapped to different chromosomes). If this is a concern,
810 please use Picard's MarkDuplicate which correctly handles these cases,
811 although a little slower.
815 Heng Li from the Sanger Institute wrote the C version of samtools. Bob
816 Handsaker from the Broad Institute implemented the BGZF library and Jue
817 Ruan from Beijing Genomics Institute wrote the RAZF library. John
818 Marshall and Petr Danecek contribute to the source code and various
819 people from the 1000 Genomes Project have contributed to the SAM format
824 Samtools website: <http://samtools.sourceforge.net>
projects / samtools.git / blob