snpQT: flexible, reproducible, and comprehensive quality control and imputation of genomic data

Implementation

snpQT was developed as a set of nine core workflow components implemented with the nextflow workflow management system⁹. Each workflow component consists of independent containerised modules, using BioContainers curated by the bioinformatics community wherever possible¹⁰. Combining independent containerised modules into workflows, and enabling multiple workflow combinations, using the nextflow architecture, enables snpQT to be a reproducible and uniquely versatile tool for the analysis of human variant genomic data. Containers improve end-user experience and enable reproducible research by automatically provisioning bioinformatics software as required and improving numerical stability¹¹. Running individual modules in independent containers also solves a common problem when installing potentially incompatible software packages, known as "dependency hell"¹². In addition, nextflow enables caching at continuous checkpoints, so users can alter thresholds without needing to rerun earlier parts of the analysis. Briefly, if a module has the same input and parameters that a previous pipeline run has already processed, then the cached work is passed to the next module in the workflow instead of recomputing new work. This means that if a user runs multiple jobs and changes parameters at a later stage in the overall pipeline, then earlier unchanged stages are skipped, saving time.

As each genomic study is unique, this requires a tailored and flexible pipeline with informative representations of intermediate quality control data. snpQT is designed to offer multiple combinations of workflows as well as modifiable threshold parameters for multiple steps (as shown in Figure 1). Workflow A runs only once, performing a local database set up, downloading and preparing reference files^13,14 and setting up specific versions of tools using Anaconda, docker or environmental modules. Workflow B has been created for the user to remap their genomic dataset from build 38 to 37 and vice versa. Workflow C performs Sample QC, including checks for missing call rate, sex discrepancies, heterozygosity, cryptic relatedness, and missing phenotypes. Workflow D performs Population Stratification: after an internal QC of the reference genome and user’s dataset, the two datasets are merged and prepared for the automatic removal of ethnic outliers using EIGENSOFT¹⁵. Principal Component Analyses are carried out before and after the outlier removal. Workflow E performs the main Variant QC, checking missing call rate, Hardy-Weinberg equilibrium deviation, minor allele frequency, missingness in case/control status, and generates covariates for GWAS, based on a user-modifiable number of Principal Components (or users may provide a covariates file). Workflow F is for pre-imputation quality control, preparing the dataset for imputation, while Workflow G performs local phasing and imputation using shapeit4¹⁶ and impute5¹⁷. Workflow H performs post-imputation QC where poorly imputed variants are removed, different categories of duplicated variant IDs are handled and the phenotypes of the dataset are updated. The workflows’ structure also allows for users to upload their data to an external imputation server, or use a different reference panel. Workflow I performs GWAS with and without adjustment of covariates (if the covariates are not provided by the user, snpQT uses the first five Principal Components generated from population stratification in Workflow D), outputting summary statistics, along with a Manhattan plot and a Q-Q plot.

As it can be challenging to choose the correct threshold for a metric, snpQT provides a "Make report" module in each of the main Workflows C, D, E, and I, that provides interactive HTML reports summarising all the plots for both before and after the chosen thresholds have been applied, enabling the user to easily inspect and check if the chosen thresholds are appropriate. Detailed summary logs and graphs are also provided throughout, depicting the total number of samples and variants in each step, for users who need easy and fast inspection of the processes, as well as for users who want a more-in-depth report prompting users towards the locations of intermediate files and logs.

Operation

snpQT is implemented in nextflow, R and Unix command line utilities. The minimum software requirements to run snpQT are Java 8, nextflow v20.10.0, and a POSIX-compatible operating system (tested on CentOS 8). The hardware requirements scale with input data and workflows: typically quality control checks require less than 16GB of RAM and 4 cores on large datasets of 40,000 individuals. However, imputation requires significant computing power - up to 50GB of RAM per chromosome per core. As well as those already listed, the following tools are used: picard (https://broadinstitute.github.io/picard/), PLINK¹⁸, PLINK2.0¹⁹, samtools²⁰, and snpflip (https://github.com/biocore-ntnu/snpflip).

The latest release of the 1,000 Genomes Project data¹³ is used as a reference panel in both VCF and processed PLINK2 formats¹⁴. A part of the population stratification and variant QC implementation was inspired by the work of 6. Optional software requirements include Docker, Anaconda, and Environment Modules which provide a simple method to install and run the underlying collection of bioinformatics software described above.

Full documentation of snpQT, including an installation guide, a Quickstart explanation of workflow combinations and commands, a complete description of workflows, and an in-depth tutorial are provided at https://snpqt.readthedocs.io/en/latest/. The following Use Case section gives examples of input and output with explanatory context, and explains all of the key parameters needed to make use of snpQT.

Installation

Before downloading and running snpQT, depending on their needs, the user should have downloaded nextflow (v20.10.0) and Docker or Anaconda. To begin installation, the repository can be cloned and set up can be initiated by running the following commands:

git clone https://github.com/nebfield/snpQT.git
cd snpQT
nextflow run main.nf

Before starting to use any of the implemented workflows, it is necessary to set up a local database of reference and auxiliary files that snpQT requires to run. Because of the large volume of reference data in imputation workflow, we have designed two types of database. The first is the core reference database, which is sufficient for build conversion, sample and variant quality control, population stratification, pre-imputation and post-imputation. The second database is required only for local imputation, and downloading the latest release of the 1,000 Genomes Project data. In either case, it is necessary to first set up the core database:

nextflow run main.nf -profile [conda/docker/modules]
--download_db core

The above command will download the core reference files, prepare them and store them in a db/ folder in the snpQT directory. On the test system this took around 1 hour to complete, but this will vary depending on the system and network. The core dataset requires ~43GB of initial storage of intermediate files, which are stored in the work/ directory. These can be removed after set up using nextflow run main.nf –download_db core && rm -r work, leaving only 19.7GB of reference files that are stored in the database directory db/.

Alternatively, an already processed .tar.gz file (17.3GB, when unzipped it requires 19.7GB of space) of the required reference data can be downloaded directly from our servers, using the following lines of code:

mkdir db/
cd db
wget https://sys-myo.com/snpQT/core.tar.gz
tar -xvf core.tar.gz

If the user is interested in local imputation then the following command should also be run (in addition to the previous commands for the core database set up):

nextflow run main.nf -profile [conda/docker/modules]
--download_db impute

This command will download the imputation reference files and store them in a snpQT/db/impute directory. Downloading and processing the imputation data may take one to two hours or more, depending on the system and network. These reference files require an additional 30GB of storage. Alternatively, an already processed .tar.gz (13GB, when unzipped it requires 15GB of space) imputation reference file can be downloaded directly from our servers, using the following lines of code, and taking around 25 minutes:

mkdir db/
cd db
wget https://sys-myo.com/snpQT/impute.tar.gz
tar -xvf impute.tar.gz

The snpQT workflows generate many intermediate files that are not shown directly to the user, and are usually not needed by the user. Nextflow stores these in the snpQT/work/ directory by default. Due to the large size to which this directory can rapidly grow, it is recommended that the work folder be deleted after setting up the reference databases, using the command shown above. The same command can also be used to regularly clear the work/ folder as needed, such as after running multiple analyses using snpQT. However, while these intermediate files are retained, snpQT will remember previous work that it has done, such that it may automatically avoid needless repetition of previously completed work if asked to run a different stage of the pipeline, or to run again with tweaked parameters, on the same input data. When the work/ folder is deleted, all work for the analysis will need to be redone.

Datasets

After the download and set up of the database, analysis can then be done. A synthetic demonstration dataset is available with snpQT, which can be used to gain familiarity with the workflows and modules, while ensuring reproducible results identical to those shown in this section. The synthetic dataset is located within the data/ folder, and consists of a .vcf.gz file and three binary plink files (.bed, .bim and .fam). The dataset contains 6,517 randomised genotypes of chromosome 1, derived from 100 female samples having balanced binary phenotypes (i.e. 51 cases vs 49 controls). The chromosome positions, alleles and SNP IDs have been updated according to human 1,000 Genomes Project data (hg37).

Human genome build conversion

The Human Genome Build Conversion workflow converts genomic files aligned in build 38 to build 37 (default mode) and vice versa, using Picard’s LiftoverVcf utility. snpQT assumes that your input genomic data are aligned to build 37, as some of the workflows are designed to accept this input. Despite that GRCh37 human reference genome is not the most recent one, it is the most frequently used build among current public reference genomic datasets (e.g. 1,000 Genomes data, Haplotype Reference Consortium panel), online imputation servers (e.g. Sanger Imputation Server) and available snp-array datasets, and for this reason snpQT has been designed to support GRCh37 (hg19) for most workflows. Hence, this workflow can be helpful for users with data aligned to build 38 to convert them to build 37, in order to run QC and population stratification. The workflow may also be helpful when a user has finished their main QC and wish to upload their data to an external imputation server that uses a reference panel aligned in b38 (i.e. TOPMed) or b37 (i.e. Haplotype Reference Concortium panel), or to run a local imputation using a reference panel which is aligned in another build.

The genomic build of the synthetic dataset (which is aligned to b37) can be converted to b38 by running the following code, with input files and options explained below:

nextflow run main.nf -profile conda -resume \
  --vcf data/toy.vcf \
  --fam data/toy.fam \
  --convert_build \
  --input_build 37 \
  --output_build 38 \
  --results convert_build_toy/

When this work is run successfully, a new folder will be created named convert_build_toy/ which contains a files/ subfolder which contains three binary plink files (the .fam file contains updated phenotype information) aligned to b38 and a .vcf.gz file for users who prefer this format for other purposes.

Main quality control

snpQT’s main QC workflow is divided into two distinct nextflow modules: Sample and Variant QC. Sample and Variant QC can be run using the parameter –-qc. The required input files are binary plink files which can be imported using the parameters –-bed, –-bim and –-fam, for .bed, .bim and .fam plink files, respectively. The main checks of Sample QC are listed below, with accompanying parameters and explanation:

The second part of the main QC is Variant QC, which is again implemented using the --qc parameter. It is considered good practice to first filter low quality samples in order to reduce the risk of removing a potentially high-risk variant during Variant QC. For this reason, the population stratification workflow (if chosen to be run by the user, as explained in the next subsection), which is essentially a Sample QC step, is designed to run between Sample QC and Variant QC, as seen in Figure 1. The Variant QC module contains the following steps:

The synthetic toy dataset does not contain sex chromosomes, so to avoid plink producing an error it is important to add --sexcheck false to skip the step which checks for sex discrepancies. Main QC on the toy dataset can be performed by running the following command:

nextflow run main.nf -profile conda -resume \
  --bed toy.bed \
  --bim toy.bim \
  --fam toy.fam \
  --qc \
  --results results_toy/ \
  --sexcheck false

On completion, the results_toy/ directory will now contain as many folders as the number of workflows that were run. Based on the example given, a results_toy/qc/ folder should be present, containing the following sub-folders and files (this structure will be similar for most other workflows):

In Figure 2 and Figure 3, we show some examples of the output plots from the Sample and Variant QC workflows for the toy dataset. Figure 2 illustrates a sample call rate histogram and a scatterplot for the toy dataset, before and after the default threshold has been applied (indicated by a red line). Figure 3 shows one of the last processes of the Variant QC workflow, where Principal Component Analysis is performed on the clean dataset both for data exploration and for generation of covariates using the first X Principal Components, which can then be used in the GWAS workflow, to account for a potential inner population sub-structure.

Figure 2.

Sample call rate for synthetic toy dataset shown as (a) a histogram and (b) a scatterplot, before and after applying the default threshold of 0.02% (red line). This synthetic randomised dataset was created for demonstration purposes, thus the sample call rate distribution may not closely resemble a real-world dataset.

Figure 3. 3D Principal Components Analysis (PCA) of the synthetic dataset.

The samples are annotated based on their phenotype (e.g. case/control status). The 3D PCA plot is available in an interactive environment incorporated into the .html reports. PCA plots are also provided in a 2D format for the first three Principal Components.

Population stratification

The aim of the Population Stratification workflow is to identify and then remove potential ethnic outliers, based on population structure, using the EBI’s latest release of a processed phased 1000 Genomes Project reference panel, aligned to human genome build 37. Population stratification is an essential step in QC analysis, since it minimises the possibility that the difference in the allele frequencies is caused by the different ancestry of the samples. The population stratification workflow requires the main QC workflow.

During this workflow, internal processing of both the 1,000 Genomes data and of the user’s dataset are performed, and then the two datasets are merged, keeping only mutual SNPs shared by both. The internal processing consists of numerous QC steps, some of which can be tailored by the user, by passing the following parameters:

When both datasets are prepared and merged, snpQT creates a racefile labeling the race of each sample. User's samples are automatically labeled as "OWN". The race label for the 1,000 Genomes data can be controlled by the --racefile parameter, using super-population labels (e.g. EUR, AFR, AMR) by default or subpopulation labels using the --racefile [super (default), sub] parameter. When the racefile and the merged dataset are ready, EIGENSOFT's smartpca software is performed for automatic outlier removal. Smartpca takes a set of parameters, which are in the form of a file. snpQT provides the option to change this file according to the users' needs using the --parfile parfile.txt parameter. Lastly, the user can choose to infer eigenvectors based on a population subset list in smartpca using the parameter --racecode ["" (default), "EUR"/"AFR"/"SAS"].

To perform population stratification on the toy dataset the user can run the following command:

nextflow run main.nf -profile conda -resume \
  --bed toy.bed \
  --bim toy.bim \
  --fam toy.fam \
  --qc \
  --pop_strat \
  --results results_toy/ \
  --sexcheck false

After successful completion, a new sub-folder pop_strat/ will be created within the –-results directory, along with the previous qc/ folder. Since the -resume parameter was used, the Sample QC processes have been cached, making the pipeline run faster. As it was mentioned above, Population Stratification runs in between Sample and Variant QC, which means now that --pop_strat is combined with --qc the Variant QC, input files have changed and therefore, the corresponding processes will run again with the updated input. Within the pop_strat/ folder, the following are included:

Figure 4 shows the Principal Component Analysis (PCA) plot for the synthetic toy dataset. The PCA topology is quite artificial here, as this is a synthetic dataset made in plink2, containing only a few thousand genotypes of chromosome 1, which are subsequently pruned, leaving a few hundred independent SNPs merged with 1000 Genomes Project data. A more natural example of PCA plots resulting from real-world data is shown in the online tutorial for the ALS dataset.

Figure 4. 3D Principal Components Analysis (PCA) of the synthetic dataset following combination with the 1000 Genomes Project data.

The samples are annotated based on their ancestry, except for the user’s data which are labelled as "OWN". Two 3D PCA plots are available in an interactive environment incorporated in the .html reports, representing before and after ethnic outlier removal using EIGENSOFT. PCA plots are also provided in a 2D format for the first three Principal Components. Since the synthetic toy dataset is artificial and contains only a few hundred of randomised independent markers, it is located a large distance from the 1000 Genome data. For this reason, the before-and-after outlier removal PCA plots are identical, so only one is shown here.

GWAS

This workflow performs both logistic and linear regression for binary and quantitative phenotypic traits, respectively, using the parameter --gwas. snpQT performs a logistic regression by default, but it is also designed to run linear regressions using the --linear parameter. These analyses can be performed with and without adjusted covariates to account for a fine-scale population structure. Covariates can be calculated at the end of the --qc workflow (preferably used with population stratification workflow) using the first X Principal Components of the generated PCA, using the --pca_covars [3 (default), 1-20] parameter; alternatively, covariates can be passed directly from the user as an argument with --covar_file [false (default), covar.txt]. The covar.txt file should follow the same format as a PLINK covariate file. The GWAS workflow requires the main QC workflow to run in advance. For a logistic regression analysis on the toy dataset, the user can run the following command (the following example assumes that the user wishes to run population stratification as well, although this is not obligatory):

nextflow run main.nf -profile conda -resume \
  --bed toy.bed \
  --bim toy.bim \
  --fam toy.fam \
  --qc \
  --pop_strat \
  --gwas \
  --results results_toy/ \
  --sexcheck false

The command above causes a total of 39 separate snpQT processes to run, and be completed in ~2 minutes on the development system. When GWAS has finished running successfully, a new gwas/ sub-folder will be created within the results_toy/ directory, along with the previous pop_strat/ and qc/ folders mentioned above. Within the gwas/ folder the following are included:

Figure 5 and Figure 6 show Q-Q and Manhattan plots, with and without covariate adjustment for the synthetic toy dataset, respectively.

Figure 5.

Q-Q (Quantile-Quantile) plots for the synthetic dataset both (a) with and (b) without covariate adjustment, with their accompanying lambda values. This synthetic randomised dataset was created for demonstration purposes using the default thresholds, thus the Q-Q distributions may not closely resemble a real-world clean dataset.

Figure 6.

Manhattan plots for the synthetic dataset both (a) with and (b) without covariate adjustment. Each dot represents a variant located based on its genomic position (x-axis) and GWAS -log(p-value) (yaxis). The synthetic dataset at this stage contains 2,312 variants located in chromosome 1.

Pre-Imputation, imputation & post-imputation

The Pre-Imputation workflow prepares a genomic dataset for imputation, including fixing issues such as flipping SNPs that are on the reverse strand, removing ambiguous and duplicated SNPs and fixing the reference allele. If the user wishes to run the Pre-imputation workflow independently (i.e. the user is not interested in local imputation) then it is required to combine the workflow with the main QC workflow (and only optionally with the Population Stratification workflow). When Pre-Imputation QC is run independently, the --pre_impute argument is used and it cannot be combined with the GWAS, Imputation and Post-Imputation workflows. To run Pre-Imputation QC on the synthetic dataset, the user can run the following command:

nextflow run main.nf -profile conda -resume \
  --bed toy.bed \
  --bim toy.bim \
  --fam toy.fam \
  --qc \
  --pre_impute \
  --results results_toy/ \
  --sexcheck false

Pre-imputation workflow creates a preImputation/files/ directory, which contains a compressed VCF and an indexed VCF files which are ready for imputation

snpQT also offers an optional Imputation workflow, using the parameter –-impute, where the user can increase the number of markers of their genomic dataset, using EBI’s latest release of the phased 1,000 Genomes Project reference panel aligned to human genome build 37. When the --impute workflow is used, the Pre-Imputation and Post-Imputation workflows are called internally before and after phasing, and local imputation takes place, accordingly. As explained above, Pre-Imputation prepares the user’s dataset for phasing using shapeit4 and local imputation with impute5. When phasing and imputation per chromosome are finished, then the Post-Imputation workflow takes as input the imputed chromosomes, and filters out all poorly imputed variants, based on Info score and MAF. The user can alter these filters using --info [0.7 (default), 0-1] and --impute_maf [0.01 (default), 0-1] parameters, respectively. The Post-Imputation workflow also annotates missing SNP IDs and handles different categories of duplicated SNPs.

To run local imputation on the synthetic toy dataset the user can run the following command:

nextflow run main.nf -profile docker -resume \
  --bed toy.bed \
  --bim toy.bim \
  --fam toy.fam \
  --qc \
  --pop_strat \
  --impute \
  --gwas \
  --results results_toy/ \
  --sexcheck false

When imputation is complete, separate preImputation/ and post_imputation/ folders are created. The post_imputation/ directory contains the following sub-folders:

Despite that the Post-Imputation workflow can be nested under the --impute workflow, it is also designed to run independently; as some users may prefer to run imputation on external online servers, or may have already imputed data and they wish to proceed only with a Post-Imputation QC. The accepted input files for Post-Imputation are a valid VCF file and an accompanying plink .fam file, which should both contain the same samples. To run Post-Imputation QC on the synthetic toy dataset, the user can run the following command:

nextflow run main.nf -profile conda -resume \
  --vcf toy.vcf \
  --fam toy.fam \
  --post_impute \
  --results results_toy/

As already explained above, the results directory contains a new post_imputation/ folder containing the same elements as in the Imputation workflow.

Lastly, snpQT provides a command-line help page which is continuously updated, using the parameter --help. We provide further information about the implemented snpQT processes, the synthetic toy dataset and the real-world ALS dataset results in the online documentation at https://snpqt.readthedocs.io/en/latest/.

Underlying data

Zenodo: snpQT reference data (Version 0.1), http://doi.org/10.5281/zenodo.4916469²⁵

Zenodo: nebfield/snpQT: v0.1.3 -Fluffy penguin, https://doi.org/10.5281/zenodo.5083093²⁴

NCBI dbGaP: Mega-GWAS ALS I. Accession number, phs000101.v5.p1: https://identifiers.org/dbgap:phs000101.v5.p1