1 .TH samtools 1 "15 March 2013" "samtools-0.1.19" "Bioinformatics tools"
4 samtools - Utilities for the Sequence Alignment/Map (SAM) format
6 bcftools - Utilities for the Binary Call Format (BCF) and VCF
9 samtools view -bt ref_list.txt -o aln.bam aln.sam.gz
11 samtools sort aln.bam aln.sorted
13 samtools index aln.sorted.bam
15 samtools idxstats aln.sorted.bam
17 samtools view aln.sorted.bam chr2:20,100,000-20,200,000
19 samtools merge out.bam in1.bam in2.bam in3.bam
21 samtools faidx ref.fasta
23 samtools pileup -vcf ref.fasta aln.sorted.bam
25 samtools mpileup -C50 -gf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam
27 samtools tview aln.sorted.bam ref.fasta
31 bcftools view in.bcf chr2:100-200 > out.vcf
33 bcftools view -Nvm0.99 in.bcf > out.vcf 2> out.afs
37 Samtools is a set of utilities that manipulate alignments in the BAM
38 format. It imports from and exports to the SAM (Sequence Alignment/Map)
39 format, does sorting, merging and indexing, and allows to retrieve reads
40 in any regions swiftly.
42 Samtools is designed to work on a stream. It regards an input file `-'
43 as the standard input (stdin) and an output file `-' as the standard
44 output (stdout). Several commands can thus be combined with Unix
45 pipes. Samtools always output warning and error messages to the standard
46 error output (stderr).
48 Samtools is also able to open a BAM (not SAM) file on a remote FTP or
49 HTTP server if the BAM file name starts with `ftp://' or `http://'.
50 Samtools checks the current working directory for the index file and
51 will download the index upon absence. Samtools does not retrieve the
52 entire alignment file unless it is asked to do so.
54 .SH SAMTOOLS COMMANDS AND OPTIONS
58 samtools view [-bchuHS] [-t in.refList] [-o output] [-f reqFlag] [-F
59 skipFlag] [-q minMapQ] [-l library] [-r readGroup] [-R rgFile] <in.bam>|<in.sam> [region1 [...]]
61 Extract/print all or sub alignments in SAM or BAM format. If no region
62 is specified, all the alignments will be printed; otherwise only
63 alignments overlapping the specified regions will be output. An
64 alignment may be given multiple times if it is overlapping several
65 regions. A region can be presented, for example, in the following
66 format: `chr2' (the whole chr2), `chr2:1000000' (region starting from
67 1,000,000bp) or `chr2:1,000,000-2,000,000' (region between 1,000,000 and
68 2,000,000bp including the end points). The coordinate is 1-based.
74 Output in the BAM format.
77 Only output alignments with all bits in INT present in the FLAG
78 field. INT can be in hex in the format of /^0x[0-9A-F]+/ [0]
81 Skip alignments with bits present in INT [0]
84 Include the header in the output.
87 Output the header only.
90 Only output reads in library STR [null]
96 Skip alignments with MAPQ smaller than INT [0]
99 Only output reads in read group STR [null]
102 Output reads in read groups listed in
107 Fraction of templates/pairs to subsample; the integer part is treated as the
108 seed for the random number generator [-1]
111 Input is in SAM. If @SQ header lines are absent, the
116 Instead of printing the alignments, only count them and print the
117 total number. All filter options, such as
122 , are taken into account.
125 This file is TAB-delimited. Each line must contain the reference name
126 and the length of the reference, one line for each distinct reference;
127 additional fields are ignored. This file also defines the order of the
128 reference sequences in sorting. If you run `samtools faidx <ref.fa>',
129 the resultant index file
136 Output uncompressed BAM. This option saves time spent on
137 compression/decomprssion and is thus preferred when the output is piped
138 to another samtools command.
153 Text alignment viewer (based on the ncurses library). In the viewer,
154 press `?' for help and press `g' to check the alignment start from a
155 region in the format like `chr10:10,000,000' or `=10,000,000' when
156 viewing the same reference sequence.
162 Output as (H)tml or (C)urses or (T)ext
165 Go directly to this position
168 Display only reads from this sample or read group
190 .IR Output file [stdout] ]
195 Generate BCF or pileup for one or multiple BAM files. Alignment records
196 are grouped by sample identifiers in @RG header lines. If sample
197 identifiers are absent, each input file is regarded as one sample.
199 In the pileup format (without
202 line represents a genomic position, consisting of chromosome name,
203 coordinate, reference base, read bases, read qualities and alignment
204 mapping qualities. Information on match, mismatch, indel, strand,
205 mapping quality and start and end of a read are all encoded at the read
206 base column. At this column, a dot stands for a match to the reference
207 base on the forward strand, a comma for a match on the reverse strand,
208 a '>' or '<' for a reference skip, `ACGTN' for a mismatch on the forward
209 strand and `acgtn' for a mismatch on the reverse strand. A pattern
210 `\\+[0-9]+[ACGTNacgtn]+' indicates there is an insertion between this
211 reference position and the next reference position. The length of the
212 insertion is given by the integer in the pattern, followed by the
213 inserted sequence. Similarly, a pattern `-[0-9]+[ACGTNacgtn]+'
214 represents a deletion from the reference. The deleted bases will be
215 presented as `*' in the following lines. Also at the read base column, a
216 symbol `^' marks the start of a read. The ASCII of the character
217 following `^' minus 33 gives the mapping quality. A symbol `$' marks the
218 end of a read segment.
224 Assume the quality is in the Illumina 1.3+ encoding.
226 Do not skip anomalous read pairs in variant calling.
229 Disable probabilistic realignment for the computation of base alignment
230 quality (BAQ). BAQ is the Phred-scaled probability of a read base being
231 misaligned. Applying this option greatly helps to reduce false SNPs
232 caused by misalignments.
235 List of input BAM files, one file per line [null]
238 Coefficient for downgrading mapping quality for reads containing
239 excessive mismatches. Given a read with a phred-scaled probability q of
240 being generated from the mapped position, the new mapping quality is
241 about sqrt((INT-q)/INT)*INT. A zero value disables this
242 functionality; if enabled, the recommended value for BWA is 50. [0]
245 At a position, read maximally
247 reads per input BAM. [250]
250 Extended BAQ computation. This option helps sensitivity especially for MNPs, but may hurt
251 specificity a little bit.
256 reference file in the FASTA format. The file can be optionally compressed by
261 BED or position list file containing a list of regions or sites where pileup or BCF should be generated [null]
264 Minimum mapping quality for an alignment to be used [0]
267 Minimum base quality for a base to be considered [13]
270 Only generate pileup in region
278 Output per-sample read depth
281 Compute genotype likelihoods and output them in the binary call format (BCF).
284 Output per-sample Phred-scaled strand bias P-value
289 except that the output is uncompressed BCF, which is preferred for piping.
291 .BR \-c ", " \-\-output =\fIFILE\fR
292 File to output bcf (vcf with \fB\-g\fR) to; default is stdout.
295 .B Options for Genotype Likelihood Computation (for -g or -u):
299 Phred-scaled gap extension sequencing error probability. Reducing
301 leads to longer indels. [20]
304 Coefficient for modeling homopolymer errors. Given an
307 run, the sequencing error of an indel of size
314 Do not perform INDEL calling
317 Skip INDEL calling if the average per-sample depth is above
322 Phred-scaled gap open sequencing error probability. Reducing
324 leads to more indel calls. [40]
327 Apply -m and -F thresholds per sample to increase sensitivity of calling.
328 By default both options are applied to reads pooled from all samples.
331 Comma dilimited list of platforms (determined by
333 from which indel candidates are obtained. It is recommended to collect
334 indel candidates from sequencing technologies that have low indel error
335 rate such as ILLUMINA. [all]
340 samtools reheader <in.header.sam> <in.bam>
342 Replace the header in
346 This command is much faster than replacing the header with a
347 BAM->SAM->BAM conversion.
351 samtools cat [-h header.sam] [-o out.bam] <in1.bam> <in2.bam> [ ... ]
353 Concatenate BAMs. The sequence dictionary of each input BAM must be identical,
354 although this command does not check this. This command uses a similar trick
357 which enables fast BAM concatenation.
361 samtools sort [-nof] [-m maxMem] <in.bam> <out.prefix>
363 Sort alignments by leftmost coordinates. File
365 will be created. This command may also create temporary files
366 .I <out.prefix>.%d.bam
367 when the whole alignment cannot be fitted into memory (controlled by
374 Output the final alignment to the standard output.
377 Sort by read names rather than by chromosomal coordinates
382 as the full output path and do not append
387 Approximately the maximum required memory. [500000000]
392 samtools merge [-nur1f] [-h inh.sam] [-R reg] <out.bam> <in1.bam> <in2.bam> [...]
394 Merge multiple sorted alignments.
395 The header reference lists of all the input BAM files, and the @SQ headers of
397 if any, must all refer to the same set of reference sequences.
398 The header reference list and (unless overridden by
404 and the headers of other files will be ignored.
410 Use zlib compression level 1 to comrpess the output
413 Force to overwrite the output file if present.
418 as `@' headers to be copied to
420 replacing any header lines that would otherwise be copied from
423 is actually in SAM format, though any alignment records it may contain
427 The input alignments are sorted by read names rather than by chromosomal
431 Merge files in the specified region indicated by
436 Attach an RG tag to each alignment. The tag value is inferred from file names.
439 Uncompressed BAM output
444 samtools index <aln.bam>
446 Index sorted alignment for fast random access. Index file
452 samtools idxstats <aln.bam>
454 Retrieve and print stats in the index file. The output is TAB delimited
455 with each line consisting of reference sequence name, sequence length, #
456 mapped reads and # unmapped reads.
460 samtools faidx <ref.fasta> [region1 [...]]
462 Index reference sequence in the FASTA format or extract subsequence from
463 indexed reference sequence. If no region is specified,
465 will index the file and create
467 on the disk. If regions are speficified, the subsequences will be
468 retrieved and printed to stdout in the FASTA format. The input file can
475 samtools fixmate <in.nameSrt.bam> <out.bam>
477 Fill in mate coordinates, ISIZE and mate related flags from a
478 name-sorted alignment.
482 samtools rmdup [-sS] <input.srt.bam> <out.bam>
484 Remove potential PCR duplicates: if multiple read pairs have identical
485 external coordinates, only retain the pair with highest mapping quality.
486 In the paired-end mode, this command
488 works with FR orientation and requires ISIZE is correctly set. It does
489 not work for unpaired reads (e.g. two ends mapped to different
490 chromosomes or orphan reads).
496 Remove duplicate for single-end reads. By default, the command works for
497 paired-end reads only.
500 Treat paired-end reads and single-end reads.
505 samtools calmd [-EeubSr] [-C capQcoef] <aln.bam> <ref.fasta>
507 Generate the MD tag. If the MD tag is already present, this command will
508 give a warning if the MD tag generated is different from the existing
509 tag. Output SAM by default.
515 When used jointly with
517 this option overwrites the original base quality.
520 Convert a the read base to = if it is identical to the aligned reference
521 base. Indel caller does not support the = bases at the moment.
524 Output uncompressed BAM
527 Output compressed BAM
530 The input is SAM with header lines
533 Coefficient to cap mapping quality of poorly mapped reads. See the
535 command for details. [0]
538 Compute the BQ tag (without -A) or cap base quality by BAQ (with -A).
541 Extended BAQ calculation. This option trades specificity for sensitivity, though the
547 samtools targetcut [-Q minBaseQ] [-i inPenalty] [-0 em0] [-1 em1] [-2 em2] [-f ref] <in.bam>
549 This command identifies target regions by examining the continuity of read depth, computes
550 haploid consensus sequences of targets and outputs a SAM with each sequence corresponding
551 to a target. When option
553 is in use, BAQ will be applied. This command is
555 designed for cutting fosmid clones from fosmid pool sequencing [Ref. Kitzman et al. (2010)].
560 samtools phase [-AF] [-k len] [-b prefix] [-q minLOD] [-Q minBaseQ] <in.bam>
562 Call and phase heterozygous SNPs.
567 Drop reads with ambiguous phase.
570 Prefix of BAM output. When this option is in use, phase-0 reads will be saved in file
574 Phase unknown reads will be randomly allocated to one of the two files. Chimeric reads
575 with switch errors will be saved in
576 .BR STR .chimeric.bam.
580 Do not attempt to fix chimeric reads.
583 Maximum length for local phasing. [13]
586 Minimum Phred-scaled LOD to call a heterozygote. [40]
589 Minimum base quality to be used in het calling. [13]
592 .SH BCFTOOLS COMMANDS AND OPTIONS
597 .RB [ \-AbFGNQSucgv ]
627 Convert between BCF and VCF, call variant candidates and estimate allele
632 .B Input/Output Options:
635 Retain all possible alternate alleles at variant sites. By default, the view
636 command discards unlikely alleles.
639 Output in the BCF format. The default is VCF.
642 Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]
645 Indicate PL is generated by r921 or before (ordering is different).
648 Suppress all individual genotype information.
651 List of sites at which information are outputted [all sites]
654 Skip sites where the REF field is not A/C/G/T
657 Output the QCALL likelihood format
660 List of samples to use. The first column in the input gives the sample names
661 and the second gives the ploidy, which can only be 1 or 2. When the 2nd column
662 is absent, the sample ploidy is assumed to be 2. In the output, the ordering of
663 samples will be identical to the one in
668 The input is VCF instead of BCF.
671 Uncompressed BCF output (force -b).
673 .B Consensus/Variant Calling Options:
676 Call variants using Bayesian inference. This option automatically invokes option
682 is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0]
685 Perform max-likelihood inference only, including estimating the site allele frequency,
686 testing Hardy-Weinberg equlibrium and testing associations with LRT.
689 Call per-sample genotypes at variant sites (force -c)
692 Ratio of INDEL-to-SNP mutation rate [0.15]
695 New model for improved multiallelic and rare-variant calling. Another
696 ALT allele is accepted if P(chi^2) of LRT exceeds the FLOAT threshold. The
697 parameter seems robust and the actual value usually does not affect the results
698 much; a good value to use is 0.99. This is the recommended calling method. [0]
701 A site is considered to be a variant if P(ref|D)<FLOAT [0.5]
704 Prior or initial allele frequency spectrum. If STR can be
708 or the file consisting of error output from a previous variant calling
712 Scaled muttion rate for variant calling [0.001]
715 Enable pair/trio calling. For trio calling, option
717 is usually needed to be applied to configure the trio members and their ordering.
718 In the file supplied to the option
720 the first sample must be the child, the second the father and the third the mother.
723 are `pair', `trioauto', `trioxd' and `trioxs', where `pair' calls differences between two input samples, and `trioxd' (`trioxs') specifies that the input
724 is from the X chromosome non-PAR regions and the child is a female (male). [null]
727 Output variant sites only (force -c)
729 .B Contrast Calling and Association Test Options:
732 Number of group-1 samples. This option is used for dividing the samples into
733 two groups for contrast SNP calling or association test.
734 When this option is in use, the following VCF INFO will be outputted:
735 PC2, PCHI2 and QCHI2. [0]
738 Number of permutations for association test (effective only with
743 Only perform permutations for P(chi^2)<FLOAT (effective only with
753 Index sorted BCF for random access.
760 .RI [ "in2.bcf " [ ... "]]]"
762 Concatenate BCF files. The input files are required to be sorted and
763 have identical samples appearing in the same order.
767 Sequence Alignment/Map (SAM) format is TAB-delimited. Apart from the header lines, which are started
768 with the `@' symbol, each alignment line consists of:
774 Col Field Description
776 1 QNAME Query template/pair NAME
778 3 RNAME Reference sequence NAME
779 4 POS 1-based leftmost POSition/coordinate of clipped sequence
780 5 MAPQ MAPping Quality (Phred-scaled)
781 6 CIAGR extended CIGAR string
782 7 MRNM Mate Reference sequence NaMe (`=' if same as RNAME)
783 8 MPOS 1-based Mate POSistion
784 9 TLEN inferred Template LENgth (insert size)
785 10 SEQ query SEQuence on the same strand as the reference
786 11 QUAL query QUALity (ASCII-33 gives the Phred base quality)
787 12+ OPT variable OPTional fields in the format TAG:VTYPE:VALUE
791 Each bit in the FLAG field is defined as:
799 0x0001 p the read is paired in sequencing
800 0x0002 P the read is mapped in a proper pair
801 0x0004 u the query sequence itself is unmapped
802 0x0008 U the mate is unmapped
803 0x0010 r strand of the query (1 for reverse)
804 0x0020 R strand of the mate
805 0x0040 1 the read is the first read in a pair
806 0x0080 2 the read is the second read in a pair
807 0x0100 s the alignment is not primary
808 0x0200 f the read fails platform/vendor quality checks
809 0x0400 d the read is either a PCR or an optical duplicate
812 where the second column gives the string representation of the FLAG field.
816 The Variant Call Format (VCF) is a TAB-delimited format with each data line consists of the following fields:
821 Col Field Description
823 1 CHROM CHROMosome name
824 2 POS the left-most POSition of the variant
825 3 ID unique variant IDentifier
826 4 REF the REFerence allele
827 5 ALT the ALTernate allele(s), separated by comma
828 6 QUAL variant/reference QUALity
829 7 FILTER FILTers applied
830 8 INFO INFOrmation related to the variant, separated by semi-colon
831 9 FORMAT FORMAT of the genotype fields, separated by colon (optional)
832 10+ SAMPLE SAMPLE genotypes and per-sample information (optional)
836 The following table gives the
838 tags used by samtools and bcftools.
844 Tag Format Description
846 AF1 double Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
847 DP int Raw read depth (without quality filtering)
848 DP4 int[4] # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
849 FQ int Consensus quality. Positive: sample genotypes different; negative: otherwise
850 MQ int Root-Mean-Square mapping quality of covering reads
851 PC2 int[2] Phred probability of AF in group1 samples being larger (,smaller) than in group2
852 PCHI2 double Posterior weighted chi^2 P-value between group1 and group2 samples
853 PV4 double[4] P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
854 QCHI2 int Phred-scaled PCHI2
855 RP int # permutations yielding a smaller PCHI2
856 CLR int Phred log ratio of genotype likelihoods with and without the trio/pair constraint
857 UGT string Most probable genotype configuration without the trio constraint
858 CGT string Most probable configuration with the trio constraint
859 VDB float Tests variant positions within reads. Intended for filtering RNA-seq artifacts around splice sites
860 RPB float Mann-Whitney rank-sum test for tail distance bias
861 HWE float Hardy-Weinberg equilibrium test, Wigginton et al., PMID: 15789306
866 Import SAM to BAM when
868 lines are present in the header:
870 samtools view -bS aln.sam > aln.bam
876 samtools faidx ref.fa
877 samtools view -bt ref.fa.fai aln.sam > aln.bam
881 is generated automatically by the
888 tag while merging sorted alignments:
890 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
891 samtools merge -rh rg.txt merged.bam ga.bam 454.bam
895 tag is determined by the file name the read is coming from. In this
908 Call SNPs and short INDELs for one diploid individual:
910 samtools mpileup -ugf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf
911 bcftools view var.raw.bcf | vcfutils.pl varFilter -D 100 > var.flt.vcf
915 option of varFilter controls the maximum read depth, which should be
916 adjusted to about twice the average read depth. One may consider to add
920 if mapping quality is overestimated for reads containing excessive
921 mismatches. Applying this option usually helps
923 but may not other mappers.
926 Generate the consensus sequence for one diploid individual:
928 samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq
931 Call somatic mutations from a pair of samples:
933 samtools mpileup -DSuf ref.fa aln.bam | bcftools view -bvcgT pair - > var.bcf
935 In the output INFO field,
937 gives the Phred-log ratio between the likelihood by treating the
938 two samples independently, and the likelihood by requiring the genotype to be identical.
941 is effectively a score measuring the confidence of somatic calls. The higher the better.
944 Call de novo and somatic mutations from a family trio:
946 samtools mpileup -DSuf ref.fa aln.bam | bcftools view -bvcgT pair -s samples.txt - > var.bcf
950 should consist of three lines specifying the member and order of samples (in the order of child-father-mother).
953 gives the Phred-log likelihood ratio with and without the trio constraint.
955 shows the most likely genotype configuration without the trio constraint, and
957 gives the most likely genotype configuration satisfying the trio constraint.
960 Phase one individual:
962 samtools calmd -AEur aln.bam ref.fa | samtools phase -b prefix - > phase.out
966 command is used to reduce false heterozygotes around INDELs.
969 Call SNPs and short indels for multiple diploid individuals:
971 samtools mpileup -P ILLUMINA -ugf ref.fa *.bam | bcftools view -bcvg - > var.raw.bcf
972 bcftools view var.raw.bcf | vcfutils.pl varFilter -D 2000 > var.flt.vcf
974 Individuals are identified from the
978 header lines. Individuals can be pooled in one alignment file; one
979 individual can also be separated into multiple files. The
981 option specifies that indel candidates should be collected only from
986 Collecting indel candidates from reads sequenced by an indel-prone
987 technology may affect the performance of indel calling.
989 Note that there is a new calling model which can be invoked by
991 bcftools view -m0.99 ...
993 which fixes some severe limitations of the default method.
995 For filtering, best results seem to be achieved by first applying the
997 filter and then applying some machine learning approach
999 vcf-annotate -f SnpGap=n
1002 Both can be found in the
1006 package (links below).
1009 Derive the allele frequency spectrum (AFS) on a list of sites from multiple individuals:
1011 samtools mpileup -Igf ref.fa *.bam > all.bcf
1012 bcftools view -bl sites.list all.bcf > sites.bcf
1013 bcftools view -cGP cond2 sites.bcf > /dev/null 2> sites.1.afs
1014 bcftools view -cGP sites.1.afs sites.bcf > /dev/null 2> sites.2.afs
1015 bcftools view -cGP sites.2.afs sites.bcf > /dev/null 2> sites.3.afs
1020 contains the list of sites with each line consisting of the reference
1021 sequence name and position. The following
1023 commands estimate AFS by EM.
1026 Dump BAQ applied alignment for other SNP callers:
1028 samtools calmd -bAr aln.bam > aln.baq.bam
1030 It adds and corrects the
1034 tags at the same time. The
1036 command also comes with the
1038 option, the same as the one in
1047 Unaligned words used in bam_import.c, bam_endian.h, bam.c and bam_aux.c.
1049 Samtools paired-end rmdup does not work for unpaired reads (e.g. orphan
1050 reads or ends mapped to different chromosomes). If this is a concern,
1051 please use Picard's MarkDuplicate which correctly handles these cases,
1052 although a little slower.
1056 Heng Li from the Sanger Institute wrote the C version of samtools. Bob
1057 Handsaker from the Broad Institute implemented the BGZF library and Jue
1058 Ruan from Beijing Genomics Institute wrote the RAZF library. John
1059 Marshall and Petr Danecek contribute to the source code and various
1060 people from the 1000 Genomes Project have contributed to the SAM format
1065 Samtools website: <http://samtools.sourceforge.net>
1067 Samtools latest source: <https://github.com/samtools/samtools>
1069 VCFtools website with stable link to VCF specification: <http://vcftools.sourceforge.net>
1071 HTSlib website: <https://github.com/samtools/htslib>
projects / samtools.git / blob