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