Compute haplotype components

MODULE STATUS

Developer: Developer's team

Scripts status: Being developed; Please don't run this module for now

Prerequisite scripts: Script 0, 01a, 01b, 01c, 01d, and 01e

Data upload method: automated upload to sftp bristol

Issue report to: Haotian Tang haotian.tang@bristol.ac.uk

Overview: Compute haplotype components (HCs), perform and compare GWASs on CpGs (unadjusted, PC-adjusted, and HC-adjusted).

Genetic principal components analysis (PCA) transforms high-dimensional genetic data into a smaller set of uncorrelated PCs that capture most of the variance. Top PCs typically reflect global ancestry, separating major ancestral groups, and later PCs capture multiple local genomic features (Atkinson et al. (2021) and Grinde et al. (2024)).

Haplotype components (HCs) capture the local, phased structure of genetic variation across the genome by using phased genetic data. This allows them to model subtle ancestry differences that genetic PCs often miss, particularly in cohorts with complex or admixed backgrounds. Theoretically, HCs offer improved control for population stratification and can reduce confounding more effectively than standard genetic PCs in downstream association and epigenetic analyses. Recent work by Yang and Lawson (2025) showed that compared with PC-corrected GWAS, HC-corrected one identifies more genome-wide significant association signals for birthplace and lifestyle-related phenotypes. HCs capture more nuanced genetic variation than PCs.

In Module 1g, we use phased VCF files located in vcf_dir (specified in your config file) to compute HCs through PBWTpaint. We select several CpG sites from EPISTRUCTURE, a set of methylation markers informative for ancestry inference. For each CpG site, three GWAS models are run: no correction, PC-corrected, and HC-corrected. The pipeline then produces Manhattan and QQ plots for each model, as well as pairwise scatter plots of the effect sizes for variants in the cis-region of each CpG site, comparing PC-corrected vs. no correction, HC-corrected vs. no correction, and HC-corrected vs. PC-corrected results.

Before you start, please ensure that ggExtra and gridExtra packages are both available in your R environment.

To run the script, use the following command:

./01g_HCs.sh -c /path/to/config/file

Module 01g is designed to run in chunks, making troubleshooting easier. Processing VCF files can take a long time (e.g., hc chunk crashes due to some error, you do not have to re-run previously successful VCF chunk).

To run the script by chunks:

# 1) VCF preprocessing and liftover
./01g_HCs.sh -c /path/to/config vcf

# 2) Haplotype component (HC) generation
./01g_HCs.sh -c /path/to/config hc

# 3) Positive-control GWAS and plots
./01g_HCs.sh -c /path/to/config gwas

Important

Please check Please check the signals from Manhattan plots, genomic inflation from QQ plots, and scatter plots in results/01/1g.

Example scatter plots (in order of hc-correction vs no-correction, pc vs no, and pc vs hs):