samtools-0.1.19/bcftools/bcftools view file.bcf1 | bcftools view

none

some

all

snps

indels

both

id

    bcftools view -Ou -s sample1,sample2 file.vcf | bcftools query -f %INFO/AC\t%INFO/AN\n

    sample1    1
    sample2    2
    sample3    2

    sample1    M
    sample2    F
    sample3    F

    ignored_column  daughterA fatherA  motherA  2
    ignored_column  sonB      fatherB  motherB  1

With the call -C alleles command, third column of the targets file must be comma-separated list of alleles, starting with the reference allele. Note that the file must be compressed and indexed. Such a file can be easily created from a VCF using:

    bcftools query -f'%CHROM\t%POS\t%REF,%ALT\n' file.vcf | bgzip -c > als.tsv.gz && tabix -s1 -b2 -e2 als.tsv.gz

    # Sample annotation file with columns CHROM, POS, STRING_TAG, NUMERIC_TAG
    1  752566  SomeString      5
    1  798959  SomeOtherString 6

    ##INFO=<ID=NUMERIC_TAG,Number=1,Type=Integer,Description="Example header line">
    ##INFO=<ID=STRING_TAG,Number=1,Type=String,Description="Yet another header line">

    bcftools annotate --set-id +'%CHROM\_%POS\_%REF\_%FIRST_ALT' file.vcf

    # Remove three fields
    bcftools annotate -x ID,INFO/DP,FORMAT/DP file.vcf.gz
    # Remove all INFO fields and all FORMAT fields except for GT and PL
    bcftools annotate -x INFO,^FORMAT/GT,FORMAT/PL file.vcf
    # Add ID, QUAL and INFO/TAG, not replacing TAG if already present
    bcftools annotate -a src.bcf -c ID,QUAL,+TAG dst.bcf
    # Carry over all INFO and FORMAT annotations except FORMAT/GT
    bcftools annotate -a src.bcf -c INFO,^FORMAT/GT dst.bcf
    # Annotate from a tab-delimited file with six columns (the fifth is ignored),
    # first indexing with tabix. The coordinates are 1-based.
    tabix -s1 -b2 -e2 annots.tab.gz
    bcftools annotate -a annots.tab.gz -h annots.hdr -c CHROM,POS,REF,ALT,-,TAG file.vcf
    # Annotate from a tab-delimited file with regions (1-based coordinates, inclusive)
    tabix -s1 -b2 -e3 annots.tab.gz
    bcftools annotate -a annots.tab.gz -h annots.hdr -c CHROM,FROM,TO,TAG input.vcf
    # Annotate from a bed file (0-based coordinates, half-closed, half-open intervals)
    bcftools annotate -a annots.bed.gz -h annots.hdr -c CHROM,FROM,TO,TAG input.vcf
    # Transfer the INFO/END tag, matching by POS,REF,ALT and ID. This example assumes
    # that INFO/END is already present in the VCF header.
    bcftools annotate -a annots.tab.gz  -c CHROM,POS,~ID,REF,ALT,INFO/END input.vcf
    # For more examples see http://samtools.github.io/bcftools/howtos/annotate.html

    X 1 60000 M 1
    X 2699521 154931043 M 1
    Y 1 59373566 M 1
    Y 1 59373566 F 0
    MT 1 16569 M 1
    MT 1 16569 F 1
    *  * *     M 2
    *  * *     F 2

    # Extract AN,AC values from an existing VCF, such 1000Genomes
    bcftools query -f'%CHROM\t%POS\t%REF\t%ALT\t%AN\t%AC\n' 1000Genomes.bcf | bgzip -c > AFs.tab.gz
    # If the tags AN,AC are not already present, use the +fill-tags plugin
    bcftools +fill-tags 1000Genomes.bcf | bcftools query -f'%CHROM\t%POS\t%REF\t%ALT\t%AN\t%AC\n' | bgzip -c > AFs.tab.gz
    tabix -s1 -b2 -e2 AFs.tab.gz
    # Create a VCF header description, here we name the tags REF_AN,REF_AC
    cat AFs.hdr
    ##INFO=<ID=REF_AN,Number=1,Type=Integer,Description="Total number of alleles in reference genotypes">
    ##INFO=<ID=REF_AC,Number=A,Type=Integer,Description="Allele count in reference genotypes for each ALT allele">
    # Now before calling, stream the raw mpileup output through `bcftools annotate` to add the frequencies
    bcftools mpileup [...] -Ou | bcftools annotate -a AFs.tab.gz -h AFs.hdr -c CHROM,POS,REF,ALT,REF_AN,REF_AC -Ou | bcftools call -mv -F REF_AN,REF_AC [...]

alleles

trio

1

2

R

A

I

LR, LA

SR, SA

1pIu, 2pIu

This option requires *-s*, unless exactly one sample is present in the VCF

    # Apply variants present in sample "NA001", output IUPAC codes for hets
    bcftools consensus -i -s NA001 -f in.fa in.vcf.gz > out.fa
    # Create consensus for one region. The fasta header lines are then expected
    # in the form ">chr:from-to".
    samtools faidx ref.fa 8:11870-11890 | bcftools consensus in.vcf.gz -o out.fa

  .gen (with --3N6 --vcf-ids)
  ---------------------------
  chr1 1:111485207_G_A rsID1 111485207 G A 0 1 0 0 1 0
  chr1 1:111494194_C_T rsID2 111494194 C T 0 1 0 0 0 1
  .gen (with --vcf-ids)
  ---------------------------
  1:111485207_G_A rsID1 111485207 G A 0 1 0 0 1 0
  1:111494194_C_T rsID2 111494194 C T 0 1 0 0 0 1
  .gen (the default)
  ------------------------------
  1:111485207_G_A 1:111485207_G_A 111485207 G A 0 1 0 0 1 0
  1:111494194_C_T 1:111494194_C_T 111494194 C T 0 1 0 0 0 1
  .sample
  -------
  ID_1 ID_2 missing
  0 0 0
  sample1 sample1 0
  sample2 sample2 0

    MaleSample    M
    FemaleSample  F

  .hap (with --vcf-ids)
  ---------------------
  1:111485207_G_A rsID1 111485207 G A 0 1 0 0
  1:111495231_A_<DEL>_111495784 rsID3 111495231 A <DEL> 0 0 1 0
  .hap (the default)
  ------------------
  1 1:111485207_G_A 111485207 G A 0 1 0 0
  1 1:111495231_A_<DEL>_111495784 111495231 A <DEL> 0 0 1 0

    MaleSample    M
    FemaleSample  F

  .hap
  -----
  0 1 0 0 1 0
  0 1 0 0 0 1
  .legend
  -------
  id position a0 a1
  1:111485207_G_A 111485207 G A
  1:111494194_C_T 111494194 C T
  .sample
  -------
  sample population group sex
  sample1 sample1 sample1 2
  sample2 sample2 sample2 2

    MaleSample    M
    FemaleSample  F

# Convert 23andme results into VCF
bcftools convert -c ID,CHROM,POS,AA -s SampleName -f 23andme-ref.fa --tsv2vcf 23andme.txt -Oz -o out.vcf.gz

    # The program looks for "CDS", "exon", "three_prime_UTR" and "five_prime_UTR" lines,
    # looks up their parent transcript (determined from the "Parent=transcript:" attribute),
    # the gene (determined from the transcript's "Parent=gene:" attribute), and the biotype
    # (the most interesting is "protein_coding").
    #
    # Attributes required for
    #   gene lines:
    #   - ID=gene:<gene_id>
    #   - biotype=<biotype>
    #   - Name=<gene_name>      [optional]
    #
    #   transcript lines:
    #   - ID=transcript:<transcript_id>
    #   - Parent=gene:<gene_id>
    #   - biotype=<biotype>
    #
    #   other lines (CDS, exon, five_prime_UTR, three_prime_UTR):
    #   - Parent=transcript:<transcript_id>
    #
    # Supported biotypes:
    #   - see the function gff_parse_biotype() in bcftools/csq.c
    1   ignored_field  gene            21  2148  .   -   .   ID=gene:GeneId;biotype=protein_coding;Name=GeneName
    1   ignored_field  transcript      21  2148  .   -   .   ID=transcript:TranscriptId;Parent=gene:GeneId;biotype=protein_coding
    1   ignored_field  three_prime_UTR 21  2054  .   -   .   Parent=transcript:TranscriptId
    1   ignored_field  exon            21  2148  .   -   .   Parent=transcript:TranscriptId
    1   ignored_field  CDS             21  2148  .   -   1   Parent=transcript:TranscriptId
    1   ignored_field  five_prime_UTR  210 2148  .   -   .   Parent=transcript:TranscriptId

a

m

r

R

s

    # Basic usage
    bcftools csq -f hs37d5.fa -g Homo_sapiens.GRCh37.82.gff3.gz in.vcf -Ob -o out.bcf
    # Extract the translated haplotype consequences. The following TBCSQ variations
    # are recognised:
    #   %TBCSQ    .. print consequences in all haplotypes in separate columns
    #   %TBCSQ{0} .. print the first haplotype only
    #   %TBCSQ{1} .. print the second haplotype only
    #   %TBCSQ{*} .. print a list of unique consequences present in either haplotype
    bcftools query -f'[%CHROM\t%POS\t%SAMPLE\t%TBCSQ\n]' out.bcf

    # Two separate VCF records at positions 2:122106101 and 2:122106102
    # change the same codon. This UV-induced C>T dinucleotide mutation
    # has been annotated fully at the position 2:122106101 with
    #   - consequence type
    #   - gene name
    #   - ensembl transcript ID
    #   - coding strand (+ fwd, - rev)
    #   - amino acid position (in the coding strand orientation)
    #   - list of corresponding VCF variants
    # The annotation at the second position gives the position of the full
    # annotation
    BCSQ=missense|CLASP1|ENST00000545861|-|1174P>1174L|122106101G>A+122106102G>A
    BCSQ=@122106101
    # A frame-restoring combination of two frameshift insertions C>CG and T>TGG
    BCSQ=@46115084
    BCSQ=inframe_insertion|COPZ2|ENST00000006101|-|18AGRGP>18AQAGGP|46115072C>CG+46115084T>TGG
    # Stop gained variant
    BCSQ=stop_gained|C2orf83|ENST00000264387|-|141W>141*|228476140C>T
    # The consequence type of a variant downstream from a stop are prefixed with *
    BCSQ=*missense|PER3|ENST00000361923|+|1028M>1028T|7890117T>C

The SNPs at positions 1 and 7 are filtered, positions 0 and 8 are not:
         0123456789
    ref  .G.GT..G..
    del  .A.G-..A..
Here the positions 1 and 6 are filtered, 0 and 7 are not:
         0123-456789
    ref  .G.G-..G..
    ins  .A.GT..A..

The second indel is filtered:
         012345678901
    ref  .GT.GT..GT..
    del  .G-.G-..G-..
And similarly here, the second is filtered:
         01 23 456 78
    ref  .A-.A-..A-..
    ins  .AT.AT..AT..

   # Check discordance of all samples from B against all sample in A
   bcftools gtcheck -g A.bcf B.bcf
   # Limit comparisons to the fiven list of samples
   bcftools gtcheck -s gt:a1,a2,a3 -s qry:b1,b2 -g A.bcf B.bcf
   # Compare only two pairs a1,b1 and a1,b2
   bcftools gtcheck -p a1,b1,a1,b2 -g A.bcf B.bcf

    bcftools isec -p dir A.vcf.gz B.vcf.gz

    bcftools isec -e'MAF<0.01' -i'dbSNP=1' -e- A.vcf.gz B.vcf.gz C.vcf.gz -n +2 -p dir

    bcftools isec -p dir -n=2 -w1 A.vcf.gz B.vcf.gz

    bcftools isec -p dir -n-1 -c all A.vcf.gz B.vcf.gz

    bcftools isec -n~1100 -c all A.vcf.gz B.vcf.gz C.vcf.gz D.vcf.gz

-m none   ..  no new multiallelics, output multiple records instead
-m snps   ..  allow multiallelic SNP records
-m indels ..  allow multiallelic indel records
-m both   ..  both SNP and indel records can be multiallelic
-m all    ..  SNP records can be merged with indel records
-m id     ..  merge by ID

By default mpileup uses heuristics to decide when to apply the BAQ algorithm. Most sequences will not be BAQ adjusted, giving a CPU time closer to --no-BAQ, but it will still be applied in regions with suspected problematic alignments. This has been tested to work well on single sample data with even allele frequency, but the reliability is unknown for multi-sample calling and for low allele frequency variants so full BAQ is still recommended in those scenarios.

    RG_ID_1
    RG_ID_2  SAMPLE_A
    RG_ID_3  SAMPLE_A
    RG_ID_4  SAMPLE_B
    RG_ID_5  FILE_1.bam  SAMPLE_A
    RG_ID_6  FILE_2.bam  SAMPLE_A
    *        FILE_3.bam  SAMPLE_C
    ?        FILE_3.bam  SAMPLE_D

FORMAT/AD   .. Allelic depth (Number=R,Type=Integer)
FORMAT/ADF  .. Allelic depths on the forward strand (Number=R,Type=Integer)
FORMAT/ADR  .. Allelic depths on the reverse strand (Number=R,Type=Integer)
FORMAT/DP   .. Number of high-quality bases (Number=1,Type=Integer)
FORMAT/SP   .. Phred-scaled strand bias P-value (Number=1,Type=Integer)
FORMAT/SCR  .. Number of soft-clipped reads (Number=1,Type=Integer)
INFO/AD     .. Total allelic depth (Number=R,Type=Integer)
INFO/ADF    .. Total allelic depths on the forward strand (Number=R,Type=Integer)
INFO/ADR    .. Total allelic depths on the reverse strand (Number=R,Type=Integer)
INFO/SCR    .. Number of soft-clipped reads (Number=1,Type=Integer)
FORMAT/DV   .. Deprecated in favor of FORMAT/AD; Number of high-quality non-reference bases, (Number=1,Type=Integer)
FORMAT/DP4  .. Deprecated in favor of FORMAT/ADF and FORMAT/ADR; Number of high-quality ref-forward, ref-reverse,
               alt-forward and alt-reverse bases (Number=4,Type=Integer)
FORMAT/DPR  .. Deprecated in favor of FORMAT/AD; Number of high-quality bases for each observed allele (Number=R,Type=Integer)
INFO/DPR    .. Deprecated in favor of INFO/AD; Number of high-quality bases for each observed allele (Number=R,Type=Integer)

    bcftools mpileup -Ou -f ref.fa aln.bam | \
    bcftools call -Ou -mv | \
    bcftools filter -s LowQual -e '%QUAL<20 || DP>100' > var.flt.vcf

    # Before atomization:
    100  CC  C,GG   1/2
    # After:
    #   bcftools norm -a .
    100         C         G      ./1
    100         CC         C      1/.
    101         C         G      ./1
    # After:
    #   bcftools norm -a '*'
    #   bcftools norm -a \*
    100         C         G,*    2/1
    100         CC         C,*    1/2
    101         C         G,*    2/1

By default, appropriate system directories are searched for installed plugins. You can override this by setting the BCFTOOLS_PLUGINS environment variable to a colon-separated list of directories to search. If BCFTOOLS_PLUGINS begins with a colon, ends with a colon, or contains adjacent colons, the system directories are also searched at that position in the list of directories.

# List options common to all plugins
bcftools plugin
# List available plugins
bcftools plugin -l
# Run a plugin
bcftools plugin counts in.vcf
# Run a plugin using the abbreviated "+" notation
bcftools +counts in.vcf
# Run a plugin from an explicit location
bcftools +/path/to/counts.so in.vcf
# The input VCF can be streamed just like in other commands
cat in.vcf | bcftools +counts
# Print usage information of plugin "dosage"
bcftools +dosage -h
# Replace missing genotypes with 0/0
bcftools +missing2ref in.vcf
# Replace missing genotypes with 0|0
bcftools +missing2ref in.vcf -- -p

// Short description used by 'bcftools plugin -l'
const char *about(void);
// Longer description used by 'bcftools +name -h'
const char *usage(void);
// Called once at startup, allows initialization of local variables.
// Return 1 to suppress normal VCF/BCF header output, -1 on critical
// errors, 0 otherwise.
int init(int argc, char **argv, bcf_hdr_t *in_hdr, bcf_hdr_t *out_hdr);
// Called for each VCF record, return NULL to suppress the output
bcf1_t *process(bcf1_t *rec);
// Called after all lines have been processed to clean up
void destroy(void);

%CHROM          The CHROM column (similarly also other columns: POS, ID, REF, ALT, QUAL, FILTER)
%END            End position of the REF allele
%END0           End position of the REF allele in 0-based coordinates
%FIRST_ALT      Alias for %ALT{0}
%FORMAT         Prints all FORMAT fields or a subset of samples with -s or -S
%GT             Genotype (e.g. 0/1)
%INFO           Prints the whole INFO column
%INFO/TAG       Any tag in the INFO column
%IUPACGT        Genotype translated to IUPAC ambiguity codes (e.g. M instead of C/A)
%LINE           Prints the whole line
%MASK           Indicates presence of the site in other files (with multiple files)
%N_PASS(expr)   Number of samples that pass the filtering expression (see *<<expressions,EXPRESSIONS>>*)
%POS0           POS in 0-based coordinates
%PBINOM(TAG)    Calculate phred-scaled binomial probability, the allele index is determined from GT
%SAMPLE         Sample name
%TAG{INT}       Curly brackets to print a subfield (e.g. INFO/TAG{1}, the indexes are 0-based)
%TBCSQ          Translated FORMAT/BCSQ. See the csq command above for explanation and examples.
%TGT            Translated genotype (e.g. C/A)
%TYPE           Variant type (REF, SNP, MNP, INDEL, BND, OTHER)
[]              Format fields must be enclosed in brackets to loop over all samples
\n              new line
\t              tab character

Everything else is printed verbatim.

# Print chromosome, position, ref allele and the first alternate allele
bcftools query -f '%CHROM  %POS  %REF  %ALT{0}\n' file.vcf.gz

# Similar to above, but use tabs instead of spaces, add sample name and genotype
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT[\t%SAMPLE=%GT]\n' file.vcf.gz

# Print FORMAT/GT fields followed by FORMAT/GT fields
bcftools query -f 'GQ:[ %GQ] \t GT:[ %GT]\n' file.vcf

# Make a BED file: chr, pos (0-based), end pos (1-based), id
bcftools query -f'%CHROM\t%POS0\t%END\t%ID\n' file.bcf

# Print only samples with alternate (non-reference) genotypes
bcftools query -f'[%CHROM:%POS %SAMPLE %GT\n]' -i'GT="alt"' file.bcf

# Print all samples at sites with at least one alternate genotype
bcftools view -i'GT="alt"' file.bcf -Ou | bcftools query -f'[%CHROM:%POS %SAMPLE %GT\n]'

# Print phred-scaled binomial probability from FORMAT/AD tag for all heterozygous genotypes
bcftools query -i'GT="het"' -f'[%CHROM:%POS %SAMPLE %GT %PBINOM(AD)\n]' file.vcf

# Print the second value of AC field if bigger than 10. Note the (unfortunate) difference in
# index subscript notation: formatting expressions (-f) uses "{}" while filtering expressions
# (-i) use "[]". This is for historic reasons and backward-compatibility.
bcftools query -f '%AC{1}\n' -i 'AC[1]>10' file.vcf.gz

Notation:
  D  = Data, AZ = autozygosity, HW = Hardy-Weinberg (non-autozygosity),
  f  = non-ref allele frequency
Emission probabilities:
  oAZ = P_i(D|AZ) = (1-f)*P(D|RR) + f*P(D|AA)
  oHW = P_i(D|HW) = (1-f)^2 * P(D|RR) + f^2 * P(D|AA) + 2*f*(1-f)*P(D|RA)
Transition probabilities:
  tAZ = P(AZ|HW)  .. from HW to AZ, the -a parameter
  tHW = P(HW|AZ)  .. from AZ to HW, the -H parameter
  ci  = P_i(C)  .. probability of cross-over at site i, from genetic map
  AZi = P_i(AZ) .. probability of site i being AZ/non-AZ, scaled so that AZi+HWi = 1
  HWi = P_i(HW)
  P_{i+1}(AZ) = oAZ * max[(1 - tAZ * ci) * AZ{i-1} , tAZ * ci * (1-AZ{i-1})]
  P_{i+1}(HW) = oHW * max[(1 - tHW * ci) * (1-AZ{i-1}) , tHW * ci * AZ{i-1}]

    bcftools query -f'%CHROM\t%POS\t%REF,%ALT\t%INFO/TAG\n' file.vcf | bgzip -c > freqs.tab.gz

    tabix -s1 -b2 -e3 file.gz

# Generate the stats
bcftools stats -s - > file.vchk

# Plot the stats
plot-vcfstats -p outdir file.vchk

# The final looks can be customized by editing the generated
# 'outdir/plot.py' script and re-running manually
cd outdir && python plot.py && pdflatex summary.tex