]> git.donarmstrong.com Git - samtools.git/blob - samtools.1
Merge branch 'master' of github.com:samtools/samtools
[samtools.git] / samtools.1
1 .TH samtools 1 "2 September 2011" "samtools-0.1.18" "Bioinformatics tools"
2 .SH NAME
3 .PP
4 samtools - Utilities for the Sequence Alignment/Map (SAM) format
5
6 bcftools - Utilities for the Binary Call Format (BCF) and VCF
7 .SH SYNOPSIS
8 .PP
9 samtools view -bt ref_list.txt -o aln.bam aln.sam.gz
10 .PP
11 samtools sort aln.bam aln.sorted
12 .PP
13 samtools index aln.sorted.bam
14 .PP
15 samtools idxstats aln.sorted.bam
16 .PP
17 samtools view aln.sorted.bam chr2:20,100,000-20,200,000
18 .PP
19 samtools merge out.bam in1.bam in2.bam in3.bam
20 .PP
21 samtools faidx ref.fasta
22 .PP
23 samtools pileup -vcf ref.fasta aln.sorted.bam
24 .PP
25 samtools mpileup -C50 -gf ref.fasta -r chr3:1,000-2,000 in1.bam in2.bam
26 .PP
27 samtools tview aln.sorted.bam ref.fasta
28 .PP
29 bcftools index in.bcf
30 .PP
31 bcftools view in.bcf chr2:100-200 > out.vcf
32 .PP
33 bcftools view -vc in.bcf > out.vcf 2> out.afs
34
35 .SH DESCRIPTION
36 .PP
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.
41
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).
47
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.
53
54 .SH SAMTOOLS COMMANDS AND OPTIONS
55
56 .TP 10
57 .B view
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 [...]]
60
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.
69
70 .B OPTIONS:
71 .RS
72 .TP 10
73 .B -b
74 Output in the BAM format.
75 .TP
76 .BI -f \ INT
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]
79 .TP
80 .BI -F \ INT
81 Skip alignments with bits present in INT [0]
82 .TP
83 .B -h
84 Include the header in the output.
85 .TP
86 .B -H
87 Output the header only.
88 .TP
89 .BI -l \ STR
90 Only output reads in library STR [null]
91 .TP
92 .BI -o \ FILE
93 Output file [stdout]
94 .TP
95 .BI -q \ INT
96 Skip alignments with MAPQ smaller than INT [0]
97 .TP
98 .BI -r \ STR
99 Only output reads in read group STR [null]
100 .TP
101 .BI -R \ FILE
102 Output reads in read groups listed in
103 .I FILE
104 [null]
105 .TP
106 .BI -s \ FLOAT
107 Fraction of templates/pairs to subsample; the integer part is treated as the
108 seed for the random number generator [-1]
109 .TP
110 .B -S
111 Input is in SAM. If @SQ header lines are absent, the
112 .B `-t'
113 option is required.
114 .TP
115 .B -c
116 Instead of printing the alignments, only count them and print the
117 total number. All filter options, such as
118 .B `-f',
119 .B `-F'
120 and
121 .B `-q'
122 , are taken into account.
123 .TP
124 .BI -t \ FILE
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
130 .I <ref.fa>.fai
131 can be used as this
132 .I <in.ref_list>
133 file.
134 .TP
135 .B -u
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.
139 .RE
140
141 .TP
142 .B tview
143 samtools tview <in.sorted.bam> [ref.fasta]
144
145 Text alignment viewer (based on the ncurses library). In the viewer,
146 press `?' for help and press `g' to check the alignment start from a
147 region in the format like `chr10:10,000,000' or `=10,000,000' when
148 viewing the same reference sequence.
149
150 .TP
151 .B mpileup
152 .B samtools mpileup
153 .RB [ \-EBug ]
154 .RB [ \-C
155 .IR capQcoef ]
156 .RB [ \-r
157 .IR reg ]
158 .RB [ \-f
159 .IR in.fa ]
160 .RB [ \-l
161 .IR list ]
162 .RB [ \-M
163 .IR capMapQ ]
164 .RB [ \-Q
165 .IR minBaseQ ]
166 .RB [ \-q
167 .IR minMapQ ]
168 .I in.bam
169 .RI [ in2.bam
170 .RI [ ... ]]
171
172 Generate BCF or pileup for one or multiple BAM files. Alignment records
173 are grouped by sample identifiers in @RG header lines. If sample
174 identifiers are absent, each input file is regarded as one sample.
175
176 In the pileup format (without
177 .BR -u or -g ),
178 each
179 line represents a genomic position, consisting of chromosome name,
180 coordinate, reference base, read bases, read qualities and alignment
181 mapping qualities. Information on match, mismatch, indel, strand,
182 mapping quality and start and end of a read are all encoded at the read
183 base column. At this column, a dot stands for a match to the reference
184 base on the forward strand, a comma for a match on the reverse strand,
185 a '>' or '<' for a reference skip, `ACGTN' for a mismatch on the forward
186 strand and `acgtn' for a mismatch on the reverse strand. A pattern
187 `\\+[0-9]+[ACGTNacgtn]+' indicates there is an insertion between this
188 reference position and the next reference position. The length of the
189 insertion is given by the integer in the pattern, followed by the
190 inserted sequence. Similarly, a pattern `-[0-9]+[ACGTNacgtn]+'
191 represents a deletion from the reference. The deleted bases will be
192 presented as `*' in the following lines. Also at the read base column, a
193 symbol `^' marks the start of a read. The ASCII of the character
194 following `^' minus 33 gives the mapping quality. A symbol `$' marks the
195 end of a read segment.
196
197 .B Input Options:
198 .RS
199 .TP 10
200 .B -6
201 Assume the quality is in the Illumina 1.3+ encoding.
202 .B -A
203 Do not skip anomalous read pairs in variant calling.
204 .TP
205 .B -B
206 Disable probabilistic realignment for the computation of base alignment
207 quality (BAQ). BAQ is the Phred-scaled probability of a read base being
208 misaligned. Applying this option greatly helps to reduce false SNPs
209 caused by misalignments.
210 .TP
211 .BI -b \ FILE
212 List of input BAM files, one file per line [null]
213 .TP
214 .BI -C \ INT
215 Coefficient for downgrading mapping quality for reads containing
216 excessive mismatches. Given a read with a phred-scaled probability q of
217 being generated from the mapped position, the new mapping quality is
218 about sqrt((INT-q)/INT)*INT. A zero value disables this
219 functionality; if enabled, the recommended value for BWA is 50. [0]
220 .TP
221 .BI -d \ INT
222 At a position, read maximally
223 .I INT
224 reads per input BAM. [250]
225 .TP
226 .B -E
227 Extended BAQ computation. This option helps sensitivity especially for MNPs, but may hurt
228 specificity a little bit.
229 .TP
230 .BI -f \ FILE
231 The
232 .BR faidx -indexed
233 reference file in the FASTA format. The file can be optionally compressed by
234 .BR razip .
235 [null]
236 .TP
237 .BI -l \ FILE
238 BED or position list file containing a list of regions or sites where pileup or BCF should be generated [null]
239 .TP
240 .BI -q \ INT
241 Minimum mapping quality for an alignment to be used [0]
242 .TP
243 .BI -Q \ INT
244 Minimum base quality for a base to be considered [13]
245 .TP
246 .BI -r \ STR
247 Only generate pileup in region
248 .I STR
249 [all sites]
250 .TP
251 .B Output Options:
252
253 .TP
254 .B -D
255 Output per-sample read depth
256 .TP
257 .B -g
258 Compute genotype likelihoods and output them in the binary call format (BCF).
259 .TP
260 .B -S
261 Output per-sample Phred-scaled strand bias P-value
262 .TP
263 .B -u
264 Similar to
265 .B -g
266 except that the output is uncompressed BCF, which is preferred for piping.
267
268 .TP
269 .B Options for Genotype Likelihood Computation (for -g or -u):
270
271 .TP
272 .BI -e \ INT
273 Phred-scaled gap extension sequencing error probability. Reducing
274 .I INT
275 leads to longer indels. [20]
276 .TP
277 .BI -h \ INT
278 Coefficient for modeling homopolymer errors. Given an
279 .IR l -long
280 homopolymer
281 run, the sequencing error of an indel of size
282 .I s
283 is modeled as
284 .IR INT * s / l .
285 [100]
286 .TP
287 .B -I
288 Do not perform INDEL calling
289 .TP
290 .BI -L \ INT
291 Skip INDEL calling if the average per-sample depth is above
292 .IR INT .
293 [250]
294 .TP
295 .BI -o \ INT
296 Phred-scaled gap open sequencing error probability. Reducing
297 .I INT
298 leads to more indel calls. [40]
299 .TP
300 .BI -P \ STR
301 Comma dilimited list of platforms (determined by
302 .BR @RG-PL )
303 from which indel candidates are obtained. It is recommended to collect
304 indel candidates from sequencing technologies that have low indel error
305 rate such as ILLUMINA. [all]
306 .RE
307
308 .TP
309 .B reheader
310 samtools reheader <in.header.sam> <in.bam>
311
312 Replace the header in
313 .I in.bam
314 with the header in
315 .I in.header.sam.
316 This command is much faster than replacing the header with a
317 BAM->SAM->BAM conversion.
318
319 .TP
320 .B cat
321 samtools cat [-h header.sam] [-o out.bam] <in1.bam> <in2.bam> [ ... ]
322
323 Concatenate BAMs. The sequence dictionary of each input BAM must be identical,
324 although this command does not check this. This command uses a similar trick
325 to
326 .B reheader
327 which enables fast BAM concatenation.
328
329 .TP
330 .B sort
331 samtools sort [-nof] [-m maxMem] <in.bam> <out.prefix>
332
333 Sort alignments by leftmost coordinates. File
334 .I <out.prefix>.bam
335 will be created. This command may also create temporary files
336 .I <out.prefix>.%d.bam
337 when the whole alignment cannot be fitted into memory (controlled by
338 option -m).
339
340 .B OPTIONS:
341 .RS
342 .TP 8
343 .B -o
344 Output the final alignment to the standard output.
345 .TP
346 .B -n
347 Sort by read names rather than by chromosomal coordinates
348 .TP
349 .B -f
350 Use
351 .I <out.prefix>
352 as the full output path and do not append
353 .I .bam
354 suffix.
355 .TP
356 .BI -m \ INT
357 Approximately the maximum required memory. [500000000]
358 .RE
359
360 .TP
361 .B merge
362 samtools merge [-nur1f] [-h inh.sam] [-R reg] <out.bam> <in1.bam> <in2.bam> [...]
363
364 Merge multiple sorted alignments.
365 The header reference lists of all the input BAM files, and the @SQ headers of
366 .IR inh.sam ,
367 if any, must all refer to the same set of reference sequences.
368 The header reference list and (unless overridden by
369 .BR -h )
370 `@' headers of
371 .I in1.bam
372 will be copied to
373 .IR out.bam ,
374 and the headers of other files will be ignored.
375
376 .B OPTIONS:
377 .RS
378 .TP 8
379 .B -1
380 Use zlib compression level 1 to comrpess the output
381 .TP
382 .B -f
383 Force to overwrite the output file if present.
384 .TP 8
385 .BI -h \ FILE
386 Use the lines of
387 .I FILE
388 as `@' headers to be copied to
389 .IR out.bam ,
390 replacing any header lines that would otherwise be copied from
391 .IR in1.bam .
392 .RI ( FILE
393 is actually in SAM format, though any alignment records it may contain
394 are ignored.)
395 .TP
396 .B -n
397 The input alignments are sorted by read names rather than by chromosomal
398 coordinates
399 .TP
400 .BI -R \ STR
401 Merge files in the specified region indicated by
402 .I STR
403 [null]
404 .TP
405 .B -r
406 Attach an RG tag to each alignment. The tag value is inferred from file names.
407 .TP
408 .B -u
409 Uncompressed BAM output
410 .RE
411
412 .TP
413 .B index
414 samtools index <aln.bam>
415
416 Index sorted alignment for fast random access. Index file
417 .I <aln.bam>.bai
418 will be created.
419
420 .TP
421 .B idxstats
422 samtools idxstats <aln.bam>
423
424 Retrieve and print stats in the index file. The output is TAB delimited
425 with each line consisting of reference sequence name, sequence length, #
426 mapped reads and # unmapped reads.
427
428 .TP
429 .B faidx
430 samtools faidx <ref.fasta> [region1 [...]]
431
432 Index reference sequence in the FASTA format or extract subsequence from
433 indexed reference sequence. If no region is specified,
434 .B faidx
435 will index the file and create
436 .I <ref.fasta>.fai
437 on the disk. If regions are speficified, the subsequences will be
438 retrieved and printed to stdout in the FASTA format. The input file can
439 be compressed in the
440 .B RAZF
441 format.
442
443 .TP
444 .B fixmate
445 samtools fixmate <in.nameSrt.bam> <out.bam>
446
447 Fill in mate coordinates, ISIZE and mate related flags from a
448 name-sorted alignment.
449
450 .TP
451 .B rmdup
452 samtools rmdup [-sS] <input.srt.bam> <out.bam>
453
454 Remove potential PCR duplicates: if multiple read pairs have identical
455 external coordinates, only retain the pair with highest mapping quality.
456 In the paired-end mode, this command
457 .B ONLY
458 works with FR orientation and requires ISIZE is correctly set. It does
459 not work for unpaired reads (e.g. two ends mapped to different
460 chromosomes or orphan reads).
461
462 .B OPTIONS:
463 .RS
464 .TP 8
465 .B -s
466 Remove duplicate for single-end reads. By default, the command works for
467 paired-end reads only.
468 .TP 8
469 .B -S
470 Treat paired-end reads and single-end reads.
471 .RE
472
473 .TP
474 .B calmd
475 samtools calmd [-EeubSr] [-C capQcoef] <aln.bam> <ref.fasta>
476
477 Generate the MD tag. If the MD tag is already present, this command will
478 give a warning if the MD tag generated is different from the existing
479 tag. Output SAM by default.
480
481 .B OPTIONS:
482 .RS
483 .TP 8
484 .B -A
485 When used jointly with
486 .B -r
487 this option overwrites the original base quality.
488 .TP 8
489 .B -e
490 Convert a the read base to = if it is identical to the aligned reference
491 base. Indel caller does not support the = bases at the moment.
492 .TP
493 .B -u
494 Output uncompressed BAM
495 .TP
496 .B -b
497 Output compressed BAM
498 .TP
499 .B -S
500 The input is SAM with header lines
501 .TP
502 .BI -C \ INT
503 Coefficient to cap mapping quality of poorly mapped reads. See the
504 .B pileup
505 command for details. [0]
506 .TP
507 .B -r
508 Compute the BQ tag (without -A) or cap base quality by BAQ (with -A).
509 .TP
510 .B -E
511 Extended BAQ calculation. This option trades specificity for sensitivity, though the
512 effect is minor.
513 .RE
514
515 .TP
516 .B targetcut
517 samtools targetcut [-Q minBaseQ] [-i inPenalty] [-0 em0] [-1 em1] [-2 em2] [-f ref] <in.bam>
518
519 This command identifies target regions by examining the continuity of read depth, computes
520 haploid consensus sequences of targets and outputs a SAM with each sequence corresponding
521 to a target. When option
522 .B -f
523 is in use, BAQ will be applied. This command is
524 .B only
525 designed for cutting fosmid clones from fosmid pool sequencing [Ref. Kitzman et al. (2010)].
526 .RE
527
528 .TP
529 .B phase
530 samtools phase [-AF] [-k len] [-b prefix] [-q minLOD] [-Q minBaseQ] <in.bam>
531
532 Call and phase heterozygous SNPs.
533 .B OPTIONS:
534 .RS
535 .TP 8
536 .B -A
537 Drop reads with ambiguous phase.
538 .TP 8
539 .BI -b \ STR
540 Prefix of BAM output. When this option is in use, phase-0 reads will be saved in file
541 .BR STR .0.bam
542 and phase-1 reads in
543 .BR STR .1.bam.
544 Phase unknown reads will be randomly allocated to one of the two files. Chimeric reads
545 with switch errors will be saved in
546 .BR STR .chimeric.bam.
547 [null]
548 .TP
549 .B -F
550 Do not attempt to fix chimeric reads.
551 .TP
552 .BI -k \ INT
553 Maximum length for local phasing. [13]
554 .TP
555 .BI -q \ INT
556 Minimum Phred-scaled LOD to call a heterozygote. [40]
557 .TP
558 .BI -Q \ INT
559 Minimum base quality to be used in het calling. [13]
560 .RE
561
562 .SH BCFTOOLS COMMANDS AND OPTIONS
563
564 .TP 10
565 .B view
566 .B bcftools view
567 .RB [ \-AbFGNQSucgv ]
568 .RB [ \-D
569 .IR seqDict ]
570 .RB [ \-l
571 .IR listLoci ]
572 .RB [ \-s
573 .IR listSample ]
574 .RB [ \-i
575 .IR gapSNPratio ]
576 .RB [ \-t
577 .IR mutRate ]
578 .RB [ \-p
579 .IR varThres ]
580 .RB [ \-P
581 .IR prior ]
582 .RB [ \-1
583 .IR nGroup1 ]
584 .RB [ \-d
585 .IR minFrac ]
586 .RB [ \-U
587 .IR nPerm ]
588 .RB [ \-X
589 .IR permThres ]
590 .RB [ \-T
591 .IR trioType ]
592 .I in.bcf
593 .RI [ region ]
594
595 Convert between BCF and VCF, call variant candidates and estimate allele
596 frequencies.
597
598 .RS
599 .TP
600 .B Input/Output Options:
601 .TP 10
602 .B -A
603 Retain all possible alternate alleles at variant sites. By default, the view
604 command discards unlikely alleles.
605 .TP 10
606 .B -b
607 Output in the BCF format. The default is VCF.
608 .TP
609 .BI -D \ FILE
610 Sequence dictionary (list of chromosome names) for VCF->BCF conversion [null]
611 .TP
612 .B -F
613 Indicate PL is generated by r921 or before (ordering is different).
614 .TP
615 .B -G
616 Suppress all individual genotype information.
617 .TP
618 .BI -l \ FILE
619 List of sites at which information are outputted [all sites]
620 .TP
621 .B -N
622 Skip sites where the REF field is not A/C/G/T
623 .TP
624 .B -Q
625 Output the QCALL likelihood format
626 .TP
627 .BI -s \ FILE
628 List of samples to use. The first column in the input gives the sample names
629 and the second gives the ploidy, which can only be 1 or 2. When the 2nd column
630 is absent, the sample ploidy is assumed to be 2. In the output, the ordering of
631 samples will be identical to the one in
632 .IR FILE .
633 [null]
634 .TP
635 .B -S
636 The input is VCF instead of BCF.
637 .TP
638 .B -u
639 Uncompressed BCF output (force -b).
640 .TP
641 .B Consensus/Variant Calling Options:
642 .TP 10
643 .B -c
644 Call variants using Bayesian inference. This option automatically invokes option
645 .BR -e .
646 .TP
647 .BI -d \ FLOAT
648 When
649 .B -v
650 is in use, skip loci where the fraction of samples covered by reads is below FLOAT. [0]
651 .TP
652 .B -e
653 Perform max-likelihood inference only, including estimating the site allele frequency,
654 testing Hardy-Weinberg equlibrium and testing associations with LRT.
655 .TP
656 .B -g
657 Call per-sample genotypes at variant sites (force -c)
658 .TP
659 .BI -i \ FLOAT
660 Ratio of INDEL-to-SNP mutation rate [0.15]
661 .TP
662 .BI -p \ FLOAT
663 A site is considered to be a variant if P(ref|D)<FLOAT [0.5]
664 .TP
665 .BI -P \ STR
666 Prior or initial allele frequency spectrum. If STR can be
667 .IR full ,
668 .IR cond2 ,
669 .I flat
670 or the file consisting of error output from a previous variant calling
671 run.
672 .TP
673 .BI -t \ FLOAT
674 Scaled muttion rate for variant calling [0.001]
675 .TP
676 .BI -T \ STR
677 Enable pair/trio calling. For trio calling, option
678 .B -s
679 is usually needed to be applied to configure the trio members and their ordering.
680 In the file supplied to the option
681 .BR -s ,
682 the first sample must be the child, the second the father and the third the mother.
683 The valid values of
684 .I STR
685 are `pair', `trioauto', `trioxd' and `trioxs', where `pair' calls differences between two input samples, and `trioxd' (`trioxs') specifies that the input
686 is from the X chromosome non-PAR regions and the child is a female (male). [null]
687 .TP
688 .B -v
689 Output variant sites only (force -c)
690 .TP
691 .B Contrast Calling and Association Test Options:
692 .TP
693 .BI -1 \ INT
694 Number of group-1 samples. This option is used for dividing the samples into
695 two groups for contrast SNP calling or association test.
696 When this option is in use, the following VCF INFO will be outputted:
697 PC2, PCHI2 and QCHI2. [0]
698 .TP
699 .BI -U \ INT
700 Number of permutations for association test (effective only with
701 .BR -1 )
702 [0]
703 .TP
704 .BI -X \ FLOAT
705 Only perform permutations for P(chi^2)<FLOAT (effective only with
706 .BR -U )
707 [0.01]
708 .RE
709
710 .TP
711 .B index
712 .B bcftools index
713 .I in.bcf
714
715 Index sorted BCF for random access.
716 .RE
717
718 .TP
719 .B cat
720 .B bcftools cat
721 .I in1.bcf
722 .RI [ "in2.bcf " [ ... "]]]"
723
724 Concatenate BCF files. The input files are required to be sorted and
725 have identical samples appearing in the same order.
726 .RE
727 .SH SAM FORMAT
728
729 Sequence Alignment/Map (SAM) format is TAB-delimited. Apart from the header lines, which are started
730 with the `@' symbol, each alignment line consists of:
731
732 .TS
733 center box;
734 cb | cb | cb
735 n | l | l .
736 Col     Field   Description
737 _
738 1       QNAME   Query template/pair NAME
739 2       FLAG    bitwise FLAG
740 3       RNAME   Reference sequence NAME
741 4       POS     1-based leftmost POSition/coordinate of clipped sequence
742 5       MAPQ    MAPping Quality (Phred-scaled)
743 6       CIAGR   extended CIGAR string
744 7       MRNM    Mate Reference sequence NaMe (`=' if same as RNAME)
745 8       MPOS    1-based Mate POSistion
746 9       TLEN    inferred Template LENgth (insert size)
747 10      SEQ     query SEQuence on the same strand as the reference
748 11      QUAL    query QUALity (ASCII-33 gives the Phred base quality)
749 12+     OPT     variable OPTional fields in the format TAG:VTYPE:VALUE
750 .TE
751
752 .PP
753 Each bit in the FLAG field is defined as:
754
755 .TS
756 center box;
757 cb | cb | cb
758 l | c | l .
759 Flag    Chr     Description
760 _
761 0x0001  p       the read is paired in sequencing
762 0x0002  P       the read is mapped in a proper pair
763 0x0004  u       the query sequence itself is unmapped
764 0x0008  U       the mate is unmapped
765 0x0010  r       strand of the query (1 for reverse)
766 0x0020  R       strand of the mate
767 0x0040  1       the read is the first read in a pair
768 0x0080  2       the read is the second read in a pair
769 0x0100  s       the alignment is not primary
770 0x0200  f       the read fails platform/vendor quality checks
771 0x0400  d       the read is either a PCR or an optical duplicate
772 .TE
773
774 where the second column gives the string representation of the FLAG field.
775
776 .SH VCF FORMAT
777
778 The Variant Call Format (VCF) is a TAB-delimited format with each data line consists of the following fields:
779 .TS
780 center box;
781 cb | cb | cb
782 n | l | l .
783 Col     Field   Description
784 _
785 1       CHROM   CHROMosome name
786 2       POS     the left-most POSition of the variant
787 3       ID      unique variant IDentifier
788 4       REF     the REFerence allele
789 5       ALT     the ALTernate allele(s), separated by comma
790 6       QUAL    variant/reference QUALity
791 7       FILTER  FILTers applied
792 8       INFO    INFOrmation related to the variant, separated by semi-colon
793 9       FORMAT  FORMAT of the genotype fields, separated by colon (optional)
794 10+     SAMPLE  SAMPLE genotypes and per-sample information (optional)
795 .TE
796
797 .PP
798 The following table gives the
799 .B INFO
800 tags used by samtools and bcftools.
801
802 .TS
803 center box;
804 cb | cb | cb
805 l | l | l .
806 Tag     Format  Description
807 _
808 AF1     double  Max-likelihood estimate of the site allele frequency (AF) of the first ALT allele
809 DP      int     Raw read depth (without quality filtering)
810 DP4     int[4]  # high-quality reference forward bases, ref reverse, alternate for and alt rev bases
811 FQ      int     Consensus quality. Positive: sample genotypes different; negative: otherwise
812 MQ      int     Root-Mean-Square mapping quality of covering reads
813 PC2     int[2]  Phred probability of AF in group1 samples being larger (,smaller) than in group2
814 PCHI2   double  Posterior weighted chi^2 P-value between group1 and group2 samples
815 PV4     double[4]       P-value for strand bias, baseQ bias, mapQ bias and tail distance bias
816 QCHI2   int     Phred-scaled PCHI2
817 RP      int     # permutations yielding a smaller PCHI2
818 CLR     int     Phred log ratio of genotype likelihoods with and without the trio/pair constraint
819 UGT     string  Most probable genotype configuration without the trio constraint
820 CGT     string  Most probable configuration with the trio constraint
821 .TE
822
823 .SH EXAMPLES
824 .IP o 2
825 Import SAM to BAM when
826 .B @SQ
827 lines are present in the header:
828
829   samtools view -bS aln.sam > aln.bam
830
831 If
832 .B @SQ
833 lines are absent:
834
835   samtools faidx ref.fa
836   samtools view -bt ref.fa.fai aln.sam > aln.bam
837
838 where
839 .I ref.fa.fai
840 is generated automatically by the
841 .B faidx
842 command.
843
844 .IP o 2
845 Attach the
846 .B RG
847 tag while merging sorted alignments:
848
849   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
850   samtools merge -rh rg.txt merged.bam ga.bam 454.bam
851
852 The value in a
853 .B RG
854 tag is determined by the file name the read is coming from. In this
855 example, in the
856 .IR merged.bam ,
857 reads from
858 .I ga.bam
859 will be attached 
860 .IR RG:Z:ga ,
861 while reads from
862 .I 454.bam
863 will be attached
864 .IR RG:Z:454 .
865
866 .IP o 2
867 Call SNPs and short INDELs for one diploid individual:
868
869   samtools mpileup -ugf ref.fa aln.bam | bcftools view -bvcg - > var.raw.bcf
870   bcftools view var.raw.bcf | vcfutils.pl varFilter -D 100 > var.flt.vcf
871
872 The
873 .B -D
874 option of varFilter controls the maximum read depth, which should be
875 adjusted to about twice the average read depth.  One may consider to add
876 .B -C50
877 to
878 .B mpileup
879 if mapping quality is overestimated for reads containing excessive
880 mismatches. Applying this option usually helps
881 .B BWA-short
882 but may not other mappers.
883
884 .IP o 2
885 Generate the consensus sequence for one diploid individual:
886
887   samtools mpileup -uf ref.fa aln.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq
888
889 .IP o 2
890 Call somatic mutations from a pair of samples:
891
892   samtools mpileup -DSuf ref.fa aln.bam | bcftools view -bvcgT pair - > var.bcf
893
894 In the output INFO field,
895 .I CLR
896 gives the Phred-log ratio between the likelihood by treating the
897 two samples independently, and the likelihood by requiring the genotype to be identical.
898 This
899 .I CLR
900 is effectively a score measuring the confidence of somatic calls. The higher the better.
901
902 .IP o 2
903 Call de novo and somatic mutations from a family trio:
904
905   samtools mpileup -DSuf ref.fa aln.bam | bcftools view -bvcgT pair -s samples.txt - > var.bcf
906
907 File
908 .I samples.txt
909 should consist of three lines specifying the member and order of samples (in the order of child-father-mother).
910 Similarly,
911 .I CLR
912 gives the Phred-log likelihood ratio with and without the trio constraint.
913 .I UGT
914 shows the most likely genotype configuration without the trio constraint, and
915 .I CGT
916 gives the most likely genotype configuration satisfying the trio constraint.
917
918 .IP o 2
919 Phase one individual:
920
921   samtools calmd -AEur aln.bam ref.fa | samtools phase -b prefix - > phase.out
922
923 The
924 .B calmd
925 command is used to reduce false heterozygotes around INDELs.
926
927 .IP o 2
928 Call SNPs and short indels for multiple diploid individuals:
929
930   samtools mpileup -P ILLUMINA -ugf ref.fa *.bam | bcftools view -bcvg - > var.raw.bcf
931   bcftools view var.raw.bcf | vcfutils.pl varFilter -D 2000 > var.flt.vcf
932
933 Individuals are identified from the
934 .B SM
935 tags in the
936 .B @RG
937 header lines. Individuals can be pooled in one alignment file; one
938 individual can also be separated into multiple files. The
939 .B -P
940 option specifies that indel candidates should be collected only from
941 read groups with the
942 .B @RG-PL
943 tag set to
944 .IR ILLUMINA .
945 Collecting indel candidates from reads sequenced by an indel-prone
946 technology may affect the performance of indel calling.
947
948 .IP o 2
949 Derive the allele frequency spectrum (AFS) on a list of sites from multiple individuals:
950
951   samtools mpileup -Igf ref.fa *.bam > all.bcf
952   bcftools view -bl sites.list all.bcf > sites.bcf
953   bcftools view -cGP cond2 sites.bcf > /dev/null 2> sites.1.afs
954   bcftools view -cGP sites.1.afs sites.bcf > /dev/null 2> sites.2.afs
955   bcftools view -cGP sites.2.afs sites.bcf > /dev/null 2> sites.3.afs
956   ......
957
958 where
959 .I sites.list
960 contains the list of sites with each line consisting of the reference
961 sequence name and position. The following
962 .B bcftools
963 commands estimate AFS by EM.
964
965 .IP o 2
966 Dump BAQ applied alignment for other SNP callers:
967
968   samtools calmd -bAr aln.bam > aln.baq.bam
969
970 It adds and corrects the
971 .B NM
972 and
973 .B MD
974 tags at the same time. The
975 .B calmd
976 command also comes with the
977 .B -C
978 option, the same as the one in
979 .B pileup
980 and
981 .BR mpileup .
982 Apply if it helps.
983
984 .SH LIMITATIONS
985 .PP
986 .IP o 2
987 Unaligned words used in bam_import.c, bam_endian.h, bam.c and bam_aux.c.
988 .IP o 2
989 Samtools paired-end rmdup does not work for unpaired reads (e.g. orphan
990 reads or ends mapped to different chromosomes). If this is a concern,
991 please use Picard's MarkDuplicate which correctly handles these cases,
992 although a little slower.
993
994 .SH AUTHOR
995 .PP
996 Heng Li from the Sanger Institute wrote the C version of samtools. Bob
997 Handsaker from the Broad Institute implemented the BGZF library and Jue
998 Ruan from Beijing Genomics Institute wrote the RAZF library. John
999 Marshall and Petr Danecek contribute to the source code and various
1000 people from the 1000 Genomes Project have contributed to the SAM format
1001 specification.
1002
1003 .SH SEE ALSO
1004 .PP
1005 Samtools website: <http://samtools.sourceforge.net>