bio-workflows-causal-genomics-pipeline

What this skill does

## 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_name` to verify parameters
- Python: `pip show <pkg>` then `help(module.function)` to check signatures
- CLI: `<tool> --version` then `<tool> --help` to 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 methods
```

## Step 0: Pre-flight - Heritability, Tissue Prioritization, Genetic Correlation

```bash
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 rg
```

**Goal:** 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

```r
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

```r
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

```r
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)

```r
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

```r
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

```r
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$cs
```

If `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

```r
library(Two