readme.Rmd

---
title: "Learning the VCF format"
output: github_document
---

```{r setup, include=FALSE}
Sys.setenv( PATH=paste0(Sys.getenv("PATH"), ":", getwd(), "/tools/bin") )
Sys.setenv( MANPATH=paste0(Sys.getenv("MANPATH"), ":", getwd(), "/tools/share/man") )
knitr::opts_chunk$set(echo = TRUE)
```

## Introduction

![Build README](https://github.com/davetang/learning_vcf_file/actions/workflows/create_readme.yml/badge.svg)
![Creating VCF examples](https://github.com/davetang/learning_vcf_file/actions/workflows/variant_call.yml/badge.svg)

Natural selection occurs under three conditions:

1. Genetic variation occurs among individuals (and this occurs mainly due to chance errors in replication)
2. The genetic variation must be heritable, i.e. passed from one generation to the next, and
3. The genetic variation results in varying fitness, i.e. individuals survive and reproduce with respect to genetic variation

The _de facto_ file format for storing genetic variation is the Variant Call Format (VCF) and was developed under the [1000 Genomes Project](https://pubmed.ncbi.nlm.nih.gov/21653522/). Currently, the Large Scale Genomics work stream of the Global Alliance for Genomics & Health (GA4GH) maintain the specification of the VCF (and other [high-throughput sequencing data formats](https://samtools.github.io/hts-specs/)). A good starting point for learning about the VCF is this [poster](http://vcftools.sourceforge.net/VCF-poster.pdf) and a portion of the poster is displayed below.

![VCF format](img/vcf_format.png)

Initially, VCFtools (and the associated scripts) was used for working with VCF files. VCFtools was also developed by the 1000 Genomes Project team but the [tool](https://github.com/vcftools/vcftools) does not seem to be actively maintained anymore. As VCF files are simply tab-delimited flat files, they are slow to process and the BCF was implemented, which is a more efficient format for data processing. A BCF file is the binary equivalent of a VCF file, akin to the SAM and BAM formats and [BCFtools](https://github.com/samtools/bcftools) is used to work with BCF (and VCF) files. BCFtools is actively maintained and therefore should be used instead of VCFtools. To learn more about BCFtools (and SAMtools), check out the paper [Twelve years of SAMtools and BCFtools](https://academic.oup.com/gigascience/article/10/2/giab008/6137722) and please cite it if you use BCFtools for your work.

Lastly, this README is created by `create_readme.sh` using GitHub Actions, which executes `readme.Rmd` and creates `README.md` and can also be viewed at <https://davetang.github.io/learning_vcf_file/>.

## Installation

The easiest way to install BCFtools is by using [Conda](https://docs.conda.io/en/latest/) (and I recommend using [Miniconda](https://docs.conda.io/en/latest/miniconda.html)).

```bash
conda install -c bioconda bcftools
```

It is also relatively straightforward to compile on Linux (if your system has all the prerequisites). The following code will install BCFtools (and HTSlib and SAMtools, which you will most likely be using as well) into the directory specified by `dir=`.

```bash
dir=$HOME/tools

ver=1.18
for tool in htslib bcftools samtools; do
   check=${tool}
   if [[ ${tool} == htslib ]]; then
      check=bgzip
   fi
   if [[ ! -e ${dir}/bin/${check} ]]; then
      url=https://github.com/samtools/${tool}/releases/download/${ver}/${tool}-${ver}.tar.bz2
      wget ${url}
      tar xjf ${tool}-${ver}.tar.bz2
      cd ${tool}-${ver}
      ./configure --prefix=${dir}
      make && make install
      cd ..
      rm -rf ${tool}-${ver}*
   fi
done
```

## Creating VCF example files

Example VCF files were generated to test the functionality of BCFtool and other VCF tools. The `aln.bt.vcf.gz`, `aln.hc.vcf.gz`, and `aln.fb.vcf.gz` VCF files were generated using a simple workflow implemented in [Groovy](http://www.groovy-lang.org/) and processed by [Bpipe](https://github.com/ssadedin/bpipe), a tool for running bioinformatics workflows. You can view the workflow in `workflow/simple/pipeline.groovy`, which carries out the following steps: (i) a random reference sequence is generated, (ii) the reference sequence is mutated, (iii) reads are derived from the mutated reference then (iv) mapped back to the original non-mutated reference (v) and finally variants are called using three separate tools: BCFtools, HaplotypeCaller, and freebayes. Additional information is available in this [blog post](http://davetang.org/muse/2015/06/04/paired-end-alignment-using-bpipe/).

To run this workflow for yourself, run the following commands:

```bash
git clone https://github.com/davetang/learning_vcf_file.git

# installs the necessary tools
cd workflow && ./install.sh

# run the workflow using Bpipe
cd simple && make
```

If you want to play around with different parameters, simply edit the variables at the start of `workflow/simple/pipeline.groovy`. The status of this workflow is checked by [GitHub Actions](https://github.com/features/actions) and on successful completion the final variant calls are copied to the `eg` folder; please refer to `.github/workflows/variant_call.yml` for more details.

## Usage

Running `bcftools` without any parameters will output the usage and the subcommands. (The `--help` option is used below to avoid an exit code of 1 since this README is generated programmatically.)

```{bash engine.opts='-l'}
bcftools --help
```

## Getting help

If BCFtools was [installed]($installation) by compiling the source code with the path to the man pages added to `MANPATH`, you can simply run `man bcftools` to access the manual pages for getting help. If BCFtools was installed using `Conda`, the `MANPATH` environment variable needs to be set [accordingly](https://github.com/conda/conda/issues/845#issuecomment-525223443). The man pages provides additional information for each subcommand and its parameters.

```{bash engine.opts='-l'}
man bcftools | grep "LIST OF COMMANDS" -B 1984
```

You can also [open an issue](https://github.com/davetang/learning_vcf_file/issues) and I will try to answer your question.

## VCF to BCF and other conversions

Use `bcftools view` or `bcftools convert` to make conversions and use `--output-type u` or `--output-type b` for uncompressed and compressed BCF, respectively.

```{bash engine.opts='-l'}
bcftools view -O u eg/aln.bt.vcf.gz > eg/aln.bt.un.bcf
bcftools view -O b eg/aln.bt.vcf.gz > eg/aln.bt.bcf

# uncompressed VCF
bcftools convert -O v eg/aln.bt.vcf.gz > eg/aln.bt.vcf
```

`bcftools convert` has more conversion options such as converting gVCF to VCF and converting VCFs from `IMPUTE2` and `SHAPEIT`.

## Viewing a VCF/BCF file

Use `bcftools view` to view a VCF, bgzipped VCF, or BCF file; BCFtools will automatically detect the format.

Viewing a VCF file.

```{bash engine.opts='-l'}
bcftools view eg/aln.bt.vcf | grep -v "^#" | head -2
```

Viewing a bgzipped VCF file.

```{bash engine.opts='-l'}
bcftools view eg/aln.bt.vcf.gz | grep -v "^#" | head -2
```

Viewing a BCF file.

```{bash engine.opts='-l'}
bcftools view eg/aln.bt.vcf.gz | grep -v "^#" | head -2
```

If you want to omit the header but do not want to type `grep -v "^#"`, you can use the `-H` option instead.

```{bash engine.opts='-l'}
bcftools view -H eg/aln.bt.vcf.gz | head -2
```

## Comparing output types

The four output types will be compared on a VCF file produced by the 1000 Genomes Project.

1. Compressed BCF (b)
2. Uncompressed BCF (u)
3. Compressed VCF (z)
4. Uncompressed VCF (v).

Download and unzip.

```{bash engine.opts='-l'}
wget --quiet https://davetang.org/file/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5_related_samples.20130502.genotypes.vcf.gz -O eg/1kgp.vcf.gz
gunzip eg/1kgp.vcf.gz
```

VCF to compressed BCF.

```{bash engine.opts='-l'}
time bcftools convert --threads 2 -O b -o eg/1kgp.bcf eg/1kgp.vcf
```

VCF to uncompressed BCF.

```{bash engine.opts='-l'}
time bcftools convert --threads 2 -O u -o eg/1kgp.un.bcf eg/1kgp.vcf
```

VCF to compressed VCF.

```{bash engine.opts='-l'}
time bcftools convert --threads 2 -O z -o eg/1kgp.vcf.gz eg/1kgp.vcf
```

File sizes

```{bash engine.opts='-l'}
du -h eg/1kgp.*
```

## VCF to PED

See my [blog post](http://davetang.org/muse/2016/07/28/vcf-to-ped/).

## VCF to BED

VCF to [Browser Extensible Data](https://www.genome.ucsc.edu/FAQ/FAQformat.html#format1) format and not [Binary PED](http://zzz.bwh.harvard.edu/plink/data.shtml#bed) format; for Binary PED, see the `plink` directory in this repo. For Browser Extensible Data use [BEDOPS](https://bedops.readthedocs.io/en/latest/index.html), specifically the [vcf2bed](https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/vcf2bed.html) tool.

Install using Conda from [Bioconda](https://anaconda.org/bioconda/bedops).

```bash
conda install -c bioconda bedops
```

Check out example file.

```{bash engine.opts='-l'}
cat eg/ex.vcf | grep -v "^#" | cut -f1-6
```

Convert to BED.

```{bash engine.opts='-l'}
vcf2bed < eg/ex.vcf | cut -f1-3
```

Note above that the deletion (rs57524763) only has 1 nucleotide but the reference should be 3 nucleotides long. Use `--deletions` to have the coordinates reflect the length of the reference sequence (and to only report deletions).

```{bash engine.opts='-l'}
vcf2bed --deletions < eg/ex.vcf | cut -f1-3
```

There is also the insertions option to report only insertions

```{bash engine.opts='-l'}
vcf2bed --insertions < eg/ex.vcf | cut -f1-3
```

To report SNVs use `snvs` but note the tool reports the MNP as a SNV and the reference length is not 2 nucleotides long.

```{bash engine.opts='-l'}
vcf2bed --snvs < eg/ex.vcf | cut -f1-3
```

## Extracting info from columns

Use the `query` subcommand to extract information from the `INFO` or `FORMAT` column.

```{bash engine.opts='-l'}
bcftools query -f 'DP=%DP\tAN=%AN\tAC=%AC\tMQ=%MQ\n' eg/aln.bt.vcf.gz | head -3
```

To extract from the `FORMAT` column, the FORMAT fields must be enclosed in square brackets, e.g. "[ %AD]".

```{bash engine.opts='-l'}
bcftools query -f 'AD=[%AD]\tGQ=[%GQ]\tAF=%AF\n' eg/ex.vcf
```

`query` can read from `STDIN` as with most of the BCFtools commands. The example below extracts information from the `INFO` column for only SNPs.

```{bash engine.opts='-l'}
bcftools view -v snps eg/aln.bt.vcf.gz | bcftools query -f 'DP=%DP\tAN=%AN\tAC=%AC\tMQ=%MQ\n' - | head -3
```

## Filtering

Many commands of BCFtools can be used for [filtering variants](https://samtools.github.io/bcftools/howtos/filtering.html) by using the `-f, --apply-filters` or the `-i, --include` and `-e, --exclude` parameters. We will download another example file to demonstrate how to perform filtering.

```{bash engine.opts='-l'}
wget --quiet https://davetang.org/file/ERR031940.merged.scored.filtered.bcf -O eg/ERR031940.bcf
```

To perform filtering on the `FILTER` column of the VCF file, use `bcftools view` and the `-f, --apply-filters` parameter.

    -f,   --apply-filters <list> require at least one of the listed FILTER strings (e.g. "PASS,.")

Keep variants with `CNN_2D_SNP_Tranche_99.90_100.00` in the `FILTER` column.

```{bash engine.opts='-l'}
bcftools view -H -f CNN_2D_SNP_Tranche_99.90_100.00 eg/ERR031940.bcf | wc -l
```

Keep variants with `CNN_2D_SNP_Tranche_99.90_100.00` and `CNN_2D_INDEL_Tranche_99.00_100.00` in the `FILTER` column.

```{bash engine.opts='-l'}
bcftools view -H -f CNN_2D_SNP_Tranche_99.90_100.00,CNN_2D_INDEL_Tranche_99.00_100.00 eg/ERR031940.bcf | wc -l
```

Keep variants that have passed filters, i.e. contain `PASS` in the `FILTER` column.

```{bash engine.opts='-l'}
bcftools view -H -f PASS eg/ERR031940.bcf | head -3
```

The `filter` command can achieve the same filtering step as above.

```{bash engine.opts='-l'}
bcftools filter -i 'FILTER="PASS"' eg/ERR031940.bcf | bcftools view -H - | head -3
```

`bcftools query` can also perform the same filtering using `-i, --include` but a format must be specified.

```{bash engine.opts='-l'}
bcftools query -i 'FILTER="PASS"' -f '%CHROM %POS %FILTER\n' eg/ERR031940.bcf | head -3
```

Which command should be used for filtering? I would recommend using `bcftools filter` for filtering since it clearly states that you are performing filtering in the command name (it is clearer to someone not familiar with BCFtools) and the `-i` and `-e` expressions are flexible enough to perform most filtering requirements. Use `bcftools query` if you want to perform filtering and modify the output. See `man bcftools` for more information about `EXPRESSIONS`, which I recommend reading.

### Filtering variant types

Use `TYPE` in the filtering expression to subset specific types of variants; types include: indel, snp, mnp, ref, bnd, other, and overlap.

SNPs.

```{bash engine.opts='-l'}
bcftools filter -i 'TYPE="snp"' eg/1kgp.bcf | bcftools view -H - | cut -f1-8 | head -2
```

Indels.

```{bash engine.opts='-l'}
bcftools filter -i 'TYPE="indel"' eg/1kgp.bcf | bcftools view -H - | cut -f1-8 | head -2
```

[Multiple Nucleotide Polymorphisms](https://genome.sph.umich.edu/wiki/Variant_classification#Definitions): The reference and alternate sequences are of the same length and have to be greater than 1 and all nucleotides in the sequences differ from one another.

```{bash engine.opts='-l'}
bcftools filter -i 'TYPE="mnps"' eg/ex.vcf | bcftools view -H - | cut -f1-8
```

Not sure what `ref` refers to but possibly structural variants. `INS:ME` refers to an insertion of a mobile element relative to the reference and the `<CN#>` refers to `Copy number allele: # copies` according to the VCF header.

```{bash engine.opts='-l'}
bcftools view -v ref eg/1kgp.bcf | bcftools view -H - | cut -f1-8 | head -3

# filter returns nothing but view does
# bcftools filter -i 'TYPE="ref"' eg/1kgp.bcf | bcftools view -H - | cut -f1-8 | head -3
```

Breakends (no variants of this class).

```{bash engine.opts='-l'}
bcftools filter -i 'TYPE="bnd"' eg/1kgp.bcf | bcftools view -H -
```

Others.

```{bash engine.opts='-l'}
bcftools filter -i 'TYPE="other"' eg/1kgp.bcf | bcftools view -H - | cut -f1-8 | head -3
```

Also not sure about `overlap` either.

```{bash engine.opts='-l'}
bcftools filter -i 'TYPE="overlap"' eg/1kgp.bcf | bcftools view -H -
```

### Filtering genotypes

Use `GT` in the filtering expression to subset specific genotypes.

```{bash engine.opts='-l'}
bcftools filter -i 'GT="1|1"' eg/1kgp.bcf | bcftools view -H - | head -3
```

Genotypes using 2nd alternate allele.

```{bash engine.opts='-l'}
bcftools filter -i 'GT~"2"' eg/1kgp.bcf | bcftools view -H - | head -3
```

You can also use abbreviations such as `GT="AA"`. See below for list of abbreviations.

```{bash engine.opts='-l'}
bcftools filter -i 'GT="AA"|GT="RA"' eg/1kgp.bcf | bcftools view -H - | head -3
```

Below are the list of genotypes you can use:

* `GT="ref"` = reference
* `GT="alt"` = alternate
* `GT="mis"` = missing genotype
* `GT="hom"` = homozygous
* `GT="het"` = heterozygous
* `GT="hap"` = haploid
* `GT="RR"` = ref-ref hom
* `GT="AA"` = alt-alt hom
* `GT="RA"` or GT="AR" = ref-alt het
* `GT="Aa"` or GT="aA" = alt-alt het
* `GT="R"` = haploid ref
* `GT="A"` = haploid alt

### Filtering INFO field/s

Use `bcftools filter` to filter out (`-e` or `--exclude`) variants. In the example below we are filtering out variants that have a depth of less than 200. We also use the `-s` parameter to name our filter and this name will be displayed in the `FILTER` column.

```{bash engine.opts='-l'}
bcftools filter -s "Depth200" -e "DP<200" eg/aln.bt.vcf.gz | grep -v "^##" | head -4
```

Use `&` to combine several other criteria.

```{bash engine.opts='-l'}
bcftools filter -s "Depth200&VDB" -e "DP<200 & VDB<0.9" eg/aln.bt.vcf.gz | grep -v "^##" | head -4
```

Note that `&` and `&&` have different meanings in the expression and make a difference if you have more than one sample in your VCF file. Using `&` requires that the conditions are satisfied within one sample and `&&` means that conditions can be satisfied across samples.

The `-e` or `-i` parameters accept an expression and we can use it to perform calculations. In the example below, variants with a QD divided by DP ratio of less than 0.3 are labelled with `QD/DP`; this ratio was created for illustrative purposes only. Note that `QD` and `DP` are prefixed with `INFO/`; this was done to explicitly state that we want the QD and DP values from the INFO field, since there is also a DP value in the FORMAT field. (This VCF file is different from the first filtering example, which only had one DP value.)

```{bash engine.opts='-l'}
bcftools filter -s "QD/DP" -e "INFO/QD / INFO/DP < 0.3" eg/aln.hc.vcf.gz | grep -v "^##" | head -4
```

## Querying

The `query` and `view` commands can be used to query a VCF file.

### Output sample names

Use `bcftools query`.

```{bash engine.opts='-l'}
bcftools query -l eg/1kgp.vcf.gz | head -5
```

### Subset sample/s from a multi-sample VCF file

Subset HG00733.

```{bash engine.opts='-l'}
bcftools view -s HG00733 eg/1kgp.vcf.gz | grep -v "^##" | head -3
```

Subset HG00124, HG00501, HG00635, HG00702, and HG00733.

```{bash engine.opts='-l'}
bcftools view -s HG00124,HG00501,HG00635,HG00702,HG00733 eg/1kgp.vcf.gz | grep -v "^##" | head -3
```

### Subset variants within a specific genomic region

Use `bcftools view` with `-r` or `-R`, which requires an index file. You can use `bcftools view` with `-t` or `-T`, which does not require an index file, but is much slower because the entire file is streamed.

```{bash engine.opts='-l'}
tabix -f eg/1kgp.bcf
time bcftools view -H -r 1:55000000-56000000 eg/1kgp.bcf | wc -l
```

`bcftools view` with `-t` streams the entire file, so is much slower.

```{bash engine.opts='-l'}
time bcftools view -H -t 1:55000000-56000000 eg/1kgp.bcf | wc -l
```

Use commas to list more than one loci.

```{bash engine.opts='-l'}
bcftools view -H -r 1:10000-50000,1:100000-200000,1:55000000-56000000 eg/1kgp.bcf | wc -l
```

Or use a file (coordinates are 1-based, which is the same as the VCF) to store regions of interest.

```{bash engine.opts='-l'}
echo -e "1\t10000\t50000\n1\t100000\t200000\n1\t55000000\t56000000" > eg/regions.txt
bcftools view -H -R eg/regions.txt eg/1kgp.bcf | wc -l
```

### Random subset of variants

Total number of variants.

```{bash engine.opts='-l'}
bcftools view -H eg/aln.bt.vcf.gz | wc -l
```

A random sample can be achieved by using a Perl one-liner but this is a very slow approach since it goes through every variant line by line. In the example below, the `srand` function sets the seed (for reproducibility) and the float controls how many variants are outputted (1%). (Note that the use of `grep -v "^#"` is only line counting purposes.)

```{bash engine.opts='-l'}
bcftools view eg/aln.bt.vcf.gz | perl -nle 'BEGIN { srand(1984) } if (/^#/){ print; next }; print if rand(1) < 0.01' | grep -v "^#" | wc -l
```

Sub-sample 1% and save as BCF file.

```{bash engine.opts='-l'}
bcftools view eg/aln.bt.vcf.gz | perl -nle 'BEGIN { srand(1984) } if (/^#/){ print; next }; print if rand(1) < 0.01' | bcftools view -O b - -o eg/aln.bt.ss.bcf
```

Sample 10%.

```{bash engine.opts='-l'}
bcftools view eg/aln.bt.vcf.gz | perl -nle 'BEGIN { srand(1984) } if (/^#/){ print; next }; print if rand(1) < 0.1' | grep -v "^#" | wc -l
```

### Split an annotation field

The [split-vep](https://samtools.github.io/bcftools/howtos/plugin.split-vep.html) plugin can be used to split a structured field. `split-vep` was written to work with VCF files created by [VEP](https://github.com/Ensembl/ensembl-vep) or `bcftools csq`.

```{bash engine.opts='-l'}
bcftools +split-vep -h || true
```

#### VEP

An [example VCF file](https://github.com/davetang/vcf_example) that was annotated with VEP is available as `eg/S1.haplotypecaller.filtered_VEP.ann.vcf.gz`. To list the annotation fields use `-l`.

```{bash engine.opts='-l'}
bcftools +split-vep -l eg/S1.haplotypecaller.filtered_VEP.ann.vcf.gz | head
```

Use `-f` to print the wanted fields in your own specified format; variants without consequences are excluded.


```{bash engine.opts='-l'}
bcftools +split-vep -f '%CHROM:%POS,%ID,%Consequence\n' eg/S1.haplotypecaller.filtered_VEP.ann.vcf.gz | head
```

Limit output to missense or more severe variants.

```{bash engine.opts='-l'}
bcftools +split-vep -f '%CHROM:%POS,%ID,%Consequence\n' -s worst:missense+ eg/S1.haplotypecaller.filtered_VEP.ann.vcf.gz | head
```

#### BCFtools csq

An [example VCF file](https://github.com/davetang/vcf_example) that was annotated with BCFtools csq is available as `eg/S1.haplotypecaller.filtered.phased.csq.vcf.gz`. The tag added by `csq` is `INFO/BCSQ`, so we need to provide this to split-vep. To list the annotation fields use `-l`.

```{bash engine.opts='-l'}
bcftools +split-vep -a BCSQ -l eg/S1.haplotypecaller.filtered.phased.csq.vcf.gz
```

Use `-f` to print the wanted fields in your own specified format; variants without consequences are excluded.

```{bash engine.opts='-l'}
bcftools +split-vep -a BCSQ -f '%CHROM:%POS,%ID,%Consequence\n' eg/S1.haplotypecaller.filtered.phased.csq.vcf.gz | head
```

The `-d` or `--duplicate` is useful to output annotations per transcript/allele on a new line.

```{bash engine.opts='-l'}
bcftools +split-vep -a BCSQ -f '%transcript,%Consequence\n' eg/S1.haplotypecaller.filtered.phased.csq.vcf.gz | head
```

Use `-d` to split.

```{bash engine.opts='-l'}
bcftools +split-vep -a BCSQ -d -f '%transcript,%Consequence\n' eg/S1.haplotypecaller.filtered.phased.csq.vcf.gz | head
```

#### SnpEff

It is possible to use the split-vep plugin with [SnpEff](https://github.com/pcingola/SnpEff), another popular variant annotation tool with some modifications to the VCF header.

SnpEff provides annotations with the `ANN` tag.

```{bash engine.opts='-l'}
bcftools view -H eg/PRJNA784038_illumina.vt.ann.vcf.gz | head -1 | cut -f8
```

We can tell `split-vep` to use the `ANN` field using the `-a` parameter but this will result in an error because the header is not defined in the expected format. (The ` || true` is added to prevent GitHub Actions from failing and you do not normally need it.)

```{bash engine.opts='-l'}
bcftools +split-vep eg/PRJNA784038_illumina.vt.ann.vcf.gz -a ANN || true
```

The [expected format](https://www.ensembl.org/info/docs/tools/vep/vep_formats.html?redirect=no) is:

```
##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position">
#CHROM  POS       ID          REF  ALT  QUAL  FILTER  INFO
21      26978790  rs75377686  T    C    .     .       CSQ=C|missense_variant|MODERATE|MRPL39|ENSG00000154719|Transcript|ENST00000419219|protein_coding|2/8||ENST00000419219.1:c.251A>G|ENSP00000404426.1:p.Asn84Ser|260|251|84
```

We can use `bcftools reheader` to modify the header for `ANN` so that is `split-vep` compatible.

```{bash engine.opts='-l'}
bcftools head eg/PRJNA784038_illumina.vt.ann.vcf.gz | perl -nle "if (/ID=ANN,/){ s/[\s\']//g; s/Functionalannotations:/Format: /; s/Annotation_Impact/Consequence/ } print" > eg/new_header.txt
bcftools reheader -h eg/new_header.txt eg/PRJNA784038_illumina.vt.ann.vcf.gz -o eg/PRJNA784038_illumina.vt.ann.vep.bcf
```

We can now run `split-vep` and list the split annotation fields.

```{bash engine.opts='-l'}
bcftools +split-vep -a ANN -l eg/PRJNA784038_illumina.vt.ann.vep.bcf
```

This variant file contains Omicron variants called using [BCFtools](https://davetang.org/muse/2022/01/26/omicron-variants/). We can use `split-vep` and `awk` to output Spike protein variants that have a `LOW` or higher consequence. The `-d` parameter is used to output consequences on a new line rather than having comma-delimited entries.

```{bash engine.opts='-l'}
bcftools +split-vep -a ANN -d -f '%CHROM:%POS\t%Gene_Name\t%Consequence\t%HGVS.c\t%HGVS.p\n' eg/PRJNA784038_illumina.vt.ann.vep.bcf | awk '$2 == "S" && $3 != "MODIFIER" {print}' | head
```

## Summarise SNPs and INDELs per sample

Use `bcftools stats` with the `-s -` parameter. The example VCF file `eg/ex.vcf` has four variants across three samples (one, two, and three).

* Sample one has two SNPs (both het) and one deletion (het)
* Sample two has two SNPs (one het and one hom alt) and insertion (het)
* Sample three has one SNP (het) and one insertion (hom alt) and one deletion (hom alt)

You can confirm the numbers from the stats output.

```{bash engine.opts='-l'}
cat eg/ex.vcf | grep -v "^#"
```

```{bash engine.opts='-l'}
bcftools stats -s - eg/ex.vcf | grep -A 4 "Per-sample counts"
```

```{bash engine.opts='-l'}
bcftools stats -s - eg/ex.vcf | grep -A 4 "Per-Sample Indels"
```

## Add AF tag to a VCF file

The allele frequency tag (AF) is missing from the VCF file generated by `BCFtools call`.

```{bash engine.opts='-l'}
bcftools view -H eg/aln.bt.vcf.gz | head -1
```

The `fill-tags` plugin can add additional tags to a VCF file including the AF tag.

```{bash engine.opts='-l'}
bcftools plugin fill-tags eg/aln.bt.vcf.gz | grep -v "^#" | head -1
```

## Add custom annotations

The `annotate` function can be used to [add additional annotations](https://samtools.github.io/bcftools/howtos/annotate.html) to a VCF file. You can use a VCF file or a tabix-indexed tab-delimited file; if a tab-delimited file is used, an additional header file is required.

To demonstrate, a tab-delimited file with incremental integers (fifth column) is created.

```{bash engine.opts='-l'}
bcftools view -H eg/aln.bt.vcf.gz | head -3 | cut -f1,2,4,5 | perl -nle '$i++; print join("\t", $_, "$i")' > eg/test_anno.tsv
cat eg/test_anno.tsv
```

`eg/test_anno.tsv` needs to be tabix-indexed with the following parameters:

* `-s` specifies the column for sequence names
* `-b` specifies the column for region start/s
* `-e` specifies the column for region end/s

```{bash engine.opts='-l'}
bgzip eg/test_anno.tsv
tabix -s1 -b2 -e2 eg/test_anno.tsv.gz
```

Create header file.

```{bash engine.opts='-l'}
echo -e "##INFO=<ID=INC,Number=1,Type=Integer,Description=\"Increment number\">" > eg/test_anno.hdr
```

Add `INC` annotation.

```{bash engine.opts='-l'}
bcftools annotate -a eg/test_anno.tsv.gz -h eg/test_anno.hdr -c CHROM,POS,REF,ALT,INC eg/aln.bt.vcf.gz | grep -v "^##" | head -4
```

To add the custom annotations to the ID column of the VCF file, specify the ID column instead of the INFO tag with the `-c` parameter. (Note that we no longer require the header file, since the ID column is a standard column in the VCF.)

```{bash engine.opts='-l'}
bcftools annotate -a eg/test_anno.tsv.gz -c CHROM,POS,REF,ALT,ID eg/aln.bt.vcf.gz | grep -v "^##" | head -4
```

In addition, `annotate` can set IDs on-the-fly using information from the VCF columns using the `--set-id` parameter.

```{bash engine.opts='-l'}
bcftools annotate --set-id '%CHROM\_%POS\_%REF\_%FIRST_ALT' eg/aln.bt.vcf.gz | grep -v "^##" | head -4
```

## Check whether the REF sequence is correct

Change the reference sequence of variant at position 151 to G.

```{bash engine.opts='-l'}
zcat eg/aln.bt.vcf.gz | perl -lane 'if($F[1] == 151){ $F[3] = G; print join("\t", @F) } else { print }' > eg/incorrect.vcf
```

Use `bcftools norm` with `-c` to check whether the REF sequence is correct.

```{bash engine.opts='-l'}
bcftools norm -f eg/test_31.fa -c w eg/incorrect.vcf > /dev/null
```

## Sorting

BCFtools has a `sort` function but the usage does not indicate the sorting order. We can create a subsetted and shuffled VCF file and check the order after sorting.

```{bash engine.opts='-l'}
bcftools view -h eg/Pfeiffer.vcf > eg/Pfeiffer_header.vcf
bcftools view eg/Pfeiffer.vcf | perl -nle 'BEGIN { srand(1984) } if (/^#/){ next }; print if rand(1) < 0.01' | shuf > eg/Pfeiffer_shuf.txt

cat eg/Pfeiffer_header.vcf eg/Pfeiffer_shuf.txt > eg/Pfeiffer_shuf.vcf

bcftools view -H eg/Pfeiffer_shuf.vcf | head
```

Sort.

```{bash engine.opts='-l'}
bcftools sort eg/Pfeiffer_shuf.vcf > eg/Pfeiffer_sorted.vcf
```

Check first 10 lines.

```{bash engine.opts='-l'}
bcftools view -H eg/Pfeiffer_sorted.vcf | head
```

Check chromosome order.

```{bash engine.opts='-l'}
bcftools view -H eg/Pfeiffer_sorted.vcf | cut -f1 | uniq
```

BCFtools sorts by lexical order (and there does not seem to be an option for another order).

If you use Linux, we can use the `sort` utility provided by your distribution instead of `bcftools sort`; we just need to separate and join the VCF header back. For this README, the following version of `sort` was used.

```{bash engine.opts='-l'}
sort --version
```

Sort using natural sort.

```{bash engine.opts='-l'}
sort -k1,1V -k2,2n eg/Pfeiffer_shuf.txt > eg/Pfeiffer_ns.txt
cat eg/Pfeiffer_header.vcf eg/Pfeiffer_ns.txt >  eg/Pfeiffer_ns.vcf
bcftools view -H eg/Pfeiffer_ns.vcf | cut -f1 | uniq
```

## Index a VCF file

An index is required for several tasks such as [subsetting variants](#subset-variants-within-a-specific-genomic-region). Using either `bcftools index` or `tabix` seems to be fine, although they generate indexes with different suffixes.

Using `bcftools index` to generate a `csi` index.

```{bash engine.opts='-l'}
bcftools index eg/aln.hc.vcf.gz
bcftools view -H -r 1000000:100-1000 eg/aln.hc.vcf.gz
rm eg/aln.hc.vcf.gz.csi
```

Using `tabix` to generate a `tbi` index.

```{bash engine.opts='-l'}
tabix eg/aln.hc.vcf.gz
bcftools view -H -r 1000000:100-1000 eg/aln.hc.vcf.gz
rm eg/aln.hc.vcf.gz.tbi
```

## Rename sample names

Use `reheader` to rename samples in a VCF using a text file containing the new sample names (one per line). Note that the default output format for `reheader` is BCF (with no option for another format).

Old sample names.

```{bash engine.opts='-l'}
bcftools query -l eg/1001genomes_snp-short-indel_only_ACGTN_5000.vcf.gz | head -5
```

Rename VCF file using new sample names in `eg/sample_name.txt`.

```{bash engine.opts='-l'}
num_sample=$(bcftools query -l eg/1001genomes_snp-short-indel_only_ACGTN_5000.vcf.gz | wc -l)
for (( i = 0 ; i < ${num_sample} ; i++ )); do echo ${i}; done > eg/sample_name.txt

bcftools reheader -s eg/sample_name.txt eg/1001genomes_snp-short-indel_only_ACGTN_5000.vcf.gz | bcftools view -h - | tail -1
```

## Remove header info

The `reheader` command can also be used to modify the header information, such as removing specific info. For example if we want to remove `GATKCommandLine`, we would first create a new header file and use it to generate a new header.

`GATKCommandLine` exists in the header before running `reheader`.

```{bash engine.opts='-l'}
bcftools view -h eg/aln.hc.vcf.gz | grep GATKCommandLine
```

`GATKCommandLine` has been removed with `reheader`.

```{bash engine.opts='-l'}
bcftools view -h eg/aln.hc.vcf.gz | grep -v GATKCommandLine > eg/new_header.txt
bcftools reheader -h eg/new_header.txt eg/aln.hc.vcf.gz | bcftools view -h -
```

## Concatenate VCF files

Concatenating is similar to joining VCF files in a row-wise manner. This only works for VCF files with the same samples that are also organised in the same order.

If we concatenate `eg/aln.bt.vcf.gz` and `eg/aln.hc.vcf.gz` we would expect that the total number of variants is the sum of both.

```{bash engine.opts='-l'}
expr $(bcftools view -H eg/aln.bt.vcf.gz | wc -l) + $(bcftools view -H eg/aln.hc.vcf.gz | wc -l)
```

Perform the concatenation with default settings and count the number of variants.

```{bash engine.opts='-l'}
bcftools concat eg/aln.bt.vcf.gz eg/aln.hc.vcf.gz  | bcftools view -H - | wc -l
```

Removing duplicates requires indexed VCF files; the `-f` parameter is used with `bcftools index` to overwrite an index if it exists.

```{bash engine.opts='-l'}
bcftools index -f eg/aln.bt.vcf.gz
bcftools index -f eg/aln.hc.vcf.gz
```

To remove duplicates (`-D`) we need to allow overlaps (`-a`).

```{bash engine.opts='-l'}
bcftools concat -a -D eg/aln.bt.vcf.gz eg/aln.hc.vcf.gz  | bcftools view -H - | wc -l
```

## Merging VCF files

Use `bcftools merge` to merge VCF files. [In this workflow](https://github.com/davetang/sars_cov_2/blob/master/.github/workflows/omicron_variants.yml), `bcftools merge` was used to create a multi-sample VCF file from individual VCF files.

```bash
bcftools merge -o PRJNA784038_illumina.vcf -O v SRR*.vcf.gz
```

## Decomposing and normalising variants

Decomposing can refer to the splitting of multi-allelic variants; we can use `bcftools norm -m` for this.

```{bash engine.opts='-l'}
bcftools view -H eg/PRJNA784038_illumina.vcf.gz | head -2 | cut -f1-5
```

Splitting.

```{bash engine.opts='-l'}
bcftools norm -m- eg/PRJNA784038_illumina.vcf.gz | grep -v "^#" | head -3 | cut -f1-5
```

Decomposing can also refer to converting MNVs into consecutive SNVs; this is achieved with `bcftools norm -a`.

```{bash engine.opts='-l'}
bcftools view -H eg/ex.vcf  | tail -1
```

Convert MNV into SNVs.

```{bash engine.opts='-l'}
bcftools norm -a eg/ex.vcf | tail -2
```

Finally, there is also left-aligning, which will be clear by viewing an example.

```{bash engine.opts='-l'}
bcftools view -H eg/aln.bt.vcf.gz | head -2
```

Left-align.

```{bash engine.opts='-l'}
bcftools norm -f eg/test_31.fa eg/aln.bt.vcf.gz | grep -v "^#" | head -2
```

Decomposing and normalising variants are very important steps when comparing VCF files.

## Comparing VCF files

BCFtools has an intersect tool (`bcftools isec`) that is useful for comparing VCF files; the tool requires an index.

```{bash engine.opts='-l'}
tabix -f eg/aln.bt.vcf.gz
tabix -f eg/aln.hc.vcf.gz

bcftools isec eg/aln.bt.vcf.gz eg/aln.hc.vcf.gz -p comp
```

Four VCF files are produced and `README.txt` explains the contents of each VCF file.

Records private to `eg/aln.bt.vcf.gz`.

```{bash engine.opts='-l'}
bcftools view -H comp/0000.vcf | wc -l
```

Records private to `eg/aln.hc.vcf.gz`.

```{bash engine.opts='-l'}
bcftools view -H comp/0001.vcf | wc -l
```

Number of overlap.

```{bash engine.opts='-l'}
bcftools view -H comp/0002.vcf | wc -l
```

[SnpSift](http://snpeff.sourceforge.net/SnpSift.html) is also useful for comparing VCF files for checking their [concordance](http://davetang.org/muse/2015/08/26/vcf-concordance/).

## Visualisation

The [vcfR](https://cran.r-project.org/web/packages/vcfR/index.html) package produces some nice plots.

```r
install.packages("vcfR")
library(vcfR)

my_vcf <- read.vcfR("~/github/learning_vcf_file/eg/Pfeiffer.vcf", verbose = FALSE)
chrom <- create.chromR(name="Pfeiffer variants", vcf=my_vcf)
chrom <- proc.chromR(chrom, verbose=TRUE)

plot(chrom)
```

![Plot chromR object](img/chrom_plot.png)

```r
chromoqc(chrom, xlim=c(860000, 900000))
```

![QC plot](img/chromoqc.png)

## Useful links

* A very useful thread on SEQanswers on learning about the VCF format: <http://seqanswers.com/forums/showthread.php?t=9345>
* Useful tutorial on VCFs files from the 1000 Genomes Project Page: <http://www.1000genomes.org/node/101>
* The author of ANNOVAR writes about VCF files: <http://annovar.openbioinformatics.org/en/latest/articles/VCF/>
* [Encoding Structural Variants in VCF](http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/VCF%20(Variant%20Call%20Format)%20version%204.0/encoding-structural-variants) version 4.0
* [BCFtools cheat sheet](https://gist.github.com/elowy01/93922762e131d7abd3c7e8e166a74a0b)

## Stargazers over time

[![Stargazers over time](https://starchart.cc/davetang/learning_vcf_file.svg)](https://starchart.cc/davetang/learning_vcf_file)