Cross-gene fine-mapping with mvSuSiE

This vignette shows how to jointly fine-map several molecular traits whose coordinates overlap the same genomic region.

Learning goals

After completing this vignette, you will be able to:

• distinguish cross-gene fine-mapping from cross-context fine-mapping;

• define a locus containing multiple eligible genes;

• run the current qtl_dataset_construct+mnm_genes workflow; and

• interpret joint variant evidence without treating correlated genes as independent discoveries.

Background and method

Nearby genes can share regulatory variants and linkage disequilibrium. Separate fine-mapping analyses may therefore assign the same regional signal to several genes without testing their effects jointly. The mnm_genes workflow uses mvSuSiE with jointSpecification="trait" to model the genes overlapping each requested locus as a multivariate outcome.

The workflow first constructs one QtlDataset. It then creates a region manifest and identifies genes whose coordinates overlap each region. A region must contain at least two joinable genes; otherwise no multivariate fit is produced. The current implementation replaces the older TAD-specific, phenotype-map and univariate-PIP prefilter procedure described by the previous vignette.

Use this route for multiple genes in one locus. Use mnm for one gene across contexts and susie_twas for a single gene-context pair. A TAD is one possible biological definition of a locus, but the workflow accepts any defensible genomic window.

Worked example

The bundled chromosome 22 example contains several genes with overlapping association windows.

Role	Example path	What it supplies
Genotype	input/finemapping/protocol_example.multigene.genotype.chr22.bed	PLINK bed/bim/fam trio
Phenotype manifest	input/finemapping/protocol_example.multigene.pheno_manifest.tsv	Gene coordinates and phenotype BED path
Covariates	input/finemapping/protocol_example.multigene.covariates.tsv	Covariates for matching samples
Region windows	input/finemapping/protocol_example.multigene.association_windows.bed	Candidate loci used to group overlapping genes

The phenotype manifest may contain many genes, but each fitted region must overlap at least two of them. Coordinates, genome build and sample identifiers must agree across inputs.

Run cross-gene fine-mapping

qtl_dataset_construct+mnm_genes constructs the shared dataset, identifies the genes in each selected locus and fits mvSuSiE jointly.

Timing: TBD

sos run pipeline/mnm_regression.ipynb qtl_dataset_construct+mnm_genes \
  --name protocol_example \
  --cwd output_mnm_genes \
  --genoFile input/finemapping/protocol_example.multigene.genotype.chr22.bed \
  --phenoFile input/finemapping/protocol_example.multigene.pheno_manifest.tsv \
  --covFile input/finemapping/protocol_example.multigene.covariates.tsv \
  --customized-association-windows input/finemapping/protocol_example.multigene.association_windows.bed \
  --region-name ENSG00000130538 \
  --transpose-covariates \
  -j 1

The command is expected to create:

• output_mnm_genes/qtl_dataset/protocol_example.qtl_dataset.rds — the harmonized input dataset.

• output_mnm_genes/multivariate_fine_mapping/protocol_example.mnm_genes_region_manifest.tsv — regions and the genes assigned to each region.

• output_mnm_genes/multivariate_fine_mapping/protocol_example.ENSG00000130538.multigene_bvsr.rds — the cross-gene QtlFineMappingResult; a region with fewer than two joinable genes instead yields no fit.

Important options are --contexts, --mvsusie-max-iter, optional --prior-twas-weights, and JSON overrides through --mvsusie-method-args.

Command reference

Check the current interface before adapting the example.

sos run pipeline/mnm_regression.ipynb -h

Results and interpretation

Cross-gene result

The output is a QtlFineMappingResult with one joint entry for the genes assigned to the region.

Component	Interpretation
Region metadata	Fitted locus and participating genes
topLoci	Variant PIPs and credible-set assignments from the joint model
mvSuSiE fit	Posterior effects across the modeled genes
Credible-set purity	Strength of localization within the regional LD structure

suppressPackageStartupMessages(library(pecotmr))

result_file <- "output_mnm_genes/multivariate_fine_mapping/protocol_example.ENSG00000130538.multigene_bvsr.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", "cs_95", "cs_95_purity", "method"), names(top_loci))
head(top_loci[, keep, drop = FALSE], 8)

A high joint PIP identifies a variant supported by the multivariate regional model; it does not by itself identify which gene mediates the signal. Interpret the variant-level evidence together with posterior effects for each gene. If several genes retain similar effects, the data may not separate shared regulation from correlated phenotypes. A broad or low-purity credible set indicates limited localization.

Limitations and common pitfalls

• Each fitted region needs at least two joinable genes; a single-gene region is not an error but cannot produce a cross-gene fit.

• Region boundaries determine which genes enter the model and must be scientifically justified.

• Nearby expression traits can be correlated for technical or biological reasons; joint association does not establish mediation.

• All phenotype and genotype files must use the same genome build and matched samples.

• The toy data demonstrate workflow structure rather than realistic power.

• Do not use the obsolete --pheno_id_map_file or --skip-analysis-pip-cutoff interface from the earlier implementation.

Next steps

For a localized signal, compare gene-specific posterior effects, annotate the highest-PIP variants and evaluate independent functional evidence. Use downstream post-processing to compare credible sets with single-gene or cross-context results rather than declaring one target gene from the joint PIP alone.