bcftools Primer

Concatenate multiple files together

If we wanted to concatenate (i.e. “stack”) multiple .vcf files together we can use bcftools concat, so long as the input files share the same fields. In this example, we’ll combine all of the chromosomes (1-22) into a single file.

The --output-type z argument specifies that the output will be compressed, and the --output flag allows us to explicitly name the resulting file:

bcftools concat ALL.chr*.clinvar.vcf.gz --output-type z --output all.clinvar.vcf.gz

NOTE: bcftools concat is not equivalent bcftools merge. For an example of the latter see below.

Select individual samples by name

bcftools view -s allows for subsetting by sample ID.

The combined all.clinvar.vcf.gz file above contains multiple samples. Here we’ll create individual compressed .vcf files for NA20536 and HG03718 samples, along with a tabix index for each file (using bcftools index -t):

bcftools view -s NA20536 all.clinvar.vcf.gz --output-type z --output NA20536.clinvar.vcf.gz
bcftools view -s HG03718 all.clinvar.vcf.gz --output-type z --output HG03718.clinvar.vcf.gz

## note: bcftools index -t is equivalent to tabix here
bcftools index -t NA20536.clinvar.vcf.gz
bcftools index -t HG03718.clinvar.vcf.gz

Filter to only include INDELs

bcftools view -v will restrict the file to specified variant types: “snps”, “indels”, “mnps”, or “other”.

We can use the command to filter the .vcf to only include INDELs:

bcftools view -v indels NA20536.clinvar.vcf.gz

Filter by rsid

With bcftools you can filter a .vcf file for certain sites by passing in a file that contains the IDs to be retained.

Assuming we have the following RSIDs in a file called snps.list²:

rs145413551
rs34610323
rs79548709
rs371163239
rs148716910
rs374704178

We can use snps.list to filter with bcftools view:

bcftools view --include ID==@snps.list NA20536.clinvar.vcf.gz

Filter by chromosome and/or position

The --regions flag takes input chromosome and/or position coordinates to filter the .vcf.

If we wanted to restrict to chromosome 5:

bcftools view --regions 5 NA20536.vcf.gz 

And if we were interested in a specific region (let’s say chromosome 10, anywhere between positions 800000:900000):

bcftools view --regions 10:800000-900000 NA20536.clinvar.vcf.gz

Format translated genotype output

bcftools query will output contents of the .vcf in text format. The contents can be specified in a string that includes fields to extract, separators, and line endings.

In this scenario, we’ll pull out the ID (RSID), chromosome, position, a translated genotype, and the “type” (SNP, INDEL, etc.) in tab-separated format:

bcftools query -f "%ID\t%CHROM\t%POS[\t%TGT]\t%TYPE\n" NA20536.clinvar.vcf.gz

Merge vcf files together

bcftools merge will combine data from multiple files.

To merge individual sample .vcf files into one:

bcftools merge NA20536.clinvar.vcf.gz HG03718.clinvar.vcf.gz --output-type z --output NA20536.HG03718.clinvar.vcf.gz

Parse genotypes for multiple samples

Given that you have a mutli-sample .vcf you can parse genotypes for each individual:

bcftools query -f '[%CHROM\t%POS\t%SAMPLE\t%TGT\n]' NA20536.HG03718.clinvar.vcf.gz

Edit chromosome names

You can rename chromosomes with bcftools annotate --rename-chrs. The command requires that you supply a tab-separated file with the desired naming convention, organized as “old\tnew”:

1\tchr1
2\tchr2
3\tchr3
4\tchr4
5\tchr5
6\tchr6
7\tchr7
8\tchr8
9\tchr9
10\tchr10
11\tchr11
12\tchr12
13\tchr13
14\tchr14
15\tchr15
16\tchr16
17\tchr17
18\tchr18
19\tchr19
20\tchr20
21\tchr21
22\tchr22

bcftools annotate --rename-chrs chromosomes.txt NA20536.clinvar.vcf.gz

View without header

To view only the results without header (i.e. remove the header) use the -H flag:

bcftools view -H NA20536.HG03718.clinvar.vcf.gz

View only header

To view only the header (i.e. extract header) use the -h flag:

bcftools view -h clinvar.vcf.gz

## ##fileformat=VCFv4.1
## ##FILTER=<ID=PASS,Description="All filters passed">
## ##fileDate=2020-02-17
## ##source=ClinVar
## ##reference=GRCh37
## ##ID=<Description="ClinVar Variation ID">
## ##INFO=<ID=AF_ESP,Number=1,Type=Float,Description="allele frequencies from GO-ESP">
## ##INFO=<ID=AF_EXAC,Number=1,Type=Float,Description="allele frequencies from ExAC">
## ##INFO=<ID=AF_TGP,Number=1,Type=Float,Description="allele frequencies from TGP">
## ##INFO=<ID=ALLELEID,Number=1,Type=Integer,Description="the ClinVar Allele ID">
## ##INFO=<ID=CLNDN,Number=.,Type=String,Description="ClinVar's preferred disease name for the concept specified by disease identifiers in CLNDISDB">
## ##INFO=<ID=CLNDNINCL,Number=.,Type=String,Description="For included Variant : ClinVar's preferred disease name for the concept specified by disease identifiers in CLNDISDB">
## ##INFO=<ID=CLNDISDB,Number=.,Type=String,Description="Tag-value pairs of disease database name and identifier, e.g. OMIM:NNNNNN">
## ##INFO=<ID=CLNDISDBINCL,Number=.,Type=String,Description="For included Variant: Tag-value pairs of disease database name and identifier, e.g. OMIM:NNNNNN">
## ##INFO=<ID=CLNHGVS,Number=.,Type=String,Description="Top-level (primary assembly, alt, or patch) HGVS expression.">
## ##INFO=<ID=CLNREVSTAT,Number=.,Type=String,Description="ClinVar review status for the Variation ID">
## ##INFO=<ID=CLNSIG,Number=.,Type=String,Description="Clinical significance for this single variant">
## ##INFO=<ID=CLNSIGCONF,Number=.,Type=String,Description="Conflicting clinical significance for this single variant">
## ##INFO=<ID=CLNSIGINCL,Number=.,Type=String,Description="Clinical significance for a haplotype or genotype that includes this variant. Reported as pairs of VariationID:clinical significance.">
## ##INFO=<ID=CLNVC,Number=1,Type=String,Description="Variant type">
## ##INFO=<ID=CLNVCSO,Number=1,Type=String,Description="Sequence Ontology id for variant type">
## ##INFO=<ID=CLNVI,Number=.,Type=String,Description="the variant's clinical sources reported as tag-value pairs of database and variant identifier">
## ##INFO=<ID=DBVARID,Number=.,Type=String,Description="nsv accessions from dbVar for the variant">
## ##INFO=<ID=GENEINFO,Number=1,Type=String,Description="Gene(s) for the variant reported as gene symbol:gene id. The gene symbol and id are delimited by a colon (:) and each pair is delimited by a vertical bar (|)">
## ##INFO=<ID=MC,Number=.,Type=String,Description="comma separated list of molecular consequence in the form of Sequence Ontology ID|molecular_consequence">
## ##INFO=<ID=ORIGIN,Number=.,Type=String,Description="Allele origin. One or more of the following values may be added: 0 - unknown; 1 - germline; 2 - somatic; 4 - inherited; 8 - paternal; 16 - maternal; 32 - de-novo; 64 - biparental; 128 - uniparental; 256 - not-tested; 512 - tested-inconclusive; 1073741824 - other">
## ##INFO=<ID=RS,Number=.,Type=String,Description="dbSNP ID (i.e. rs number)">
## ##INFO=<ID=SSR,Number=1,Type=Integer,Description="Variant Suspect Reason Codes. One or more of the following values may be added: 0 - unspecified, 1 - Paralog, 2 - byEST, 4 - oldAlign, 8 - Para_EST, 16 - 1kg_failed, 1024 - other">
## ##contig=<ID=1>
## ##contig=<ID=2>
## ##contig=<ID=3>
## ##contig=<ID=4>
## ##contig=<ID=5>
## ##contig=<ID=6>
## ##contig=<ID=7>
## ##contig=<ID=8>
## ##contig=<ID=9>
## ##contig=<ID=10>
## ##contig=<ID=11>
## ##contig=<ID=12>
## ##contig=<ID=13>
## ##contig=<ID=14>
## ##contig=<ID=15>
## ##contig=<ID=16>
## ##contig=<ID=17>
## ##contig=<ID=18>
## ##contig=<ID=19>
## ##contig=<ID=20>
## ##contig=<ID=21>
## ##contig=<ID=22>
## ##contig=<ID=X>
## ##contig=<ID=Y>
## ##contig=<ID=MT>
## ##bcftools_viewVersion=1.10.2-27-g9d66868+htslib-1.10.2-33-g1bbcd02
## ##bcftools_viewCommand=view -h clinvar.vcf.gz; Date=Fri Feb 28 19:06:40 2020
## #CHROM   POS ID  REF ALT QUAL    FILTER  INFO