Multi-trait colocalization using ColocBoost

Multi-trait colocalization across molecular contexts and against GWAS, using ColocBoost.

Overview

Two signals at the same locus - a QTL and a GWAS hit, or QTLs in two cell types - may share a causal variant or merely sit in the same LD block. Colocalization tries to tell those apart. Classical pairwise methods assume one causal variant per trait and compare traits two at a time, which loses power when a region has several independent signals or when the shared signal is weak in any single pair.

ColocBoost treats it as a multi-task problem instead: a gradient boosting framework that couples traits as it selects causal variants, so evidence that is weak in isolation can still support a shared signal across many contexts. It scales to hundreds of traits and allows multiple causal variants per region (Cao et al., 2025).

Requires fine-mapping outputs (.susie.rds files) from mnm_regression or rss_analysis.

ColocBoost performs colocalization analyses (Cao et al., 2025). In short, ColocBoost performs colocalization while accounting for multiple causal variants within a genomic region of interest and can scale to hundreds of traits. GWAS summary statistics may be included or omitted. Required inputs are individual xQTL data from the same cohort (multiple phenotypes, same genotype). Linkage disequilibrium reference data is required if GWAS summary statistics are used.

When to run it. After QTL association testing, when you have several contexts to compare, or a GWAS to integrate against. Colocalization is integration, not fine-mapping: run mnm_regression if what you want is credible sets for one trait.

Input

Both entry points read per-region fine-mapping outputs (.susie.rds) produced by mnm_regression or rss_analysis. Individual-level genotype data must be in PLINK bed format; summary statistics come in through the GWAS and LD meta-data files instead.

• --name colocboost_xqtl (tag used in output filenames and as the study name)

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

Individual-level xQTL inputs.

• --genoFile input/colocboost/example.chr22.bed (genotype data in PLINK bed format; either one file for the whole genome or a list of per-chromosome files)

• --phenoFile input/colocboost/pheno_manifest_multicontext.tsv (a manifest of phenotype files per region, one row per region and context, in UCSC bed.gz format with a bed.gz.tbi index)

• --covFile input/colocboost/example_covariates.tsv (covariate files matching the phenotype list)

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

Region selection.

• --customized-association-windows input/colocboost/association_windows.bed (association windows for variants, cis or trans. Without it a fixed-size window of --cis-window (default 1000000) is used around each region.)

• --region-name ENSG00000130538 (restrict the analysis to named regions. If dropped, the union of all regions across all phenotype region lists is analyzed.)

• --region-list (a file of regions to analyze; the last column must be the region name)

Summary-statistics inputs (GWAS integration).

• --gwas-meta-data (GWAS summary statistics meta-data file)

• --ld-meta-data (LD reference meta-data file. See rss_analysis for the expected formats.)

Analysis mode.

• --xqtl-coloc (colocalize across the xQTL contexts; on by default)

• --separate-gwas / --no-separate-gwas (treat each GWAS as its own outcome, or pool them. --separate-gwas is the default; the xQTL-only example passes --no-separate-gwas.)

• --no-indel (drop indel variants. Applies to both individual-level genotype loading and summary-statistics QC; indels are kept by default.)

Output

• <cwd>/qtl_dataset/<study>.qtl_dataset.rds (colocboost_1) (the assembled QTL dataset for the study, built once and reused by the colocalization steps. Not present in the example output.)

• <cwd>/colocboost/<name>.colocboost_manifest.tsv (colocboost_2) (manifest listing the per-region result files produced for this run. Not present in the example output.)

• <cwd>/colocboost/<name>.<region_id>.colocboost.rds (colocboost_3) (the fitted ColocBoost model, one RDS per analysis region – per gene of interest, and per GWAS when GWAS integration is used)

List of 1
$ ENSG00000283047: List of 7  <colocboost>
 ..$ ucos_summary: NULL
 ..$ vpa: NULL
 ..$ data_info: List of 6
 .. ..$ n_outcomes: integer [1:1] 1
 .. ..$ n_variables: integer [1:1] 184
 .. ..$ outcome_info: List of 4  <data.frame>
 .. .. ..$ outcome_names: ...
 .. .. ..$ sample_size: ...
 .. .. ..$ is_sumstats: ...
 .. .. ..$ is_focal: ...
 .. ..$ variables: character [1:184] chr22:10685239:C:T chr22:10761227:G:A chr22:11285239:C:T chr22:11601805:G:A chr22:11791698:G:T ...
 .. ..$ coef: List of 1
 .. .. ..$ trait_A_ENSG00000283047: ...
 .. ..$ z: List of 1
 .. .. ..$ trait_A_ENSG00000283047: ...
 ..$ model_info: List of 9
 .. ..$ model_coveraged: logical [1:1] TRUE
 .. ..$ n_updates: integer [1:1] 352
 .. ..$ profile_loglik: numeric [1:352] 1 0.973943746514459 0.964750650763505 0.955751398118438 0.946916816050365 ...
 .. ..$ outcome_profile_loglik: List of 1
 .. .. ..$ trait_A_ENSG00000283047: ...
 .. ..$ outcome_proximity_obj: List of 1
 .. .. ..$ trait_A_ENSG00000283047: ...
 .. ..$ outcome_coupled_best_update_obj: List of 1
 .. .. ..$ trait_A_ENSG00000283047: ...
 .. ..$ outcome_model_coveraged: List of 1
 .. .. ..$ trait_A_ENSG00000283047: ...
 .. ..$ outcome_n_updates: integer [1:1] 352
 .. ..$ jk_star: matrix [351, 1] 125 125 125 131 125 ...
 ..$ ucos_details: NULL
 ..$ region_info: List of 3
 .. ..$ region_coord: List of 3  <data.frame>
 .. .. ..$ chrom: ...
 .. .. ..$ start: ...
 .. .. ..$ end: ...
 .. ..$ grange: List of 3  <data.frame>
 .. .. ..$ chrom: ...
 .. .. ..$ start: ...
 .. .. ..$ end: ...
 .. ..$ region_name: character [1:2] ENSG00000283047 ENSG00000283047
 ..$ computing_time: List of 3
 .. ..$ Loading: difftime [1:1] 10.3424470424652  <difftime>
 .. ..$ QC: difftime [1:1] 0.0479960441589355  <difftime>
 .. ..$ Analysis: List of 3
 .. .. ..$ xqtl_coloc: ...
 .. .. ..$ joint_gwas: ...
 .. .. ..$ separate_gwas: ...

The example files on disk are named <name>.<region_id>.cb_xqtl.rds, not .colocboost.rds as the step declares. The code is authoritative; the example data predates the current naming. In this example ucos_summary and vpa are NULL because a single outcome yields no colocalization events.

Each RDS contains:

• ucos_summary – a summary table of the colocalization events.

• vpa – the variable colocalized probability for each variant (the probability of a variant being colocalized with at least one trait).

• data_info – information on the input data.

• model_info – information on the fitted ColocBoost model.

• ucos_details – trait-specific (uncolocalized) effects information.

• region_info – information on the analysis region.

Minimal Working Example

The phenotype and covariate examples are duplicated below to show that when there are multiple phenotypes for the same genotype it is possible to use this pipeline to analyze all of them (more than two is accepted as well).

Here --region-name focuses the analysis on 3 genes. In practice if this parameter is dropped, the union of all regions in all phenotype region lists will be analyzed. It is possible for some of the regions there are no genotype data, in which case the pipeline will output RDS files with a warning message to indicate the lack of genotype data to analyze.

Note: Suggested output naming convention is cohort_modality, eg ROSMAP_snRNA_pseudobulk.

ColocBoost with xQTL only

This is the minimal working example and runs end-to-end on the toy gene-expression data, producing a *.cb_xqtl.rds colocalization result for the region. The same command also accepts peaks or any other molecular phenotype - just swap in the phenotype .bed.gz/manifest.

Timing: TBD (on toy dataset)

sos run pipeline/colocboost.ipynb colocboost \
    --name colocboost_xqtl --cwd output/colocboost_xqtl \
    --genoFile input/colocboost/example.chr22.bed \
    --phenoFile input/colocboost/pheno_manifest_multicontext.tsv \
    --covFile input/colocboost/example_covariates.tsv --transpose-covariates \
    --customized-association-windows input/colocboost/association_windows.bed \
    --region-name ENSG00000130538 \
    --no-separate-gwas --xqtl-coloc -j1

ColocBoost with GWAS

Adds GWAS summary statistics (via --gwas-meta-data) and an LD reference (via --ld-meta-data). The xQTL ColocBoost completes, but the separate GWAS-xQTL colocalization is still under development with the installed package (see the note in the command cell).

Timing: TBD (on toy dataset)

sos run pipeline/colocboost.ipynb colocboost \
    --name colocboost_gwas --cwd output/colocboost_gwas \
    --genoFile input/colocboost/example.chr22.bed \
    --phenoFile input/colocboost/pheno_manifest_multicontext.tsv \
    --covFile input/colocboost/example_covariates.tsv --transpose-covariates \
    --customized-association-windows input/colocboost/association_windows.bed \
    --gwas-meta-data input/colocboost/gwas_meta.txt \
    --ld-meta-data input/ld_reference/protocol_example.ld_meta_file.tsv \
    --region-name ENSG00000130538 \
    --separate-gwas --xqtl-coloc -j1

ColocBoost with multiple LD references (multi-ancestry)

Uses a per-study LD reference declared in the GWAS meta-data. It shares the separate-GWAS code path with Step 2 and is therefore still under development on the toy dataset.

Timing: TBD (on toy dataset)

sos run pipeline/colocboost.ipynb colocboost \
    --name colocboost_multi_ld --cwd output/colocboost_multi_ld \
    --genoFile input/colocboost/example.chr22.bed \
    --phenoFile input/colocboost/pheno_manifest_multicontext.tsv \
    --covFile input/colocboost/example_covariates.tsv --transpose-covariates \
    --customized-association-windows input/colocboost/association_windows.bed \
    --gwas-meta-data input/colocboost/gwas_meta.txt \
    --ld-meta-data input/ld_reference/protocol_example.ld_meta_file.tsv \
    --region-name ENSG00000130538 \
    --separate-gwas --xqtl-coloc -j1

Command Interface

sos run pipeline/colocboost.ipynb -h

usage: sos run pipeline/colocboost.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:
  colocboost

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)
                        PLINK1 genotype: pass the .bed file (.bim/.fam are read
                        automatically).
  --phenoFile  paths

                        QtlDataset phenotype manifest TSV. Columns: ID, #chr,
                        start, end, path, cond, and an optional cov_path (per-
                        context covariates). The first path is used. Replaces
                        the legacy list-of-phenotype-BEDs input.
  --phenoIDFile  paths()

                        Legacy phenotype ID-mapping file(s); accepted for CLI
                        compatibility but not consumed by the manifest-based
                        flow (the manifest carries the IDs).
  --covFile  paths()

                        Uniform covariate file applied across all contexts
                        (QTLtools-format rows=covariates supported via
                        --transpose-covariates). Per-context covariates instead
                        come from the manifest cov_path column.
  --gwas-meta-data . (as path)
                        GWAS summary-statistics meta TSV. Columns: study_id,
                        chrom, file_path; optional n_sample, n_case, n_control,
                        ld_meta_data, column_mapping_file.
  --ld-meta-data . (as path)
                        Default LD-reference meta (#chr, start, end, path). Per-
                        study overrides come from the gwas meta ld_meta_data
                        column.
  --region-list . (as path)
                        Optional region list whose LAST column lists region/gene
                        IDs to analyze.
  --region-name  (as list)
                        Region/gene IDs to analyze (space separated on the CLI).
  --customized-association-windows . (as path)
                        Optional per-gene association windows (BED: #chr start
                        end ID) overriding the default cis_window for the GWAS
                        extraction window.
  --cis-window 1000000 (as int)
                        cis window (bp) on each side of the gene for fine-
                        mapping and the default GWAS extraction window.
  --[no-]save-data (default to False)
                        Save the per-gene colocboost input bundle (currently
                        unused; reserved).
  --study study
                        QtlDataset study label and covariate handling (legacy
                        CLI retained).
  --[no-]transpose-covariates (default to False)
  --genotype-covariates . (as path)
  --imiss 1.0 (as float)
                        Remove a variant if it has more than imiss fraction
                        missing genotypes.
  --maf 0.0025 (as float)
                        MAF and variance-of-X cutoffs.
  --xvar-cutoff 0.0 (as float)
  --mac 5 (as int)
                        MAC cutoff, on top of the MAF cutoff.
  --[no-]indel (default to True)
                        Keep indels when True.
  --keep-samples . (as path)
                        Only keep a subset of samples / variants
                        (whitespace/newline-separated IDs).
  --keep-variants . (as path)
  --skip-analysis-pip-cutoff  (as list)
                        Per-context single-effect skip cutoff passed to
                        colocboost.R (--pip-cutoff-to-skip): a scalar applied to
                        every context, or one or more context=value items
                        (negative -> data-driven 3/ncol(X); empty -> 0).
  --skip-sumstats-analysis-pip-cutoff -1.0 (as float)
                        Summary-stats single-effect skip cutoff passed to
                        summaryStatsQc via gwas_sumstats_construct.R.
  --[no-]xqtl-coloc (default to True)
                        ColocBoost variants to run.
  --[no-]joint-gwas (default to False)
  --[no-]separate-gwas (default to True)
  --qc-method none
                        Summary-stats LD-mismatch QC method (none | slalom |
                        dentist).
  --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 of commands to run per job.
  --walltime 1h
  --mem 20G
  --numThreads 1 (as int)

Sections
  colocboost_1:
  colocboost_2:
  colocboost_3:

Workflow implementation

[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")
# PLINK1 genotype: pass the .bed file (.bim/.fam are read automatically).
parameter: genoFile = path
# QtlDataset phenotype manifest TSV. Columns: ID, #chr, start, end, path,
# cond, and an optional cov_path (per-context covariates). The first path is
# used. Replaces the legacy list-of-phenotype-BEDs input.
parameter: phenoFile = paths
# Legacy phenotype ID-mapping file(s); accepted for CLI compatibility but not
# consumed by the manifest-based flow (the manifest carries the IDs).
parameter: phenoIDFile = paths()
# Uniform covariate file applied across all contexts (QTLtools-format
# rows=covariates supported via --transpose-covariates). Per-context
# covariates instead come from the manifest cov_path column.
parameter: covFile = paths()
# GWAS summary-statistics meta TSV. Columns: study_id, chrom, file_path;
# optional n_sample, n_case, n_control, ld_meta_data, column_mapping_file.
parameter: gwas_meta_data = path()
# Default LD-reference meta (#chr, start, end, path). Per-study overrides come
# from the gwas meta ld_meta_data column.
parameter: ld_meta_data = path()
# Optional region list whose LAST column lists region/gene IDs to analyze.
parameter: region_list = path()
# Region/gene IDs to analyze (space separated on the CLI).
parameter: region_name = []
# Optional per-gene association windows (BED: #chr start end ID) overriding the
# default cis_window for the GWAS extraction window.
parameter: customized_association_windows = path()
# cis window (bp) on each side of the gene for fine-mapping and the default
# GWAS extraction window.
parameter: cis_window = 1000000
# Save the per-gene colocboost input bundle (currently unused; reserved).
parameter: save_data = False
# QtlDataset study label and covariate handling (legacy CLI retained).
parameter: study = name
parameter: transpose_covariates = False
parameter: genotype_covariates = path('.')
# Remove a variant if it has more than imiss fraction missing genotypes.
parameter: imiss = 1.0
# MAF and variance-of-X cutoffs.
parameter: maf = 0.0025
parameter: xvar_cutoff = 0.0
# MAC cutoff, on top of the MAF cutoff.
parameter: mac = 5
# Keep indels when True.
parameter: indel = True
# Only keep a subset of samples / variants (whitespace/newline-separated IDs).
parameter: keep_samples = path()
parameter: keep_variants = path()
# Per-context single-effect skip cutoff passed to colocboost.R
# (--pip-cutoff-to-skip): a scalar applied to every context, or one or more
# context=value items (negative -> data-driven 3/ncol(X); empty -> 0).
parameter: skip_analysis_pip_cutoff = []
# Summary-stats single-effect skip cutoff passed to summaryStatsQc via
# gwas_sumstats_construct.R.
parameter: skip_sumstats_analysis_pip_cutoff = -1.0
# ColocBoost variants to run.
parameter: xqtl_coloc = True
parameter: joint_gwas = False
parameter: separate_gwas = True
# Summary-stats LD-mismatch QC method (none | slalom | dentist).
parameter: qc_method = "none"
# 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 of commands to run per job.
parameter: job_size = 200
parameter: walltime = "1h"
parameter: mem = "20G"
parameter: numThreads = 1

[colocboost_1]
# Build one pecotmr::QtlDataset from the phenotype manifest + shared genotype
# and serialize to RDS (mirrors mnm_regression's qtl_dataset_construct). No
# fan-out: the per-gene [colocboost_3] step loads this single RDS and selects
# the focal gene/contexts at analysis time. Per-context covariates come from
# the manifest cov_path column; a uniform --covFile (or --genotype-covariates)
# applies across all contexts. ColocBoost runs on unscaled residuals
# (--no-scale-residuals), matching the legacy pipeline.
cov_arg = covFile[0] if (len(covFile) > 0 and covFile[0].is_file()) else genotype_covariates
input: None
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, entrypoint = entrypoint
    Rscript ${modular_script_dir}/pecotmr_integration/qtl_dataset_construct.R \
        --study ${study} \
        --genotype-prefix ${genoFile:n} \
        --phenotype-manifest ${phenoFile[0]} \
        --genotype-covariates ${cov_arg if cov_arg.is_file() else '""'} \
        ${'--transpose-covariates' if transpose_covariates else ''} \
        --maf-cutoff ${maf} \
        --xvar-cutoff ${xvar_cutoff} \
        --mac-cutoff ${mac} \
        --imiss-cutoff ${imiss} \
        ${('--keep-samples ' + str(keep_samples)) if keep_samples.is_file() else ''} \
        ${('--keep-variants ' + str(keep_variants)) if keep_variants.is_file() else ''} \
        ${'' if indel else '--drop-indel'} \
        --no-scale-residuals \
        --output ${_output}

[colocboost_2]
# Resolve per-gene analysis units (gene -> ld_block + grouped GWAS sources)
# into one manifest TSV via colocboost_manifest.R, so the next step can fan out
# over its rows with inline csv.DictReader (no notebook-local Python parsing).
input: None
output: f"{cwd:a}/colocboost/{name}.colocboost_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/colocboost_manifest.R \
        --pheno-manifest ${phenoFile[0]} \
        --cis-window ${cis_window} \
        ${('--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 ''} \
        ${('--gwas-meta ' + str(gwas_meta_data)) if gwas_meta_data.is_file() else ''} \
        ${('--ld-meta ' + str(ld_meta_data)) if ld_meta_data.is_file() else ''} \
        --output ${_output}

ColocBoost analysis

[colocboost_3]
# Per-gene ColocBoost over the pre-built QtlDataset. Fans out over the manifest
# rows; per gene (1) builds a region GwasSumStats via gwas_sumstats_construct.R
# when the gene has GWAS studies, then (2) runs colocboost.R (xQTL /
# joint-GWAS / separate-GWAS variants) into ONE combined RDS. Replaces the
# legacy load_multitask_regional_data + colocboost_analysis_pipeline path.
import csv
manifest = f"{cwd:a}/colocboost/{name}.colocboost_manifest.tsv"
jobs = list(csv.DictReader(open(manifest), delimiter='\t'))
stop_if(len(jobs) == 0, "colocboost: empty manifest; check --region-name / phenotype manifest.")
input: f"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds", for_each = "jobs"
output: f"{cwd:a}/colocboost/{name}.{_jobs['region_id']}.colocboost.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, entrypoint = entrypoint
    set -e
    gene=${_jobs['gene_id']}
    studies="${_jobs['studies']}"
    gss_arg=""
    if [ -n "$studies" ]; then
        gss="${_output:nn}.gwas_sumstats.rds"
        Rscript ${modular_script_dir}/pecotmr_integration/gwas_sumstats_construct.R \
            --study "$studies" \
            --gwas-tsv "${_jobs['gwas_tsvs']}" \
            --ld-block "${_jobs['ld_block']}" \
            --ld-meta "${_jobs['ld_meta']}" \
            --n-case "${_jobs['n_cases']}" \
            --n-control "${_jobs['n_controls']}" \
            --pip-cutoff-to-skip ${skip_sumstats_analysis_pip_cutoff} \
            --qc-method ${qc_method} \
            ${('--column-mapping "' + _jobs['column_mappings'] + '"') if _jobs['column_mappings'] else ''} \
            --output "$gss"
        gss_arg="--gwas-sumstats $gss"
    fi
    Rscript ${modular_script_dir}/pecotmr_integration/colocboost.R \
        --qtl-dataset ${_input} \
        --gene-id $gene \
        --cis-window ${cis_window} \
        $gss_arg \
        ${'' if xqtl_coloc else '--no-xqtl-coloc'} \
        ${'--joint-gwas' if joint_gwas else ''} \
        ${'--separate-gwas' if separate_gwas else ''} \
        --pip-cutoff-to-skip "${",".join([str(x) for x in skip_analysis_pip_cutoff])}" \
        --output ${_output}