Prepare aligned UKB phenotypes, plan REGENIE and PLINK2 commands, and run an ICD-10 PheWAS inside RAP.
Scope and privacy boundary
This module prepares the phenotype-side inputs for an end-to-end UK Biobank GWAS followed by PheWAS. It does not bundle UK Biobank source data, cache participant-level records, download genetic data, or submit a DNAnexus job. Participant-level files are only written after an explicit call to ukb_write_gwas_phewas_phenotype(), which defaults to requiring RAP and a path below /mnt/project.
The interface is additive: the existing disease, survival, extraction, and analysis functions are unchanged. Their outputs can be passed directly to the new builder.
The literature screen was deliberately used for tool prioritization rather than as proof that one method is universally superior. A snapshot on 2026-07-29 found the following OpenAlex citation counts for the original method papers: PLINK 2, 14,183; BOLT-LMM, 1,808; SAIGE, 1,588; REGENIE, 1,582; the PheWAS R package, 532; PHESANT, 251; and DeepPheWAS, 23. A separate rough Europe PMC title/abstract/full-text term search combining each tool name with “UK Biobank” returned 4,839, 917, 671, 669, 1,981, 166, and 18 records, respectively. These are discovery statistics, not validated counts of analyses that actually used each program; ambiguous names such as “PheWAS” especially inflate term-search results. REGENIE was selected over BOLT-LMM for the first implementation because it is the current official RAP example path, supports binary and quantitative multi-trait analysis, and has maintained upstream guidance for UKB-scale runs. SAIGE remains important when case-control imbalance is the dominant concern.
Use ukb_gwas_phewas_tools() to inspect the same support decisions in R.
1. Build the disease phenotype with existing UKBAnalytica functions
The example below studies incident hypertension. Prevalent cases are encoded as missing by build_survival_dataset(), so they are not silently used as controls. If the scientific question is lifetime/ever hypertension instead, construct that explicit 0/1 definition before calling the new builder.
library(UKBAnalytica)# raw_dt is an approved UKB table already available inside the RAP project.hypertension_def <-get_predefined_diseases()["Hypertension"]analysis_dt <-build_survival_dataset(dt = raw_dt,disease_definitions = hypertension_def,prevalent_sources =c("ICD10", "ICD9", "Self-report", "Death"),outcome_sources =c("ICD10", "ICD9", "Death"),primary_disease ="Hypertension",baseline_col ="p53_i0",censor_date =as.Date("2023-10-31"),show_flow =TRUE)icd10_long <-parse_icd10_diagnoses(raw_dt)
outcome_status is the GWAS trait in this single-disease example: 0 is an eligible non-case, 1 is an incident case, and NA is not analysed for that phenotype.
2. Create aligned GWAS and PheWAS inputs
pc_cols <-intersect(paste0("p22009_a", 1:10), names(analysis_dt))gwas_covariates <-intersect(c("p31", "p21022", "p22000", pc_cols),names(analysis_dt))# Read IID values from the actual genotype .sample/.psam file used in the run.# Supplying them here makes the intersection explicit and auditable.imputed_sample_ids <- genotype_sample_table$IIDphenotype <-build_gwas_phewas_phenotype(data = analysis_dt,phenotype_cols ="outcome_status",diagnoses = icd10_long,covariates = gwas_covariates,categorical_covariates =intersect(c("p31", "p22000"), gwas_covariates),id_col ="eid",diagnosis_id_col ="eid",diagnosis_code_col ="icd10_code",diagnosis_date_col ="diag_date",genotype_ids = imputed_sample_ids,trait_types =c(outcome_status ="binary"),covariate_missing ="complete_case",sex_col ="p31",sample_qc =ukb_gwas_qc_profile("none"))phenotypephenotype$trait_summaryphenotype$qc_flowphenotype$covariate_missing
The default QC profile does not impose an ancestry restriction. This is intentional: ancestry definition, multi-ancestry analysis, relatedness policy, and transferability checks belong in the study protocol. The following profile reproduces the four sample filters in the DNAnexus tutorial and is available only when it is scientifically appropriate:
The planner generates direct REGENIE commands for an existing Swiss Army Knife or controlled RAP execution environment. It does not submit a job.
regenie_plan <-ukb_plan_regenie(files = files,bed_prefix ="/mnt/project/Bulk/Genotype/array/ukb_array_qc",step1_extract ="/mnt/project/Analysis/GWAS_PheWAS/qc/array_step1.snplist",step2_bgen =sprintf("/mnt/project/Bulk/Imputation/UKB imputation from genotype/ukb_c%d_b0_v3.bgen",1:22 ),step2_sample ="/mnt/project/Bulk/Imputation/UKB imputation from genotype/ukb22828_c1_b0_v3_s487395.sample",step2_extract ="/mnt/project/Analysis/GWAS_PheWAS/qc/imputed_qc.snplist",output_prefix ="/mnt/project/Analysis/GWAS_PheWAS/results/hypertension",threads =16,lowmem =TRUE,firth =TRUE,step1_args =c("--cv", "10", "--keep-l0"),step2_args =c("--minMAC", "20", "--minINFO", "0.8", "--gz"))print(regenie_plan)# Dry run: returns the plan and performs no computation.ukb_run_regenie(regenie_plan)# Explicit execution only inside an environment where the reviewed REGENIE# binary and all inputs are already available.# ukb_run_regenie(regenie_plan, execute = TRUE, require_rap = TRUE)
step1_args and step2_args preserve additional token boundaries while the structured wrapper continues to own inputs, phenotype selection, trait type, prediction files, core association settings, and output paths. Passing a managed flag such as --out through these arguments is rejected instead of silently creating an ambiguous command.
Access to the rest of REGENIE
For survival analysis, gene-based testing, conditional analysis, interactions, LD calculation, split level-0 jobs, or a newer official option that does not fit the structured wrapper, use the raw-token planner:
The generic planner validates that every entry is a non-empty token and quotes tokens before execution. It deliberately does not maintain a second copy of REGENIE’s evolving option table or validate cross-option semantics. Check the official REGENIE option reference and record regenie --version; an older binary can reject options such as --t2e, which was added in REGENIE 4.0.
For the managed DNAnexus REGENIE app, use the generated phenotype, covariate, keep, and variant-list files as versioned app inputs. The direct runner is not a replacement for DNAnexus job orchestration or project billing controls.
4. Export lead-variant dosages with PLINK2
After GWAS QC, significance testing, and protocol-defined LD clumping, put the selected variant IDs in a one-column file. The PLINK2 planner follows the official RAP example’s --bgen ... ref-first --extract ... --export A path.
convert_gwas_datatype() provides a validated high-level interface for common genotype format conversions. It follows the same safety contract: dry-run by default, explicit RAP execution, and no automatic software installation.
Format conversion is not automatically lossless. PLINK2 collapses raw GEN/BGEN genotype probabilities to dosages on import, BED primarily stores hard calls, and VCF dosage output requires an explicit vcf_dosage_export mode. Prefer PGEN when reference alleles, phase, dosage, or multiallelic variants need to be retained, and verify the output before deleting or retiring a source fileset.
Access to the rest of PLINK2
The dedicated wrappers expose the parameters needed for their specific tasks. They do not try to duplicate PLINK2’s complete, frequently updated flag surface as hundreds of R formal arguments.
Additional compatible filters and resource options can be passed as tokens:
This raw-argument route can express current PLINK2 flags, but UKBAnalytica cannot validate every flag combination. The installed PLINK2 version remains the authority and will reject unsupported or incompatible arguments. Review the complete PLINK2 flag index and record plink2 --version with the analysis provenance.
5. Run ICD-10 PheWAS
The optional upstream PheWAS dependency is not installed automatically:
# Run once in a reviewed RAP R environment if PheWAS is not already present.# remotes::install_github("PheWAS/PheWAS")
Read and merge the .raw files produced by the PLINK2 plan, retaining IID and dosage columns:
ukb_run_phewas() uses the upstream ICD-10 phecode map, preserves the full eligible population when creating controls, applies phecode exclusion rules, and can apply sex-specific restrictions when sex_col was supplied to the builder.
Python and basilisk
The initial supported path does not need Python: REGENIE and PLINK2 are external binaries already commonly run through RAP environments, and PheWAS is an R package. Therefore this module does not create a hidden Conda environment. If a later backend has a genuine Python API and a stable pinned dependency set, basilisk can be added for that backend without changing this interface.
Interpretation safeguards
Freeze phenotype definitions, covariates, ancestry strategy, QC thresholds, genetic build, and primary/secondary analyses before the final scan.
Do not interpret discovery counts as independent replication.
Harmonize variant IDs and alleles before clumping or dosage extraction.
Use association methods appropriate to case-control imbalance and rare variants; consider the planned SAIGE backend when imbalance is severe.
Report all sample attrition using qc_flow, covariate missingness, and the DNAnexus job/version provenance.