VCF to aGDS (vcf2agds) Overview

Description of VCF to aGDS conversion and annotation on the UK Biobank RAP

Due to the number of VCFs and their file sizes on the DNAnexus platform, it may be necessary to reduce the amount of data contained within VCFs if wanting to merge all VCFs associated with a chromosome on a workstation for subsequent format conversion and annotation. The following suite of applets are provided here that may help facilitate this process. The examples provided below focus on the UK Biobank 500k WGS data.

Workflow Overview

Step 1: Trimming down data in the VCFs

See vcf_trimmer for an applet that removes fields and performs required filtering through bcftools. This applet takes in a list of files as an input to process. This applet could be used across multiple jobs submitted on the DNAnexus platform which would require unique lists of VCFs to be split. Note, left-alignment and normalisation not presently set.

Example: Processing WGS (500k) VCFs for chromosome 17 over 500 jobs:

Generating the input file lists:

If you wanted to process the UK Biobank 500k WGS data release for chromosome 17 then vcf_trimmer would require a list of VCFs associated with that chromosome. Furthermore, you may also want to submit 500 jobs whereby the list of VCFs associated with chromosome 17 is split into 500 unique and equally sized VCF lists. Running the get_vcf_file_list_by_chr_and_split.sh bash script in this repo will produce all the relevant input files for chromosome 17 (chr17_vcf_list_[1:500]) as well as for all other chromosomes for a total of 151,561 VCFs currently available.

bash get_vcf_file_list_by_chr_and_split.sh

This script has been provided as a starting point that can be amended to suit your requirements. Note, the file lists generated by this bash script will need to be subsequently uploaded to a project folder on the DNAnexus platform.

Submitting the jobs:

Once you have a set of 500 files listing the VCFs uploaded to a project folder on the RAP, you can run the vcf_trimmer with each of the files generated above. For example, if we wanted to run vcf_trimmer on chromosome 17 to remove FORMAT fields except GT and keep only variants with AAscore >= 0.5, then the command would look something like this:

for i in {1..500}
do
  dx run /path/to/vcf_trimmer \
    -ivcf_file_list=/path/to/chr17_vcf_list_${i} \
    -ifile_label=trimmed \
    -ioutput_dir=/path/to/output/dir \
    -iqc_thresholds="INFO/AAScore>=0.5" \
    -ifields_to_remove="FORMAT/FT,FORMAT/AD,FORMAT/MD,FORMAT/DP,FORMAT/RA,FORMAT/PP,FORMAT/GQ,FORMAT/PL" \
    -y
done

Step 2: Merging VCFs

See vcf_merger. Another applet wrapper for bcftools. This applet also takes in a list of VCF files. Continuing the example above, to obtain a list of files in block order of a directory on the RAP for chromosome 17 after running vcf_trimmer, type:

dx ls "/path/to/output/dir/ukb23374_c17_b*_v1_trimmed.vcf.gz" | sort -t"b" -k3.1 -n | awk -v d="/path/to/output/dir" '{print d $0}' > chr17_vcfs_to_merge

The output file should be subsequently uploaded to the RAP for processing with vcf_merger.

Step 3: Converting VCF to SeqArray GDS for subsequent annotation in the FAVOR database

See vcf2gds. This applet comes with an R library that will be unpacked during runtime and used for data conversion. The SeqArray Genomic Data Structure (GDS) format is designed as a storage-efficient high-performance data format for WGS variant calls, utilizing a hierarchical structure to support efficient data access.

Step 4: Annotating GDS to aGDS for subsequent association analysis in the STAARpipeline

See favorannotator. This applet pulls the docker image that will be executed during runtime and used for annotation. The annotated GDS (aGDS) format further extends the capabilities of SeqArray GDS by incorporating multi-faceted variant annotations within an all-in-one file to facilitate a range of functionally-informed downstream analyses.