Fine Mapping Coronary Artery Disease Loci in UK Biobank’s Whole Genome Sequencing Data by Optimising Usage of the Cloud-Based Research Analysis Platform

RAPpoet: A Driver-Worker Approach to Running Jobs on the UKB RAP

As of writing, the web UI of RAP and available tutorials are tailored towards processing UKB’s genomic data sequentially or in small batches. While this was well suited to earlier and smaller versions of the UKB datasets (e.g., array imputed data), with 145,207 pVCFs this approach has become impractical.

To address this, we have designed a scalable orchestration engine called RAP parallelisation orchestration engine template (RAPpoet), to streamline our pipeline across the UKB pVCFs using the DNAnexus CLI (Figure 1).

• Download figure
• Open in new tab

Figure 1. A Driver-Worker Approach (RAPpoet) for Managing the Configuration of Requested Instances on UK Biobank RAP and for Parallelising the Processing of Files on Cloud Instances.

The driver script, which remains on the user’s local machine, includes configurable environmental variables, such as the type of instance requested and number of batches to process in each instance. Meanwhile, the worker script, located in the RAP project environment, contains the actual commands to be executed on the instance that can be set by the user depending on their requirements.

RAPpoet utilises two scripts: the ‘driver’ and the ‘worker’. The ‘driver’ script sits within the user’s local compute and includes commands for configuring the instance environment, uploading required files, and initiating the execution of the ‘worker’ scripts on requested instances. The ‘worker’ script is uploaded onto each instance and defines the processes to be run over the uploaded files.

This setup enables the parallelisation of tasks as the ‘worker’ script can execute the defined process concurrently via the xargs tool, optimising resource utilisation by getting a single instance to handle multiple files. This also decreases the number of instances needed making it easier to manage and monitor jobs.

Additionally, RAPpoet streamlines the workflow by centralising the control and coordination in the ‘driver’ script. This improves the overall manageability of jobs by consolidating file setup and instance configuration within the ‘driver’ script. Essential variables, such as instance types, file lists, batch sizes, and the number of files run in parallel, are defined, ensuring consistency across all instances. Furthermore, centralising control to the ‘driver’ script mitigates potential issues related to file duplication and ensures a more organised and streamlined execution.

Finally, RAPpoet facilitates efficient multi-tasking. While the RAP’s web UI typically executes one command at a time on uploaded files, our approach allows the user to define a more comprehensive process within the ‘worker’ script, such as running multiple bcftool commands in a single job. This capability also enables the execution of multiple steps calling on different tools within a single job, further enhancing the flexibility and efficiency of each task.

Machine Learning Methods Show Greater Sensitivity Than Logistic Regression

To demonstrate how RAPpoet streamlines an association analysis, a pipeline was run from quality control (QC) to fine mapping.

The significance of the well-established 9p21.3 (CDKN2B-AS1) risk locus [7,8] was confirmed with single-SNV LR (sLR) (Figure 2, Supplementary Table S2). The unprecedented density of the UKB WGS dataset further enabled fine mapping to identify the most likely causal variant. The same 13 SNVs were identified in the first credible sets of both FINEMAP [9] and SuSie [10] (Figure 3; Supplementary Tables S3 and S4). Of these, rs10757274 has been linked to CDKN2B-AS1 expression [11,12], is part of the core risk haplotype region [13], and has been associated with CAD in several populations [14,15], suggesting that rs10757274 is the most likely causative SNV in the risk locus.

• Download figure
• Open in new tab

Figure 2. Miami Plot of Chromosome 9 for single-SNV Logistic Regression (Top) and VariantSpark (Bottom).

The locus 9p21.3 reached genome-wide significance in the logistic regression analysis and is the top association from VariantSpark. Note that the Gini score (bottom y-axis) is used to measure variable importance and while it is analogous to the traditional P-values of logistic regression to rank associations relative to other associations, it does not include multiple testing and does not have guiding significance thresholds.

• Download figure
• Open in new tab

Figure 3. Scatterplot of SNVs in the 9p21.3 locus with Posterior Inclusion Probability (PIP) (y-axis) against Chromosome 9 Position (x-axis) for FINEMAP (A) and SuSie (B).

The first credible set of both approaches include the same 13 SNVs which are highlighted in dark purple and annotated with rsIDs.

The same region and rs10757274 were also identified using the ML GWAS platforms VariantSpark [16] (Figure 2; Supplementary Table S5) and REGENIE [17] (Supplementary Table S6).

To further test the differences between sLR and ML on this richer data resource, the same association pipeline was also applied to chromosome 21. Figure 4 shows that sLR did not find any significant genome-wide SNVs (Supplementary Table S7). However, REGENIE found one significant genome-wide SNV, rs28451064 (Supplementary Table S8) in the known CAD locus 21q22.11 (KCNE2) [18] and VariantSpark also identifies this SNV as the top association (Supplementary Table S9), demonstrating that these ML-based approaches are more sensitive than sLR. The rs28451064 SNV is intergenic between long non-coding RNA (lncRNA), LINC00310 and the KCNE2 gene and has been previously associated with CAD as a putative functional SNV [19,20].

• Download figure
• Open in new tab

Figure 4. Miami Plot of Chromosome 21 for Logistic Regression (Top) and REGENIE (Bottom).

The locus 21q22.11, which has been previously associated with CAD, is only showing marginal association in the single-SNV logistic regression analysis. Conversely, a single SNV in the locus reached genome-wide significance in the REGENIE analysis.

Taken together, these results highlight that the era of ML in genomics could be ushered in by the unprecedented data volumes, both in cohort size and density, that UKB-style mega-biobanks offer, potentially leading to novel findings that traditional statistical approaches may overlook.

Optimising Cloud Compute and Parallelisation is Crucial to Reducing Cost and Runtime on RAP

On the RAP, analyses are carried out on AWS Elastic Compute Cloud (EC2) instances with different options for storage, memory capacity, and number of cores. In this study, we evaluated the impact of different compute architectures and of running files in parallel versus sequentially on cost and runtime, particularly for the requisite quality control (QC) steps.

Running QC step one (Figure 5) on a single chromosome 21 pVCF of 10.5GB, the average pVCF size, using the RAP’s Web UI and Workflow approach took 30 minutes. In contrast, using RAPpoet to parallelise across multiple files, the average runtime for one pVCF was 1.75 minutes (Supplementary Table S10), a 94% reduction in runtime compared to the sequential approach.

• Download figure
• Open in new tab

Figure 5. Schematic of Pipeline for UK Biobank 500K Whole Genome Sequencing Dataset on RAP.

RAPpoet was used to facilitate this pipeline on the UKB RAP.

Furthermore, by optimising compute instances, the cost associated with QC step one of the pipeline (Figure 5) was reduced by 44% (£0.029 per file vs £0.052) between chromosome 9 and 21 while keeping runtime constant. Initially, the instance type mem1_ssd1_v2_x72 was selected resulting in an average runtime of 1.75 minutes and an average cost of £0.052 per pVCF (Supplementary Table S10). However, by switching to the smaller instance mem2_ssd1_v2_x48 type (i.e., half the number of vCPUs), the average runtime and cost per pVCF for QC step one on chromosome 9 was 1.80 minutes and £0.029 (Supplementary Table S11).

Other settings that impact cost and runtimes such as job priority and file I/O were also benchmarked (Supplementary Material S1). Notably, the dxFUSE I/O filesystem limits the number of files that can be processed concurrently, severely hindering the level of parallelisation of RAPpoet and requiring a looping system within RAPpoet to process the 6,106 pVCFs from chromosome 9.

Cohort Determination

Participants were determined to have coronary artery disease (CAD) using International Classification of Disease (ICD)-9, ICD-10, and OPCS Classification of Interventions and Procedures (OPCS)-4 codes (i.e., 410, 411, 413, 414, I20-25, Z951, Z955, K40-46, K49, K50, and K75) from their linked electronic health records. The healthy controls were participants who reported to be non-smoking with a BMI < 30 and had no diagnoses of other cardiovascular diseases and known comorbidities, such as epilepsy and arrhythmias.

Sample Quality Control (QC)

Sample QC measures that were provided by the UKB[24] (provided in square brackets are their FieldIDs from the UKB AMS) were also used to filter our cohort. Samples were retained that (i) had matching reported and genetic sex [22001], (ii) were not carriers of full or mosaic sex chromosome aneuploidies [22019], (iii) were not related to other UKB participants to the third degree [22021], (iv) were genetically identified as white British [22006]. Further GWAS QC thresholds (i.e., call rates > 98% and +/-3 standard deviations of calculated heterozygosity rate mean) were included to filter out low sample quality.

After sample QC, there were 24,954 CAD cases with 25,858 controls. While there were more healthy controls available (n = 182,772), the ratio of CAD cases and controls were kept relatively even to control for inflation of test statistics when using PLINK2[3] for single-SNV logistic regression (sLR) and to prevent convergence issues with REGENIE[17].

Variant Quality Control (QC)

The following thresholds on GraphTyper metrics were used to ensure that variants called in the pVCFs were of high-quality: AAScore > 0.15, ABHet > 0.175, ABHom < 0.9, PASS_ratio > 0.05, QD > 6, and QUAL >= 10. Structural variants such as indels were also removed and multi-allelic variants were split for ease in downstream analysis.

Further GWAS QC thresholds were used to make certain only reliable variants were included in subsequent analysis. This included Hardy-Weinberg equilibrium P < 1 x 10^-6, minor allele frequency (MAF) < 1 x 10^-5, and a call rate > 98%.

After variant QC, there were 6,392,685 single nucleotide variants (SNVs) in chromosome 9 and 1,958,444 SNVs in chromosome 21. Due to the available sample size after QC, we added a further MAF 0.1% filter to remove ultra-rare SNVs before running sLR, REGENIE, fine mapping and VariantSpark. At this more stringent MAF threshold, there were 597,370 SNVs and 187,146 SNVs in chromosomes 9 and 21 respectively.

Single-SNV Logistic Regression (sLR) with PLINK

A sLR with CAD status as response and the QC SNVs were modelled using PLINK2. The covariates included in the model were age, sex, BMI, and the first 10 principal components. The principal components were calculated using the UKB imputed genotypes to get genome-wide coverage. Genome-wide significance was the accepted P < 5 x 10^-8 while marginal significance was P < 1 x 10^-5.

Whole-Genome Regression with REGENIE[17]

Whole genome regression, employing the REGENIE method, was used as an alternative approach for testing SNV associations with the same sets of covariate and phenotype data as above. The REGENIE run was divided into two steps. Initially, ridge regression predictors (j) are constructed for each block of the genome by partitioning it into consecutive blocks of B SNVs. These predictors are then combined through cross-validation to yield a single optimal predictor for estimating phenotype values. Subsequently, this predictor is divided into 23 chromosome-specific predictions for a leave-one-chromosome-out (LOCO) approach which is effective in preserving the power to detect associations by preventing proximal contamination[25,26]. The second step of REGENIE tests the association of each SNV by conditioning on the LOCO predictors, which are used as covariates to account for underlying background polygenic effects. The fast Firth logistic regression method provided by

REGENIE was used for association testing, which improves the accuracy of logistic regression models by minimising bias in maximum-likelihood estimates through penalised likelihood. This approach ensures well-calibrated Type 1 errors and provides reliable estimates of SNV effect sizes and standard errors, especially when rare variants are being tested in the SNV set.

While both steps were available as applets in the RAP tool library, the initial step was conducted on our internal high-performance computing (HPC) environment using the previously downloaded imputed array genotype dataset to mitigate costs on the RAP. This approach was taken assuming that the predicted phenotype values would consistently capture the background genome-wide polygenic effect on the phenotype. Specifically, the genome-wide imputed array SNVs (with MAF > 0.01, missingness < 0.02, Het > 3, and GBR individuals) were employed to develop the regression predictors. The parameters included B = 500, j = 5 (default shrinkage parameter value = [0.01 0.25 0.5 0.75 0.99]), and 5-fold cross-validation to select the most effective predictors.

The LOCO predictions were then exported from the HPC to the RAP for the second step where the WGS QC’d SNVs on chromosomes 9 and 21 were tested using the fast Firth approximation implementation. Genome-wide significance was the accepted P < 5 x 10^-8 while marginal significance was P < 1 x 10^-5.

Machine Learning Analysis with VariantSpark

VariantSpark is an Apache Spark-based application for building genome-wide random forest models and is not yet included in the RAP tool library. Following the guidelines available at https://documentation.dnanexus.com/developer/apps/developing-spark-apps, VariantSpark’s executable was created as a custom applet to run on RAP’s Spark cluster cloud architecture within the project space. To resolve conflicts between JAR file versions, VariantSpark needed to be compiled within the startup script rather than at applet build time, sacrificing some initialisation performance to ensure compatibility.

Smaller instance sizes are unable to handle the Spark application, which necessitates the use of more compute intensive and expensive instances for worker nodes. Because RAP does not currently differentiate between worker nodes and the primary node when specifying instance types, this necessarily increases the cost associated with the now underutilised driver node. As such, specifications allowing for the sizes of worker and driver nodes to be optimised separately by matching node capabilities more closely to their operational roles will result in better resource management and cost-efficiency.

The hyperparameters of the random forest models for both chromosomes 9 and 21 were no_of_tree (number of trees in forest) = 10K, min_node_size (minimum number of samples in node to be processed) = 10K, and mtry_fraction (fraction of variants evaluated at each node of a tree) = 0.1. The importance of a variable in the model was assessed using average Gini impurity measure.

The VariantSpark analysis of chromosome 9 was run on eight on-demand ‘mem2_ssd2_v2_x48’ instances which completed in 13 hours and 16 minutes, costing £179.31. In comparison, the VariantSpark analysis of chromosome 21 used six on-demand ‘mem2_ssd2_v2_x32’ instances which took 8 hours and 19 minutes to complete and cost £56.14.

Annotating Association Results with Annovar[27]

Annovar was used to annotate association results from sLR, REGENIE, and VariantSpark. Namely, mapping SNVs to dbSNP rsIDs, genes and genetic function from UCSC’s RefSeq, and previous studies from the GWAS Catalog.

Fine Mapping with PolyFun (SuSIE and FineMap)

A set of 5,179 SNVs at the associated CDKN2B locus, including a 500Kb flanking region, was subjected to functionally informed fine mapping using PolyFun[28] incorporating both the SuSie[10] and FINEMAP[9] methods. The approach adhered to the standard PolyFun workflow, integrating prior causal probabilities based on per-SNV heritability and summary LD information pre-calculated using individuals of European ancestry from the UK Biobank cohorts.

As the pre-calculated information was based on the hg19/GRCh37 reference, liftOver[29] was used to map the sLR summary statistics from hg38/GRCh38 to hg19/GRCh37. Following this, variants without pre-computed values were filtered out, leaving 4,843 remaining variants. For both methods, a maximum of 10 causal SNVs was specified for the locus. The locus spanned four overlapping LD regions, with variants assigned to region-specific causal sets.

Note that in the absence of a PolyFun applet on the RAP, the summary statistics file generated via sLR was downloaded from the RAP, and fine mapping was performed on an internal HPC environment.