DistPCA is a distributed out-of-core C++ framework for tera-scale genomic Principal Component Analysis (PCA), designed to scale efficiently across both single- and multi-node computing systems. Built on top of Message Passing Interface (MPI), it employs a hybrid multi-level parallelism scheme combining multiprocessing, OpenMP multithreading, SIMD vectorization, and double buffering across all three stages of the PCA pipeline (I/O, data preprocessing, numerical method). Evaluated on datasets reaching up to 11 TB, DistPCA achieves speedups of up to 58.2× and over 98% reduction in wall-clock time, while maintaining parallel efficiency above 82% and preserving the accuracy of the recovered principal components (PCs). For a detailed description of the framework and experimental evaluation, please refer to our preprint.
- Prerequisites & Installation
- Usage
- Datasets
- Performance Evaluation
- Reproducibility
- File Structure
- Planned Features
- Citation
- Acknowledgments
Clone the repository:
git clone https://github.com/CEID-HPCLAB/DistPCA.git
cd DistPCAInstall Intel MKL (Base Toolkit) and Intel MPI + OpenMP (HPC Toolkit), which provides the mpicxx and mpicc wrappers:
sudo wget -qO- https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB | gpg --dearmor | sudo tee /usr/share/keyrings/oneapi-archive-keyring.gpg > /dev/null
echo "deb [signed-by=/usr/share/keyrings/oneapi-archive-keyring.gpg] https://apt.repos.intel.com/oneapi all main" | sudo tee /etc/apt/sources.list.d/oneAPI.list
sudo apt update
sudo apt install intel-basekit intel-hpckitThen, initialize the environment and build:
source /opt/intel/oneapi/setvars.sh
make # compile
make clean # remove build artifacts before rebuildingThe executable will be available at build/DistPCA.exe.
Important
Before running DistPCA, make sure the Intel oneAPI environment is initialized with source /opt/intel/oneapi/setvars.sh
Set the number of OpenMP threads per MPI process:
export OMP_NUM_THREADS=<num_threads>After compilation, run DistPCA from the repository root:
mpirun -np <num_processes> ./build/DistPCA.exe \
-bfile <file_path> \
-nsv <nsv> \
-nrhs <nrhs> \
-power <num_power_iterations> \
-crit <convergence_criterion> \
-tol <convergence_tolerance> \
-bsize <block_size> \
-miter <max_iterations> \
-verbose <verbose> \
-fwrite <save_output> \
-fullSVD <full_svd>| Parameter | Description |
|---|---|
-np |
Number of MPI processes (Mandatory) |
-bfile |
Path to the input .bed dataset file (e.g., ../example/ToyHapmap) (Mandatory) |
-nsv |
Number of sought principal components (Default: 10) |
-nrhs |
Dimension of the target subspace (Default: 2 * nsv) |
-power |
Number of power iterations to perform (Default: 1) |
-crit |
Convergence criterion (Default: 2) |
-tol |
Convergence tolerance (Default: 1e-3) |
-bsize |
Total number of SNPs per block (Default: 100) |
-miter |
Maximum iterations to run if convergence criterion is taking longer to achieve (Default: 100) |
-verbose |
Logging level. If set to 2, detailed convergence info is printed (Default: 1) |
-fwrite |
Boolean flag. If set to 1, stores the singular values and singular vectors (Default: 0) |
-fullSVD |
Boolean flag. If set to 1, computes full SVD using LAPACKE (only if dataset fits in RAM) (Default: 0) |
Note
DistPCA supports three convergence criteria. The first is the trace-based criterion, which monitors the relative change of the total explained variance (trace) between successive iterations. The second is the individual eigenvalue criterion, which checks the relative change of each singular value and requires all components to satisfy the specified tolerance. The third is the Mean Explained Variance (MEV) criterion, which assesses subspace convergence by measuring the average squared cosine similarity between successive eigenvector estimates. By default, the MEV criterion is used.
Note
DistPCA supports three MPI-based parallelism schemes for computing the sought PCs. The first scheme, implemented in SubspaceIteration_MPI, is an in-core method used when each MPI process can fully load its assigned portion of the dataset into RAM. The second scheme supports out-of-core computation of PCs and uses three levels of parallelism (multiprocessing, OpenMP multithreading, and SIMD vectorization). It is accessible through BlockSubspaceIter_MPI_OOC. The third scheme is the full DistPCA implementation presented in the paper and extends the second scheme by additionally supporting compute–transfer overlap using a double-buffering strategy. It is implemented in BlockSubspaceIter_MPI_OOC_double_buffering.
The datasets used in this research consist of three real-world and three synthetic datasets. Real-world datasets require preprocessing with PLINK, which can be installed as follows:
# Download and install PLINK
wget -q https://s3.amazonaws.com/plink1-assets/plink_linux_x86_64_20231211.zip
unzip -q plink_linux_x86_64_20231211.zip -d plink_tmp
sudo mv plink_tmp/plink /usr/local/bin/
# Clean up intermediate files
rm -rf plink_tmp
rm -f plink_linux_x86_64_20231211.zipThe three real-world datasets used in this work are the 1000 Genomes Project, the Simons Genome Diversity Project (SGDP), and the Human Genome Diversity Project (HGDP). Their dimensions after preprocessing are summarized below, together with the commands required for downloading and preprocessing each dataset.
| Dataset | Individuals | SNPs |
|---|---|---|
| 1000 Genomes | 2,490 | 1,664,505 |
| SGDP | 345 | 694,659 |
| HGDP | 942 | 133,594 |
1000 Genomes Dataset
# Download from figshare
wget -qO 1000G.zip "https://api.figshare.com/v2/articles/9208979/download"
unzip -q 1000G.zip
rm -f 1000G_phase3_common_norel.fam2
unzip -q 1000G_phase3_common_norel.zip
# Population ancestry panel (used to color populations in Figure 6)
wget -q https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel
mv integrated_call_samples_v3.20130502.ALL.panel ./docs/results/accuracy
# Preprocessing
plink --bfile 1000G_phase3_common_norel --maf 0.01 --make-bed --out 1000G.qc
plink --bfile 1000G.qc --indep-pairwise 1000 50 0.2 --out 1000G.qc.prune
plink --bfile 1000G.qc --extract 1000G.qc.prune.prune.in --make-bed --out 1000G.qc.pruned
mv 1000G.qc.pruned.{bed,bim,fam} ./scripts/experiments
# Clean up intermediate files
rm -f 1000G.zip 1000G_phase3_common_norel.{zip,bed,bim,fam,log,fam2} 1000G.qc.* 1000G.qc.pruned.log 1000G-phase3-common-norel.RNote
The 1000 Genomes dataset can also be downloaded by running fetch_1000G.sh, located in scripts/data/.
Simons Genome Diversity Project (SGDP) Dataset
# Download from Reich Lab
URL="https://sharehost.hms.harvard.edu/genetics/reich_lab/sgdp/variant_set/cteam_extended.v4.maf0.1perc"
wget -q "${URL}.bed"
wget -q "${URL}.bim.zip"
wget -q "${URL}.fam"
unzip -q cteam_extended.v4.maf0.1perc.bim.zip
rm cteam_extended.v4.maf0.1perc.bim.zip
# Preprocessing
plink --bfile cteam_extended.v4.maf0.1perc --maf 0.01 --make-bed --out sgdp.qc
plink --bfile sgdp.qc --indep-pairwise 1000 50 0.2 --out sgdp.qc.prune
plink --bfile sgdp.qc --extract sgdp.qc.prune.prune.in --make-bed --out sgdp.qc.pruned
mv sgdp.qc.pruned.{bed,bim,fam} ./scripts/experiments
# Clean up intermediate files
rm -f cteam_extended.v4.maf0.1perc.{bed,bim,fam,log}
rm -f sgdp.qc.{bed,bim,fam,log,hh,nosex}
rm -f sgdp.qc.prune.*
rm -f sgdp.qc.pruned.{log,hh,nosex}Note
The SGDP dataset can also be downloaded by running fetch_SGDP.sh, located in scripts/data/.
Human Genome Diversity Project (HGDP) Dataset
# Download from Reich Lab
wget -q https://reichdata.hms.harvard.edu/pub/datasets/humanOrigins/Harvard_HGDP-CEPH.tgz
tar -xzf Harvard_HGDP-CEPH.tgz
# Preprocessing
plink --file Harvard_HGDP-CEPH/all_snp --make-bed --out hgdp
plink --bfile hgdp --maf 0.01 --make-bed --out hgdp.qc
plink --bfile hgdp.qc --indep-pairwise 1000 50 0.2 --out hgdp.qc.prune
plink --bfile hgdp.qc --extract hgdp.qc.prune.prune.in --make-bed --out hgdp.qc.pruned
mv hgdp.qc.pruned.{bed,bim,fam} ./scripts/experiments/
# Clean up intermediate files
rm -f hgdp.{bed,bim,fam,log,hh}
rm -f hgdp.qc.{bed,bim,fam,log,hh}
rm -f hgdp.qc.prune.{prune.in,prune.out,log,hh}
rm -f hgdp.qc.pruned.{log,hh}
rm -rf Harvard_HGDP-CEPH*Note
The ΗGDP dataset can also be downloaded by running fetch_ΗGDP.sh, located in scripts/data/.
Warning
After running the commands above, or the corresponding .sh script for each dataset, the resulting PLINK binary files (.bed, .bim, .fam) will be stored in scripts/experiments/.
The three synthetic datasets used in this work are:
| Dataset | Individuals | SNPs |
|---|---|---|
| 50K Genomes | 50,000 | 6,000,000 |
| 500K Genomes | 500,000 | 3,000,000 |
| 1M Genomes | 1,000,000 | 1,000,000 |
Synthetic datasets can be generated using DataSimulator. First, install the required dependencies and build:
sudo apt-get install libboost-all-dev libgsl-dev
git clone https://github.com/eugeniamaria/DataSimulator.git
cd DataSimulator && makeThen, run:
./GeneticDataSimulator -npop [int] -nregions [int] -nindividuals [int] -nSNP [int] -minfreq [double] -txtoutput [int] -filename [char]This generates two output files: output_file.map (SNP information) and output_file.ped (individual genotypes).
Warning
Generating large datasets (e.g., 1M individuals × 1M SNPs) may require substantial disk space due to the size of the resulting text-based .ped file. For large-scale datasets, it is recommended to generate the data in parts and merge the resulting files into PLINK binary format (.bed, .bim, .fam) using PLINK, rather than generating a single large .ped file.
The performance of DistPCA was evaluated on three synthetic datasets and three publicly available real-world datasets of varying sizes, as described in Section Datasets. Throughout all experiments, the underlying Randomized Subspace Iteration (RSI) method targets the leading
The experiments were conducted on the ARIS supercomputer, a national Greek HPC cluster facility, using four thin compute nodes. Each thin node is partitioned into eight Non-Uniform Memory Access (NUMA) domains and is configured as follows:
| Component | Details |
|---|---|
| CPU | Dual-socket AMD EPYC 7742 (128 cores, 2.25 GHz) |
| RAM | 512 GB (restricted to 64 GB per node for all experiments) |
| Filesystem | GPFS |
A detailed overview of the ARIS infrastructure is available here.
Note
MPI ranks were distributed across NUMA domains, with OpenMP threads pinned to cores within each domain and fixed to 8 per rank throughout all experiments. Hyperthreading was disabled and MKL routines were accessed via Intel oneAPI (v2025.0.1).
Note
To further evaluate the scalability of DistPCA across different computing environments, additional experiments were conducted on ATHENA, a CPU server with two nodes, each equipped with a dual-socket Intel Xeon Gold 6430 CPU (32 cores, 2.1 GHz) and 126 GB of RAM. Unless otherwise specified, the reported results were obtained on the ARIS supercomputer.
DistPCA demonstrates near-linear scalability, achieving speedups of up to 58.2× and over 98% reduction in wall-clock time, while maintaining parallel efficiency above 82% across all evaluated scenarios. As shown in the figures, the SGDP and HGDP datasets are omitted, as they complete in under 5 seconds even with 8 MPI ranks.
Figure 1: Runtime performance on the ARIS supercomputer
Figure 2: Strong scaling speedup (left) and parallel efficiency (right)
These performance gains are achieved while preserving the accuracy of the recovered PCs, as illustrated in the following figures.
Figure 3
Left: Entry-wise relative error of the 10 leading eigenvectors computed by DistPCA for the 1000 Genomes dataset, compared to the eigenvectors returned by the full-rank SVD
Right: Projection of the samples of the 1000 Genomes dataset on the top two left singular vectors, as computed by DistPCA. Samples are grouped into five populations: AFR, AMR, EAS, EUR, and SAS
As observed from the following table, DistPCA consistently outperforms PCAone [1, 2], the current state-of-the-art solution for large-scale genomic PCA, across all datasets.
| Dataset | PCAone | DistPCA | Speedup | Reduction % |
|---|---|---|---|---|
| 1000 Genomes | 173s | 47s | 3.68x | 72.8% |
| 50K Genomes | 9.1h | 7.8h | 1.17x | 14.1% |
| 500K Genomes | 12.1h | 2.3h | 5.26x | 78.5% |
| 1M Genomes | 7.9h | 2.6h | 3.04x | 67.9% |
Important
PCAone was employed with the window-based Randomized SVD (RSVD) method proposed in [1] to compute the 20 leading PCs, while the same stopping criterion was used for both frameworks (--memory argument for each dataset to ensure that its window-based RSVD method processed blocks of 100 SNPs. All other PCAone parameters were set to their default values.
All precomputed results from the conducted experiments are available here. To regenerate the figures directly from these outputs, run:
cd scripts/plots/
python3 -m venv venv
source venv/bin/activate
pip install -r requirements.txt
python3 runtime.py # Runtime performance (Figure 3 in the paper)
python3 speedup.py # Strong scaling speedup (Figure 4 in the paper)
python3 par_efficiency.py # Parallel efficiency
python3 rel_error.py # Entry-wise relative error of eigenvectors (Figure 5 in the paper)
python3 pop_structure.py # Population structure (PC1 vs PC2) (Figure 6 in the paper)To reproduce the reported runtime results from scratch, first follow the Datasets Section to download and preprocess the real-world datasets and generate the synthetic ones. Once ready, make sure that the .bed, .bim, and .fam files for each dataset are stored under scripts/experiments/, and run:
# Move dataset files to scripts/experiments/ if not already there
mv <dataset>.bed <dataset>.bim <dataset>.fam scripts/experiments/
cd scripts/experiments/
bash run_1000G.sh
bash run_50K.sh
bash run_500K.sh
bash run_1M.shNote
After executing the above scripts, the wall-clock time results for each worker configuration will be stored in the docs/results/runtime/ directory in separate .txt files, one per dataset.
By running the corresponding Python plot scripts located under scripts/plots/ (runtime.py, speedup.py), Figures 3 and 4 of the paper can be generated.
Important
A key parameter of the proposed hybrid multi-level parallelism scheme is the block size used to partition the dataset, as this determines the number of SNPs processed by each MPI rank at a time and, thus, the parallel I/O performance. The optimal size depends on several factors, including dataset size, RAM, I/O bandwidth, Last-Level Cache (LLC), and the number of compute nodes. In practice, the block size parameter should be selected such that the I/O workload generated by the workers matches the storage system’s capabilities (e.g., bandwidth and latency) and maximizes the LLC hit ratio.
Warning
Experiments were conducted on the ARIS supercomputer using four thin compute nodes. Wall-clock time results may exhibit slight variations depending on cluster infrastructure, node availability, and storage system.
To reproduce the reported accuracy results from scratch, after downloading the 1000 Genomes dataset and moving it to scripts/experiments/ (see Section Datasets), run:
cd scripts/experiments/
bash run_1000G_accuracy.shNote
After executing the above script, the eigenvalues and corresponding eigenvectors computed by DistPCA will be stored in the docs/results/accuracy/ directory in separate .txt files. For reference, eigenvalues and eigenvectors computed via full SVD using LAPACKE are also stored in the same directory in a separate file.
By running the corresponding Python plotting scripts located under scripts/plots/ (rel_error.py, pop_structure.py), Figures 5 and 6 of the paper can be generated.
Caution
The execution of run_1000G_accuracy.sh includes the in-core computation of PCs via full SVD, which requires the entire dataset to be loaded into main memory in uncompressed form. For computing the PCs using LAPACKE's dgesvd, at least 105 GB of available RAM is required.
To reproduce the reported results for the PCAone framework (Table 3 in the paper), first run the setup.sh script from scripts/experiments/PCAone/:
cd scripts/experiments/PCAone/
bash setup.shThe script clones the PCAone repository and builds the PCAone framework using Intel oneAPI.
Then, use the run.sh script, located in scripts/experiments/PCAone/, to run PCAone and compute the leading PCs for the 1000, 50K, 500K, 1M Genomes datasets:
cd scripts/experiments/PCAone/
# 1000 Genomes dataset
bash run.sh 1000
# 50K Genomes dataset
bash run.sh 50K
# 500K Genomes dataset
bash run.sh 500K
# 1M Genomes dataset
bash run.sh 1MTip
Upon completion for each dataset, PCAone reports the elapsed runtime in seconds to stdout.
Note
The script is configured to run PCAone with the parameters described in the Comparison with PCAone subsection. Specifically, PCAone uses the window-based RSVD method proposed in [1] to compute the 20 leading PCs, with
Important
When run.sh is executed without a dataset argument, PCAone automatically runs across all four datasets.
DistPCA/
├── docs/
│ ├── figures/ # Generated figures for the paper
│ └── results/ # Precomputed experimental results
│ ├── runtime/ # Runtime performance results (Figures 3 and 4 of the paper)
│ └── accuracy/ # Evaluation results for computed PCs (Figures 5 and 6 of the paper)
│
├── scripts/
│ ├── data/ # Scripts for downloading and preprocessing real-world datasets
│ ├── plots/ # Scripts to reproduce all figures
│ └── experiments/ # Scripts for running all experiments
│ └── PCAone/ # Scripts for installing, building, and reproducing PCAone experiments
│
├── src/ # Core implementation of DistPCA
├── example/ # Toy dataset for testing and demonstration
│
├── Makefile
- Improve the API documentation
- Provide a Python API
- Support additional genetic data formats (e.g., PLINK2 binary fileset)
- Integrate multi-GPU support
- Support alternative out-of-core methods for PCs approximation
If you find DistPCA useful for your research, please cite:
@article{mermigkis2026distpca,
title = {DistPCA: Tera-Scale Genomic PCA via Out-of-Core Distributed Parallelism},
author = {Mermigkis, Georgios and Sofotasios, Argiris and Kontopoulou, Eugenia-Maria and Gallopoulos, Efstratios and Hadjidoukas, Panagiotis},
journal = {bioRxiv},
year = {2026},
doi = {10.64898/2026.05.15.725487},
url = {https://www.biorxiv.org/content/10.64898/2026.05.15.725487v1}
}This work was supported by computational time granted from the National Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National HPC facility - ARIS - under project ID pa260203distpca.