VCF import

VariantGrid reads VCF files and inserts the records in a SQL database. A VCF file can have 2M rows with 100+ samples (400M reference/alt genotypes) and some labs have thousands of VCFs a year.

The difficult part is doing this at speed while ensuring each variant has a single record, so that you can see every sample that ever contained that variant.

The image below shows an exome import. Colors represent different steps, with lines from start/end time of each task, with a new vertical row started if there is any overlap (to show parallel processing).

Everything is not perfectly parallel, we need to finish some stages before others can start.

For instance you need to normalize variants and insert previously unknown variants before you can annotate them, or insert genotypes against the Variant IDs. Here are descriptions of some steps:

Create Data from VCF Header

• Stores header text
• Creates VCF and sample database objects
• Find what format/info fields map to AD/DP/GQ/PL (different for some VCFs)
• Converts VCF filters into single character codes. Writes filter descriptions and mappings to single character to DB
• Stores importer version used, eg: PythonKnownVariantsImporter (v.4). Git: 7cbb9ae2876168b24924e5b97a785e4925a9d2f0. Uses cyvcf2 (v.0.8.9)

Preprocess VCF

Implemented in upload/vcf/vcf_preprocess.py. Everything is assembled into a single shell pipeline (each stage's stdout piped to the next stage's stdin) so the file is only read once:

• vcf_clean_and_filter (management command) - handle Contigs, replace the header with a cleaned one, clean the file, track what's discarded
• bcftools norm - --multiallelics=- splits rows with multiple alts into multiple ref/alt rows, and normalises so variants from multiple files "match up". --check-ref=s sets the ref from the fasta (eg replacing "N"), and --old-rec-tag records the pre-normalisation record in an INFO field
• bcftools view --no-header - strip the header (it gets re-attached to each split file below)
• vcf_clean_alts (management command) - now that multi-allelics have been split, filter out bad alts while keeping the rest of the record
• GNU split - split the VCF into chunks (settings.VCF_IMPORT_FILE_SPLIT_ROWS) so they can be processed in parallel. The --filter re-attaches the cleaned header and bgzips each chunk

We deliberately don't deduplicate here: bcftools norm --rm-dup drops --old-rec-tag, so we'd lose the normalisation provenance (bcftools issue #2225).

Historical note: this used to be vt (decompose / normalize / uniq). NORM_TOOL_VT survives as an enum value in upload/models/models.py so old imports still report the tool they actually used.

Any variants lost or modified in the above process are tracked with the following database records.

• VCFSkippedContigs - contig name/num records skipped by vcf_clean_and_filter
• ModifiedImportedVariant - populated from the bcftools --old-rec-tag INFO field (ModifiedImportedVariant.BCFTOOLS_OLD_VARIANT_TAG), whose format is CHR|POS|REF|ALT|USED_ALT_IDX. This records the original multi-allelic / un-normalised representation of each row.

There is not yet any way to search or see these via the GUI (but they are being stored)

Create Unknown Loci and Variants

Variant primary keys are looked up in Postgres via snpdb.variant_pk_lookup.VariantPKLookup. It batches lookups by contig + position and matches on a hash built in SQL (position_ref_id_alt_id_svlen), so a batch of coordinates resolves to existing Variant ids in a couple of queries.

If a variant/locus is unknown, we write it to a file and launch jobs to insert them. Unknown insert jobs are run one at a time on the "variant_id_single_worker" queue to avoid race conditions (inserting duplicate records).

• Reads the unknown variants csv and re-checks whether the variant now exists (the unknown variant generation tasks run in parallel, and it may have been inserted by another VCF since the file was written)
• Writes remaining unknowns to CSV, inserts into Postgres (Locus/Variants assigned an auto integer primary key)

Historical note: the original implementation (RedisVariantPKLookup) kept variant hash → PK in Redis, which cost a few gigs of RAM. It was removed in August 2021 — Redis is still used for the Django cache and as the Celery result backend, but no longer holds variant state.

CheckStartAnnotationTask

Waits for the tasks in 'Unknown Variants' stage to finish.

• If at least one "Create Unknown Loci and Variants" task was run during this import, run the annotation scheduler, which dumps any unannotated variants to disk and runs the annotation pipeline.

Schedule Parallel VCF Processing tasks

Waits for all 'Unknown Variants' stage tasks to finish, so that we know every variant in the VCF file has been inserted into the database.

Launch a task for each sub-vcf file made from GNU split.

Process VCF File

• Jobs run in parallel on sub-vcfs, uses CyVCF2
• Batch queries Postgres (via VariantPKLookup) to map from variant coordinate -> integer primary key
• Writes out CSVs of CohortVariantZygosityCounts, launch tasks to insert them
• Writes out VCF filter records

CohortGenotype SQL COPY

Inserts records generated by previous step. These have information from multiple samples packed together in a single database row, so we can perform fast multi-sample zygosity queries.

Gemini uses this technique, but we store our data as Postgres arrays instead of binary blobs so we can do everything in SQL.

VCFLocusFilter SQL COPY

• Inserts VCF filters generated by "Process VCF File"
• As these are 1 per VCF row, we store these per Locus (not Variant)
• They are single character codes, generated in "Create Data from VCF Header", mappings which look like:

{'filter_code': '$',
 'description': 'Truth sensitivity tranche level for SNP model at VQS Lod: -1.4325 <= x < 2.3577',
 'filter_id': 'VQSRTrancheSNP99.00to99.90'}

We only store non-PASS filters. The final records are quite small - an integer (locus ID) and 1 character per filter. We can quickly perform filters by using Postgres SQL Regex.

GenotypeVCFCheckAnnotationTask

Check whether every variant in the VCF has been annotated. This runs after all above steps are finished (data insertion stage).

• Creates a database record UploadedVCFPendingAnnotation for the Uploaded VCF

Runs attempt_schedule_annotation_stage_steps which:

• Finds the lowest unannotated variant (for latest annotation version) in the database
• Compare lowest unannotated variant id to the maximum variant ID inserted for the VCF (Each "Process VCF File" task updates UploadedVCF with highest record it used)
• If all variants are annotated, close this step, otherwise leave it open.

Every time the annotation scheduler finishes, it looks for unfinished UploadedVCFPendingAnnotation objects then calls attempt_schedule_annotation_stage_steps()

UpdateVariantZygosityCountsTask

This runs after all data is inserted. Keeps track of the total amount of het/hom samples that share a variant in the database. Updates from the CohortVariantZygosityCount het/hom counts.

Works like Reference Counting, counts from CohortVariantZygosityCount het/hom counts are subtracted when the VCF/samples are deleted.

SampleLocusCountsTask

This runs after all data is inserted. Counts how many alts per locus a sample had, eg:

Locus Count	Count
1 (HOM)	1102186
2 (HET)	110271
3	4324
4	953
5	245

This is useful for looking for germline sample contamination.

Calculate Sample Stats

Collects information about variants in a sample. Runs after annotation is completed, as it counts eg ClinVar pathological variants, genes that have an OMIM phenotype, Snpeff impacts etc for both all variants and those with filter=PASS.

These can be viewed on the Sample and VCF pages, and are also used to lookup the node counts in the analysis rather than calculate them each time.

ImportGenotypeVCFSuccessTask

Runs after all other steps are completed. Sets VCF and samples to "SUCCESS" and sends a message/email saying the data is ready.