
Bio Workflows Causal Genomics Pipeline
- 1 installs
- 1.1k repo stars
- Updated July 25, 2026
- gptomics/bioskills
Run a post-GWAS causal inference pipeline triangulating Mendelian randomization, colocalization, fine-mapping, TWAS, and drug-target MR from summary statistics.
About
Orchestrates an end-to-end post-GWAS causal inference workflow spanning heritability partitioning, MR with CHP-aware sensitivity, colocalization, SuSiE/FOCUS fine-mapping, mediation, TWAS, and effector-gene prioritization. A developer uses it when triangulating causal evidence, prioritizing tissues, or nominating drug targets from GWAS summary statistics.
- Multi-method MR triangulation with STROBE-MR compliance
- Extensive QC checkpoints across each inference stage
Bio Workflows Causal Genomics Pipeline by the numbers
- 1 all-time installs (skills.sh)
- Ranked #1,803 of 2,064 Data Science & ML skills by installs in the Skillselion catalog
- Data as of Aug 5, 2026 (Skillselion catalog sync)
npx skills add https://github.com/gptomics/bioskills --skill bio-workflows-causal-genomics-pipelineAdd your badge
Show developers this skill is listed on Skillselion. Paste this into your README.
| Installs | 1 |
|---|---|
| repo stars | ★ 1.1k |
| Last updated | July 25, 2026 |
| Repository | gptomics/bioskills ↗ |
What it does
Run a post-GWAS causal inference pipeline triangulating Mendelian randomization, colocalization, fine-mapping, TWAS, and drug-target MR from summary statistics.
Files
Version Compatibility
Reference examples tested with: TwoSampleMR 0.5+, MR-PRESSO 1.0+, coloc 5.2+, susieR 0.12+, MendelianRandomization 0.9+, ldsc 1.0.1 (python3 fork), MetaXcan 0.7+, pyfocus 0.6+, MAGMA 1.10+, MRlap 0.0.3+, cause 1.2+, lhcMR 0.0.0.9000+, HDL 1.4+, LAVA 0.1+, GenomicSEM 0.0.5+.
Before using code patterns, verify installed versions match. If versions differ:
- R:
packageVersion('<pkg>')then?function_nameto verify parameters - Python:
pip show <pkg>thenhelp(module.function)to check signatures - CLI:
<tool> --versionthen<tool> --helpto confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
Causal Genomics Pipeline
"Run post-GWAS causal inference from summary statistics" -> Orchestrate heritability partitioning and tissue prioritization, genetic-correlation diagnostics, instrument selection, Mendelian randomization with CHP-aware sensitivity, colocalization, fine-mapping (SuSiE / FOCUS), mediation, TWAS triangulation, cis-pQTL drug-target MR, effector-gene prioritization, and (optionally) GenomicSEM common-factor GWAS to triangulate causal evidence and nominate publication-grade causal exposures, genes, and mechanisms.
Pipeline Overview
GWAS Summary Statistics (exposure + outcome)
|
v
[0. Pre-flight: h2 + tissue prioritization + rg diagnostic]
LDSC / S-LDSC baseline-LD / Finucane 2018 cell-type
Cross-trait LDSC / HDL / LAVA --> if abs(rg) > 0.3 then CHP-aware MR required
|
v
[1. Instrument Selection] -----> LD clumping, F-stat filtering, Steiger pre-filter
|
v
[2. Mendelian Randomization] --> IVW, MR-Egger, Weighted Median/Mode, MR-RAPS
|
+--> [3. Sensitivity] -------> MR-PRESSO, Egger intercept (Isq), leave-one-out, Steiger
|
+--> [3b. CHP-aware MR] -----> CAUSE (delta_elpd), LHC-MR posterior (if rg > 0.3)
|
v
[4. Colocalization] -----------> coloc.abf / coloc.susie / HyPrColoc / SMR-HEIDI
|
v
[5. Fine-Mapping] -------------> SuSiE rss + estimate_s_rss / FINEMAP-inf / PolyFun
|
v
[6. Mediation Analysis] -------> Network MR / MVMR / CMAverse 4-way
|
+--> [7. TWAS triangulation] -> FUSION / S-PrediXcan / FOCUS PIP >= 0.8
|
+--> [8. Cis-pQTL drug-target MR] -> UKB-PPP / deCODE + cross-platform replication
|
v
[9. Effector-gene prioritization] -> Open Targets L2G + PoPS + cS2G + coloc + TWAS (>= 3 of 6 evidence)
|
v
[10. (optional) GenomicSEM common-factor GWAS] -> factor model + Q_SNP
|
v
Triangulated causal-evidence summary across methodsStep 0: Pre-flight - Heritability, Tissue Prioritization, Genetic Correlation
ldsc.py --h2 trait.sumstats.gz --ref-ld-chr eur_w_ld_chr/ --w-ld-chr eur_w_ld_chr/ --out trait.h2
ldsc.py --h2 trait.sumstats.gz --ref-ld-chr baselineLD.,cts_v8. --w-ld-chr weights. --ref-ld-chr-cts Multi_tissue_gene_expr.ldcts --out trait.cts
ldsc.py --rg trait1.sumstats.gz,trait2.sumstats.gz --ref-ld-chr eur_w_ld_chr/ --w-ld-chr eur_w_ld_chr/ --out rgGoal: Confirm heritable signal, pick the right tissue for TWAS / V2G, and detect shared heritable confounding that mandates CHP-aware MR (CAUSE / LHC-MR).
Reconciliation: S-LDSC mean chi-squared > 1.02 with h2 SE < 0.02 and intercept ratio < 0.3 is required. Cell-type prioritization with coefficient_p < 0.05 / N_tissues nominates the tissue for downstream TWAS weights and ABC enhancer-gene priors. If cross-trait LDSC abs(rg) > 0.3 (and HDL sample-overlap < 5%), Step 3b becomes mandatory. See causal-genomics/heritability-partitioning and causal-genomics/genetic-correlation.
Step 1: Instrument Selection
library(TwoSampleMR)
exposure_dat <- read_exposure_data(filename = 'exposure_gwas.tsv', sep = '\t',
snp_col = 'SNP', beta_col = 'BETA', se_col = 'SE',
effect_allele_col = 'A1', other_allele_col = 'A2',
eaf_col = 'EAF', pval_col = 'P')
exposure_dat <- subset(exposure_dat, pval.exposure < 5e-8)
exposure_dat <- clump_data(exposure_dat, clump_r2 = 0.001, clump_kb = 10000)
exposure_dat$F_stat <- (exposure_dat$beta.exposure / exposure_dat$se.exposure)^2
exposure_dat <- subset(exposure_dat, F_stat >= 10)
exposure_dat <- subset(exposure_dat, !(eaf.exposure > 0.42 & eaf.exposure < 0.58 & substr(effect_allele.exposure,1,1) %in% c('A','T') & substr(other_allele.exposure,1,1) %in% c('A','T')))For cis-MR (drug target) use clump_r2 = 0.1 within +/- 500 kb of the gene. Use 5e-9 if M > 5M variants tested.
Step 2: Mendelian Randomization
outcome_dat <- read_outcome_data(filename = 'outcome_gwas.tsv', sep = '\t',
snp_col = 'SNP', beta_col = 'BETA', se_col = 'SE',
effect_allele_col = 'A1', other_allele_col = 'A2',
eaf_col = 'EAF', pval_col = 'P')
dat <- harmonise_data(exposure_dat, outcome_dat)
mr_results <- mr(dat, method_list = c('mr_ivw', 'mr_egger_regression',
'mr_weighted_median', 'mr_weighted_mode'))Concordance across IVW, Egger, weighted median, and weighted mode is the headline causal claim. See causal-genomics/mendelian-randomization.
Step 3: Sensitivity Analysis
library(MRPRESSO)
presso <- mr_presso(BetaOutcome = 'beta.outcome', BetaExposure = 'beta.exposure',
SdOutcome = 'se.outcome', SdExposure = 'se.exposure',
OUTLIERtest = TRUE, DISTORTIONtest = TRUE, data = dat,
NbDistribution = 5000, SignifThreshold = 0.05)
egger_int <- mr_pleiotropy_test(dat)
isq <- Isq(dat$beta.exposure, dat$se.exposure)
het <- mr_heterogeneity(dat)
loo <- mr_leaveoneout(dat)
steiger <- directionality_test(dat)Isq >= 0.9 is required for the MR-Egger NOME assumption; below that, run SIMEX correction or drop Egger. See causal-genomics/pleiotropy-detection.
Step 3b: CHP-aware MR (when rg > 0.3 or shared confounder suspected)
library(cause)
library(lhcMR)
cause_fit <- cause(X = cause_dat, variants = top_vars, param_ests = params)
elpd <- summary(cause_fit)$elpd
lhc_fit <- lhc_mr(trait.df = sumstats_df, LD.filepath = ld_path, rho.filepath = rho_path,
nStep = 2, SP_single = 3, SP_pair = 50, run_local = FALSE)CAUSE reports delta_elpd of sharing-vs-causal model; z > 1.96 favors true causation over CHP. LHC-MR jointly estimates causal effect and confounder effect via likelihood; the 95% credible interval excluding zero is the causal-effect verdict. Required when LDSC abs(rg) > 0.3. See causal-genomics/pleiotropy-detection.
Step 4: Colocalization
library(coloc)
d1 <- list(beta = exposure_locus$BETA, varbeta = exposure_locus$SE^2,
snp = exposure_locus$SNP, position = exposure_locus$BP,
type = 'quant', N = exposure_n, MAF = exposure_locus$EAF)
d2 <- list(beta = outcome_locus$BETA, varbeta = outcome_locus$SE^2,
snp = outcome_locus$SNP, position = outcome_locus$BP,
type = 'cc', N = outcome_n, s = case_fraction, MAF = outcome_locus$EAF)
result <- coloc.abf(d1, d2, p1 = 1e-4, p2 = 1e-4, p12 = 1e-5)
sens <- coloc::sensitivity(result, rule = 'H4 > 0.7')PP.H4 >= 0.7 for triangulation, >= 0.8 for publication, >= 0.95 for industry-grade target packages. Always sweep p12 over 1e-6 to 5e-5; conclusions must be stable across the sweep. For allelic heterogeneity use coloc.susie. See causal-genomics/colocalization-analysis.
Step 5: Fine-Mapping with SuSiE
library(susieR)
R <- as.matrix(read.csv('ld_matrix.csv', row.names = 1))
diag_s <- estimate_s_rss(z = locus_stats$BETA / locus_stats$SE, R = R, n = sample_size)
fitted <- susie_rss(bhat = locus_stats$BETA, shat = locus_stats$SE,
R = R, n = sample_size, L = 10, coverage = 0.95, min_abs_corr = 0.5)
cs <- fitted$sets$csIf estimate_s_rss lambda > 0.05 the external LD is mismatched; rerun with in-sample LD or use SuSiE-inf / FINEMAP-inf. min_abs_corr >= 0.5 (r-squared >= 0.25) is the purity threshold for retaining a credible set. For HLA use L = 20-30. See causal-genomics/fine-mapping.
Step 6: Mediation Analysis
library(TwoSampleMR)
mv_exposures <- mv_extract_exposures(c('ieu-a-2', 'ieu-a-1089'))
mv_outcome <- mv_extract_outcome(mv_exposures$SNP, 'ieu-a-7')
mvdat <- mv_harmonise_data(mv_exposures, mv_outcome)
mvmr_result <- mv_multiple(mvdat)Indirect effect = total - direct. For molecular mediators (expression, methylation, protein), prefer two-step MR with cis-instruments at the mediator. Run Imai sensitivity (rho_crit) or mediational E-value > 2. See causal-genomics/mediation-analysis.
Step 7: TWAS Triangulation
python MetaXcan/SPrediXcan.py --model_db_path mashr_Whole_Blood.db \
--covariance mashr_Whole_Blood.txt.gz --gwas_file gwas.txt.gz \
--snp_column SNP --effect_allele_column A1 --non_effect_allele_column A2 \
--beta_column BETA --se_column SE --pvalue_column P \
--output_file twas.csv
focus finemap gwas.sumstats.gz LD.bcor mashr.db --locations EUR.locations.bed --out twas.focusGoal: Nominate gene-level causal hits and prune LD-induced TWAS false positives.
Tissue is picked from Step 0 stratified LDSC; Bonferroni at 0.05 / N_tissues. FOCUS PIP >= 0.8 retains a single candidate causal gene per region; without FOCUS, co-regulated TWAS hits cannot be distinguished. Cross-reference TWAS hits with coloc.susie PP.H4 and cis-eQTL MR for triangulation. See causal-genomics/transcriptome-wide-association.
Step 8: Cis-pQTL Drug-Target MR
library(TwoSampleMR)
pqtl_dat <- extract_instruments(outcomes = 'prot-a-XXX', p1 = 5e-8, clump = TRUE,
r2 = 0.1, kb = 1000)
pqtl_dat <- subset(pqtl_dat, chr == target_chr & pos > target_tss - 500000 & pos < target_tss + 500000)
out_dat <- extract_outcome_data(snps = pqtl_dat$SNP, outcomes = c('ieu-a-7', 'ieu-b-31'))
dat <- harmonise_data(pqtl_dat, out_dat)
mr_results <- mr(dat, method_list = c('mr_ivw', 'mr_wald_ratio'))Goal: Mimic pharmacological inhibition of a drug target via cis-pQTL and triangulate with coloc.
Cross-platform replication on Olink (UKB-PPP) and SomaScan (deCODE) is mandatory; the two platforms are concordant for ~60% of proteins and discordant calls are platform artifacts. Run pheWAS for on-target adverse effects, PAV-excluded sensitivity, and coloc.susie PP.H4 >= 0.8 at the cis-pQTL locus. See causal-genomics/proteome-mr-drug-target.
Step 9: Effector-Gene Prioritization
magma --bfile g1000_eur --pval gwas.tsv N=N --gene-annot genes.annot --out trait
pops_features.py --features pops.features.txt --gene_annot genes.txt --out trait.features
pops.py --features trait.features --gene_results trait.genes.raw --out trait.popsGoal: Map each fine-mapped credible set to a candidate effector gene by integrating six evidence streams.
Integrate: (1) Open Targets L2G (Mountjoy 2021), (2) PoPS similarity score (Weeks 2023), (3) cS2G combined SNP-to-gene (Gazal 2022), (4) coloc.susie PP.H4 with eQTL/pQTL, (5) FOCUS TWAS PIP, (6) ABC / ENCODE-rE2G enhancer-gene linking. Require >= 3 of 6 concordant evidence streams for high-confidence claim. L2G and PoPS disagree by design (different feature regimes); report both. See causal-genomics/effector-gene-prioritization.
Step 10 (optional): GenomicSEM Common-Factor GWAS
library(GenomicSEM)
ldsc_output <- ldsc(traits = c('t1.sumstats.gz','t2.sumstats.gz','t3.sumstats.gz'),
sample.prev = c(NA, NA, NA), population.prev = c(NA, NA, NA),
ld = ld_path, wld = ld_path, trait.names = c('t1','t2','t3'))
model <- 'F1 =~ NA*t1 + t2 + t3\nF1 ~~ 1*F1'
fit <- usermodel(ldsc_output, model = model, estimation = 'DWLS')
factor_gwas <- userGWAS(covstruc = ldsc_output, SNPs = sumstats_combined,
model = paste0(model, '\nF1 ~ SNP\nt1 + t2 + t3 ~ 0*SNP'),
estimation = 'DWLS', sub = c('F1~SNP'))Heywood cases (negative residual variance) require fixing residuals positive or dropping the indicator. Verify CFI > 0.95 and RMSEA < 0.06. Q_SNP p > 0.05 confirms factor-level (not trait-specific) signal. See causal-genomics/genomic-sem.
Parameter Recommendations
| Step | Parameter | Recommendation |
|---|---|---|
| Instruments | p-value | 5e-8 standard; 5e-9 if M > 5M variants tested |
| Instruments | F-statistic | >= 10 two-sample; >= 20 one-sample |
| Instruments | clump_r2 | 0.001 polygenic; 0.1 cis-MR |
| Instruments | clump_kb | 10000 (10 Mb) |
| Coloc | p12 prior | 1e-5 standard; 5e-6 conservative; 1e-6 trans-eQTL |
| Coloc | PP.H4 | >= 0.7 triangulation; >= 0.8 publication; >= 0.95 industry |
| SuSiE | L | 10 default; 20-30 HLA |
| SuSiE | coverage | 0.95 standard; 0.9 if N < 1000 |
| SuSiE | min_abs_corr | 0.5 default (r-squared >= 0.25) |
| MR-PRESSO | NbDistribution | 1000 exploratory; >= 5000 publication; >= 10000 stringent |
| TWAS | FOCUS PIP | >= 0.8 candidate causal gene |
| TWAS | tissue Bonferroni | 0.05 / N_tissues (~2.5e-4 for 200 tissues) |
| LDSC | mean chi-squared | > 1.02 for h2 interpretability |
| LDSC | rg trigger for CHP-MR | abs(rg) > 0.3 |
| LDSC | HDL sample overlap | < 5% |
Troubleshooting
| Issue | Likely Cause | Solution |
|---|---|---|
| No instruments | Underpowered GWAS | Relax p-value to 5e-6 with caution |
| Weak instruments (F < 10) | Small effect SNPs | Drop weak instruments; use better-powered GWAS |
| Inconsistent MR methods | Pleiotropy | Check MR-PRESSO outliers; use weighted median |
| Egger intercept p < 0.05 | Directional pleiotropy | Report Egger estimate; check Isq >= 0.9 |
| PP.H3 > PP.H4 | Different causal variants | Use coloc.susie for allelic heterogeneity |
| No credible sets | LD matrix issues | Check estimate_s_rss lambda; use in-sample LD |
| Steiger reverse direction | Reverse causation | Run bidirectional MR; pre-filter on Steiger |
| MR-PRESSO blind to apparent confounder | Correlated horizontal pleiotropy | Run CAUSE or LHC-MR (see pleiotropy-detection) |
| TWAS hit at gene-dense locus | LD-induced false positive | Run FOCUS fine-mapping; require PIP >= 0.8 |
| Cis-pQTL MR positive, replication fails | Olink/SomaScan platform discordance | Cross-platform replication; PAV-excluded sensitivity |
| L2G + PoPS disagree | Different feature regimes | Report both; require concordance for high-confidence claim |
| Q_SNP heterogeneity in common-factor GWAS | Factor mis-specification | userGWAS with per-trait paths; report Q_pval |
Related Skills
causal-genomics/mendelian-randomization - IVW, Egger, MR-RAPS, MVMR causal-genomics/colocalization-analysis - coloc.abf, coloc.susie, HyPrColoc, SMR-HEIDI causal-genomics/fine-mapping - SuSiE rss, FINEMAP-inf, PolyFun, SuSiEx causal-genomics/pleiotropy-detection - MR-PRESSO, CAUSE, LHC-MR, contamination-mixture causal-genomics/mediation-analysis - Two-step MR, MVMR, CMAverse 4-way, HIMA causal-genomics/transcriptome-wide-association - FUSION, S-PrediXcan, FOCUS, UTMOST causal-genomics/heritability-partitioning - LDSC, S-LDSC, LDAK, HDL, HESS causal-genomics/proteome-mr-drug-target - UKB-PPP, deCODE, cis-pQTL MR, pheWAS causal-genomics/effector-gene-prioritization - L2G, PoPS, cS2G, MAGMA, FLAMES causal-genomics/genetic-correlation - Cross-trait LDSC, HDL, LAVA, Popcorn causal-genomics/genomic-sem - Common-factor GWAS, Q_SNP, MTAG reconciliation population-genetics/association-testing - Upstream GWAS methods atac-seq/enhancer-gene-linking - ABC / ENCODE-rE2G priors for effector-gene step single-cell/preprocessing - scRNA / scATAC tissue priors for stratified LDSC
# Reference: TwoSampleMR 0.5+, coloc 5.2+, susieR 0.12+ | Verify API if version differs
# Complete post-GWAS causal inference pipeline
# Mendelian randomization -> sensitivity analysis -> colocalization -> fine-mapping
library(TwoSampleMR)
library(MRPRESSO)
library(coloc)
library(susieR)
# Configuration
# Example: LDL cholesterol (exposure) -> coronary heart disease (outcome)
# Use GWAS summary statistics from public repositories:
# - IEU OpenGWAS: https://gwas.mrcieu.ac.uk/
# - GWAS Catalog: https://www.ebi.ac.uk/gwas/
EXPOSURE_FILE <- 'exposure_gwas.tsv'
OUTCOME_FILE <- 'outcome_gwas.tsv'
EXPOSURE_N <- 100000
OUTCOME_N <- 50000
CASE_FRACTION <- 0.3
# ============================================================
# Step 1: Instrument Selection
# ============================================================
cat('=== Step 1: Instrument Selection ===\n')
exposure_dat <- read_exposure_data(
filename = EXPOSURE_FILE,
sep = '\t',
snp_col = 'SNP',
beta_col = 'BETA',
se_col = 'SE',
effect_allele_col = 'A1',
other_allele_col = 'A2',
eaf_col = 'EAF',
pval_col = 'P'
)
# p < 5e-8: Standard GWAS significance; use 5e-6 for underpowered exposures
exposure_dat <- subset(exposure_dat, pval.exposure < 5e-8)
cat(sprintf('Genome-wide significant SNPs: %d\n', nrow(exposure_dat)))
# LD clumping for independent instruments
# clump_r2 = 0.001: Strict independence (standard for MR)
# clump_kb = 10000: 10 Mb window
exposure_dat <- clump_data(exposure_dat, clump_r2 = 0.001, clump_kb = 10000)
cat(sprintf('After LD clumping: %d independent instruments\n', nrow(exposure_dat)))
# F-statistic > 10 to avoid weak instrument bias
exposure_dat$F_stat <- (exposure_dat$beta.exposure / exposure_dat$se.exposure)^2
cat(sprintf('F-statistic range: %.1f - %.1f (mean: %.1f)\n',
min(exposure_dat$F_stat), max(exposure_dat$F_stat), mean(exposure_dat$F_stat)))
weak <- sum(exposure_dat$F_stat < 10)
if (weak > 0) {
cat(sprintf('WARNING: Removing %d weak instruments (F < 10)\n', weak))
exposure_dat <- subset(exposure_dat, F_stat >= 10)
}
cat(sprintf('Final instrument count: %d\n', nrow(exposure_dat)))
# ============================================================
# Step 2: Mendelian Randomization
# ============================================================
cat('\n=== Step 2: Mendelian Randomization ===\n')
outcome_dat <- read_outcome_data(
filename = OUTCOME_FILE,
sep = '\t',
snp_col = 'SNP',
beta_col = 'BETA',
se_col = 'SE',
effect_allele_col = 'A1',
other_allele_col = 'A2',
eaf_col = 'EAF',
pval_col = 'P'
)
dat <- harmonise_data(exposure_dat, outcome_dat)
cat(sprintf('Harmonised SNPs: %d\n', nrow(subset(dat, mr_keep == TRUE))))
# Run multiple MR methods for triangulation
# IVW: Primary (assumes balanced pleiotropy)
# Egger: Allows directional pleiotropy (needs >= 3 SNPs)
# Weighted median: Robust to up to 50% invalid instruments
# Weighted mode: Most robust, least power
mr_results <- mr(dat, method_list = c(
'mr_ivw',
'mr_egger_regression',
'mr_weighted_median',
'mr_weighted_mode'
))
cat('\nMR Results:\n')
print(mr_results[, c('method', 'nsnp', 'b', 'se', 'pval')])
# QC: Check consistency across methods
directions <- sign(mr_results$b)
if (length(unique(directions)) == 1) {
cat('OK: All methods agree on direction of effect\n')
} else {
cat('WARNING: Inconsistent effect directions - investigate pleiotropy\n')
}
# ============================================================
# Step 3: Sensitivity Analysis
# ============================================================
cat('\n=== Step 3: Sensitivity Analysis ===\n')
# MR-PRESSO: detect and correct for outlier instruments
# NbDistribution=3000: Standard; use 10000 for publication
presso <- mr_presso(
BetaOutcome = 'beta.outcome',
BetaExposure = 'beta.exposure',
SdOutcome = 'se.outcome',
SdExposure = 'se.exposure',
OUTLIERtest = TRUE,
DISTORTIONtest = TRUE,
data = dat,
NbDistribution = 3000,
SignifThreshold = 0.05
)
presso_global_p <- presso$`MR-PRESSO results`$`Global Test`$Pvalue
cat(sprintf('MR-PRESSO global test p: %.4f', presso_global_p))
if (presso_global_p < 0.05) {
cat(' (outliers detected)\n')
n_outliers <- sum(presso$`MR-PRESSO results`$`Outlier Test`$Pvalue < 0.05, na.rm = TRUE)
cat(sprintf('Outlier instruments: %d\n', n_outliers))
} else {
cat(' (no significant outliers)\n')
}
# MR-Egger intercept for directional pleiotropy
# p > 0.05: No evidence of directional pleiotropy (good)
egger_int <- mr_pleiotropy_test(dat)
cat(sprintf('Egger intercept: %.4f, p: %.4f', egger_int$egger_intercept, egger_int$pval))
if (egger_int$pval < 0.05) {
cat(' (directional pleiotropy present!)\n')
} else {
cat(' (no directional pleiotropy)\n')
}
# Cochran Q for heterogeneity
het <- mr_heterogeneity(dat)
ivw_het <- het[het$method == 'Inverse variance weighted', ]
cat(sprintf('Cochran Q p: %.4f\n', ivw_het$Q_pval))
# Steiger directionality test
steiger <- directionality_test(dat)
cat(sprintf('Steiger correct direction: %s (p: %.2e)\n',
steiger$correct_causal_direction, steiger$steiger_pval))
if (!steiger$correct_causal_direction) {
cat('WARNING: Steiger suggests reverse causation. Consider bidirectional MR.\n')
}
# Leave-one-out analysis
loo <- mr_leaveoneout(dat)
loo_range <- range(loo$b[loo$SNP != 'All'])
cat(sprintf('Leave-one-out estimate range: [%.3f, %.3f]\n', loo_range[1], loo_range[2]))
# ============================================================
# Step 4: Colocalization
# ============================================================
cat('\n=== Step 4: Colocalization ===\n')
# Read full summary statistics for the top locus
# Identify top MR locus from strongest instrument
top_snp <- exposure_dat[which.max(exposure_dat$F_stat), ]
cat(sprintf('Top locus: %s (chr%s)\n', top_snp$SNP, top_snp$chr.exposure))
# Extract +/- 500 kb around top SNP for colocalization
exposure_full <- read.delim(EXPOSURE_FILE)
outcome_full <- read.delim(OUTCOME_FILE)
locus_chr <- top_snp$chr.exposure
locus_pos <- top_snp$pos.exposure
# 500 kb window: Standard for colocalization
window <- 500000
exposure_locus <- subset(exposure_full, CHR == locus_chr & BP >= (locus_pos - window) & BP <= (locus_pos + window))
outcome_locus <- subset(outcome_full, CHR == locus_chr & BP >= (locus_pos - window) & BP <= (locus_pos + window))
common_snps <- intersect(exposure_locus$SNP, outcome_locus$SNP)
exposure_locus <- exposure_locus[match(common_snps, exposure_locus$SNP), ]
outcome_locus <- outcome_locus[match(common_snps, outcome_locus$SNP), ]
cat(sprintf('SNPs in locus: %d\n', length(common_snps)))
d1 <- list(
beta = exposure_locus$BETA, varbeta = exposure_locus$SE^2,
snp = exposure_locus$SNP, position = exposure_locus$BP,
type = 'quant', N = EXPOSURE_N, MAF = exposure_locus$EAF
)
d2 <- list(
beta = outcome_locus$BETA, varbeta = outcome_locus$SE^2,
snp = outcome_locus$SNP, position = outcome_locus$BP,
type = 'cc', N = OUTCOME_N, s = CASE_FRACTION, MAF = outcome_locus$EAF
)
# Priors: p1=1e-4, p2=1e-4, p12=1e-5 (standard)
# Use p12=5e-6 for more conservative colocalization
coloc_result <- coloc.abf(d1, d2, p1 = 1e-4, p2 = 1e-4, p12 = 1e-5)
cat('\nColocalization posterior probabilities:\n')
cat(sprintf(' PP.H0 (neither): %.3f\n', coloc_result$summary['PP.H0.abf']))
cat(sprintf(' PP.H1 (exposure only): %.3f\n', coloc_result$summary['PP.H1.abf']))
cat(sprintf(' PP.H2 (outcome only): %.3f\n', coloc_result$summary['PP.H2.abf']))
cat(sprintf(' PP.H3 (both, different variants): %.3f\n', coloc_result$summary['PP.H3.abf']))
cat(sprintf(' PP.H4 (shared causal variant): %.3f\n', coloc_result$summary['PP.H4.abf']))
# PP.H4 > 0.8: Strong evidence for shared causal variant
# PP.H4 > 0.5: Suggestive
pp4 <- coloc_result$summary['PP.H4.abf']
if (pp4 > 0.8) {
cat('Strong colocalization evidence\n')
} else if (pp4 > 0.5) {
cat('Suggestive colocalization\n')
} else {
cat('Weak colocalization. Traits may have different causal variants at this locus.\n')
}
# ============================================================
# Step 5: Fine-Mapping with SuSiE
# ============================================================
cat('\n=== Step 5: Fine-Mapping ===\n')
# Fine-mapping requires an LD matrix from a matched reference panel
# In practice: compute from 1000 Genomes or UK Biobank
# LD_MATRIX_FILE <- 'ld_matrix.csv'
# R <- as.matrix(read.csv(LD_MATRIX_FILE, row.names = 1))
# If LD matrix is available:
# L=10: Maximum number of causal variants (standard)
# coverage=0.95: 95% credible set
# fitted <- susie_rss(
# bhat = exposure_locus$BETA,
# shat = exposure_locus$SE,
# R = R,
# n = EXPOSURE_N,
# L = 10,
# coverage = 0.95
# )
#
# cs <- fitted$sets$cs
# for (i in seq_along(cs)) {
# snps_in_cs <- exposure_locus$SNP[cs[[i]]]
# pip <- fitted$pip[cs[[i]]]
# cat(sprintf('Credible set %d: %d SNPs, min PIP = %.3f\n', i, length(snps_in_cs), min(pip)))
# }
#
# high_pip <- exposure_locus$SNP[fitted$pip > 0.5]
# cat(sprintf('Causal SNPs (PIP > 0.5): %s\n', paste(high_pip, collapse = ', ')))
cat('Fine-mapping requires LD matrix from matched reference panel.\n')
cat('Compute from 1000 Genomes or in-sample genotypes, then run susie_rss().\n')
# ============================================================
# Summary
# ============================================================
cat('\n=== Causal Evidence Summary ===\n')
ivw <- mr_results[mr_results$method == 'Inverse variance weighted', ]
cat(sprintf('IVW estimate: beta = %.3f (95%% CI: %.3f to %.3f), p = %.2e\n',
ivw$b, ivw$b - 1.96 * ivw$se, ivw$b + 1.96 * ivw$se, ivw$pval))
cat(sprintf('Direction consistency: %s\n', ifelse(length(unique(sign(mr_results$b))) == 1, 'YES', 'NO')))
cat(sprintf('Pleiotropy (Egger intercept p): %.4f\n', egger_int$pval))
cat(sprintf('Outliers (MR-PRESSO p): %.4f\n', presso_global_p))
cat(sprintf('Correct direction (Steiger): %s\n', steiger$correct_causal_direction))
cat(sprintf('Colocalization (PP.H4): %.3f\n', pp4))
cat('\nPipeline complete.\n')
Causal Genomics Pipeline - Usage Guide
Overview
Complete post-GWAS causal inference workflow that triangulates evidence across heritability partitioning, genetic correlation, Mendelian randomization (with CHP-aware sensitivity via CAUSE / LHC-MR), colocalization, fine-mapping (SuSiE / FOCUS), mediation, TWAS, cis-pQTL drug-target MR, effector-gene prioritization (L2G / PoPS / cS2G), and (optionally) GenomicSEM common-factor GWAS. Suitable for nominating causal exposures, picking the right tissue for downstream analyses, mapping a lead SNP to a candidate effector gene, de-risking or nominating a drug target, and producing a STROBE-MR-compliant publication-grade evidence battery from GWAS summary statistics.
Prerequisites
install.packages('remotes')
remotes::install_github('MRCIEU/TwoSampleMR')
remotes::install_github('rondolab/MR-PRESSO')
remotes::install_github('jean997/cause')
remotes::install_github('zhenin/HDL/HDL')
remotes::install_github('josefin-werme/LAVA')
remotes::install_github('jrs95/MRlap')
remotes::install_github('LizaDarrous/lhcMR')
install.packages('coloc')
install.packages('susieR')
install.packages('MendelianRandomization')
install.packages('GenomicSEM', repos = c('https://genomic-sem.github.io/repo', 'https://cran.r-project.org'))git clone https://github.com/bulik/ldsc.git
conda env create -n ldsc -f ldsc/environment.yml
git clone https://github.com/hakyimlab/MetaXcan.git
pip install pyfocus
git clone https://ctg.cncr.nl/software/MAGMA/magma_v1.10.zip
pip install popsInput data:
- GWAS summary statistics for exposure and outcome (TSV with SNP, BETA, SE, A1, A2, EAF, P, N columns)
- Pre-computed LDSC weights and baseline-LD annotation (eur_w_ld_chr, baselineLD)
- TWAS prediction weights (PredictDB MASHR models or FUSION elastic net) matched to the tissue prioritized in Step 0
- LD matrix for fine-mapping loci (in-sample preferred; reference-panel acceptable with estimate_s_rss check)
- Cis-pQTL summary stats for drug-target MR (UKB-PPP, deCODE, Fenland, INTERVAL)
- Open Targets L2G features and PoPS gene-feature table for effector-gene step
Quick Start
Tell your AI agent what you want to do:
- "Run the full post-GWAS causal-inference pipeline on my GWAS summary stats"
- "Estimate heritability, prioritize tissues, then run MR with CHP-aware sensitivity"
- "Triangulate TWAS, colocalization, and cis-eQTL MR to nominate the causal gene"
- "Run cis-pQTL drug-target MR with Olink and SomaScan cross-platform replication"
- "Prioritize an effector gene from a GWAS lead locus using L2G, PoPS, and ABC"
- "Test a common-factor GWAS model across my correlated traits via GenomicSEM"
Example Prompts
Full Pipeline
"I have GWAS summary statistics for LDL cholesterol and coronary heart disease. Run the complete causal-inference pipeline including heritability partitioning, tissue prioritization, MR with CHP-aware sensitivity, coloc, fine-mapping, and effector-gene prioritization."
"Perform MR, sensitivity analysis, and colocalization to test if BMI causally affects type 2 diabetes, with CAUSE if rg > 0.3."
Heritability and Tissue Prioritization
"Estimate SNP heritability with LDSC, partition with the baseline-LD model, and prioritize tissues via Finucane 2018 cell-type LDSC for my schizophrenia GWAS."
"Run stratified LDSC on my GWAS to pick the right GTEx tissue for downstream TWAS."
Genetic Correlation
"Compute cross-trait LDSC rg between my exposure and outcome to decide whether CAUSE or LHC-MR is required."
"Run LAVA local genetic correlation across the genome to find regions of shared signal."
MR and Sensitivity
"Select strong instruments and check F-statistics for my exposure GWAS."
"Run MR-PRESSO with NbDistribution = 10000 to detect outlier instruments."
"Run CAUSE because LDSC rg > 0.3 between my traits and CHP is plausible."
Colocalization and Fine-Mapping
"Colocalize my GWAS locus with whole-blood eQTL using coloc.susie because the region has multiple signals."
"Fine-map my GWAS locus with SuSiE rss, run estimate_s_rss for LD diagnostics, and report 95% credible sets with min_abs_corr >= 0.5."
TWAS Triangulation
"Run S-PrediXcan in whole-blood MASHR models and fine-map TWAS hits with FOCUS at PIP >= 0.8."
"Triangulate TWAS, coloc.susie, and cis-eQTL MR to nominate a single causal gene per locus."
Drug-Target MR
"Run cis-pQTL MR for IL6R using UKB-PPP and replicate in deCODE SomaScan with PAV-excluded sensitivity."
"Screen on-target adverse phenotypes with pheWAS at the cis-pQTL for my drug target."
Effector-Gene Prioritization
"Map this GWAS lead variant to a candidate effector gene using Open Targets L2G, PoPS, cS2G, coloc, and ABC."
"Reconcile L2G and PoPS calls at my locus and report concordance across 6 evidence streams."
Mediation
"Test if CRP mediates the causal effect of BMI on CHD using two-step MR or MVMR with Imai sensitivity."
"Run CMAverse 4-way decomposition for exposure-mediator interaction."
GenomicSEM
"Fit a common-factor GenomicSEM model across my 5 psychiatric GWAS and run a common-factor GWAS with Q_SNP heterogeneity testing."
"Compare GenomicSEM common-factor GWAS against MTAG for the same trait set."
What the Agent Will Do
1. Estimate SNP heritability with LDSC and run partitioned S-LDSC for functional categories 2. Prioritize tissues via Finucane 2018 cell-type-specific LDSC (informs Step 7 TWAS tissue choice) 3. Compute cross-trait LDSC rg; if abs(rg) > 0.3, mandate CHP-aware MR (CAUSE / LHC-MR) in Step 3b 4. Select genome-wide-significant independent instruments with LD clumping, F-statistic, palindromic filters, and Steiger pre-filter 5. Run primary MR across IVW, MR-Egger, weighted median, and weighted mode for concordance 6. Run sensitivity battery: MR-PRESSO, Egger intercept (with Isq >= 0.9 for NOME), leave-one-out, Cochran Q, Steiger directionality 7. If rg > 0.3, run CAUSE delta_elpd and LHC-MR posterior to distinguish causation from CHP 8. Colocalize the locus with coloc.abf or coloc.susie; sweep p12 prior; require PP.H4 >= 0.7 (triangulation) or 0.8 (publication) 9. Fine-map with susie_rss; run estimate_s_rss for LD diagnostics; require min_abs_corr >= 0.5 10. Mediation: two-step MR or MVMR with Imai rho_crit / mediational E-value sensitivity 11. TWAS: run S-PrediXcan in the prioritized tissue, fine-map with FOCUS PIP >= 0.8 12. Cis-pQTL drug-target MR: replicate across Olink and SomaScan platforms; run pheWAS and PAV sensitivity 13. Effector-gene prioritization: integrate L2G, PoPS, cS2G, coloc, TWAS, and ABC; require >= 3 of 6 concordant streams 14. Optional GenomicSEM: fit common-factor model, verify CFI > 0.95 and RMSEA < 0.06, run userGWAS with Q_SNP 15. Triangulate evidence and produce STROBE-MR-compliant summary
Tips
- Always use multiple MR methods; if directions disagree, investigate pleiotropy and run CAUSE / LHC-MR
- F-statistic >= 10 is minimum for two-sample MR; use >= 20 for one-sample MR to avoid weak-instrument bias toward the observational estimate
- Filter palindromic SNPs with MAF near 0.5 (0.42 to 0.58); they cannot be reliably harmonised
- MR-PRESSO global test p > 0.05 means no significant horizontal pleiotropy; below that, report the outlier-corrected estimate, not the raw IVW
- Isq < 0.9 invalidates the MR-Egger NOME assumption; run SIMEX or drop Egger from the headline
- Cross-trait LDSC abs(rg) > 0.3 triggers CHP-aware sensitivity (CAUSE delta_elpd z > 1.96 or LHC-MR posterior excluding zero); HDL is biased when sample overlap is above 5%
- PP.H4 >= 0.7 is the triangulation threshold, >= 0.8 publication-grade, >= 0.95 industry-grade
- Always sweep coloc.abf p12 over 1e-6 to 5e-5; conclusions must be stable across the sweep
- Use coloc.susie (not coloc.abf) for loci with allelic heterogeneity (multiple independent signals)
- Fine-mapping requires accurate LD; check estimate_s_rss lambda < 0.05 before trusting external reference-panel LD
- min_abs_corr >= 0.5 (r-squared >= 0.25) is the SuSiE purity threshold; lower sets are unreliable
- For TWAS, pick the tissue from Step 0 stratified LDSC; do not run all 49 GTEx tissues by default
- FOCUS PIP >= 0.8 is required to retain a single candidate causal gene per region; without it, co-regulated TWAS hits cannot be distinguished
- Cis-pQTL MR requires cross-platform replication on Olink (UKB-PPP) and SomaScan (deCODE); ~40% of proteins are platform-discordant
- For effector-gene prioritization, require >= 3 of 6 concordant evidence streams (L2G, PoPS, cS2G, coloc, TWAS, ABC); L2G and PoPS disagree by design and both should be reported
- Steiger pre-filtering of instruments (variance explained in exposure > variance explained in outcome) reduces reverse-causation bias
- GenomicSEM common-factor GWAS requires CFI > 0.95 and RMSEA < 0.06; Q_SNP p > 0.05 confirms factor-level (not trait-specific) signal
- For mediation, run Imai sensitivity (rho_crit > 0.3 robust) or mediational E-value > 2; molecular mediators are best tested via two-step MR with cis-instruments
Related Skills
causal-genomics/mendelian-randomization - IVW, Egger, MR-RAPS, MVMR causal-genomics/colocalization-analysis - coloc.abf, coloc.susie, HyPrColoc, SMR-HEIDI causal-genomics/fine-mapping - SuSiE rss, FINEMAP-inf, PolyFun, SuSiEx causal-genomics/pleiotropy-detection - MR-PRESSO, CAUSE, LHC-MR, contamination-mixture causal-genomics/mediation-analysis - Two-step MR, MVMR, CMAverse 4-way, HIMA causal-genomics/transcriptome-wide-association - FUSION, S-PrediXcan, FOCUS, UTMOST causal-genomics/heritability-partitioning - LDSC, S-LDSC, LDAK, HDL, HESS causal-genomics/proteome-mr-drug-target - UKB-PPP, deCODE, cis-pQTL MR, pheWAS causal-genomics/effector-gene-prioritization - L2G, PoPS, cS2G, MAGMA, FLAMES causal-genomics/genetic-correlation - Cross-trait LDSC, HDL, LAVA, Popcorn causal-genomics/genomic-sem - Common-factor GWAS, Q_SNP, MTAG reconciliation population-genetics/association-testing - Upstream GWAS methods atac-seq/enhancer-gene-linking - ABC / ENCODE-rE2G priors for effector-gene step single-cell/preprocessing - scRNA / scATAC tissue priors for stratified LDSC