xQTL-GWAS pairwise enrichment and colocalization

Estimates a global enrichment between xQTL and GWAS SuSiE fine-mapping signals, then colocalises them on the overlapping regions.

Overview

This workflow processes fine-mapping results for xQTL, generated by susie_twas in the mnm_regression.ipynb notebook for cis xQTL, and GWAS fine-mapping results produced by susie_rss in the rss_analysis.ipynb notebook. It is designed to perform enrichment and colocalization analysis, particularly when fine-mapping results originate from different regions in the case of cis-xQTL and GWAS. The pipeline is capable to integrate and analyze data across these distinct regions. Originally tailored for cis-xQTL and GWAS integration, this pipeline can be applied to other pairwise integrations. An example of such application is in trans analysis, where the fine-mapped regions might be identical between trans-xQTL and GWAS, representing a special case of this broader implementation.

Fine-mapping gives credible sets per trait; it does not say whether an xQTL signal and a GWAS signal are the same signal. This step takes both sets of fine-mapping results and asks that question directly, using the full posterior rather than a single lead variant per locus, which is what makes it robust when a region contains more than one causal variant.

Enrichment

1. Identify GWAS Blocks: Select GWAS blocks with top loci from Bellenguez data and using a single variant regression method. Locate the original filenames (original_data) for these blocks and map the analysis regions to identify overlapping gene analysis regions and corresponding QTL files with top loci table (condition_top_loci).

2. Contexts in xQTL Meta-data: Find contexts within the xQTL meta-data that include top loci results for each gene. Retrieve the corresponding original xQTL files (original_data) by referencing the context meta-data.

3. Execute xQTL GWAS Enrichment: Perform xqtl_gwas_enrichment analysis to generate enrichment parameters (a0, a1) and calculate priors (p1, p2, p12).

Coloc

1. xQTL Meta-data Contexts: Within the xQTL meta-data, identify contexts that include top loci results for each gene, indicated by condition_top_loci. Retrieve the corresponding original xQTL files (original_data) by referring back to the context meta-data. (Optional) Subset the QTL table with specifies gene list.

2. Original GWAS Files: Identify GWAS blocks that contain overlapping top loci variants for each gene (block_top_loci) in above file. Utilize the GWAS list to locate the original filenames (original_data) for the identified GWAS block candidates.

3. Enrichment and Coloc Analysis: Automatically execute xqtl_gwas_enrichment to enable enrichment analysis, generating priors for subsequent coloc analysis through logic defined in enrichment. Then, apply susie_coloc for coloc analysis on the selected xQTL and GWAS original files, analyzing each gene under each identified condition.

When to run it. After susie_twas (in mnm_regression) for the xQTL side and susie_rss (in rss_analysis) for the GWAS side. Both must be fine-mapped first - this step consumes results, it does not produce them.

Input

Both workflows read lists of SuSiE fine-mapping output objects in RDS format, of class(susie) in R, described by two meta-files.

• --xqtl-meta-data input/susie_enloc_data/protocol_example.enloc.xqtl_meta.tsv (one row per xQTL analysis region. Columns chr, start, end, original_data, conditions_top_loci, block_top_loci, where original_data is an RDS file, conditions_top_loci shows which contexts have a top-loci table (potential signals), and block_top_loci gives the blocks whose overlapped top-loci variants coincide with xQTL data. For enrichment this comes from fine_mapping_post_processing/cis_results_export; for coloc from fine_mapping_post_processing/overlap_qtl_gwas. The difference is the last two columns, block_top_loci and final_combined_data, which document overlapped block information at variant level.)

#chr	start	end	region_id	TSS	original_data	combined_data	combined_data_sumstats	...
chr1	20520000	24080000	ENSG00000142798	21937309	protocol_example.enloc.MiGA_eQTL.ENSG00000142798.univariate_susie_twas_weights.rds,protocol_example.enloc.KNIGHT_pQTL.ENSG00000142798.univariate_susie_twas_weights.rds	protocol_example.enloc.qtl.ENSG00000142798.cis_results_db.export.rds	protocol_example.enloc.qtl.ENSG00000142798.cis_results_db.export_sumstats.rds	...

• --gwas-meta-data input/susie_enloc_data/protocol_example.enloc.gwas_meta.tsv (one row per GWAS analysis region. Columns chr, start, end, original_data, block_top_loci and onwards, where original_data is an RDS file. Output of fine_mapping_post_processing/gwas_results_export.)

#chr	start	end	region_id	TSS	original_data	combined_data	combined_data_sumstats	...
chr1	17351816	20110062	chr1_17351816_20110062	NA	protocol_example.enloc.RSS_QC_RAISS_imputed.chr1_17351816_20110062.univariate_susie_rss.rds	protocol_example.enloc.gwas.chr1_17351816_20110062.cis_results_db.export.rds	protocol_example.enloc.gwas.chr1_17351816_20110062.cis_results_db.export_sumstats.rds	...
chr1	20110062	22020160	chr1_20110062_22020160	NA	protocol_example.enloc.RSS_QC_RAISS_imputed.chr1_20110062_22020160.univariate_susie_rss.rds	protocol_example.enloc.gwas.chr1_20110062_22020160.cis_results_db.export.rds	protocol_example.enloc.gwas.chr1_20110062_22020160.cis_results_db.export_sumstats.rds	...

• --context-meta input/susie_enloc_data/protocol_example.enloc.context_meta.tsv (maps each analysis name to its cohort and the contexts it contains, so the raw data for each context can be identified)

analysis_name	cohort	context
KNIGHT_pQTL	KNIGHT	Knight_pQTL_brain,Knight_eQTL_brain
MSBB_eQTL	MSBB	BM_10_MSBB_eQTL,BM_22_MSBB_eQTL,BM_36_MSBB_eQTL,BM_44_MSBB_eQTL

• --qtl-path input/susie_enloc_data (directory holding the xQTL RDS objects named in the xQTL meta-file)

• --gwas-path input/susie_enloc_data (directory holding the GWAS RDS objects named in the GWAS meta-file)

The -obj flags are not paths: each names the nested keys to follow inside the RDS object to reach the value the pipeline needs.

• --xqtl-finemapping-obj preset_variants_result susie_result_trimmed (path within the xQTL object to the SuSiE fit)

• --xqtl-varname-obj preset_variants_result variant_names (path within the xQTL object to the variant names)

• --xqtl-region-obj region_info grange (path within the xQTL object to the region definition)

• --gwas-finemapping-obj AD_Bellenguez_2022 RSS_QC_RAISS_imputed susie_result_trimmed (path within the GWAS object to the SuSiE fit)

• --gwas-varname-obj AD_Bellenguez_2022 RSS_QC_RAISS_imputed variant_names (path within the GWAS object to the variant names)

susie_coloc takes two further options:

• --ld-meta-file-path input/ld_reference/protocol_example.ld_meta_file.tsv (LD reference index used by coloc postprocessing. Columns chr, start, end, path. When supplied, coloc credible sets are also reported.)

#chr	start	end	path
chr1	16103	2888443	chr1/protocol_example.LD.chr1
chr1	2888443	4320284	chr1/protocol_example.LD.chr1

• --skip-enrich (omit the enrichment prior. Without it, coloc reads <cwd>/<name>.<context>.enrichment.rds to tune its prior, so xqtl_gwas_enrichment must have been run first.)

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

Output

• <cwd>/<name>.<context>.enrichment.rds (xqtl_gwas_enrichment) (one global enrichment estimate per GWAS and xQTL-context pair, combining all input data for that context)

enrich <- readRDS("output/xqtl_gwas_enrichment/protocol_example.enloc.protocol_example.enloc.Knight_eQTL_brain.enrichment.rds")

str(enrich)
# data.frame:  1 obs. of  6 variables:
# $ gwasStudy         : chr "AD_Bellenguez_2022"
# $ qtlStudy          : chr "protocol_example.enloc.protocol_example.enloc"
# $ qtlContext        : chr "Knight_eQTL_brain"
# $ enrichment        : num NA
# $ enrichmentSe      : num NA
# $ enrichmentLogOdds : num NA

Note on toy data: one row per gwas/qtl-context pair. Enrichment values are NA here because a single gene x single GWAS block does not give the underlying logistic-regression enrichment model enough data to fit; this is expected for the toy MWE, not a pipeline error. Reliable enrichment estimates require many gene x GWAS-region pairs.

• <cwd>/susie_coloc/<name>.<context>@<gene>.coloc.rds (susie_coloc) (colocalisation results for the regions of interest, one file per context and gene. When --ld-meta-file-path is given, coloc credible sets are reported too.)

res <- readRDS("output/susie_coloc/susie_coloc/protocol_example.enloc.Knight_eQTL_brain@ENSG00000142798.coloc.rds")

length(res)   # 1
names(res)    # [1] "ENSG00000142798"

str(res, max.level = 2)
# List of 1
#  $ ENSG00000142798:List of 2
#   ..$ (unnamed)      : chr "No coloc results due to the absence of a GWAS log Bayes factor matrix filtered by prior tolerance."
#   ..$ analysis_region: chr "chr1_20110062_22020160"

res$ENSG00000142798$analysis_region   # "chr1_20110062_22020160"

The output is a named list with one entry per gene. Each gene entry contains:

• [[1]] (unnamed): the colocalization result — either a data.frame of PP values (positive result) or a character string explaining why coloc was skipped.

• $analysis_region: the genomic block ID (e.g. "chr1_20110062_22020160") that was analyzed.

When colocalization runs successfully, res[[gene]][[1]] is a data.frame with one row per SuSiE effect pair (xQTL effect × GWAS effect) and columns:

Column	Description
nsnps	Number of variants in the shared analysis window
PP.H0.abf	Posterior probability: neither trait has a causal variant in this region
PP.H1.abf	Posterior probability: only xQTL has a causal variant
PP.H2.abf	Posterior probability: only GWAS has a causal variant
PP.H3.abf	Posterior probability: both have causal variants, but different ones
PP.H4.abf	Posterior probability: both traits share the same causal variant (colocalization)
hit1	Top variant ID from the xQTL credible set
hit2	Top variant ID from the GWAS credible set

PP.H0 + PP.H1 + PP.H2 + PP.H3 + PP.H4 = 1. PP.H4.abf > 0.8 is the standard threshold for strong colocalization evidence.

Note on toy data: this example returns a message string instead of a data.frame because the toy GWAS fine-mapping signal in this region is too weak (max PIP ≈ 0.04, no credible set). This is an expected negative result for a toy dataset, not a pipeline error.

Minimal Working Example

The two workflows are independent as run here. xqtl_gwas_enrichment estimates the global enrichment and susie_coloc colocalises, and the command below passes --skip-enrich so coloc does not read the enrichment estimate. On real data, dropping that flag lets coloc pick up <name>.<context>.enrichment.rds and use it to tune its prior; SoS will not run enrichment for you, so run it first if you want that.

Global xQTL-GWAS enrichment

xqtl_gwas_enrichment estimates one enrichment value per GWAS and xQTL-context pair, combining all input data for that context.

Timing: TBD (on toy dataset)

sos run pipeline/SuSiE_enloc.ipynb xqtl_gwas_enrichment \
    --gwas-meta-data input/susie_enloc_data/protocol_example.enloc.gwas_meta.tsv \
    --xqtl-meta-data input/susie_enloc_data/protocol_example.enloc.xqtl_meta.tsv \
    --xqtl-finemapping-obj preset_variants_result susie_result_trimmed \
    --xqtl-varname-obj preset_variants_result variant_names \
    --gwas-finemapping-obj AD_Bellenguez_2022 RSS_QC_RAISS_imputed susie_result_trimmed \
    --gwas-varname-obj  AD_Bellenguez_2022 RSS_QC_RAISS_imputed variant_names \
    --xqtl-region-obj   region_info grange \
    --qtl-path input/susie_enloc_data \
    --gwas-path input/susie_enloc_data \
    --context-meta input/susie_enloc_data/protocol_example.enloc.context_meta.tsv \
    --cwd output/xqtl_gwas_enrichment

Pairwise colocalisation

susie_coloc runs colocalisation on the overlapping regions, analysing each gene under each identified condition.

Timing: TBD (on toy dataset)

sos run pipeline/SuSiE_enloc.ipynb susie_coloc \
    --gwas-meta-data input/susie_enloc_data/protocol_example.enloc.gwas_meta.tsv \
    --xqtl-meta-data input/susie_enloc_data/protocol_example.enloc.xqtl_meta.tsv \
    --xqtl-finemapping-obj preset_variants_result susie_result_trimmed \
    --xqtl-varname-obj preset_variants_result variant_names \
    --gwas-finemapping-obj AD_Bellenguez_2022 RSS_QC_RAISS_imputed susie_result_trimmed \
    --gwas-varname-obj  AD_Bellenguez_2022 RSS_QC_RAISS_imputed variant_names \
    --xqtl-region-obj   region_info grange \
    --qtl-path input/susie_enloc_data \
    --gwas-path input/susie_enloc_data \
    --context-meta input/susie_enloc_data/protocol_example.enloc.context_meta.tsv \
    --ld-meta-file-path input/ld_reference/protocol_example.ld_meta_file.tsv \
    --skip-enrich \
    --cwd output/susie_coloc

Command Interface

sos run pipeline/SuSiE_enloc.ipynb -h

usage: sos run <path>/xqtl-protocol/code/pecotmr_integration/SuSiE_enloc.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:
  get_analysis_regions
  get_overlapped_analysis_regions
  xqtl_gwas_enrichment
  susie_coloc

Global Workflow Options:
  --cwd output (as path)
                        Workdir
  --gwas-meta-data VAL (as path, required)
                        A list of file paths for fine-mapped GWAS results.
  --xqtl-meta-data VAL (as path, required)
                        A list of file paths for fine-mapped xQTL results.
  --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
  --name  f"{xqtl_meta_data:bnn}.{gwas_meta_data:bnn}"

                        It is required to input the name of the analysis
  --container ''
  --job-size 200 (as int)
                        For cluster jobs, number commands to run per job
  --walltime 5m
                        Wall clock time expected
  --mem 16G
                        Memory expected: quite large for enrichment analysis but small for xQTL colocalization
  --numThreads 1 (as int)
                        Number of threads
  --xqtl-finemapping-obj  (as list)
                        Optional table name in xQTL RDS files to get fine mapping results (eg 'susie_result_trimmed ').
  --xqtl-varname-obj  (as list)
                        Optional table name in xQTL RDS files to get variable names (eg ' variant_names ').
  --gwas-finemapping-obj  (as list)
                        Optional table name in GWAS RDS files to get fine mapping results (eg 'AD_Bellenguez_2022 single_effect_regression
                        susie_result_trimmed ').
  --gwas-varname-obj  (as list)
                        Optional table name in GWAS RDS files to get variable names(eg 'AD_Bellenguez_2022 single_effect_regression
                        variant_names ').
  --xqtl-region-obj  (as list)
                        Optional table name in xQTL RDS files to get region info (eg 'susie_result_trimmed region_info grange ').
  --gwas-region-obj  (as list)
                        Optional table name in GWAS RDS files to get region info (eg 'AD_Bellenguez_2022 single_effect_regression region_info
                        grange ').
  --gwas-path ''
                        Directory path for GWAS orignal finemapping results
  --qtl-path ''
                        Directory path for xQTL orignal finemapping results
  --context-meta . (as path)
                        a meta file showing the context and corresponding analysis_name
  --[no-]minp (default to False)
                        use conditions_top_loci_minp column or not, which can reduce a lot computing resource by pick top 1 isoform

Sections
  get_analysis_regions: get overlapped regions by gene and block region for enrichment analysis
  get_overlapped_analysis_regions: get overlapped regions from flatten table, to get the regions with overlapped variants and select those
                        regions for coloc analysis
  xqtl_gwas_enrichment:
  susie_coloc:
    Workflow Options:
      --[no-]skip-enrich (default to False)
                        depends: sos_step("xqtl_gwas_enrichment") #changed skip enrichment or not
      --[no-]filter-lbf-cs (default to False)
                        filter lbf by cs as default in susie coloc. default is False means filtering by V > prior_tol
      --ld-meta-file-path . (as path)
                        ld reference path for postprocessing

Workflow implementation

[global]
# Workdir
parameter: cwd = path("output")
parameter: modular_script_dir = path('code/script')
# A list of file paths for fine-mapped GWAS results. 
parameter: gwas_meta_data = path
# A list of file paths for fine-mapped xQTL results. 
parameter: xqtl_meta_data = path
# Optional: if a region list is provide the enrichment analysis will be focused on provided region. 
# The LAST column of this list will contain the ID of regions to focus on
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 = []
# It is required to input the name of the analysis
parameter: name = f"{xqtl_meta_data:bnn}.{gwas_meta_data:bnn}"
parameter: container = ""
# For cluster jobs, number commands to run per job
parameter: job_size = 200
# Wall clock time expected
parameter: walltime = "5m"
# Memory expected: quite large for enrichment analysis but small for xQTL colocalization
parameter: mem = "16G"
# Number of threads
parameter: numThreads = 1
#Optional table name in xQTL RDS files to get fine mapping results (eg 'susie_result_trimmed ').
parameter: xqtl_finemapping_obj = []
#Optional table name in xQTL RDS files to get variable names (eg ' variant_names ').
parameter: xqtl_varname_obj = []
#Optional table name in GWAS RDS files to get fine mapping results (eg 'AD_Bellenguez_2022 single_effect_regression susie_result_trimmed ').
parameter: gwas_finemapping_obj = []
#Optional table name in GWAS RDS files to get variable names(eg 'AD_Bellenguez_2022 single_effect_regression variant_names ').
parameter: gwas_varname_obj = []
#Optional table name in xQTL RDS files to get region info (eg 'susie_result_trimmed region_info grange ').
parameter: xqtl_region_obj = []
#Optional table name in GWAS RDS files to get region info (eg 'AD_Bellenguez_2022 single_effect_regression region_info grange ').
parameter: gwas_region_obj = []
#Directory path for GWAS orignal finemapping results 
parameter: gwas_path = ''
#Directory path for xQTL orignal finemapping results 
parameter: qtl_path = ''
# a meta file showing the context and corresponding analysis_name
parameter: context_meta = 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
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 = []
# use conditions_top_loci_minp column or not, which can reduce a lot computing resource by pick top 1 isoform
parameter: minp = False

[get_analysis_regions: shared = {'enrichment_rows': "list(__import__('csv').DictReader(open(enrichment_manifest), delimiter=chr(9)))"}]
# Resolve the per-condition xQTL x GWAS enrichment units into a manifest TSV via
# enloc_manifest.R -- the coordinate-overlap pairing + per-condition study
# routing that the legacy get_analysis_regions did in-notebook with pandas. The
# parsed rows are shared as enrichment_rows, so [xqtl_gwas_enrichment]
# auto-triggers this step via depends: sos_variable -- the notebook is invoked
# exactly as before (e.g. sos run SuSiE_enloc.ipynb xqtl_gwas_enrichment).
input: None
enrichment_manifest = f"{cwd:a}/get_analysis_regions/{name}.enrichment_manifest.tsv"
bash: expand = "${ }", container = container
    mkdir -p ${cwd:a}/get_analysis_regions
    Rscript ${modular_script_dir}/pecotmr_integration/enloc_manifest.R --mode enrichment \
        --xqtl-meta ${xqtl_meta_data} --gwas-meta ${gwas_meta_data} \
        ${("--context-meta " + str(context_meta)) if context_meta.is_file() else ""} \
        ${("--qtl-path " + str(qtl_path)) if qtl_path else ""} \
        ${("--gwas-path " + str(gwas_path)) if gwas_path else ""} \
        ${("--region-list " + str(region_list)) if region_list.is_file() else ""} \
        ${("--region-name " + ",".join(region_name)) if len(region_name) > 0 else ""} \
        ${("--gwas-finemapping-obj " + " ".join([str(x) for x in gwas_finemapping_obj])) if len(gwas_finemapping_obj) > 0 else ""} \
        ${"--minp" if minp else ""} \
        --output ${enrichment_manifest}

[get_overlapped_analysis_regions: shared = {'coloc_rows': "list(__import__('csv').DictReader(open(coloc_manifest), delimiter=chr(9)))"}]
# Resolve the per-(condition,region) xQTL x GWAS coloc units into a manifest TSV
# via enloc_manifest.R --mode coloc -- the block_top_loci variant-overlap pairing
# + per-condition study routing that the legacy get_overlapped_analysis_regions
# did in-notebook with pandas. Shared as coloc_rows, so [susie_coloc]
# auto-triggers this step (sos run SuSiE_enloc.ipynb susie_coloc) unchanged.
input: None
coloc_manifest = f"{cwd:a}/get_overlapped_analysis_regions/{name}.coloc_manifest.tsv"
bash: expand = "${ }", container = container
    mkdir -p ${cwd:a}/get_overlapped_analysis_regions
    Rscript ${modular_script_dir}/pecotmr_integration/enloc_manifest.R --mode coloc \
        --xqtl-meta ${xqtl_meta_data} --gwas-meta ${gwas_meta_data} \
        ${("--context-meta " + str(context_meta)) if context_meta.is_file() else ""} \
        ${("--qtl-path " + str(qtl_path)) if qtl_path else ""} \
        ${("--gwas-path " + str(gwas_path)) if gwas_path else ""} \
        ${("--region-list " + str(region_list)) if region_list.is_file() else ""} \
        ${("--region-name " + ",".join(region_name)) if len(region_name) > 0 else ""} \
        ${("--gwas-finemapping-obj " + " ".join([str(x) for x in gwas_finemapping_obj])) if len(gwas_finemapping_obj) > 0 else ""} \
        ${"--minp" if minp else ""} \
        --output ${coloc_manifest}

[xqtl_gwas_enrichment]
depends: sos_variable("enrichment_rows")
jobs = enrichment_rows
stop_if(len(jobs) == 0, f'No files left for analysis')
_qtl_groups = [j['qtl_files'].split(',') for j in jobs]
input: [f for grp in _qtl_groups for f in grp], group_by = lambda x: _qtl_groups, group_with = "jobs"
output: f'{cwd:a}/{name}.{_jobs["unit_id"]}.enrichment.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = '${ }', stdout = f"{_output:n}.stdout", stderr = f"{_output:n}.stderr", container = container
    set -e
    # The xQTL and GWAS sides are already S4 FineMappingResult collections
    # (fine_mapping output); qtl_enrichment.R combines its multi-file inputs
    # itself, so the legacy legacy_enloc_finemap_convert.R step is gone.
    Rscript ${modular_script_dir}/pecotmr_integration/qtl_enrichment.R \
        --qtl-fine-mapping ${" ".join([str(x) for x in _input])} \
        --gwas-fine-mapping ${_jobs["gwas_files"].replace(",", " ")} \
        --ncore ${numThreads} \
        --output ${_output}

[susie_coloc]
depends: sos_variable("coloc_rows")
# depends: sos_step("xqtl_gwas_enrichment") #changed
stop_if(len(coloc_rows) == 0, f'No files left for analysis')
# skip enrichment or not
parameter: skip_enrich=False
# filter lbf by cs as default in susie coloc. default is False means filtering by V > prior_tol
parameter: filter_lbf_cs = False
# ld reference path for postprocessing 
parameter: ld_meta_file_path = path()
jobs = coloc_rows
_qtl_groups = [j['qtl_files'].split(',') for j in jobs]
input: [f for grp in _qtl_groups for f in grp], group_by = lambda x: _qtl_groups, group_with = "jobs"
output: f'{cwd:a}/{step_name}/{name}.{_jobs["unit_id"]}.coloc.rds'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output:bn}'
bash: expand = '${ }', stdout = f"{_output:n}.stdout", stderr = f"{_output:n}.stderr", container = container
    set -e
    # xQTL + GWAS sides are already S4 FineMappingResult collections; coloc.R
    # combines its multi-file inputs itself (no legacy conversion). The enrichment
    # file (context before '@' in unit_id) tunes colocPipeline's prior.
    Rscript ${modular_script_dir}/pecotmr_integration/coloc.R \
        --qtl-fine-mapping ${" ".join([str(x) for x in _input])} \
        --gwas-input ${_jobs["gwas_files"].replace(",", " ")} \
        ${('--filter-lbf-cs' if filter_lbf_cs else '')} \
        ${('--enrichment '+cwd.absolute().joinpath(name + "." + _jobs["unit_id"].split("@")[0] + ".enrichment.rds").as_posix()) if not skip_enrich else ''} \
        --output ${_output}