Advanced regression models for association analysis with individual-level data

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

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/qtl_mini/example_covariates.tsv \
    --customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
    --region-name ENSG00000130538 --transpose-covariates --save-data -j1

sos run pipeline/mnm_regression.ipynb qtl_dataset_construct+mnm \
    --name protocol_example --cwd output_mnm \
    --genoFile tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \
    --phenoFile tests/fixtures/qtl_mini/protocol_example.pheno_manifest_context.tsv \
    --covFile tests/fixtures/qtl_mini/example_covariates.tsv \
    --customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
    --transpose-covariates --save-data -j1

sos run pipeline/mnm_regression.ipynb susie_twas \
    --name protocol_example \
    --cwd output_10peaks \
    --genoFile tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \
    --phenoFile tests/fixtures/qtl_mini/pheno_manifest.tsv \
    --covFile tests/fixtures/covariate_hidden_factor/covariates.tsv \
    --customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
    -j1

sos run pipeline/mnm_regression.ipynb qtl_dataset_construct+mnm \
    --name protocol_example \
    --cwd output_10peaks \
    --genoFile tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \
    --phenoFile tests/fixtures/qtl_mini/pheno_manifest.tsv \
    --covFile tests/fixtures/covariate_hidden_factor/covariates.tsv \
    --customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
    -j1

sos run pipeline/mnm_regression.ipynb mnm \
    --name protocol_example --cwd output_mnm3 \
    --genoFile tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \
    --phenoFile tests/fixtures/qtl_mini/protocol_example.pheno_manifest_context.tsv \
    --covFile tests/fixtures/qtl_mini/example_covariates.tsv \
    --customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
    --save-data -j1

sos run pipeline/mnm_regression.ipynb mnm_genes \
    --name protocol_example --cwd output_mnmgenes \
    --genoFile tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \
    --phenoFile tests/fixtures/qtl_mini/protocol_example.pheno_manifest_context.tsv \
    --covFile tests/fixtures/qtl_mini/example_covariates.tsv \
    --customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
    --pheno_id_map_file tests/fixtures/qtl_mini/pheno_id_map.tsv \
    --keep-samples tests/fixtures/qtl_mini/keep_samples.txt \
    --save-data -j1

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

sos run pipeline/mnm_regression.ipynb mvfsusie \
    --name protocol_example --cwd output/mvfsusie \
    --genoFile tests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed \
    --phenoFile tests/fixtures/qtl_mini/pheno_manifest.tsv \
    --covFile tests/fixtures/covariate_hidden_factor/covariates.tsv \
    --customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
    --save-data -j1

sos run pipeline/mnm_regression.ipynb -h

[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 protocol_example.pheno_manifest_context.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]
# 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})