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 tooutput.)
Individual-level xQTL inputs.
--genoFiletests/fixtures/qtl_mini/protocol_example.genotype.chr22.bed(genotype data in PLINKbedformat; either one file for the whole genome or a list of per-chromosome files)--phenoFiletests/fixtures/qtl_mini/protocol_example.pheno_manifest_context.tsv(a manifest of phenotype files per region, one row per region and context, in UCSCbed.gzformat with abed.gz.tbiindex)--covFiletests/fixtures/qtl_mini/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-windowstests/fixtures/qtl_mini/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 Fine-mapping with GWAS Summary Statistics 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-gwasis 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 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 --transpose-covariates \
--customized-association-windows tests/fixtures/qtl_mini/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 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 --transpose-covariates \
--customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
--gwas-meta-data tests/fixtures/qtl_mini/gwas_meta.txt \
--ld-meta-data tests/fixtures/ld_reference/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 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 --transpose-covariates \
--customized-association-windows tests/fixtures/qtl_mini/association_windows.bed \
--gwas-meta-data tests/fixtures/qtl_mini/gwas_meta.txt \
--ld-meta-data tests/fixtures/ld_reference/ld_meta_file.tsv \
--region-name ENSG00000130538 \
--separate-gwas --xqtl-coloc -j1Command Interface¶
sos run pipeline/colocboost.ipynb -husage: 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
# RNG seed forwarded to colocboost.R (stochastic gradient boosting) for reproducible fits.
parameter: seed = 999
# 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 ''} \
--seed ${seed} \
--pip-cutoff-to-skip "${",".join([str(x) for x in skip_analysis_pip_cutoff])}" \
--output ${_output}
- Cao, X., Sun, H., Feng, R., Mazumder, R., Najar, C. F. B. A., Li, Y. I., de Jager, P. L., Bennett, D., Dey, K. K., & Wang, G. (2025). Integrative multi-omics QTL colocalization maps regulatory architecture in aging human brain. 10.1101/2025.04.17.25326042