This vignette shows how to jointly fine-map multiple related molecular measurements in one locus with fSuSiE.
Learning goals¶
After completing this vignette, you will be able to:
decide when fSuSiE is appropriate instead of a single-trait SuSiE analysis;
organize genotype, phenotype, covariate and region files for a joint fit;
run the current
qtl_dataset_construct+fsusieworkflow; andinspect posterior inclusion probabilities and credible-set assignments without printing the entire result object.
Background and method¶
fSuSiE is intended for a locus measured through multiple related functional traits, such as correlated chromatin features or other molecular measurements. It fits those measurements jointly so that shared evidence can help identify variants associated with the multivariate functional signal. A single gene-expression trait should instead use the ordinary susie_twas route.
The workflow first builds a QtlDataset, then defines the requested region and calls fineMappingPipeline() with the fsusie method. The main settings control the prior (--prior), maximum prior scale (--max-scale), credible-set purity threshold (--min-purity) and post-processing method (--post-processing). Setting --susie-top-pc above zero additionally fine-maps selected phenotype principal components with univariate SuSiE; it is optional and is not part of the default fSuSiE fit.
The bundled files below are a small multi-context molecular-trait example used to demonstrate the current interface. They are not a biological epigenomic benchmark, so substantive analyses should replace them with related functional measurements from the same samples and locus.
Worked example¶
The analysis unit is one locus (ENSG00000130538) measured in two contexts for the same individuals.
| Role | Example path | Required structure |
|---|---|---|
| Genotype | tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed | PLINK bed/bim/fam trio with samples matching the phenotypes |
| Phenotype manifest | input/finemapping/protocol_example.protocol_example.pheno_manifest_context.tsv | One row per context, with coordinates, trait identifier and phenotype BED path |
| Covariates | tests/fixtures/covariate_hidden_factor/covariates.tsv | Covariates by sample; this example is transposed with --transpose-covariates |
| Association windows | tests/fixtures/qtl_mini/association_windows.bed | Chromosome, start, end and region identifier |
All four inputs must use compatible genome coordinates and sample identifiers. Multiple phenotype rows must represent measurements that are scientifically meaningful to analyze jointly.
Run fSuSiE¶
qtl_dataset_construct+fsusie builds the regional dataset and then performs the joint fSuSiE fit.
Timing: TBD
sos run pipeline/mnm_regression.ipynb qtl_dataset_construct+fsusie \
--name protocol_example \
--cwd output_fsusie \
--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 \
-j 1The command is expected to create:
output_fsusie/qtl_dataset/protocol_example.qtl_dataset.rds— harmonized genotype, phenotype and covariate data.output_fsusie/fsusie/protocol_example.fsusie_region_manifest.tsv— the locus selected for the fit.output_fsusie/fsusie/protocol_example.ENSG00000130538.fsusie.rds— theQtlFineMappingResultcontaining the fSuSiE fit and variant-level summaries.
Key optional settings are:
| Option | Purpose |
|---|---|
--prior | Prior family used by fSuSiE; default mixture_normal |
--max-scale | Upper scale bound for the effect-size prior |
--min-purity | Minimum credible-set purity retained in post-processing |
--post-processing | fSuSiE post-processing rule |
--susie-top-pc | Number of phenotype PCs to fine-map separately; default 0 |
Command reference¶
Show all current workflows and options before adapting the example to a new dataset.
sos run pipeline/mnm_regression.ipynb -hResults and interpretation¶
fSuSiE result¶
The result is a QtlFineMappingResult. Its table-level metadata identifies the region and fitted method; each entry stores the fit plus a topLoci table with one row per variant.
| Level | Contents |
|---|---|
| Result table | Context, trait or region identifiers, method and one entry per fit |
entry | Method-specific fSuSiE fit and attached summaries |
topLoci | Variant identifiers, PIPs, posterior effects and credible-set labels |
| Credible-set columns | Membership at the requested coverage and the corresponding purity |
suppressPackageStartupMessages(library(pecotmr))
result_file <- "output_fsusie/fsusie/protocol_example.ENSG00000130538.fsusie.rds"
result <- readRDS(result_file)
cat("Class:", class(result), "\n")
print(result[, setdiff(names(result), "entry"), drop = FALSE])
entry <- result$entry[[1]]
top_loci <- getTopLoci(entry)
top_loci <- top_loci[order(top_loci$pip, decreasing = TRUE), ]
keep <- intersect(
c("variant_id", "pip", "posterior_mean", "posterior_sd",
"cs_95", "cs_95_purity", "method"),
names(top_loci)
)
head(top_loci[, keep, drop = FALSE], 8)Interpret variants by PIP first, then use credible-set membership and purity to judge localization. A high-PIP variant in a compact, high-purity credible set is more strongly localized than a similarly ranked variant in a large or low-purity set. Lack of a retained credible set in a toy dataset is not evidence that the locus has no functional effect; small sample size, weak signal, phenotype definition and LD can all prevent localization.
Limitations and common pitfalls¶
fSuSiE requires at least two joinable traits or contexts. Use ordinary SuSiE for a single trait.
The traits should describe related measurements in the same locus and samples; combining unrelated outcomes only because they share coordinates is not a valid joint analysis.
Sample identifiers must agree across genotype, phenotype and covariate files.
Chromosome naming and genome build must match across all inputs.
The example data demonstrate file structure and command composition, not realistic epigenomic power.
If the regional fit fails before producing the final RDS, inspect the per-region
.stderrfile rather than interpreting the intermediate dataset as a completed analysis.
Next steps¶
For a retained credible set, annotate the highest-PIP variants and compare their effects across the functional measurements represented in the joint fit. The result can then be passed to the appropriate post-processing or biological-validation workflow; do not treat fine-mapping probabilities alone as functional validation.