PGSXplorer: an integrated nextflow pipeline for comprehensive quality control and polygenic score model development

Pipeline development

PGSXplorer was developed using Nextflow (v24.04.3) (Tommaso et al., 2017) and the DSL2 language and executed within a Docker container. Dockerization was employed to ensure a consistent and reproducible environment through encapsulation of all dependencies and software required for the analysis. A suite of specialized tools for PGS model development were integrated into the pipeline, including PLINK (v1.9) (Chang et al., 2015), PRSice-2 (Choi & O’Reilly, 2019), LD-Pred2 (Privé, Arbel & Vilhjálmsson, 2020), Lassosum2 (Privé et al., 2021), MegaPRS (Zhang et al., 2021), SBayesR-C (Zheng et al., 2024), PRS-CSx (Ruan et al., 2022), and MUSSEL (Jin et al., 2023).

QC workflow of the pipeline

The QC workflow of our pipeline begins with the initial GWAS QC. This module automates key QC processes for GWAS summary statistics to ensure data integrity and reliability in downstream analyses. The QC module performs the following steps:

1. MAF filtering: Variants with MAF values below 0.01 are excluded to remove rare variants that may cause noise in the analysis.

2. INFO filtering: Variants with INFO (imputation score) values below 0.8 are removed to enhance genotyping accuracy and minimize potential bias.

3. Duplicate SNP removal: Duplicate SNPs are systematically identified and removed to avoid redundant or conflicting data.

Following these steps, the filtered GWAS summary statistics are reconstructed and prepared for subsequent stages of the workflow. These QC procedures follow standard guidelines outlined by a previous study (Choi, Mak & O’Reilly, 2020).

Our pipeline was also designed to optimize the QC of the user-provided genotype data through a series of automated steps: (i) filtering missing SNPs, (ii) filtering missing individuals, (iii) minor allele frequency (MAF) filtering, (iv) Hardy–Weinberg equilibrium (HWE) filtering, (v) relatedness check, (vi) heterozygosity assessment, and (vii) removal of duplicate SNPs. These QC steps are streamlined and are executed using PLINK (v1.9) with the parameters --geno, --mind, --maf, --hwe, --rel-cutoff, and --het (Chang et al., 2015). Additionally, customized scripts written in R (v4.1.0) were employed to visualize heterozygosity, HWE, relatedness, and MAF distributions of the given data. The specific parameters and filtering criteria are detailed below:

1. Filtering missing SNPs: To preserve the integrity and accuracy of the genotype data, the plink --geno parameter was used with a threshold value to remove SNPs with a certain percentage of missing genotypes. A value of 0.02 was determined as the default for PGSXplorer, eliminating variants with more than 2% missing SNP data across the samples (Turner et al., 2011).

2. Filtering missing individuals: After filtering missing SNPs, individuals with a high rate of missing genotypes are filtered to maintain data quality. In this step of the pipeline, the plink --mind parameter with a value of 0.02 was set as the default value (Turner et al., 2011).

3. Filtering by MAF: The MAF threshold was determined to filter out rare variants that may introduce noise or bias into the analysis. Generally, values of 0.01 to 0.05 are used for this filtering step (Pavan et al., 2020). In this step, the default value was set to 0.05 using the plink --maf parameter.

4. Filtering by HWE: Deviations from HWE were evaluated to identify genotyping errors and ensure data quality. In our study, we applied HWE filtering in a two-step process using PLINK. The first filtering step applied a strict HWE threshold of 1e−6 to the control group. This step ensured removal of SNPs that deviated significantly from HWE among controls. We then applied a second HWE threshold of 1e−10 to the case group. This less stringent step only targeted SNPs that showed extreme deviations from HWE in the case data, as the stringent threshold had already been applied to controls (Marees et al., 2018).

5. Relatedness checking: In population studies, the maximum degree of relatedness between any pair of individuals is typically expected to be less than that of second-degree relatives (Turner et al., 2011). To address this, PGSXplorer identifies and filters possible sample mix-ups or family relationships that may introduce bias downstream analyses. In this step, a default value of 0.1875 was defined for the plink --rel-cutoff parameter.

6. Heterozygosity assessment: Monitoring heterozygosity levels is essential for identifying potential contamination or issues with genotyping data quality. In the pipeline, heterozygosity filtering is performed based on +/− 3 Standard Deviations (SD) from the mean.

7. Removal of duplicate SNPs: To identify and remove duplicate SNPs from given genotype dataset, we first listed the duplicate SNPs using the following command: awk ‘{print $2}’ <input_bim_file> | sort | uniq -d > <output_duplicate_snps_list>

This command extracts SNP identifiers from the .bim file, sorts them, identifies duplicates, saves them to a list file. Then PLINK –exclude command was used to remove these duplicate SNPs from dataset. This process ensured that final dataset was free of duplicate SNPs, improving the quality and reliability of our downstream analyses.

Integration of QC, phasing and imputation steps

Following the QC steps of the workflow, the filtered datasets in PLINK format are automatically converted to VCF format for further processing. These VCF files are then phased using Eagle (v2.4.1) (Loh et al., 2016) with reference files tailored to the GRCh38 genome build (details regarding the reference files are available on PGSXplorer’s GitHub page). After phasing, imputation is performed using the same GRCh38-compatible reference datasets with Beagle (v5.4) (Browning, Zhou & Browning, 2018). Subsequent to imputation, additional QC is performed based on imputation scores to ensure data accuracy. The final files are then automatically converted back to PLINK format and made ready for target ancestry inference and PGS calculations.

To include the QC steps, phasing, imputation, and visualizations in the PGSXplorer, the inputs and outputs of each step were defined and assigned to different channels, with separate modules created for each process as shown in Fig. 1. They were carefully designed to filter out low-quality data, thereby enhancing the reliability of the downstream analyses. The initial input, provided in PLINK (bed, bim, fam) or VCF formats, completes the target and GWAS QC steps, phasing, and imputation in an automated and optimized manner. This process prepares the data for PGS calculations. The default parameters used in these steps are determined according to a previous study (Marees et al., 2018), but users can provide the desired values as parameters.

Target ancestry inference

Detecting ancestry components in genomic data is a critical step in PGS model development. In this study, the Fastmixture (Santander, Martinez & Meisner, 2024) tool was integrated into the pipeline to perform target ancestry inference. Fastmixture utilizes probabilistic modeling approaches to efficiently and accurately determine the proportions of individuals belonging to different ancestral groups by processing genotype data (Santander, Martinez & Meisner, 2024). The resulting outputs of the target ancestry inference step were used in downstream analyses to account for population structure.

PGS modeling of genomic data

After completing the QC steps, PGSXplorer integrates four well-known PGS algorithms —PLINK, PRSice-2, LD-Pred2 (both grid and auto), and Lassosum2 along with MegaPRS and SBayesR-C—to generate robust PGS models from GWAS summary statistics. Each of these tools was chosen for its distinct role and contribution to the development of PGS models. Their selection in this study was based on their widespread recognition and specialized capabilities in the field: (i) PLINK is a well-established tool in genetic association studies, widely recognized for its efficiency in processing large-scale genetic data. Its ability to filter SNPs meeting certain p-value thresholds significantly contributes to PGS calculation (Purcell et al., 2007), (ii) PRSice-2 is highly regarded for its versatility in handling large cohorts, requiring users to have a solid understanding of bioinformatics. It enables the estimation of disease risk based on genetic variants and provides flexibility in PGS model development by allowing adjustment of p-value thresholds to suit specific research goals (Choi & O’Reilly, 2019), and (iii) LD-Pred2 enhances the accuracy of PGS models by incorporating linkage disequilibrium (LD) information, a critical factor in capturing the genetic architecture of complex traits. By leveraging LD, LD-Pred2 improves the predictive power for identifying genetic variants associated with complex traits (Privé, Arbel & Vilhjálmsson, 2020). (iv) Lassosum2 estimates PGS using GWAS summary statistics alone and has demonstrated consistent improvements in prediction accuracy, particularly when modeling multiple PGS derived from various parameters. It is a valuable tool within a reference-standardized framework, especially for its ability to handle high-dimensional genetic data and improve trait prediction (Pain et al., 2020; Privé, Arbel & Vilhjálmsson, 2020). (v) MegaPRS accurately estimates the effect of genetic variants on phenotypic traits by incorporating complex inheritance models. Using the BLD-LDAK model, it considers LD, MAF and functional characteristics of SNPs, which increases accuracy and enables better modeling of genetic variance. Standing out for its computational efficiency, MegaPRS offers an effective and flexible solution for genetic predictions across different populations and traits (Zhang et al., 2021). (vi) SBayesR-C calculates PGS using GWAS summary statistics and LD matrices, modeling genetic effect sparsity through a finite mixture of normal distributions. It handles large-scale genetic data efficiently, shows competitive prediction accuracy compared to methods such as LD-Pred2 and Lassosum, and allows integration with functional annotations to increase its power. This makes it a valuable tool in the field of genetic epidemiology and individualized medicine (Zheng et al., 2024).

Multi-ancestry PGS tools increase the accuracy of PGS models by exploiting genetic variation across different populations, allowing for more precise modeling of allele frequencies and LD patterns (Chen et al., 2014; Ruan et al., 2022). This approach helps overcome the limitations of traditional methods that rely heavily on European-ancestry data and leads to better estimates for non-European populations (Ge et al., 2022; Shim et al., 2023). To improve the accuracy and comprehensiveness of PGS estimates, PGSXplorer includes PRS-CSx and MUSSEL, both of which significantly enhance the pipeline’s capabilities. (vii) PRS-CSx improves PGS predictions by integrating GWAS summary statistics from multiple populations, which is crucial for accurately estimating disease risk across diverse genetic backgrounds. This tool addresses the common problem of reduced predictive power in non-European populations using a cross-population approach (Ruan et al., 2022). (viii) MUSSEL further strengthens the pipeline by applying advanced statistical techniques such as clustering, thresholding, empirical Bayes, which are particularly effective in optimizing PGS across different ancestries (Jin et al., 2023). Together, these tools overcome the limitations of single ancestor models, making the PGSXplorer workflow for diverse global populations, and thus increasing its value in personalized medicine and risk assessment.

Integration of PGS tools

Incorporating PGS models into PGSXplorer required carefully structuring the inputs and outputs for each step, which were subsequently assigned to distinct Nextflow channels to ensure seamless data flow throughout the pipeline. As illustrated in Fig. 2, each modeling process is encapsulated within a separate module, enhancing modularity and flexibility of our pipeline’s design. Tools that model PGS using GWAS data from a single population are categorized as Single PGS, whereas those that incorporate data from at least two distinct populations are defined as Multi PGS.

Generation and preparation of synthetic data

To validate the functionality of PGSXplorer, synthetic genotyping data were generated using the HAPNEST (Wharrie et al., 2023). HAPNEST facilitates the creation of synthetic genomic datasets representing various ethnicities, making it ideal for testing and validation purposes. Specifically designed to support genomic research, HAPNEST utilizes containerization through Docker or Singularity, ensuring reproducibility and ease of use. Key features of HAPNEST include the ability to fetch diverse reference datasets customize parameters for specific research needs, and standardize the software environment via containerization.

In our study, the sizes of the synthetic datasets were chosen to present different population sizes to test the performance and efficiency of our pipeline under various conditions. The datasets included 500, 1,000 and 10,000 individuals of European (EUR) origin, and 3,000 and 10,000 individuals of East Asian (EAS) origin. The datasets of 500 EUR, 1,000 EUR, and 3,000 EAS individuals were labeled as T1, T2, and T3, respectively. Using these datasets, both the data processing capacity and computation times of the pipeline were evaluated. Synthetic genomic data were generated for all chromosomes, ensuring comprehensive coverage. Different populations were simulated by modifying polygenicity and genotype proportion values in the HAPNEST configuration file. The commands used to generate these datasets with HAPNEST are as follows:

Genotype data generation: docker run -v /HAPNEST/data:/data -it sophiewharrie/intervene-synthetic-data generate_geno 16 /data/config.yaml

Phenotype data generation: docker run -v /HAPNEST/data:/data -it sophiewharrie/intervene-synthetic-data generate_pheno data/config.yaml

The configuration files used in this process, which detail the parameters and population structures, have been shared on the PGSXplorer GitHub page (https://github.com/tutkuyaras/PGSXplorer), along with the datasets themselves. These synthetic datasets were also utilized to calculate GWAS summary statistics using PLINK2 (Chang et al., 2015). Logistic regression models were employed with the command: plink2 --bfile target --pheno phenotype_file.txt --glm hide-covar --covar target.eigenvec --ci 0.95 --out gwas_sumstat.

Systematic archiving and visualization of genomic quality control

PGSXplorer is designed to systematically archive the outputs generated at each stage of the QC process in dedicated directories. This approach ensures that the data from all seven QC steps are comprehensively recorded and filtered according to user-defined parameters, enhancing transparency and facilitating downstream analyses. The outputs, ranging from initial data filtration to final QC reports, are clearly organized and readily accessible for further review.

In addition to this structured archiving, the QC module also includes automated graphical representations of key metrics. As shown in Fig. 3, distributions of heterozygosity, HWE, inbreeding coefficients (Pi-hat or IBD), and MAF are visualized and automatically generated for user inspection. QC graph results for T1 and T3 are also given in Figs. S1 and S2, respectively. These visual outputs offer critical insights into the quality of genomic data, enabling rapid identification of potential issues such as population stratification or genotyping errors. By combining progressive data archiving with real-time graphical analysis, the QC process becomes both comprehensive and user-friendly, ensuring high-quality data is prepared for subsequent genomic analyses.

The SNP and individual information eliminated according to the parameters used in the QC steps are presented to the user. Table 1 shows the number of SNPs eliminated during QC for the three datasets tested with PGSXplorer. Since a synthetic dataset was used and the GRCh38 rsID (Reference SNP cluster ID) list provided by HAPNEST is common, the initial number of variants is identical. However, the number of eliminated variants and the final number of variants remaining differ across populations. The default parameters of PGSXplorer were applied during QC steps, but users can modify these parameters to suit their specific research objectives and study requirements.

Target ancestry inference

Fastmixture, integrated into the pipeline for target ancestry inference, offers an efficient method for analyzing the population structure of genetic data (Santander, Martinez & Meisner, 2024). The Q file generated by this module provides users with ancestry proportions for individuals, indicating the probabilities of each individual belonging to different populations. The .p file, on the other hand, contains allele frequencies specific to each population for genetic variants. This file is particularly useful for examining genetic differences between populations and identifying population-specific genetic patterns.

PGS modeling

The genomic data generated during the QC steps serve as input to the PGS modules, enabling customized analyses. Users can specify which models to execute by using the command nextflow run main.nf–help, ensuring that only the desired modules are and unnecessary computations are avoided. By default, all tools are set to run automatically, except for MUSSEL, which is disabled due to its high computational demands. MUSSEL requires substantial processing power and is designed to be executed only on servers, providing users with the flexibility to opt in for resource-intensive analyses.

As outputs of the PLINK algorithm, the first of the single PGS models, individual-based PGS were generated across seven different p-value thresholds. The resulting scores were stored separately in the “outputs” folder. Table S1 presents the AUC and R² values calculated for the PGS from three different target datasets during the experimentation and optimization phases of PGSXplorer.

The PRSice-2 module generates visual outputs that summarize the performance of PGS models. Figure 4 illustrates the bar plot and high-resolution plot generated based on p-value thresholds and PGS model fit. Along with these visual outputs, the pipeline provides several additional result files, including a list of the best-fitting PGS values determined by the regression model (.best), the number of SNPs used for PGS calculation at each p-value threshold, along with the corresponding R² and p-values (.prsice), a summary of the model (.summary), and detailed log files. For the LD-Pred2 method, two approaches—auto and grid—have been integrated into the pipeline, offering graphical outputs such as heritability (h²) and p-value plots for the auto model, and GLM-Z score plots for the grid model, as shown in Fig. 5. Additionally, files containing individual principal components and PGS are included in the outputs folder.

Lassosum2, which leverages penalized regression to account for LD, provides an efficient and scalable approach for large datasets. As part of the output, graphical representations, including GLM-Z score graph and post-processed (Postp) plots, are provided to visualize the performance of the Lassosum2 model, as shown in Fig. 6 (the plot for T3 is presented in Fig. S3).

MegaPRS and SBayesR-C offer complementary functionalities for PGS modeling, each providing detailed outputs tailored to user needs. MegaPRS supports various statistical models, with BayesR set as the default, allowing users to select the most appropriate model for their data. The outputs include summary files, parameters, correlations, and the model that delivers the best result, all stored in the output folder. Similarly, SBayesR-C generates comprehensive results, including SNP weights and other key metrics for PGS calculation. The output file contains information such as SNP IDs, effective alleles, combined effects at the genotype scale (BETA), probabilities of causation (PIP), and effects at the last iteration (BETAlast).

PRS-CSx, a cross-population PGS modeling approach, calculates PGSs separately for each of the two or more GWAS datasets used. If the --meta parameter is applied, the output includes a meta PGS file that combines SNP effect sizes across populations using inverse-variance weighted meta-analysis of the population-specific posterior effect size estimates. All generated files are saved in the outputs folder. Figure 7 shows the distribution of PRS-CSx module results computed with PGSXplorer, using the T2 (EUR-1,000) and T3 (EAS-3,000) datasets. Figure S4 presents similar results for the T1 (EUR-500) dataset. The integrated MUSSEL tool generates several important outputs. It produces population-specific PGS files with scores calculated for each population based on their genetic datasets and summary statistics, as well as meta-PGS files created using the inverse variance weighted meta-analysis method to combine scores across populations. Additionally, it includes files with SNP-level posterior effect size estimates for each population and meta-analysis results. Configuration files detailing the parameters used in the analysis, such as selected SNPs and population settings, are also provided to the user.

Computational performance metrics of PGSXplorer

Computational metrics of the PGSXplorer tool, such as CPU utilization and execution time, are provided to users through Nextflow parameters: -with-report pipeline_report.html and -with-timeline timeline.html. These metrics offer detailed insights into the computational performance of the pipeline. We present these metrics for three synthetic datasets of varying sizes—T1 (EUR-500), T2 (EUR-1,000), and T3 (EAS-3,000)—used to demonstrate the proper functioning of PGSXplorer. Analyses were conducted on chromosomes 1 and 2, and detailed metrics, including runtime and CPU utilization for each dataset, are provided in Table S2. This comprehensive reporting provides users with a clear understanding of the pipeline’s resource requirements across various data sizes and scenarios.