Genotype VCF File Quality Control - FunGen-xQTL Consortium

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 -- tests/fixtures/vcf_qc/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 for qc): 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 for qc): FASTA for left-normalization and the REF check.
--cwd (default output): work directory for the analysis.
--modular-script-dir (default code/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 toy protocol_example.genotype.* VCFs, whose FORMAT carries only GT -- with no DP, GQ or AD those 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 (default 0.2): maximum per-variant missingness.
--DP-snp (default 10) and --DP-indel (default 10): minimum read depth per genotype.
--GQ (default 20): minimum genotype quality.
--AB-snp (default 0.15) and --AB-indel (default 0.2): allele-balance thresholds for heterozygous calls.
--hwe-filter (default 0.0): HWE p-value cutoff. The parameter help calls 0.0 “no HWE filter”, but the script applies HWE > <cutoff> literally, so 0.0 still drops any site whose HWE p-value is exactly 0 or whose HWE tag 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) -- from qc_1: left-normalized, variant-ID-assigned, dbSNP-annotated VCF.
<geno>.leftnorm.bcftools_qc.vcf.gz (+ .tbi) -- from qc_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_tstv and <geno>.leftnorm.bcftools_qc.novel_variant.snipsift_tstv -- from qc_3: TS/TV and Hom/Het summaries, split by whether the variant carries a dbSNP rsID.
<geno>.leftnorm.bcftools_qc.known_variant_sumstats and <geno>.leftnorm.bcftools_qc.novel_variant_sumstats -- from qc_3: the full bcftools stats report per split, with the SN record counts and the per-sample PSC block.
a .stdout / .stderr log 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 output/vcf_qc/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 tests/fixtures/vcf_qc/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 tests/fixtures/vcf_qc/protocol_example.genotype.chr22.vcf.gz \
    --dbsnp-variants <path/to/00-All.add_chr.variants.gz> \
    --reference-genome <path/to/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 tests/fixtures/vcf_qc/protocol_example.genotype.vcf_list.txt \
    --dbsnp-variants <path/to/00-All.add_chr.variants.gz> \
    --reference-genome <path/to/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta> \
    --cwd output/vcf_qc \
    --skip_vcf_header_filtering True \
    -j 2

Step 3b — Quality Control on data with DP/GQ/AD tags (default path)¶

The genotype VCFs above are genotype-only (GT), so QC is run with --skip_vcf_header_filtering True to bypass the depth/quality filters. To demonstrate the default QC path (which filters on read depth DP, genotype quality GQ, and allele balance from AD), the protocol ships protocol_example.genotype.withfmt.chr22.vcf.gz. Its DP/GQ/AD/AB values are synthetic placeholders (seed=22), not real measurements, and exist only to exercise the default filters.

Timing: TBD (on toy dataset)

sos run pipeline/VCF_QC.ipynb qc \
    --genoFile output/vcf_qc/protocol_example.genotype.withfmt.chr22.vcf.gz \
    --dbsnp-variants <path/to/00-All.add_chr.variants.gz> \
    --reference-genome <path/to/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta> \
    --cwd output/vcf_qc_withfmt \
    -j 2

Producing the following results (toy data, per chromosome):

In the toy protocol_example data almost all variants are novel (absent from the dbSNP 00-All set), so the known_variant Ts/Tv files are typically empty and the informative ratios come from the novel_variant files.
Example total Ts/Tv for novel variants: chr1 = 1.807, chr2 = 2.012, chr22 = 2.490.

The total Ts/Tv is read from the last (Total) column of the snpsift Ts/Tv table.

grep Ts/Tv output/vcf_qc/protocol_example.genotype.chr1.leftnorm.novel_variant.snipsift_tstv | rev | cut -d',' -f1 | rev

For known variants: in this toy dataset the known_variant Ts/Tv files are usually empty because the example variants are not present in dbSNP.

grep Ts/Tv output/vcf_qc/protocol_example.genotype.chr1.leftnorm.known_variant.snipsift_tstv | rev | cut -d',' -f1 | rev

To inspect novel-variant Ts/Tv on another chromosome (e.g. chr22):

grep Ts/Tv output/vcf_qc/protocol_example.genotype.chr22.leftnorm.novel_variant.snipsift_tstv | rev | cut -d',' -f1 | rev

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}

References¶

Carson, A. R., Smith, E. N., Matsui, H., Brækkan, S. K., Jepsen, K., Hansen, J.-B., & Frazer, K. A. (2014). Effective filtering strategies to improve data quality from population-based whole exome sequencing studies. BMC Bioinformatics, 15(1). 10.1186/1471-2105-15-125
Lek, M., Karczewski, K. J., Minikel, E. V., Samocha, K. E., Banks, E., Fennell, T., O’Donnell-Luria, A. H., Ware, J. S., Hill, A. J., Cummings, B. B., Tukiainen, T., Birnbaum, D. P., Kosmicki, J. A., Duncan, L. E., Estrada, K., Zhao, F., Zou, J., Pierce-Hoffman, E., Berghout, J., … MacArthur, D. G. (2016). Analysis of protein-coding genetic variation in 60,706 humans. Nature, 536(7616), 285–291. 10.1038/nature19057
Szustakowski, J. D., Balasubramanian, S., Kvikstad, E., Khalid, S., Bronson, P. G., Sasson, A., Wong, E., Liu, D., Wade Davis, J., Haefliger, C., Katrina Loomis, A., Mikkilineni, R., Noh, H. J., Wadhawan, S., Bai, X., Hawes, A., Krasheninina, O., Ulloa, R., Lopez, A. E., … Ye, Z. (2021). Advancing human genetics research and drug discovery through exome sequencing of the UK Biobank. Nature Genetics, 53(7), 942–948. 10.1038/s41588-021-00885-0