Advanced regression models for association analysis with individual-level data

Advanced regression models for association analysis with individual-level data#

Fits univariate and multi-context fine-mapping models to individual-level genotype and molecular phenotype data, and derives TWAS weights.

Overview#

An association scan tells you a region matters; it does not tell you which variant in it is responsible, because variants in linkage disequilibrium carry nearly the same signal.

SuSiE reframes the question as variable selection: it fits a sum of single effects and returns credible sets - small groups of variants, each likely to contain one causal variant - together with posterior inclusion probabilities (Wang et al., 2020).

When the same locus is measured across several contexts - tissues, cell types, conditions - fine-mapping each separately throws away the shared structure. mvSuSiE learns the patterns of sharing from the data and uses them to sharpen credible sets (Zou et al., 2026).

The same fit also yields TWAS weights, the per-variant effects used later to predict expression into a GWAS cohort - which is why fine-mapping and weights come from one command rather than two, with skip_twas_weights defaulting to False. In practice the pipeline runs per region: susie_twas performs univariate fine-mapping and trains weights with several methods, writing one *.univariate_twas_weights.rds per region; run it with --save-data so the combined fine-mapping object is kept. mnm then consumes that output through a --fine-mapping-meta table and produces the multi-context ensemble weights, one *.multicontext_bvsr.rds per gene.

Two input modes take identical commands: a bulk RNA-seq matrix keyed by Ensembl gene ID with association windows at TSS +/- 500 kb, or a list of peaks or any other molecular phenotype - swap the phenotype BED and its manifest. mnm_genes (multi-gene), fsusie/mvfsusie (functional) and mvsusie (multivariate) are supported as well (Wang et al., 2020; Zou et al., 2026).

When to run it. After QTL association testing, with individual-level genotypes, phenotypes and covariates in hand. From summary statistics only, use rss_analysis.

Requires processed genotype (PLINK/VCF) and phenotype files from data_preprocessing. For TWAS, also requires the LD reference generated by rss_ld_sketch.

Input#

Gene-expression mode (input/colocboost/) - the example below:

  • --genoFile input/colocboost/example.chr22.bed (PLINK1 binary genotype set. Pass the .bed; the .bim and .fam are read automatically. chr22, 1,995 variants and 49 samples.)

  • --phenoFile input/colocboost/pheno_manifest_multicontext.tsv (phenotype list, one row per region and context. Columns #chr start end ID path cond cov_path: path is the bgzipped and tabixed phenotype matrix, cond names the context, cov_path the per-context covariates. 32 rows in the toy file.)

#chr	start	end	ID	path	cond	cov_path
chr22	10939387	10961338	ENSG00000283047	input_geneexpr/example_geneexpr.bed.gz	context1	input_geneexpr/example_covariates.tsv
chr22	10939387	10961338	ENSG00000283047	input_geneexpr/example_geneexpr_ctx2.bed.gz	context2	input_geneexpr/example_covariates.tsv

  • --covFile input/colocboost/example_covariates.tsv (covariates, one row per covariate and one column per sample; passed with --transpose-covariates.)

ID	SAMPLE_001	SAMPLE_002	SAMPLE_003	...
sex	1	1	1	...
age	67.06	82.92	87.76	...

  • --customized-association-windows input/colocboost/association_windows.bed (analysis window per region, columns #chr start end ID: the gene TSS +/- 500 kb, clamped at the chromosome edges. 16 regions.)

#chr	start	end	ID
chr22	10439387	11461338	ENSG00000283047
chr22	15028191	16029139	ENSG00000130538

Peak mode (input/) - same interface, alternative input:

  • --genoFile input/genotype/protocol_example.genotype.chr22.bed (same PLINK1 layout, a different cohort: chr22, 1,995 variants and 1,287 samples with S#### identifiers.)

  • --phenoFile input/phenotype/protocol_example.pheno_manifest.tsv (peak list, 5 columns #chr start end ID path; one row per peak, path pointing at the per-peak bgzipped matrix.)

#chr	start	end	ID	path
chr22	10662152	10664678	C22P107555	input/proteomics/protocol_example.peaks_split/C22P107555.bed.gz
chr22	16094161	16096165	C22P107557	input/proteomics/protocol_example.peaks_split/C22P107557.bed.gz

  • --covFile input/covariate/protocol_example.covariates.tsv (same shape as the gene-expression covariates, keyed by #id.)

#id	SAMPLE_001	SAMPLE_002	SAMPLE_003	...
sex	1	1	1	...
age	90.97	80.24	83.9	...

  • --customized-association-windows input/finemapping/protocol_example.association_windows.bed (cis window per peak, same four columns and no header row.)

chr22	0	19959512	C22P107555
chr22	0	19959512	C22P107557

Key requirement: the 4th-column ID of the association-windows file must match the 4th-column ID of the phenotype list, otherwise the window for that region is not found and the region is skipped.

The key is that the 4th column ID should match with the 4th column ID in the phenotype list. Otherwise the association-window to analyze will not be found.

About indels#

Option --no-indel will remove indel from analysis.

These apply across the workflows:

  • --cwd output/mnm (working directory all outputs are written under. Defaults to output.)

  • --name protocol_example (tag used in every output filename)

  • --customized-association-windows input/colocboost/association_windows.bed (association windows per region. Without it a fixed window of --cis-window is used.)

  • --cis-window 1000000 (half-width of the fixed association window used when no window file is given)

  • --region-name ENSG00000283047 (restrict the run to named regions; the union of all regions is used if omitted)

  • --save-data (keep the assembled per-region dataset alongside the fine-mapping result)

  • --transpose-covariates (set when the covariate file has samples in columns rather than rows)

  • --keep-samples (a file of sample IDs to retain before fine-mapping)

  • --max-cv-variants 5000 (cap on the variants carried into cross-validated TWAS weights; -1 for no cap)

  • --pheno-id-map-file (maps phenotype IDs to the IDs used in the region list, when the two differ)

  • --fine-mapping-meta output/mnm/fine_mapping_meta_mnm.tsv (the manifest of per-region fine-mapping results, consumed by the multi-gene and fSuSiE workflows)

Output#

qtl_dataset_construct

  • output/mnm/qtl_dataset/protocol_example.qtl_dataset.rds (the assembled QtlDataset: a genotype handle plus per-region phenotype and covariate data, reused by every method downstream.)

S4 object  <QtlDataset>
$ keepIndel: logical [1:1] TRUE
$ study: character [1:1] protocol_example
$ genotypes: S4 object  <GenotypeHandle>
 ..$ chromPaths: character [1:0] 
 ..$ path: character [1:1] input/colocboost/example.chr22
 ..$ format: character [1:1] plink1
 ..$ snpInfo: List of 5  <data.frame>
 .. ..$ SNP: character [1:1995] chr22:10414272:A:T chr22:10416111:A:T chr22:10416669:A:T chr22:10420511:A:T chr22:10424460:A:T ...
 .. ..$ CHR: character [1:1995] 22 22 22 22 22 ...
 .. ..$ BP: integer [1:1995] 10414272 10416111 10416669 10420511 10424460 ...
 .. ..$ A1: character [1:1995] A A A A A ...
... (truncated)

susie_twas

  • output/mnm/fine_mapping/protocol_example.ENSG00000130538.univariate_bvsr.rds (from susie_twas, the SuSiE fine-mapping result)

QtlFineMappingResult: S4 DataFrame, 2 rows x 5 columns
$ study: character [1:2] protocol_example protocol_example
$ context: character [1:2] context1 context2
$ trait: character [1:2] ENSG00000130538 ENSG00000130538
$ method: character [1:2] susie susie
$ entry: S4 object  <SimpleList>
 ..$ listData: List of 2
 .. ..$ : S4 object  <FineMappingEntry>
 .. .. ..$ variantIds: ...
 .. .. ..$ susieFit: ...
 .. .. ..$ topLoci: ...
 .. .. ..$ cvResult: ...
... (truncated)

  • output/mnm/twas_weights/protocol_example.ENSG00000130538.univariate_twas_weights.rds (from susie_twas, per-variant TWAS weights from each method, with cross-validation performance)

TwasWeights: S4 DataFrame, 22 rows x 5 columns
$ study: character [1:22] protocol_example protocol_example protocol_example protocol_example protocol_example ...
$ context: character [1:22] context1 context1 context1 context1 context1 ...
$ trait: character [1:22] ENSG00000130538 ENSG00000130538 ENSG00000130538 ENSG00000130538 ENSG00000130538 ...
$ method: character [1:22] mrash susie susie_inf enet lasso ...
$ entry: S4 object  <SimpleList>
 ..$ listData: List of 22
 .. ..$ : S4 object  <TwasWeightsEntry>
 .. .. ..$ variantIds: ...
 .. .. ..$ weights: ...
 .. .. ..$ fits: ...
 .. .. ..$ cvPerformance: ...
... (truncated)

mnm

  • multivariate_fine_mapping/{name}.{region}.multicontext_bvsr.rds (the multi-context fit, one per region.)

QtlFineMappingResult: S4 DataFrame, 2 rows x 5 columns
$ study: character [1:2] protocol_example protocol_example
$ context: character [1:2] context1 context2
$ trait: character [1:2] ENSG00000130538 ENSG00000130538
$ method: character [1:2] mvsusie mvsusie
$ entry: S4 object  <SimpleList>
 ..$ listData: List of 2
 .. ..$ : S4 object  <FineMappingEntry>
 .. .. ..$ variantIds: ...
 .. .. ..$ susieFit: ...
 .. .. ..$ topLoci: ...
 .. .. ..$ cvResult: ...
 .. ..$ : S4 object  <FineMappingEntry>
... (truncated)

mnm_genes

  • multivariate_fine_mapping/{name}.mnm_genes_region_manifest.tsv

  • multivariate_fine_mapping/{name}.{region_id}.multigene_bvsr.rds

fsusie

  • fsusie/{name}.fsusie_region_manifest.tsv

  • fsusie/{name}.{region_id}.fsusie.rds

mvfsusie

  • {name}.{region}.mvfsusie_{prior}.rds

library(pecotmr)
w <- readRDS('output/mnm/multivariate_fine_mapping/protocol_example.ENSG00000130538.multicontext_bvsr.rds')
w$context                      # contexts fitted, one row each
w$method                       # e.g. "mvsusie"
str(w$entry[[1]], max.level = 1)   # variantIds / susieFit / topLoci / cvResult

Minimal Working Example#

Univariate SuSiE fine-mapping and TWAS weights (susie_twas) - gene expression#

This is the main example: cis xQTL fine-mapping on gene expression. Two steps, pointed at the input/colocboost/ files. Use the 1-gene phenotype list for the fastest run; swap in pheno_manifest.tsv / association_windows.bed to analyze all 16 in-range genes.

Timing: ~3 min (on toy dataset)

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

Step 2. Multi-context fine-mapping and ensemble TWAS weights#

This reads those weights and adds the multi-context ensemble.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_regression.ipynb mnm \
    --name protocol_example --cwd output_mnm \
    --genoFile input/colocboost/example.chr22.bed \
    --phenoFile input/colocboost/pheno_manifest_multicontext.tsv \
    --covFile input/colocboost/example_covariates.tsv \
    --customized-association-windows input/colocboost/association_windows.bed \
    --fine-mapping-meta input/colocboost/fine_mapping_meta.tsv \
    --transpose-covariates --save-data -j1

Univariate SuSiE fine-mapping and TWAS weights (susie_twas) - peak#

The same pipeline also accepts peaks (chromatin/epigenomic) or any other molecular phenotype - just swap in the phenotype .bed.gz and its manifest; the command is otherwise identical to Mode A. Example below uses the 10 toy chr22 peaks.

Step 1. SuSiE + TWAS weights#

Timing: TBD (on toy dataset)

sos run pipeline/mnm_regression.ipynb susie_twas \
    --name protocol_example \
    --cwd output_10peaks \
    --genoFile input/genotype/protocol_example.genotype.chr22.bed \
    --phenoFile input/phenotype/protocol_example.pheno_manifest.tsv \
    --covFile input/covariate/protocol_example.covariates.tsv \
    --customized-association-windows input/finemapping/protocol_example.association_windows.bed \
    -j1

Step 2. Multi-context fine-mapping and ensemble TWAS weights#

Timing: TBD (on toy dataset)

sos run pipeline/mnm_regression.ipynb mnm \
    --name protocol_example \
    --cwd output_10peaks \
    --genoFile input/genotype/protocol_example.genotype.chr22.bed \
    --phenoFile input/phenotype/protocol_example.pheno_manifest.tsv \
    --covFile input/covariate/protocol_example.covariates.tsv \
    --customized-association-windows input/finemapping/protocol_example.association_windows.bed \
    --fine-mapping-meta output_10peaks/fine_mapping_meta.tsv \
    -j1

Multivariate analysis with mvSuSiE and mr.mash (mnm)#

Fine-maps multiple molecular phenotypes jointly across a shared window using mvSuSiE with mr.mash priors. This step needs multiple contexts (e.g. the same genes measured under two or more conditions). The toy set ships a single context, so this MWE adds a second, synthetic context (example/colocboost_ctx2.bed.gz, derived from context 1 with added noise — demo data only, not a real measurement) and a multi-context phenotype manifest. With that in place the mnm step now runs end-to-end on the toy data and produces multicontext_bvsr.rds, multicontext_bvsr.rds, and multicontext_data.rds per gene.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_regression.ipynb mnm \
    --name protocol_example --cwd output_mnm3 \
    --genoFile input/colocboost/example.chr22.bed \
    --phenoFile input/colocboost/pheno_manifest_multicontext.tsv \
    --covFile input/colocboost/example_covariates.tsv \
    --customized-association-windows input/colocboost/association_windows.bed \
    --save-data -j1

Multi-gene / trans fine-mapping (mnm_genes)#

Multivariate fine-mapping treating each gene as a separate condition across a shared window. On the toy set every region returns a NULL multigene_bvsr (no multi-gene/trans signal), which is the correct, honest result for these independent single-gene regions.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_regression.ipynb mnm_genes \
    --name protocol_example --cwd output_mnmgenes \
    --genoFile input/colocboost/example.chr22.bed \
    --phenoFile input/colocboost/pheno_manifest_multicontext.tsv \
    --covFile input/colocboost/example_covariates.tsv \
    --customized-association-windows input/colocboost/association_windows.bed \
    --pheno_id_map_file input/colocboost/pheno_id_map.tsv \
    --fine_mapping_meta input/colocboost/fine_mapping_meta.tsv \
    --keep-samples input/colocboost/keep_samples.txt \
    --save-data -j1

Functional regression fSuSiE for epigenomic QTL (fsusie)#

Functional SuSiE for fine-mapping a functional phenotype across genomic positions in a region. This workflow produces a *.fsusie_mixture_normal_TI__top_pc_weights.rds weight file.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_regression.ipynb fsusie \
    --name protocol_example --cwd output/fsusie \
    --genoFile input/genotype/protocol_example.genotype.chr22.bed \
    --phenoFile input/phenotype/protocol_example.pheno_manifest.tsv \
    --covFile input/covariate/protocol_example.covariates.tsv \
    --customized-association-windows input/finemapping/protocol_example.association_windows.bed \
    --cis-window 0 --max-cv-variants 5000 --save-data -j1

Functional regression fSuSiE with another modality (mvfsusie)#

Multivariate functional SuSiE (mvfSuSiE). Still under development for the toy dataset.

Timing: TBD (on toy dataset)

sos run pipeline/mnm_regression.ipynb mvfsusie \
    --name protocol_example --cwd output/mvfsusie \
    --genoFile input/genotype/protocol_example.genotype.chr22.bed \
    --phenoFile input/phenotype/protocol_example.pheno_manifest.tsv \
    --covFile input/covariate/protocol_example.covariates.tsv \
    --customized-association-windows input/finemapping/protocol_example.association_windows.bed \
    --save-data -j1

Command Interface#

sos run pipeline/mnm_regression.ipynb -h

usage: sos run pipeline/mnm_regression.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  qtl_dataset_construct
  susie_twas
  mnm
  mnm_genes
  fsusie
  mvfsusie

Global Workflow Options:
  --modular-script-dir code/script (as path)
  --name VAL (as str, required)
                        It is required to input the name of the analysis
  --cwd output (as path)
  --genoFile VAL (as path, required)
                        A list of file paths for genotype data, or the genotype
                        data itself.
  --phenoFile  paths

                        One or multiple lists of file paths for phenotype data.
  --phenoIDFile  paths()

                        One or multiple lists of file paths for phenotype ID
                        mapping file. The first column should be the original
                        ID, the 2nd column should be the ID to be mapped to.
  --covFile  paths

                        Covariate file path
  --region-list . (as path)
                        Optional: if a region list is provide the analysis will
                        be focused on provided region. The LAST column of this
                        list will contain the ID of regions to focus on
                        Otherwise, all regions with both genotype and phenotype
                        files will be analyzed
  --region-name  (as list)
                        Optional: if a region name is provided the analysis
                        would be focused on the union of provides region list
                        and region names
  --keep-samples . (as path)
                        Only focus on a subset of samples
  --keep-variants . (as path)
                        Only focus on a subset of variants
  --customized-association-windows . (as path)
                        An optional list documenting the custom association
                        window for each region to analyze, with four column,
                        chr, start, end, region ID (eg gene ID). If this list is
                        not provided, the default `window` parameter (see below)
                        will be used.
  --cis-window -1 (as int)
                        Specify the cis window for the up and downstream radius
                        to analyze around the region of interest in units of bp
                        When this is set to negative, we will rely on using
                        customized_association_windows
  --[no-]save-data (default to False)
                        save data object or not
  --phenotype-names  [f'{x:bn}' for x in phenoFile]

                        Name of phenotypes
  --[no-]trans-analysis (default to False)
                        And indicator whether it is trans-analysis, ie, not
                        using phenotypic coordinate information
  --seed 999 (as int)
  --init-L 8 (as int)
                        association analysis paramters initial number of single
                        effects for SuSiE
  --L 30 (as int)
                        maximum number of single effects to use for SuSiE
  --[no-]estimate-residual-variance (default to True)
  --imiss 1.0 (as float)
                        remove a variant if it has more than imiss missing
                        individual level data
  --maf 0.0025 (as float)
                        MAF and variance of X cutoff
  --xvar-cutoff 0.0 (as float)
  --mac 5 (as int)
                        MAC cutoff, on top of MAF cutoff
  --[no-]indel (default to True)
                        Remove indels if indel = False
  --pip-cutoff 0.025 (as float)
  --coverage 0.95 0.7 0.5 (as list)
  --skip-analysis-pip-cutoff  (as list)
                        If this value is not 0, then an initial single effect
                        analysis will be performed to determine if follow up
                        analysis will be continued or to simply return NULL If
                        this is negative we use a default way to determine this
                        cutoff which is conservative but still useful
  --[no-]skip-fine-mapping (default to False)
                        Skip fine-mapping
  --[no-]skip-twas-weights (default to False)
                        Skip TWAS weights computation
  --twas-cv-folds 5 (as int)
                        Perform K folds valiation CV for TWAS Set it to zero if
                        this is to be skipped
  --twas-cv-threads  twas_cv_folds

  --max-cv-variants -1 (as int)
                        maximum number of variants to consider for CV We will
                        randomly pick a subset of it for CV purpose We can set
                        it to eg 8000 to save computational burden althought may
                        risk overfitting for methods comparison purpose When set
                        to -1 we don't use this feature
  --ld-reference-meta-file . (as path)
  --container ''
                        Analysis environment settings
  --entrypoint  ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""

  --job-size 200 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 1h
                        Wall clock time expected
  --mem 20G
                        Memory expected
  --numThreads 1 (as int)
                        Number of threads

Sections
  qtl_dataset_construct:
    Workflow Options:
      --study study
                        Build one pecotmr::QtlDataset from the phenotype
                        manifest + shared genotype and serialize to RDS,
                        replacing the legacy get_analysis_regions /
                        regional_data loader. No fan-out: downstream per-gene
                        fine-mapping / TWAS load this single RDS and select
                        contexts/genes at analysis time.  Manifest schema
                        (qtl_dataset_construct.R): columns ID, path, cond and an
                        optional cov_path; extra columns (e.g. #chr/start/end)
                        are ignored. The toy MWE pheno_manifest_multicontext.tsv
                        already matches -- a single-context analysis just points
                        at a one-`cond` manifest (or selects the context
                        downstream). QTLtools-format covariate TSVs (rows =
                        covariates) need --transpose-covariates.
      --[no-]transpose-covariates (default to False)
                        Set for QTLtools-format covariate TSVs (rows =
                        covariates, cols = samples); required for the toy MWE
                        covariates.
      --genotype-covariates . (as path)
                        Optional genotype-derived covariates applied uniformly
                        across contexts; per-context covariates come from the
                        manifest cov_path column.
      --maf-cutoff 0.0 (as float)
                        Construct-time variant/sample QC (stored on the
                        QtlDataset; cutoff 0 = off).
      --mac-cutoff 0.0 (as float)
      --imiss-cutoff 0.0 (as float)
      --[no-]drop-indel (default to False)
                        Drop indel variants (default keeps them, matching
                        QtlDataset()).
  susie_twas:
    Workflow Options:
      --fine-mapping-methods susie
      --twas-methods 'susie,mrash,enet,lasso,mcp,scad,l0learn,bayes_r,bayes_c'
      --fine-mapping-coverage 0.95 (as float)
      --contexts ''
                        Comma-separated context names to restrict both passes
                        to; empty = all contexts.
      --pip-cutoff-to-skip 0.0 (as float)
                        Fine-mapping SER pre-screen (0 = off, <0 = adaptive
                        3/nVariants).
      --twas-maf-cutoff 0.01 (as float)
                        TWAS weight-learning knobs. TWAS learns weights applied
                        out-of-sample; keep an explicit MAF / X-variance cutoff
                        for the TWAS pass (passed via the shared --maf-
                        cutoff/--xvar-cutoff).
      --twas-xvar-cutoff 0.01 (as float)
      --cv-folds 5 (as int)
      --cv-threads 1 (as int)
      --maf-cutoff -1.0 (as float)
                        Opt-in per-analysis overrides of the QtlDataset's
                        construct-time filters, applied to both passes.
                        Sentinels: a negative cutoff / drop_indel False / "."
                        path = leave the dataset's stored value untouched.
      --mac-cutoff -1.0 (as float)
      --imiss-cutoff -1.0 (as float)
      --[no-]drop-indel (default to False)
      --fine-mapping-method-args ''
                        Optional JSON of per-method kwargs forwarded to the
                        pipelines, e.g. '{"susie":{"L":10}}' for fine-mapping or
                        '{"lasso":{"nfolds":10}}' for TWAS.
      --twas-method-args ''
  mnm:
    Workflow Options:
      --mvsusie-max-iter 200 (as int)
      --pip-cutoff-to-skip 0.0 (as float)
      --contexts ''
                        Comma-separated context names to restrict to; empty =
                        all contexts.
      --prior-twas-weights . (as path)
                        Optional TwasWeights RDS (preceding mr.mash run)
                        supplying the data-driven prior.
      --prior-weights-min 1e-10 (as float)
      --mvsusie-method-args ''
                        Optional JSON of extra mvsusie kwargs (overrides the
                        max_iter default), e.g. '{"mvsusie":{"max_iter":500}}'.
  mnm_genes_1:
  mnm_genes_2:
    Workflow Options:
      --mvsusie-max-iter 200 (as int)
                        Cross-gene (multi-trait) mvSuSiE per locus: region mode
                        + jointSpecification="trait" joins the genes whose
                        coordinates overlap the locus. Optional mr.mash data-
                        driven prior via --prior-twas-weights (canonical
                        otherwise). Replaces the legacy multigene combine +
                        multigene_udr machinery.
      --pip-cutoff-to-skip 0.0 (as float)
      --contexts ''
      --prior-twas-weights . (as path)
      --prior-weights-min 1e-10 (as float)
      --mvsusie-method-args ''
  fsusie_1:
  fsusie_2:
    Workflow Options:
      --prior 'mixture_normal'
                        Functional SuSiE per locus (fsusieR::susiF) over the
                        pre-built QtlDataset. fsusie-specific knobs ride on
                        --method-args; --susie-top-pc>0 additionally fine-maps
                        each context's top principal components with univariate
                        SuSiE (usePCA; ports the legacy fsusie.R
                        susie_on_top_pc). Replaces fsusie.R.
      --max-scale 10 (as int)
      --min-purity 0.5 (as float)
      --post-processing TI
      --max-SNP-EM 100 (as int)
      --[no-]small-sample-correction (default to False)
      --susie-top-pc 0 (as int)
                        >0 -> also fine-map the top `susie_top_pc` PCs per
                        context with univariate SuSiE.
      --pip-cutoff-to-skip 0.0 (as float)
      --contexts ''
      --fsusie-method-args ''
                        Optional JSON overriding the assembled fsusie kwargs,
                        e.g. '{"fsusie":{"prior":"mixture_normal_per_scale"}}'.
  mvfsusie:
    Workflow Options:
      --prior 'mixture_normal_per_scale'
                        prior can be either of ["mixture_normal",
                        "mixture_normal_per_scale"]
      --max-SNP-EM 1000 (as int)

Workflow implementation#

The runnable SoS workflow sections below are carried over verbatim from the original mnm_regression.ipynb. Sections still under active development are marked as WIP placeholders and are not part of the supported MWE.

[global]
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# It is required to input the name of the analysis
parameter: name = str
parameter: cwd = path("output")
# A list of file paths for genotype data, or the genotype data itself. 
parameter: genoFile = path
# One or multiple lists of file paths for phenotype data.
parameter: phenoFile = paths
# One or multiple lists of file paths for phenotype ID mapping file. The first column should be the original ID, the 2nd column should be the ID to be mapped to.
parameter: phenoIDFile = paths()
# Covariate file path
parameter: covFile = paths
# Optional: if a region list is provide the analysis will be focused on provided region. 
# The LAST column of this list will contain the ID of regions to focus on
# Otherwise, all regions with both genotype and phenotype files will be analyzed
parameter: region_list = path()
# Optional: if a region name is provided 
# the analysis would be focused on the union of provides region list and region names
parameter: region_name = []
# Only focus on a subset of samples
parameter: keep_samples = path()
# Only focus on a subset of variants
parameter: keep_variants = path()
# An optional list documenting the custom association window for each region to analyze, with four column, chr, start, end, region ID (eg gene ID).
# If this list is not provided, the default `window` parameter (see below) will be used.
parameter: customized_association_windows = path()
# Specify the cis window for the up and downstream radius to analyze around the region of interest in units of bp
# When this is set to negative, we will rely on using customized_association_windows
parameter: cis_window = -1
# save data object or not
parameter: save_data = False
# Name of phenotypes
parameter: phenotype_names = [f'{x:bn}' for x in phenoFile]
# And indicator whether it is trans-analysis, ie, not using phenotypic coordinate information
parameter: trans_analysis = False
parameter: seed = 999

# association analysis paramters
# initial number of single effects for SuSiE
parameter: init_L = 8
# maximum number of single effects to use for SuSiE
parameter: L = 30
parameter: estimate_residual_variance = True
# remove a variant if it has more than imiss missing individual level data
parameter: imiss = 1.0
# MAF and variance of X cutoff
parameter: maf = 0.0025
parameter: xvar_cutoff = 0.0
# MAC cutoff, on top of MAF cutoff
parameter: mac = 5
# Remove indels if indel = False
parameter: indel = True
parameter: pip_cutoff = 0.025
parameter: coverage = [0.95, 0.7, 0.5]
# If this value is not 0, then an initial single effect analysis will be performed 
# to determine if follow up analysis will be continued or to simply return NULL
# If this is negative we use a default way to determine this cutoff which is conservative but still useful
parameter: skip_analysis_pip_cutoff = []
# Skip fine-mapping
parameter: skip_fine_mapping = False
# Skip TWAS weights computation
parameter: skip_twas_weights = False
# Perform K folds valiation CV for TWAS
# Set it to zero if this is to be skipped
parameter: twas_cv_folds = 5
parameter: twas_cv_threads = twas_cv_folds
# maximum number of variants to consider for CV
# We will randomly pick a subset of it for CV purpose
# We can set it to eg 8000 to save computational burden althought may risk overfitting for methods comparison purpose
# When set to -1 we don't use this feature
parameter: max_cv_variants = -1
parameter: ld_reference_meta_file = path()
# Analysis environment settings
parameter: container = ""
import re
parameter: entrypoint= ('micromamba run -a "" -n' + ' ' + re.sub(r'(_apptainer:latest|_docker:latest|\.sif)$', '', container.split('/')[-1])) if container else ""
# For cluster jobs, number commands to run per job
parameter: job_size = 200
# Wall clock time expected
parameter: walltime = "1h"
# Memory expected
parameter: mem = "20G"
# Number of threads
parameter: numThreads = 1

if len(phenoFile) != len(covFile):
    raise ValueError("Number of input phenotypes files must match that of covariates files")
if len(phenoFile) != len(phenotype_names):
    raise ValueError("Number of input phenotypes files must match the number of phenotype names")
if len(phenoIDFile) > 0 and len(phenoFile) != len(phenoIDFile):
    raise ValueError("Number of input phenotypes files must match the number of phenotype ID mapping files")

if len(skip_analysis_pip_cutoff) == 0:
    skip_analysis_pip_cutoff = [0.0] * len(phenoFile)
if len(skip_analysis_pip_cutoff) == 1:
    skip_analysis_pip_cutoff = skip_analysis_pip_cutoff * len(phenoFile)
if len(skip_analysis_pip_cutoff) != len(phenoFile):
    raise ValueError(f"``skip_analysis_pip_cutoff`` should have either length 1 or length the same as phenotype files ({len(phenoFile)} in this case)")

# make it into an R List string
skip_analysis_pip_cutoff = [f"'{y}'={x}" for x,y in zip(skip_analysis_pip_cutoff, phenotype_names)]

[qtl_dataset_construct]
# Build one pecotmr::QtlDataset from the phenotype manifest + shared genotype and
# serialize to RDS, replacing the legacy get_analysis_regions / regional_data
# loader. No fan-out: downstream per-gene fine-mapping / TWAS load this single
# RDS and select contexts/genes at analysis time.
#
# Manifest schema (qtl_dataset_construct.R): columns ID, path, cond and an
# optional cov_path; extra columns (e.g. #chr/start/end) are ignored. The toy
# MWE pheno_manifest_multicontext.tsv already matches -- a single-context
# analysis just points at a one-`cond` manifest (or selects the context
# downstream). QTLtools-format covariate TSVs (rows = covariates) need
# --transpose-covariates.
parameter: study = name
# Set for QTLtools-format covariate TSVs (rows = covariates, cols = samples);
# required for the toy MWE covariates.
parameter: transpose_covariates = False
# Optional genotype-derived covariates applied uniformly across contexts;
# per-context covariates come from the manifest cov_path column.
parameter: genotype_covariates = path('.')
# Construct-time variant/sample QC (stored on the QtlDataset; cutoff 0 = off).
parameter: maf_cutoff = 0.0
parameter: xvar_cutoff = 0.0
parameter: mac_cutoff = 0.0
parameter: imiss_cutoff = 0.0
# Drop indel variants (default keeps them, matching QtlDataset()).
parameter: drop_indel = False
# Optional whitespace-delimited ID files to restrict samples / variants ("." = none).
parameter: keep_samples = "."
parameter: keep_variants = "."
output: f"{cwd:a}/qtl_dataset/{study}.qtl_dataset.rds"
task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f"{step_name}_{_output:bn}"
bash: expand = '${ }', stderr = f"{_output}.stderr", stdout = f"{_output}.stdout", container = container
    Rscript ${modular_script_dir}/pecotmr_integration/qtl_dataset_construct.R \
        --study ${study} \
        --genotype-prefix ${genoFile:n} \
        --phenotype-manifest ${phenoFile[0]} \
        --genotype-covariates ${genotype_covariates if genotype_covariates.is_file() else '""'} \
        ${'--transpose-covariates' if transpose_covariates else ''} \
        --maf-cutoff ${maf_cutoff} \
        --xvar-cutoff ${xvar_cutoff} \
        --mac-cutoff ${mac_cutoff} \
        --imiss-cutoff ${imiss_cutoff} \
        ${'--drop-indel' if drop_indel else ''} \
        ${('--keep-samples ' + str(keep_samples)) if keep_samples != '.' else ''} \
        ${('--keep-variants ' + str(keep_variants)) if keep_variants != '.' else ''} \
        --output ${_output}

[susie_twas]
# Univariate SuSiE fine-mapping + TWAS weights over the pre-built QtlDataset
# (from the qtl_dataset_construct step), one gene per fan-out unit. The gene(s)
# to analyze are passed on the command line via --region-name; the QtlDataset
# RDS already holds the genotype + per-context phenotypes, so this step neither
# re-reads the phenotype manifest nor rebuilds regional_data. Fine-mapping runs
# first (pecotmr_integration/fine_mapping.R -> fineMappingPipeline); the TWAS
# weights pass then reuses its SuSiE fits via --fine-mapping-result
# (pecotmr_integration/twas_weights.R -> twasWeightsPipeline). Replaces the
# legacy regional_data + susie_twas.R path.
parameter: cis_window = 1000000
parameter: fine_mapping_methods = "susie"
parameter: twas_methods = "susie,mrash,enet,lasso,mcp,scad,l0learn,bayes_r,bayes_c"
parameter: fine_mapping_coverage = 0.95
# Comma-separated context names to restrict both passes to; empty = all contexts.
parameter: contexts = ""
# Fine-mapping SER pre-screen (0 = off, <0 = adaptive 3/nVariants).
parameter: pip_cutoff_to_skip = 0.0
# TWAS weight-learning knobs.
# TWAS learns weights applied out-of-sample; keep an explicit MAF / X-variance
# cutoff for the TWAS pass (passed via the shared --maf-cutoff/--xvar-cutoff).
parameter: twas_maf_cutoff = 0.01
parameter: twas_xvar_cutoff = 0.01
parameter: max_cv_variants = 5000
parameter: cv_folds = 5
parameter: cv_threads = 1
# Reproducibility: integer RNG seed; -1 = unset.
parameter: seed = -1
# Opt-in per-analysis overrides of the QtlDataset's construct-time filters,
# applied to both passes. Sentinels: a negative cutoff / drop_indel False /
# "." path = leave the dataset's stored value untouched.
parameter: maf_cutoff = -1.0
parameter: mac_cutoff = -1.0
parameter: xvar_cutoff = -1.0
parameter: imiss_cutoff = -1.0
parameter: drop_indel = False
parameter: keep_samples = "."
parameter: keep_variants = "."
# Optional JSON of per-method kwargs forwarded to the pipelines, e.g.
# '{"susie":{"L":10}}' for fine-mapping or '{"lasso":{"nfolds":10}}' for TWAS.
parameter: fine_mapping_method_args = ""
parameter: twas_method_args = ""
fail_if(len(region_name) == 0, "susie_twas: pass --region-name <gene-id> [<gene-id> ...] to choose which gene(s) to analyze.")
qtl_dataset = path(f"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds")
# Assemble the shared opt-in worker flags once (command construction only).
seed_arg = f"--seed {seed}" if seed >= 0 else ""
# Non-MAF/xvar filters are shared by both passes (empty = use the dataset slots).
filter_overrides = " ".join(filter(None, [
    f"--mac-cutoff {mac_cutoff}"       if mac_cutoff   >= 0 else "",
    f"--imiss-cutoff {imiss_cutoff}"   if imiss_cutoff >= 0 else "",
    "--drop-indel"                     if drop_indel        else "",
    f"--keep-samples {keep_samples}"   if keep_samples  != "." else "",
    f"--keep-variants {keep_variants}" if keep_variants != "." else "",
]))
# MAF / X-variance ride the shared --maf-cutoff/--xvar-cutoff, handled per pass:
# fine-mapping opt-in only (default off = dataset slot); TWAS always explicit
# (its portability cutoff) unless a general override is set, which wins for both.
fm_maf_xvar = " ".join(filter(None, [
    f"--maf-cutoff {maf_cutoff}"   if maf_cutoff  >= 0 else "",
    f"--xvar-cutoff {xvar_cutoff}" if xvar_cutoff >= 0 else "",
]))
twas_maf_xvar = " ".join([
    f"--maf-cutoff {maf_cutoff if maf_cutoff >= 0 else twas_maf_cutoff}",
    f"--xvar-cutoff {xvar_cutoff if xvar_cutoff >= 0 else twas_xvar_cutoff}",
])
input: qtl_dataset, for_each = "region_name"
output: fine_mapping = f"{cwd:a}/fine_mapping/{name}.{_region_name}.univariate_bvsr.rds",
        twas_weights = f"{cwd:a}/twas_weights/{name}.{_region_name}.univariate_twas_weights.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f"{step_name}_{_output[0]:bnn}"
bash: expand = "${ }", stderr = f"{_output[0]:nn}.susie_twas.stderr", stdout = f"{_output[0]:nn}.susie_twas.stdout", container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \
        --qtl-dataset ${_input} \
        --gene-id ${_region_name} \
        --cis-window ${cis_window} \
        --methods ${fine_mapping_methods} \
        --coverage ${fine_mapping_coverage} \
        --pip-cutoff-to-skip ${pip_cutoff_to_skip} \
        ${seed_arg} ${filter_overrides} ${fm_maf_xvar} \
        ${("--contexts " + contexts) if contexts else ""} \
        ${("--method-args '" + fine_mapping_method_args + "'") if fine_mapping_method_args else ""} \
        --output ${_output["fine_mapping"]}

    Rscript ${modular_script_dir}/pecotmr_integration/twas_weights.R \
        --qtl-dataset ${_input} \
        --gene-id ${_region_name} \
        --cis-window ${cis_window} \
        --methods ${twas_methods} \
        --max-cv-variants ${max_cv_variants} \
        --cv-folds ${cv_folds} \
        --cv-threads ${cv_threads} \
        --fine-mapping-result ${_output["fine_mapping"]} \
        ${seed_arg} ${filter_overrides} ${twas_maf_xvar} \
        ${("--contexts " + contexts) if contexts else ""} \
        ${("--method-args '" + twas_method_args + "'") if twas_method_args else ""} \
        --output ${_output["twas_weights"]}

[mnm]
# Multivariate cross-context fine-mapping with mvSuSiE over the pre-built
# QtlDataset (from qtl_dataset_construct), one gene per fan-out unit. Each gene's
# contexts are fit jointly (jointSpecification="context"). The canonical mixture
# prior is used by default; pass --prior-twas-weights (a TwasWeights RDS from a
# preceding `twas_weights.R --methods mrmash` run) for the mr.mash data-driven
# reweighted prior. Replaces the legacy mnm.R + multivariate_analysis_pipeline.
parameter: cis_window = 1000000
parameter: mvsusie_max_iter = 200
parameter: pip_cutoff_to_skip = 0.0
# Comma-separated context names to restrict to; empty = all contexts.
parameter: contexts = ""
# Optional TwasWeights RDS (preceding mr.mash run) supplying the data-driven prior.
parameter: prior_twas_weights = path()
parameter: prior_weights_min = 1e-10
# Optional JSON of extra mvsusie kwargs (overrides the max_iter default), e.g. '{"mvsusie":{"max_iter":500}}'.
parameter: mvsusie_method_args = ""
fail_if(len(region_name) == 0, "mnm: pass --region-name <gene-id> [<gene-id> ...] to choose which gene(s) to analyze.")
qtl_dataset = path(f"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds")
mvsusie_args = mvsusie_method_args if mvsusie_method_args else f'{{"mvsusie":{{"max_iter":{mvsusie_max_iter}}}}}'
input: qtl_dataset, for_each = "region_name"
output: f"{cwd:a}/multivariate_fine_mapping/{name}.{_region_name}.multicontext_bvsr.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f"{step_name}_{_output:bn}"
bash: expand = "${ }", stderr = f"{_output:n}.stderr", stdout = f"{_output:n}.stdout", container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \
        --qtl-dataset ${_input} \
        --gene-id ${_region_name} \
        --cis-window ${cis_window} \
        --methods mvsusie \
        --joint-specification context \
        --coverage ${coverage[0]} \
        ${("--secondary-coverage " + ",".join([str(x) for x in coverage[1:]])) if len(coverage) > 1 else ""} \
        --pip-cutoff ${pip_cutoff} \
        --pip-cutoff-to-skip ${pip_cutoff_to_skip} \
        --seed ${seed} \
        ${("--contexts " + contexts) if contexts else ""} \
        --method-args '${mvsusie_args}' \
        ${("--twas-weights " + str(prior_twas_weights)) if prior_twas_weights.is_file() else ""} \
        --data-driven-prior-weights-cutoff ${prior_weights_min} \
        --output ${_output}

[mnm_genes_1]
# Resolve per-locus regions (region_id -> chr:start-end) into a manifest via
# region_manifest.R for the cross-gene (jointSpecification="trait") mvSuSiE fit.
# Loci come from --customized-association-windows (or --region-list); restrict
# with --region-name. Falls back to per-gene cis windows from the manifest.
input: None
output: f"{cwd:a}/multivariate_fine_mapping/{name}.mnm_genes_region_manifest.tsv"
task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f"{step_name}_{_output:bn}"
bash: expand = "${ }", stderr = f"{_output}.stderr", stdout = f"{_output}.stdout", container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/region_manifest.R \
        --pheno-manifest ${phenoFile[0]} \
        --cis-window ${cis_window if cis_window > 0 else 1000000} \
        ${("--customized-association-windows " + str(customized_association_windows)) if customized_association_windows.is_file() else ""} \
        ${("--region-name " + ",".join(region_name)) if len(region_name) > 0 else ""} \
        ${("--region-list " + str(region_list)) if region_list.is_file() else ""} \
        --output ${_output}

[mnm_genes_2]
# Cross-gene (multi-trait) mvSuSiE per locus: region mode +
# jointSpecification="trait" joins the genes whose coordinates overlap the
# locus. Optional mr.mash data-driven prior via --prior-twas-weights (canonical
# otherwise). Replaces the legacy multigene combine + multigene_udr machinery.
parameter: mvsusie_max_iter = 200
parameter: pip_cutoff_to_skip = 0.0
parameter: contexts = ""
parameter: prior_twas_weights = path()
parameter: prior_weights_min = 1e-10
parameter: mvsusie_method_args = ""
import csv
manifest = path(f"{cwd:a}/multivariate_fine_mapping/{name}.mnm_genes_region_manifest.tsv")
jobs = list(csv.DictReader(open(manifest), delimiter='\t'))
stop_if(len(jobs) == 0, "mnm_genes: empty region manifest; check --region-name / --customized-association-windows.")
qtl_dataset = path(f"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds")
mvsusie_args = mvsusie_method_args if mvsusie_method_args else f'{{"mvsusie":{{"max_iter":{mvsusie_max_iter}}}}}'
input: qtl_dataset, for_each = "jobs"
output: f"{cwd:a}/multivariate_fine_mapping/{name}.{_jobs['region_id']}.multigene_bvsr.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f"{step_name}_{_output:bn}"
bash: expand = "${ }", stderr = f"{_output:n}.stderr", stdout = f"{_output:n}.stdout", container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \
        --qtl-dataset ${_input} \
        --region ${_jobs['ld_block']} \
        --methods mvsusie \
        --joint-specification trait \
        --coverage ${coverage[0]} \
        ${("--secondary-coverage " + ",".join([str(x) for x in coverage[1:]])) if len(coverage) > 1 else ""} \
        --pip-cutoff ${pip_cutoff} \
        --pip-cutoff-to-skip ${pip_cutoff_to_skip} \
        --seed ${seed} \
        ${("--contexts " + contexts) if contexts else ""} \
        --method-args '${mvsusie_args}' \
        ${("--twas-weights " + str(prior_twas_weights)) if prior_twas_weights.is_file() else ""} \
        --data-driven-prior-weights-cutoff ${prior_weights_min} \
        --output ${_output}

[fsusie_1]
# Resolve per-locus regions into a manifest via region_manifest.R for the
# functional SuSiE fit (one locus = one fsusie unit; fsusieR::susiF jointly fits
# the multi-trait function within each context). Loci come from
# --customized-association-windows (or --region-list); restrict with --region-name.
input: None
output: f"{cwd:a}/fsusie/{name}.fsusie_region_manifest.tsv"
task: trunk_workers = 1, trunk_size = 1, walltime = walltime, mem = mem, cores = numThreads, tags = f"{step_name}_{_output:bn}"
bash: expand = "${ }", stderr = f"{_output}.stderr", stdout = f"{_output}.stdout", container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/region_manifest.R \
        --pheno-manifest ${phenoFile[0]} \
        --cis-window ${cis_window if cis_window > 0 else 1000000} \
        ${("--customized-association-windows " + str(customized_association_windows)) if customized_association_windows.is_file() else ""} \
        ${("--region-name " + ",".join(region_name)) if len(region_name) > 0 else ""} \
        ${("--region-list " + str(region_list)) if region_list.is_file() else ""} \
        --output ${_output}

[fsusie_2]
# Functional SuSiE per locus (fsusieR::susiF) over the pre-built QtlDataset.
# fsusie-specific knobs ride on --method-args; --susie-top-pc>0 additionally
# fine-maps each context's top principal components with univariate SuSiE
# (usePCA; ports the legacy fsusie.R susie_on_top_pc). Replaces fsusie.R.
parameter: prior = "mixture_normal"
parameter: max_scale = 10
parameter: min_purity = 0.5
parameter: post_processing = "TI"
parameter: max_SNP_EM = 100
parameter: small_sample_correction = False
# >0 -> also fine-map the top `susie_top_pc` PCs per context with univariate SuSiE.
parameter: susie_top_pc = 0
parameter: pip_cutoff_to_skip = 0.0
parameter: contexts = ""
# Optional JSON overriding the assembled fsusie kwargs, e.g. '{"fsusie":{"prior":"mixture_normal_per_scale"}}'.
parameter: fsusie_method_args = ""
import csv, json
manifest = path(f"{cwd:a}/fsusie/{name}.fsusie_region_manifest.tsv")
jobs = list(csv.DictReader(open(manifest), delimiter='\t'))
stop_if(len(jobs) == 0, "fsusie: empty region manifest; check --region-name / --customized-association-windows.")
qtl_dataset = path(f"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds")
fsusie_args = fsusie_method_args if fsusie_method_args else json.dumps({"fsusie": {"prior": prior, "max_scale": max_scale, "min_purity": min_purity, "post_processing": post_processing, "max_SNP_EM": max_SNP_EM, "cor_small": bool(small_sample_correction)}})
input: qtl_dataset, for_each = "jobs"
output: f"{cwd:a}/fsusie/{name}.{_jobs['region_id']}.fsusie.rds"
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f"{step_name}_{_output:bn}"
bash: expand = "${ }", stderr = f"{_output:n}.stderr", stdout = f"{_output:n}.stdout", container = container, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \
        --qtl-dataset ${_input} \
        --region ${_jobs['ld_block']} \
        --methods fsusie \
        --coverage ${coverage[0]} \
        ${("--secondary-coverage " + ",".join([str(x) for x in coverage[1:]])) if len(coverage) > 1 else ""} \
        --pip-cutoff ${pip_cutoff} \
        --pip-cutoff-to-skip ${pip_cutoff_to_skip} \
        --seed ${seed} \
        ${("--contexts " + contexts) if contexts else ""} \
        --method-args '${fsusie_args}' \
        ${("--use-pca --n-pcs " + str(susie_top_pc)) if susie_top_pc > 0 else ""} \
        --output ${_output}

mvfsusie (WIP placeholder)#

This section is still under development in the original pipeline and is included here as a placeholder only. It is not part of the runnable MWE.

[mvfsusie]
# prior can be either of ["mixture_normal", "mixture_normal_per_scale"]
parameter: prior  = "mixture_normal_per_scale"
parameter: max_SNP_EM = 1000
depends: sos_variable("regional_data")
# Check if both 'data' and 'meta_info' are empty lists
stop_if(len(regional_data['data']) == 0, f'Either genotype or phenotype data are not available for region {", ".join(region_name)}.')

meta_info = regional_data['meta_info']
input: regional_data["data"], group_by = lambda x: regional_data["data"], group_with = "meta_info"
output: f'{cwd:a}/{step_name[:-2]}/{name}.{_meta_info[0]}.mvfsusie_{prior}.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
R: expand = '${ }', stdout = f"{_output:n}.stdout", stderr = f"{_output:n}.stderr", container = container, entrypoint = entrypoint
    # Load regional association data
    fdat = load_regional_association_data(genotype = ${_input[0]:anr},
                                          phenotype = c(${",".join(['"%s"' % x.absolute() for x in _input[1::2]])}),
                                          covariate = c(${",".join(['"%s"' % x.absolute() for x in _input[2::2]])}),
                                          region = ${'"%s:%s-%s"' % (_meta_info[1], _meta_info[2], _meta_info[3])},
                                          maf_cutoff = ${maf},
                                          mac_cutoff = ${mac},
                                          imiss_cutoff = ${imiss})
    # Fine-mapping with mvfSuSiE
    library("mvf.susie.alpha")
    Y = map(fdat$residual_Y, ~left_join(fdat$X[,1]%>%as.data.frame%>%rownames_to_column("rowname"), .x%>%t%>%as.data.frame%>%rownames_to_column("rowname") , by = "rowname")%>%select(-2)%>%column_to_rownames("rowname")%>%as.matrix )
    fitted <- multfsusie(Y_f = list(Y[[1]],Y[[3]]), 
                         Y_u = Reduce(cbind, Y[[2]]),
                         pos = list(pos1 =fdat$phenotype_coordiates[[1]], pos2 = fdat$phenotype_coordiates[[3]]),
                         X=X,
                         L=${max_L},
                         data.format="list_df")
    saveRDS(fitted, ${_output:ar})

Troubleshooting#

Step

Problem

Possible Reason

Solution

susie_twas

A region is skipped / no window found

4th-column ID in the association-windows file does not match the phenotype list ID

Make the ID columns match exactly.

susie_twas

susie_twas: pass --region-name <gene-id>

No region selected

Pass --region-name <gene-id> (e.g. ENSG00000130538) or a --region-list.

susie_twas

ERROR: One of the local workers has been killed

Out-of-memory: too many regions in one process (-j4, or all 10 serially)

Use -j1; re-run heavy regions one at a time with a single-region phenotype list.

mnm

Target unavailable: ...qtl_dataset.rds

qtl_dataset_construct not run first

Chain it: qtl_dataset_construct+susie_twas, then mnm.

mnm

--phenotype-manifest missing required column(s): cond

Single-context manifest used

Use the multi-context manifest (pheno_manifest_multicontext.tsv).

both

No shared samples between phenotype and phenotype-covariate

Covariates are in QTLtools (transposed) layout

Add --transpose-covariates.

both

Sample IDs do not overlap between genotype/covariates and phenotype

Genotype .fam IDs differ from the phenotype matrix sample names

Rename one side so they match (e.g. the gene-expression mode renames genotype/covariate IDs to SAMPLE_001..060).