Genotype VCF File Quality Control#
This implements some recommendations from UK Biobank on sequence data quality control.
Overview#
Genotype problems are cheap to fix here and expensive to find later. Chromosome naming
differs between references (1 versus chr1), the same variant can be represented several
ways, and low-depth or low-quality calls look like real genotypes to every tool downstream.
This step normalises naming, annotates variants against dbSNP so identifiers are consistent
across cohorts, and applies the depth and quality filters, so everything after it works from
one convention.
A major challenge in biomedical research is the quality control (QC) of sequencing data. False positive variant calls can hinder the ability to detect disease associated variants or introduce spurious associations, therefore the need for a rigorous QC. This module focuses on QC after the variant calling stage and requires project Variant Calling Format (pVCF) as input files. Default thresholds are defined for genotype and variant-level hard filtering based on recommendations from the UK Biobank team and a thorough review of the literature [cf. Carson et al. BMC Bioinformatics (2014),cf. Lek et al. Nature (2016),cf. Szustakowski et al. Nature Genetics (2021)]. Bcftools performs the QC steps. Multi-allelic sites are first split into bi-allelic records. An optional workflow keeps only bi-allelic sites in the data. Variants are then annotated based on dbSNP data. Genotypes are kept if they have a Genotype Depth (DP) >= 10 and a Genotype Quality (GQ) >= 20. Variants are included if at least one sample has an allelic balance (AB) >= 0.15 for Single Nucleotide Variants (SNVs) and AB>=0.2 for indels, variant missigness is below 20% and, when an HWE cutoff is supplied, the Hardy-Weinberg Equilibrium p-value exceeds it. Allele balance is calculated for heterozygotes as the number of bases supporting the least-represented allele over the total number of base observations. Output summary statistics, such as transistion/transversion ratios (TS/TV ratio) are calculated to determine the effectiveness of QC.
When to run it. First, on the raw VCF, before formatting, PCA or association testing. A meta-analysis across cohorts is where naming and identifier mismatches surface, and by then it is far harder to trace.
Input#
--genoFile input/genotype/protocol_example.genotype.chr*.vcf.gz(+.tbi): the per-chromosome pVCFs to QC, from upstream genotype calling. A manifest may be given instead –input/genotype/protocol_example.genotype.vcf_list.txt, two columns, chromosome number and VCF path – which lets the chromosomes run in parallel. This is the only required data input; everything else is reference material or a threshold.--dbsnp-variants input/reference_data/00-All.add_chr.variants.gz(+.tbi) (required forqc): chr-renamed dbSNP table,CHROM/POS/ID/REF/ALT, used to split variants into known and novel.--reference-genome input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta(+.fai) (required forqc): FASTA for left-normalization and the REF check.--cwd(defaultoutput): work directory for the analysis.--modular-script-dir(defaultcode/script): directory holding the modular wrapper scripts.--remove-samples/--keep-samples(both default.): optional sample lists applied before QC.--skip-vcf-header-filtering(default off): skip the header-based contig filtering.--gt-only-vcf-qc(default off): treat the input as genotype-only, so the depth and quality filters below are not applied. Required for the toyprotocol_example.genotype.*VCFs, whose FORMAT carries onlyGT– with noDP,GQorADthose filters have nothing to read.
Variant preprocessing (qc_1):
--bi-allelic(default off): keep only bi-allelic sites.--snp-only(default off): drop indels.
Variant- and genotype-level filters (qc_2):
--geno-filter(default0.2): maximum per-variant missingness.--DP-snp(default10) and--DP-indel(default10): minimum read depth per genotype.--GQ(default20): minimum genotype quality.--AB-snp(default0.15) and--AB-indel(default0.2): allele-balance thresholds for heterozygous calls.--hwe-filter(default0.0): HWE p-value cutoff. The parameter help calls0.0“no HWE filter”, but the script appliesHWE > <cutoff>literally, so0.0still drops any site whose HWE p-value is exactly 0 or whoseHWEtag is missing.
Cluster resources: --numThreads (default 1), --job-size (default 1), --walltime (default
24h), --mem (default 60G). A whole-genome VCF needs the full 60G; qc_1 and qc_2 process
roughly 14G of data per hour, so scale --walltime with input size.
On choosing the filter thresholds#
For WES data the UK Biobank recommendation is DP>10 for both SNPs and indels. For WGS the depth
thresholds can be relaxed, or GQ relied on alone – --DP-snp 1 --DP-indel 1 effectively disables
the depth filter.
The allele-balance filters keep a site if at least one sample passes the threshold. Allele balance is computed for heterozygotes as the number of bases supporting the least-represented allele over the total number of base observations.
Filtering is done with bcftools; this
cheatsheet is a useful
reference.
Output#
<geno>.leftnorm.vcf.gz(+.tbi) – fromqc_1: left-normalized, variant-ID-assigned, dbSNP-annotated VCF.<geno>.leftnorm.bcftools_qc.vcf.gz(+.tbi) – fromqc_2: the same VCF after the missingness, depth, quality and allele-balance filters. This is the file downstream modules read.<geno>.leftnorm.bcftools_qc.known_variant.snipsift_tstvand<geno>.leftnorm.bcftools_qc.novel_variant.snipsift_tstv– fromqc_3: TS/TV and Hom/Het summaries, split by whether the variant carries a dbSNP rsID.<geno>.leftnorm.bcftools_qc.known_variant_sumstatsand<geno>.leftnorm.bcftools_qc.novel_variant_sumstats– fromqc_3: the fullbcftools statsreport per split, with theSNrecord counts and the per-samplePSCblock.a
.stdout/.stderrlog and the generated shell script beside each target.
All paths are relative to --cwd, and <geno> is the --genoFile stem, e.g.
protocol_example.genotype.chr22.
With --gt-only-vcf-qc the intermediate and final names carry an extra .gt_only segment, e.g.
<geno>.leftnorm.gt_only.bcftools_qc.vcf.gz.
No example output is committed for the QC stages. On the toy chr22 data every
bcftools_qc.*_sumstats file reports number of records: 0, with or without --gt-only-vcf-qc,
while the input to that step holds 18,455 variants of which bcftools view -c1 retains 4,259. The
qc_2 filters are discarding everything and the step still exits successfully, so the committed
example output is empty and no dump can be shown here yet.
bcftools stats command provides useful summary statistics including TS/TV ratio, which is routinely used as a quality measure of variant calls. With dbSNP based annotation of novel and known variants, bcftools can compute TS/TV for novel and known variants at variant level, and at sample level. It should be noted that variant level TS/TV does not take sample genotype into consideration – it simply counts the TS and TV event for observed SNPs in the data. Other tools, such as snpsift, implements variant level TS/TV by counting TS and TV events in sample genotypes and compute the ratio after summing up TS and TV across all samples. See here some discussions on this issue. TS/TV is reported before and after QC, but users should be aware of the difference when interpreting the results.
Minimal Working Example#
The steps below run on the toy protocol_example genotype data. Steps 1 and 2 are optional helpers; Step 3 is the main QC workflow.
Step 1. Rename Chromosomes (optional)#
Add a chr prefix to chromosome names so they match the reference fasta. Timing: ~40 min on full data.
The toy genotypes are already chr-prefixed, so to demonstrate rename_chrs the protocol ships a small derived input protocol_example.genotype.rawchr.chr22.vcf.gz whose chromosome is named 22 (no chr). Running rename_chrs adds the prefix back, writing protocol_example.genotype.rawchr.add_chr.vcf.gz next to the input.
Timing: TBD (on toy dataset)
sos run pipeline/VCF_QC.ipynb rename_chrs \
--genoFile input/genotype/protocol_example.genotype.rawchr.chr22.vcf.gz \
--cwd output/vcf_qc
Extracts CHROM/POS/ID/REF/ALT from a VCF into the compact .variants.gz dbSNP table used by the QC step. It is run here on a toy genotype VCF to produce protocol_example.genotype.chr22.variants.gz.
Annotates against the dbSNP database in VCF format.
Timing: TBD (on toy dataset)
sos run pipeline/VCF_QC.ipynb dbsnp_annotate \
--genoFile input/genotype/protocol_example.genotype.chr22.vcf.gz \
--cwd output/vcf_qc
Step 2. Merge per-chromosome files (handled elsewhere)#
Converting the QC-d VCFs to PLINK format and merging them is done by genotype_formatting.ipynb
(vcf_to_plink, then merge_plink), not by this module. See that notebook for the commands and its
own --name option.
Step 3. Quality Control#
The main QC workflow: split multi-allelic sites, left-normalize indels, assign variant IDs, and annotate known/novel variants against dbSNP.
Run on a single chromosome (chr22) below, or in parallel over all chromosomes.
Timing: TBD (on toy dataset)
sos run pipeline/VCF_QC.ipynb qc \
--genoFile input/genotype/protocol_example.genotype.chr22.vcf.gz \
--dbsnp-variants input/reference_data/00-All.add_chr.variants.gz \
--reference-genome input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
--cwd output/vcf_qc \
--skip_vcf_header_filtering True
To run in parallel for all genotype data listed in protocol_example.genotype.vcf_list.txt:
Note on --skip_vcf_header_filtering True: the toy protocol_example.genotype.* VCFs contain only the GT (genotype) FORMAT field. The default depth/quality filtering relies on the DP FORMAT tag, which is absent here, so --skip_vcf_header_filtering True runs left-normalization, variant-ID assignment and dbSNP annotation while skipping the DP-based filtering — the appropriate path for genotype-only VCFs. To see the default path (which does apply the DP/GQ/AD filters), use the protocol_example.genotype.withfmt.chr22.vcf.gz file in Step 3b below, which carries those tags.
Timing: TBD (on toy dataset)
sos run pipeline/VCF_QC.ipynb qc \
--genoFile input/genotype/protocol_example.genotype.vcf_list.txt \
--dbsnp-variants input/reference_data/00-All.add_chr.variants.gz \
--reference-genome input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
--cwd output/vcf_qc \
--skip_vcf_header_filtering True \
-j 2
Command Interface#
sos run pipeline/VCF_QC.ipynb -h
usage: sos run pipeline/VCF_QC.ipynb
[workflow_name | -t targets] [options] [workflow_options]
workflow_name: Single or combined workflows defined in this script
targets: One or more targets to generate
options: Single-hyphen sos parameters (see "sos run -h" for details)
workflow_options: Double-hyphen workflow-specific parameters
Workflows:
rename_chrs
dbsnp_annotate
qc
Global Workflow Options:
--modular-script-dir code/script (as path)
--genoFile paths
input can either be 1 vcf genoFile, or a list of vcf
genoFile.
--remove-samples . (as path)
The path to the file that contains the list of samples
to remove (format FID, IID)
--keep-samples . (as path)
The path to the file that contains the list of samples
to keep (format FID, IID)
--cwd output (as path)
Workdir
--numThreads 1 (as int)
Number of threads
--job-size 1 (as int)
For cluster jobs, number commands to run per job
--walltime 24h
Walltime
--mem 60G
Usually a whole genome VCF.gz file has the size of
200+GB, after testing, a minimum of 60GB of mem is
requried.
--[no-]gt-only-vcf-qc (default to False)
Whether to run GT-only VCF QC for inputs that do not
carry per-sample DP/GQ/AD FORMAT fields.
--[no-]skip-vcf-header-filtering (default to False)
Deprecated alias for GT-only VCF QC.
Sections
rename_chrs:
dbsnp_annotate:
qc_1: Handel multi-allelic sites, left normalization of indels
and add variant ID
Workflow Options:
--dbsnp-variants VAL (as path, required)
Path to dbSNP variants generated previously
--reference-genome VAL (as path, required)
Path to fasta file for HG reference genome, eg
GRCh38_full_analysis_set_plus_decoy_hla.fa
--[no-]bi-allelic (default to False)
--[no-]snp-only (default to False)
qc_2: genotype QC
Workflow Options:
--geno-filter 0.2 (as float)
Maximum missingess per-variant, default to 0.2
--DP-snp 10 (as int)
Sample level QC - read depth (DP) to filter out SNPs
below this value Default to 10, with WES data in mind
But for WGS, setting it to 2 may be fine considering the
WGS may have low DP but the GQ filter should be good
enough
--GQ 20 (as int)
Sample level QC - genotype quality (GQ) of specific
sample. This measure tells you how confident we are that
the genotype we assigned to a particular sample is
correct
--DP-indel 10 (as int)
Sample level QC - read depth (DP) to filter out indels
below this value
--AB-snp 0.15 (as float)
Allele balance for snps
--AB-indel 0.2 (as float)
Allele balance for indels
--hwe-filter 0.0 (as float)
HWE filter, default to 0.0 which means no HWE filter is
applied
qc_3:
Workflow implementation#
[global]
parameter: modular_script_dir = path('code/script') # override with --modular-script-dir
# input can either be 1 vcf genoFile, or a list of vcf genoFile.
parameter: genoFile = paths
# The path to the file that contains the list of samples to remove (format FID, IID)
parameter: remove_samples = path('.')
# The path to the file that contains the list of samples to keep (format FID, IID)
parameter: keep_samples = path('.')
# Workdir
parameter: cwd = path("output")
# Number of threads
parameter: numThreads = 1
# For cluster jobs, number commands to run per job
parameter: job_size = 1
# Walltime
parameter: walltime = '24h'
# Usually a whole genome VCF.gz file has the size of 200+GB, after testing, a minimum of 60GB of mem is requried.
parameter: mem = '60G'
import re
# Whether to run GT-only VCF QC for inputs that do not carry per-sample DP/GQ/AD FORMAT fields.
parameter: gt_only_vcf_qc = False
# Deprecated alias for GT-only VCF QC.
parameter: skip_vcf_header_filtering = False
gt_only_vcf_qc = bool(gt_only_vcf_qc or skip_vcf_header_filtering)
vcf_qc_mode_suffix = '.gt_only' if gt_only_vcf_qc else ''
# use this function to edit memory string for PLINK input
from sos.utils import expand_size
cwd = path(f"{cwd:a}")
import os
def get_genotype_file(geno_file_paths):
#
def valid_geno_file(x):
suffixes = path(x).suffixes
if suffixes[-1] == '.bed':
return True
elif suffixes[-1] == '.vcf':
return True
elif len(suffixes)>1 and ''.join(suffixes[-2:]) == ".vcf.gz" :
return True
return False
#
def complete_geno_path(x, geno_file):
if not valid_geno_file(x):
raise ValueError(f"Genotype file {x} should be VCF (end with .vcf.gz) or PLINK bed file (end with .bed)")
if not os.path.isfile(x):
# relative path
if not os.path.isfile(f'{geno_file:ad}/' + x):
raise ValueError(f"Cannot find genotype file {x}")
else:
x = f'{geno_file:ad}/' + x
return x
#
def format_chrom(chrom):
if chrom.startswith('chr'):
chrom = chrom[3:]
return chrom
# Inputs are either VCF or bed, or a vector of them
if len(geno_file_paths) > 1:
if all([valid_geno_file(x) for x in geno_file_paths]):
return paths(geno_file_paths)
else:
raise ValueError(f"Invalid input {geno_file_paths}")
# Input is one genotype file or text list of genotype files
geno_file = geno_file_paths[0]
if valid_geno_file(geno_file):
return paths(geno_file)
else:
units = [x.strip().split() for x in open(geno_file).readlines() if x.strip() and not x.strip().startswith('#')]
if all([len(x) == 1 for x in units]):
return paths([complete_geno_path(x[0], geno_file) for x in units])
elif all([len(x) == 2 for x in units]):
genos = dict([(format_chrom(x[0]), path(complete_geno_path(x[1], geno_file))) for x in units])
else:
raise ValueError(f"{geno_file} should contain one column of file names, or two columns of chrom number and corresponding file name")
return genos
genoFile = get_genotype_file(genoFile)
Annotate known and novel variants#
You can download the known variant reference from this link.
For a detailed explanation of the procedure and its rationale, please refer to this post.
[rename_chrs: provides = '{genoFile:nn}.add_chr.vcf.gz']
# This file can be downloaded from https://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606_b150_GRCh38p7/VCF/00-All.vcf.gz.
input: genoFile
output: f'{_input:nn}.add_chr.vcf.gz'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "${ }", stderr = f'{_output:nn}.stderr', stdout = f'{_output:nn}.stdout'
bash ${modular_script_dir}/data_preprocessing/genotype/VCF_QC.sh rename_chrs \
--genoFile "${_input}" \
--output "${_output}"
[dbsnp_annotate]
input: genoFile
output: f"{cwd}/{_input:bnn}.variants.gz"
task: trunk_workers = 1, trunk_size=job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
bash ${modular_script_dir}/data_preprocessing/genotype/VCF_QC.sh dbsnp_annotate \
--genoFile "${_input}" \
--output "${_output}"
Genotype QC#
This step handles multi-allelic sites and annotate variants to known and novel. An RS ID is added to variants in dbSNP. Variants without rsID are considered novel variants. For every hour it can produce ~14Gb of data, please set the –walltime parameter according to the size of your input files.
# Handel multi-allelic sites, left normalization of indels and add variant ID
[qc_1 (variant preprocessing)]
# Path to dbSNP variants generated previously
parameter: dbsnp_variants = path
# Path to fasta file for HG reference genome, eg GRCh38_full_analysis_set_plus_decoy_hla.fa
parameter: reference_genome = path
parameter: bi_allelic = False
parameter: snp_only = False
input: genoFile, group_by = 1
output: f'{cwd}/{_input:bnn}.{"leftnorm" if not bi_allelic else "biallelic"}{".snp" if snp_only else ""}{vcf_qc_mode_suffix}.vcf.gz'
task: trunk_workers = 1, trunk_size=job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
bash ${modular_script_dir}/data_preprocessing/genotype/VCF_QC.sh qc \
--cwd "${cwd}" \
--genoFile "${_input}" \
--output "${_output}" \
--dbsnp-variants "${dbsnp_variants}" \
--reference-genome "${reference_genome}" \
--bi-allelic "${bi_allelic}" \
--snp-only "${snp_only}" \
--gt-only-vcf-qc "${gt_only_vcf_qc}" \
--numThreads ${numThreads}
This step filter variants based on FILTER PASS, DP and QC, fraction of missing genotypes (all samples), and on HWE, for snps and indels. It will also remove monomorphic sites – using bcftools view -c1.
# genotype QC
[qc_2 (variant level QC)]
# Maximum missingess per-variant, default to 0.2
parameter: geno_filter = 0.2
# Sample level QC - read depth (DP) to filter out SNPs below this value
# Default to 10, with WES data in mind
# But for WGS, setting it to 2 may be fine considering the WGS may have low DP but the GQ filter should be good enough
parameter: DP_snp = 10
# Sample level QC - genotype quality (GQ) of specific sample. This measure tells you how confident we are that the genotype we assigned to a particular sample is correct
parameter: GQ = 20
# Sample level QC - read depth (DP) to filter out indels below this value
parameter: DP_indel = 10
# Allele balance for snps
parameter: AB_snp = 0.15
# Allele balance for indels
parameter: AB_indel = 0.2
# HWE filter, default to 0.0 which means no HWE filter is applied
parameter: hwe_filter = 0.0
output: f"{_input:nn}.bcftools_qc.vcf.gz"
task: trunk_workers = 1, trunk_size=job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand= "${ }", stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout'
bash ${modular_script_dir}/data_preprocessing/genotype/VCF_QC.sh qc_2 \
--cwd "${cwd}" \
--genoFile "${_input}" \
--output "${_output}" \
--gt-only-vcf-qc "${gt_only_vcf_qc}" \
--geno-filter ${geno_filter} \
--DP-snp ${DP_snp} \
--GQ ${GQ} \
--DP-indel ${DP_indel} \
--AB-snp ${AB_snp} \
--AB-indel ${AB_indel} \
--hwe-filter ${hwe_filter} \
--numThreads ${numThreads}
[qc_3 (genotype data summary statistics)]
input: output_from('qc_1'), output_from('qc_2'), group_by = 1
output: f"{cwd}/{_input:bnn}.novel_variant_sumstats",
f"{cwd}/{_input:bnn}.known_variant_sumstats",
f"{cwd}/{_input:bnn}.novel_variant.snipsift_tstv",
f"{cwd}/{_input:bnn}.known_variant.snipsift_tstv"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand= "${ }", stderr = f'{_output[0]:n}.stderr', stdout = f'{_output[0]:n}.stdout'
bash ${modular_script_dir}/data_preprocessing/genotype/VCF_QC.sh qc_3 \
--cwd "${cwd}" \
--genoFile "${_input}" \
--novel-sumstats "${_output[0]}" \
--known-sumstats "${_output[1]}" \
--novel-tstv "${_output[2]}" \
--known-tstv "${_output[3]}" \
--numThreads ${numThreads}