Univariate fine-mapping and TWAS with SuSiE

Univariate fine-mapping and TWAS with SuSiE#

Fine-map cis-xQTL signals and estimate cross-validated transcriptome-wide association study (TWAS) weights for one molecular trait in one or more contexts.

Learning goals#

After completing this vignette, readers should be able to:

  1. define the analysis unit and cis-association window used by susie_twas;

  2. distinguish SuSiE fine-mapping outputs from TWAS prediction weights;

  3. compare prediction methods using out-of-fold performance; and

  4. identify credible sets, posterior inclusion probabilities, and common input failures.

Background and method#

qtl_dataset_construct+susie_twas in mnm_regression.ipynb builds a QtlDataset from genotype, phenotype, and covariate files. For each region, it residualizes genotype and phenotype matrices, fits SuSiE, and trains ten TWAS prediction models (mrash, susie, susie_inf, enet, lasso, mcp, scad, l0learn, bayes_r, and bayes_c) plus a cross-validation-based ensemble.

SuSiE represents the regional effect as a sum of single effects and reports posterior inclusion probabilities (PIPs) and credible sets. In this example, the maximum number of effects is L = 5. TWAS weights solve a different problem: they predict the molecular trait from regional genotypes. Fine-mapping evidence and predictive performance should therefore be interpreted separately.

The optional mnm workflow pools cross-validated predictions across contexts and estimates a constrained method-by-context mixture, producing multi-context TWAS weights.

Worked example#

The example fine-maps one gene in two molecular-trait contexts and estimates within-context TWAS weights.

Input files#

File

Purpose

input/finemapping/protocol_example.genotype.chr22.bed with matching .bim and .fam

PLINK genotype data

input/finemapping/protocol_example.pheno_manifest_multicontext.tsv

Maps each gene and context to its indexed phenotype BED and covariate file

input/finemapping/protocol_example.gene_expression.context1.bed.gz with .tbi

Context 1 molecular-trait values

input/finemapping/protocol_example.gene_expression.context2.bed.gz with .tbi

Context 2 molecular-trait values

input/finemapping/protocol_example.covariates.tsv

Covariates for the two contexts

input/finemapping/protocol_example.association_windows.bed

Cis-association window for each gene

Run univariate fine-mapping and estimate TWAS weights#

This command constructs the analysis dataset, fits SuSiE separately in each context, and estimates TWAS weights. --save-data retains the residualized matrices for analyses that reuse the fitted dataset.

sos run pipeline/mnm_regression.ipynb qtl_dataset_construct+susie_twas \
  --name protocol_example \
  --cwd output_uni \
  --genoFile input/finemapping/protocol_example.genotype.chr22.bed \
  --phenoFile input/finemapping/protocol_example.pheno_manifest_multicontext.tsv \
  --covFile input/finemapping/protocol_example.covariates.tsv \
  --customized-association-windows input/finemapping/protocol_example.association_windows.bed \
  --region-name ENSG00000130538 \
  --transpose-covariates \
  --save-data \
  -j1

The worked example produces:

  • output_uni/fine_mapping/protocol_example.ENSG00000130538.univariate_bvsr.rds: context-specific SuSiE fine-mapping results.

  • output_uni/twas_weights/protocol_example.ENSG00000130538.univariate_twas_weights.rds: prediction weights and cross-validation results.

  • output_uni/fine_mapping_meta.tsv: manifest linking regions to the generated fine-mapping and weight files.

Optional: adapt the analysis to peak-based traits#

Peak- or CpG-based traits use the same fine-mapping workflow, but first require a phenotype manifest and biologically appropriate association windows. For peak data, the helper below splits a multi-sample BED into indexed region files and assigns each feature to its smallest containing extended TAD. Replace the placeholders with your peak BED, extended-TAD reference, and output directory, then use the generated manifest and windows in the main command.

bash code/misc/preprocess/prepare_peak_inputs.sh \
  <raw_peaks.bed.gz> \
  <extended_TADB.bed> \
  <out_dir>

# Produces:
# <out_dir>/pheno_manifest.tsv
# <out_dir>/association_windows.bed
# <out_dir>/peaks_split/<ID>.bed.gz{,.tbi}

Parameters used in the example#

Parameter

Why it matters here

--phenoFile

Supplies the context-specific molecular traits through the phenotype manifest.

--customized-association-windows

Defines the variants considered for each gene instead of using a fixed cis window.

--region-name

Restricts this demonstration to ENSG00000130538.

--transpose-covariates

Matches the orientation of the supplied covariate table.

--save-data

Saves residualized matrices that can be reused downstream.

-j1

Runs one job for a reproducible, easy-to-debug example.

See the module link or command reference below for all remaining options and defaults.

Command reference#

sos run pipeline/mnm_regression.ipynb -h

Results and interpretation#

Fine-mapping result#

  • output_uni/fine_mapping/protocol_example.ENSG00000130538.univariate_bvsr.rds

The RDS contains a QtlFineMappingResult:

Level

Access

Contents

Result

attr(result, "listData")

Study, context, trait, method, and context-specific entries

Context entry

entries$entry[[i]]

Fine-mapping information for one context

topLoci

attr(entry, "topLoci")

Variant IDs, PIPs, posterior effects, credible-set membership, and purity

susieFit

attr(entry, "susieFit")

Fitted SuSiE model and retained credible sets

cvResult

attr(entry, "cvResult")

Cross-validation information used during weight estimation

 
suppressPackageStartupMessages(library(pecotmr))

result <- readRDS(
  "output_uni/fine_mapping/protocol_example.ENSG00000130538.univariate_bvsr.rds"
)
entries <- attr(result, "listData")
entry <- entries$entry[[which(entries$context == "context1")]]

top_loci <- attr(entry, "topLoci")
head(
  top_loci[order(top_loci$pip, decreasing = TRUE),
           c("variant_id", "pip", "cs_95", "cs_95_purity")],
  5
)
# variant_id             pip   cs_95   cs_95_purity
# chr22:16445263:A:T   0.032   susie_0       0
# chr22:16418860:A:T   0.031   susie_0       0
# chr22:16214125:A:T   0.027   susie_0       0

fit <- attr(entry, "susieFit")
length(fit$sets$cs)
# 0 retained credible sets

The largest PIPs in context 1 are approximately 0.03, every displayed candidate-set label has purity 0, and SuSiE retains no credible set after filtering. This toy example therefore does not support a localized fine-mapped signal. It demonstrates how to inspect PIPs and credible-set diagnostics, not a biological discovery.

TWAS prediction result#

  • output_uni/twas_weights/protocol_example.ENSG00000130538.univariate_twas_weights.rds

The RDS contains a TwasWeights object. attr(tw, "listData") gives the study, context, trait, method, and one entry per context-method combination.

Entry slot

Contents

variantIds

Variants aligned with the weight vector

weights

Per-variant prediction weights

fits

Fitted method-specific prediction model

cvResult

CV metrics for individual methods; methodCoef and methodPerformance for the ensemble

standardized

Whether genotype and phenotype inputs were standardized

dataType

Optional molecular-trait type metadata

 
suppressPackageStartupMessages(library(pecotmr))

tw <- readRDS(
  "output_uni/twas_weights/protocol_example.ENSG00000130538.univariate_twas_weights.rds"
)
entries <- attr(tw, "listData")

# Inspect the SuSiE predictor for context 1.
i <- which(entries$context == "context1" & entries$method == "susie")
susie_entry <- entries$entry[[i]]
round(attr(susie_entry, "cvResult")$metrics[
  c("corr", "rsq", "adj_rsq", "pval")
], 4)
# corr = 0.9798; rsq = 0.9599; adj_rsq = 0.9591; pval < 1e-4

# Compare cross-validated R2 where a method returned CV metrics.
i <- which(entries$context == "context1")
method_rsq <- sapply(entries$entry[i], function(entry) {
  cv <- attr(entry, "cvResult")
  if (!is.list(cv) || is.null(cv$metrics) || !is.numeric(cv$metrics)) {
    return(NA_real_)
  }
  cv$metrics["rsq"]
})
setNames(round(method_rsq, 3), entries$method[i])

# Inspect nonzero ensemble coefficients.
ensemble <- entries$entry[[
  which(entries$context == "context1" & entries$method == "ensemble")
]]
coef <- attr(ensemble, "cvResult")$methodCoef
round(coef[coef > 0.001], 3)
# bayes_r = 0.955; enet = 0.045

For context 1, the SuSiE predictor has cross-validated R2 = 0.960. The ensemble assigns approximately 0.955 of its coefficient weight to bayes_r and 0.045 to enet. These values show how to assess prediction performance and ensemble composition. Because the example contains only 49 samples, they should not be treated as stable estimates of out-of-sample accuracy.

Fine-mapping and TWAS answer different questions. Fine-mapping evaluates which variants are supported within the regional model; TWAS evaluates how well cis-genotypes predict the molecular trait. Before downstream TWAS, confirm allele alignment and validate prediction performance in an adequately sized dataset. A TWAS association does not by itself establish mediation or causality.

Limitations and common pitfalls#

  • Sample-ID drift. The phenotype BED’s header columns (5+), the bfile FID/IID, and the covariate TSV header must all be drawn from the same set. Sos does not warn on missing samples — it silently drops them and may report zero overlap.

  • Empty regions. If a region has no variants in its association window (e.g. the manifest points at a window with no bfile coverage) sos marks it completed with an empty output. Verify the expected file count in the INFO: susie_twas output: line.

  • --cis-window vs --customized-association-windows. Either/or. If both are passed the customized windows win. Pass --cis-window -1 to force the customized path explicitly.

Next steps#

Use the fine-mapping and TWAS mini-protocol to choose among univariate, multivariate, functional, multi-gene, and summary-statistic routes. Refer to mnm_regression.ipynb for the complete command interface and parameter defaults.

After validating credible sets and cross-validated prediction performance, the TWAS weights can enter the GWAS-integration workflows. Fine-mapping PIPs and TWAS associations answer different questions and should not be treated as interchangeable evidence.