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:
define the analysis unit and cis-association window used by
susie_twas;distinguish SuSiE fine-mapping outputs from TWAS prediction weights;
compare prediction methods using out-of-fold performance; and
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 |
|---|---|
tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed with matching .bim and .fam | PLINK genotype data |
input/finemapping/protocol_example.protocol_example.pheno_manifest_context.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 |
tests/fixtures/covariate_hidden_factor/covariates.tsv | Covariates for the two contexts |
tests/fixtures/qtl_mini/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 tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \
--phenoFile tests/fixtures/qtl_mini/protocol_example.pheno_manifest_context.tsv \
--covFile tests/fixtures/covariate_hidden_factor/covariates.tsv \
--customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
--region-name ENSG00000130538 \
--transpose-covariates \
--save-data \
-j1The 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 setsThe 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.045For 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-windowvs--customized-association-windows. Either/or. If both are passed the customized windows win. Pass--cis-window -1to 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.