TWAS, cTWAS and MR#
Runs transcriptome-wide association analysis, quantile TWAS, and cTWAS fine-mapping to identify which gene and SNP signals are likely to be causal.
Overview#
A TWAS scan tests each gene’s genetically predicted expression against a trait, but a
significant gene is not necessarily a causal one: nearby variants with direct effects on
the trait, and the predicted expression of neighbouring genes, are correlated with the
gene’s own eQTLs and act as confounders. cTWAS addresses this by fine-mapping genes and
variants jointly within a region, so a gene is credited only for signal that its expression
explains beyond the surrounding variants and genes, and reports a posterior inclusion
probability rather than a p-value (Zhao et al., 2024).
finemapCtwasRegions) and offers four workflows:
This module provides software implementations for transcriptome-wide association analysis (TWAS), Quantile TWAS, and variant selection that yields sparse signals for cTWAS (causal TWAS) analysis following the multi-group cTWAS method of Qian et al. (2024+). It additionally performs Mendelian Randomization using fine-mapping instrumental variables (IV) as described in Zhang et al. (2020) for “causal” effect estimation and model validation. The unit of analysis is a single gene-trait pair.
This notebook runs TWAS followed by cTWAS fine-mapping on a toy chr22 dataset: starting from pre-computed SuSiE-TWAS prediction weights for one gene, it tests each molecular context for association with a GWAS trait, keeps the imputable models selected by cross-validation, and then jointly fine-maps genes and SNPs within an LD block to identify which signals are likely to be directly causal rather than driven by correlation. Quantile TWAS extends traditional TWAS by testing genetic effects at different quantiles of the trait distribution.
This procedure is a continuation of the SuSiE-TWAS workflow: it assumes that xQTL fine-mapping has been performed and molecular-trait prediction weights pre-computed (to be used for TWAS). Cross-validation of TWAS weights is optional but highly recommended.
Prerequisites: TWAS weights (protocol_example.twas_weights.rds) from mnm_regression, plus GWAS summary statistics and an LD matrix for the region of interest.
Timing (toy chr22 dataset): twas ~50 sec; ctwas ~42 sec.
MR can be run for candidate genes afterwards, limited to genes with cTWAS significance and a strong instrumental variable, using fine-mapped xQTL with GWAS data; multiple IVs aggregate by fixed-effect meta-analysis, and results with severe exclusion-restriction violations are excluded.
When to run it. After you have TWAS weights (from mnm_regression) and GWAS summary
statistics. cTWAS consumes weights, it does not produce them.
Input#
--gwas_meta_datainput/twas/protocol_example.twas.gwas_meta.tsv(one row per GWAS study. Columnsstudy_id,chrom,file_path,sample_size;file_pathpoints at the summary statistics for that chromosome.)
study_id chrom file_path sample_size
protocol_example_twas_chr22 22 protocol_example.twas.gwas_sumstats.chr22.tsv.gz 200000
--xqtl_meta_datainput/twas/protocol_example.twas.xqtl_meta.tsv(one row per gene or region with its weight file. Columns#chr,start,end,region_id,TSS,original_data,contexts.)
#chr start end region_id TSS original_data contexts
chr22 10000000 19000000 ENSG00000130538 15528191 input/twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.univariate_twas_weights.rds bulk_rnaseq
--ld_meta_datainput/ld_reference/protocol_example.ld_meta_file.tsv(LD blocks and where the matrices live. Columns#chr,start,end,path. Produced byld_reference_generation.)
#chr start end path
chr1 16103 2888443 chr1/protocol_example.LD.chr1
chr1 2888443 4320284 chr1/protocol_example.LD.chr1
--regionsinput/twas/protocol_example.twas.LD_blocks.chr22.bed(the LD blocks to analyze. Columnschr,start,stop.)
chr start stop
chr22 10000000 19000000
--xqtl_type_tableinput/twas/protocol_example.twas.data_type_table.txt(maps each xQTL context to its modality, used when grouping contexts for cTWAS)
context type
bulk_rnaseq eQTL
--cwd output/twas(working directory all outputs are written under)--name protocol_example(tag used in every output filename)--region-name chr22_10000000_19000000(optional; restrict to named regions)
Model and filtering options:
--ld_reference_sample_size(sample size of the LD reference panel, used to scale the LD matrices)--prior_var_structure(how the cTWAS prior variance is shared across groups)--rsq_cutoff(minimum cross-validation r-squared for a weight model to be used)--rsq_pval_cutoff(maximum cross-validation p-value for a weight model to be used)--twas_weight_cutoff(drop weights below this magnitude before assembling cTWAS inputs)--thin(fraction of SNPs retained when thinning for parameter estimation)--run_param_est(estimate cTWAS group priors)--run_finemapping(run cTWAS fine-mapping after the priors are estimated)--skip_assembly(reuse an existingctwas_inputs.rdsinstead of rebuilding it)
File formats:
GWAS summary statistics are tab-delimited and tabix-indexed by
chrom/pos; the first four columns arechrom,pos,A1,A2. For MR,effect_allele_frequencyand a sample-size column are also required.The GWAS meta-file may carry an optional
column_mapping_filecolumn pointing at a YAML that renames columns: requiredchrom,pos,A1,A2,z(orbetahat/sebetahat); optionaln,var_y.chrom 0means genome-wide.The LD reference meta-file gives the LD-matrix path as a comma-separated
matrix,bimpair, plus the genome build. See ld_reference_generation.xQTL weights are an RDS database keyed region to context to weight matrix. Weights from LASSO, Elastic Net and mr.ash are taken as-is for QTL variants overlapping GWAS variants; SuSiE weights are adjusted to match GWAS variants exactly.
Use
--comment_string "#"if the summary statistics carry comment lines; no comment character is assumed by default.--cs_min_cor(minimum correlation for a cTWAS credible set; default 0)--min_pip_cutoff(minimum PIP for a variant to be reported; default 0)--max_num_variants(cap on variants considered per region; default unlimited)--multi_group(use the multi-group cTWAS model instead of a single shared prior)
Output#
output/twas/twas/<name>.<region_id>.<gene>.twas.rds(twas) (per-gene TWAS result, one file per gene in the region)output/twas/twas/<name>.<region_id>.twas.tsv.gz(twas, side file) (the TWAS table for the whole region. Columnsgwas_study,chrom,start,end,block,gene,TSS,context,is_imputable,method,is_selected_method,rsq_adj_cv,pval_cv,twas_z,twas_pval.is_selected_methodmarks the best-performing model per gene-context pair – the one with the highest cross-validation r-square.)
chr molecular_id TSS start end context gwas_study method is_imputable
22 ENSG00000130538 15528191 10000000 19000000 bulk_rnaseq protocol_example_twas_chr22 bayes_c TRUE
22 ENSG00000130538 15528191 10000000 19000000 bulk_rnaseq protocol_example_twas_chr22 bayes_r TRUE
output/twas/twas/<name>.<region_id>.mr_result.tsv.gz(twas, side file) (Mendelian Randomization results for the same region)
gene_name num_CS num_IV cpip meta_eff se_meta_eff meta_pval Q
ENSG00000130538
<cwd>/ctwas/<name>.ctwas_inputs.rds(ctwas_1) (assembled cTWAS inputs: weights, LD and GWAS z-scores for the regions to be fine-mapped)<cwd>/ctwas/<name>.ctwas_est.rds(ctwas_2) (estimated cTWAS group prior parameters)<cwd>/ctwas/<name>.ctwas_finemap.rds(ctwas_3) (cTWAS fine-mapping result: posterior inclusion probabilities for genes and SNPs in each LD block)
Only the .twas.rds files are declared in the twas step output: statement. The .twas.tsv.gz, .mr_result.tsv.gz, .twas_data.rds, .gwas_sumstats.rds and .twas_weights.rds files, and the cor_cache/ directory, are written alongside but not tracked by SoS.
Minimal Working Example#
The three methods below are independent: run whichever you need. All use the chr22 example data in this repository.
Standard TWAS scan#
twas extracts the GWAS z-scores and LD matrix for each region, processes the SuSiE-TWAS weights, and runs the association test for each molecular context, keeping the best cross-validated model per gene. Genes that clear neither --rsq_cutoff nor --rsq_pval_cutoff in any method are dropped as non-imputable.
Timing: TBD (on toy dataset)
sos run pipeline/twas_ctwas.ipynb twas \
--cwd output --name protocol_example \
--gwas_meta_data input/twas/protocol_example.twas.gwas_meta.tsv \
--xqtl_meta_data input/twas/protocol_example.twas.xqtl_meta.tsv \
--ld_meta_data input/ld_reference/protocol_example.ld_meta_file.tsv \
--ld_reference_sample_size 17000 \
--regions input/twas/protocol_example.twas.LD_blocks.chr22.bed \
--xqtl_type_table input/twas/protocol_example.twas.data_type_table.txt \
--rsq_pval_cutoff 0.05 --rsq_cutoff 0.01 \
--region-name chr22_10000000_19000000
cTWAS fine-mapping#
ctwas combines the selected TWAS models, gene and SNP z-scores and the LD reference into a per-region input, estimates the global group prior across regions, then fine-maps each region for posterior inclusion probabilities over genes and SNPs. It is one workflow run three times, switched by --run_param_est and --run_finemapping; --skip_assembly reuses the assembled inputs rather than rebuilding them.
Requires the cTWAS R package: remotes::install_github("xinhe-lab/ctwas", ref = "multigroup").
-n cannot preview this chain: ctwas_1 reads a manifest that get_analysis_regions only writes at run time, so a dry run stops at the first step.
Step 1.: Assemble region data#
Timing: TBD (on toy dataset)
sos run pipeline/twas_ctwas.ipynb ctwas \
--cwd output --name protocol_example \
--thin 1 --prior_var_structure shared_all \
--gwas_meta_data input/twas/protocol_example.twas.gwas_meta.tsv \
--xqtl_meta_data input/twas/protocol_example.twas.xqtl_meta.tsv \
--ld_meta_data input/ld_reference/protocol_example.ld_meta_file.tsv \
--regions input/twas/protocol_example.twas.LD_blocks.chr22.bed \
--twas_weight_cutoff 0 \
--region-name chr22_10000000_19000000
Step 2.: Estimate global parameters#
Timing: TBD (on toy dataset)
sos run pipeline/twas_ctwas.ipynb ctwas \
--run_param_est --skip_assembly --thin 1 \
--prior_var_structure shared_all \
--cwd output --name protocol_example \
--gwas_meta_data input/twas/protocol_example.twas.gwas_meta.tsv \
--xqtl_meta_data input/twas/protocol_example.twas.xqtl_meta.tsv \
--ld_meta_data input/ld_reference/protocol_example.ld_meta_file.tsv \
--regions input/twas/protocol_example.twas.LD_blocks.chr22.bed
Step 3.: Fine-map causal genes and SNPs#
Timing: TBD (on toy dataset)
sos run pipeline/twas_ctwas.ipynb ctwas \
--run_finemapping --skip_assembly \
--prior_var_structure shared_all \
--cwd output --name protocol_example \
--gwas_meta_data input/twas/protocol_example.twas.gwas_meta.tsv \
--xqtl_meta_data input/twas/protocol_example.twas.xqtl_meta.tsv \
--ld_meta_data input/ld_reference/protocol_example.ld_meta_file.tsv \
--regions input/twas/protocol_example.twas.LD_blocks.chr22.bed \
--region-name chr22_10000000_19000000
Quantile TWAS#
Use pre-computed TWAS weights for quantile-specific testing: for each quantile level, cluster and integrate by fixed and dynamic region groups, extract the relevant GWAS z-scores and LD matrix for the region, and perform quantile region-specific association tests, identifying variants whose effects vary across different quantile regions of the phenotype distribution. Use --region-name (formatted as chr_start_stop) to focus on specific blocks, or --region to analyze a selected list of regions.
Timing: TBD (on toy dataset)
sos run pipeline/twas_ctwas.ipynb quantile_twas \
--cwd output --name protocol_example \
--gwas_meta_data input/twas/protocol_example.twas.gwas_meta.tsv \
--xqtl_meta_data input/twas/protocol_example.twas.xqtl_meta.tsv \
--ld_meta_data input/ld_reference/protocol_example.ld_meta_file.tsv \
--ld_reference_sample_size 17000 \
--regions input/twas/protocol_example.twas.LD_blocks.chr22.bed \
--xqtl_type_table input/twas/protocol_example.twas.data_type_table.txt \
--region-name chr22_10000000_19000000
Command Interface#
sos run pipeline/twas_ctwas.ipynb -h
usage: sos run pipeline/twas_ctwas.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
twas
ctwas
quantile_twas
Global Workflow Options:
--cwd output (as path)
--modular-script-dir code/script (as path)
--gwas-meta-data . (as path)
--xqtl-meta-data . (as path)
--ld-meta-data . (as path)
--ld-reference-sample-size 17000 (as int)
--xqtl-type-table ''
--gwas-name (as list)
--gwas-data (as list)
--column-mapping (as list)
--regions . (as path)
--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:bn}.{gwas_meta_data:bn}"
--container ''
--job-size 100 (as int)
--walltime 5m
--mem 8G
--numThreads 1 (as int)
--name-suffix ''
name suffix add to end of the name variable for the
output files
--gwas-study (as list)
optional parameter in ctwas to only perform gwas
specific analysis
Sections
get_analysis_regions:
twas:
Workflow Options:
--coverage 'cs_coverage_0.95'
Legacy CLI retained. Wired to the pecotmr selection
knobs: rsq_cutoff, rsq_pval_cutoff, rsq_option,
rsq_pval_option (CV weight selection). MR is deferred,
so mr_pval_cutoff and the remaining legacy params are
kept declared for CLI stability but are not consumed by
the S4 path.
--rsq-cutoff 0.01 (as float)
Thresholds for rsq and CV p-value for imputability/best-
model selection
--rsq-pval-cutoff 0.05 (as float)
--mr-pval-cutoff 0.05 (as float)
--[no-]save-ctwas-data (default to True)
--[no-]save-mr-result (default to True)
--rsq-option rsq
--rsq-pval-option adj_rsq_pval pval (as list)
--batch-load-memory 500 (as int)
load by batches if memory resource is limited, default
to load all at once
--event-filter-rules . (as path)
--comment-string NULL
--[no-]rename-column (default to False)
ctwas_1:
Workflow Options:
--chrom ''
chromosome to assemble (integer or 'chrN'); default: the
lone gene-bearing chrom
--twas-weight-cutoff 0.0 (as float)
weight pre-filters passed to assembleCtwasInputs
--max-num-variants Inf
--cs-min-cor 0.0 (as float)
--min-pip-cutoff 0.0 (as float)
--thin 1.0 (as float)
declared for CLI stability; consumed downstream / not by
assembleCtwasInputs
--maxSNP 20000 (as int)
--[no-]multi-group (default to True)
--[no-]skip-assembly (default to False)
--twas-weights (as list)
TwasWeights source (option c): empty -> the upstream
twas step's per-gene weights
({cwd}/twas/{name}.{region}.{gene}.twas_weights.rds);
override with a list of prebuilt TwasWeights RDS (e.g.
the blessed ctwas weights extracted by the blessed S4
ctwas weights) when the toy twas weights carry no
signal.
ctwas_2:
Workflow Options:
--[no-]run-param-est (default to False)
cTWAS step 2 (estimate priors): estCtwasParam via
ctwas_est.R. Reads the assembled inputs from [ctwas_1]
on disk so it can be run independently. --fallback-to-
prefit recovers the prefit estimate when the accurate EM
cannot be estimated (e.g. the single-gene MWE),
mirroring the legacy ctwas_2 workaround.
--[no-]skip-assembly (default to False)
--thin 1.0 (as float)
--prior-var-structure 'shared_all'
--niter 50 (as int)
--[no-]multi-group (default to True)
declared for CLI stability; not consumed by the S4 path
ctwas_3:
Workflow Options:
--[no-]run-finemapping (default to False)
cTWAS step 3 (screen + finemap): screenCtwasRegions +
finemapCtwasRegions via ctwas_finemap.R. Reads the
estimated params from [ctwas_2] on disk so it can be run
independently. Emits the per-gene cTWAS fine-mapping
result.
--L 5 (as int)
--min-nonSNP-PIP 0.5 (as float)
--[no-]merge-regions (default to False)
Boundary-gene region merging (legacy ctwas_3
merge_regions, default-off): when True, after fine-
mapping, each high-PIP boundary gene's adjacent LD
blocks are merged and re-fine-mapped via
mergeCtwasBoundaryRegions (require-in-CS, like the
legacy susie_pip>0.5 & !is.na(cs) selection); maxSNP
caps SNPs per merged region.
--maxSNP 20000 (as int)
--[no-]keep-snps (default to True)
retain the SNP background as a dedicated
study=context="SNP" CtwasResult row; default True
matches the legacy pipeline, which kept SNP rows in the
finemap output
--thin 1.0 (as float)
declared for CLI stability; not consumed by the S4 path
--max-iter 0 (as int)
--prior-var-structure 'shared_all'
--subset-context (as list)
--[no-]multi-group (default to True)
--p-diff-thresh 5e-08 (as float)
--alias NULL
quantile_twas:
Workflow Options:
--[no-]save-ctwas-data (default to True)
--[no-]save-mr-result (default to False)
Workflow implementation#
The following cells define the SoS workflow steps ([global], [get_analysis_regions], [twas], [ctwas_1], [ctwas_2], [ctwas_3], [quantile_twas]) invoked by the commands in the steps above.
[global]
parameter: cwd = path("output/")
parameter: modular_script_dir = path('code/script')
parameter: gwas_meta_data = path()
parameter: xqtl_meta_data = path()
parameter: ld_meta_data = path()
parameter: ld_reference_sample_size = 17000
parameter: xqtl_type_table = ''
parameter: gwas_name = []
parameter: gwas_data = []
parameter: column_mapping = []
parameter: regions = 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 = []
parameter: name = f"{xqtl_meta_data:bn}.{gwas_meta_data:bn}"
parameter: container = ''
parameter: job_size = 100
parameter: walltime = "5m"
parameter: mem = "8G"
parameter: numThreads = 1
# name suffix add to end of the name variable for the output files
parameter: name_suffix = ""
# optional parameter in ctwas to only perform gwas specific analysis
parameter: gwas_study=[]
[get_analysis_regions: shared = {'twas_manifest_rows': "list(__import__('csv').DictReader(open(twas_manifest_file), delimiter=chr(9)))"}]
# Resolve the per-region TWAS analysis units into a manifest TSV via
# twas_manifest.R -- the gene<->region TSS-overlap binning, GWAS study/chrom
# resolution, and meta-relative path adaptation that the legacy
# extract_regional_data did in-notebook with pandas. The parsed rows are shared
# as twas_manifest_rows, so [twas] / [ctwas] auto-trigger this step via
# depends: sos_variable -- the notebook is invoked exactly as before
# (e.g. sos run twas_ctwas.ipynb twas). Gene -> home-LD-block placement is NOT
# done here: pecotmr's assembleCtwasInputs places each gene into its block from
# the weight `region` provenance.
input: None
twas_manifest_file = f"{cwd:a}/get_analysis_regions/{name}.twas_manifest.tsv"
bash: expand = "${ }", container = container
mkdir -p ${cwd:a}/get_analysis_regions
Rscript ${modular_script_dir}/pecotmr_integration/twas_manifest.R \
--gwas-meta ${gwas_meta_data} \
--xqtl-meta ${xqtl_meta_data} \
${("--regions " + str(regions)) if regions.is_file() else ""} \
${("--region-name " + ",".join(region_name)) if len(region_name) > 0 else ""} \
${("--gwas-name " + " ".join([str(x) for x in gwas_name])) if len(gwas_name) > 0 else ""} \
${("--gwas-data " + " ".join([str(x) for x in gwas_data])) if len(gwas_data) > 0 else ""} \
${("--column-mapping " + " ".join([str(x) for x in column_mapping])) if len(column_mapping) > 0 else ""} \
--output ${twas_manifest_file}
[twas]
# Per-region TWAS-Z + Mendelian Randomization, bridged to the pecotmr S4
# wrappers (no inline analysis R; no Python helpers). Restores the legacy
# per-region fan-out: one task per analysis region, looping the region's genes
# inside. Per region it (1) builds one multi-study GwasSumStats over the region
# via gwas_sumstats_construct.R (its LD sketch spans whatever LD-reference
# blocks the region overlaps), then per gene (2) converts the legacy
# per-gene S4 TwasWeights (from twas_weights.ipynb) and
# runs twas.R (causalInferencePipeline) on the per-gene S4 TwasWeights
# to emit a per-gene .twas.rds. MR is deferred (no fine-mapping result is
# passed), so only TWAS Z + p-value are produced.
# Fan-out consumes the get_analysis_regions manifest via its shared
# twas_manifest_rows, which auto-triggers that step -- so this runs exactly as
# before, e.g. sos run twas_ctwas.ipynb twas.
depends: sos_variable('twas_manifest_rows')
# Legacy CLI retained. Wired to the pecotmr selection knobs: rsq_cutoff,
# rsq_pval_cutoff, rsq_option, rsq_pval_option (CV weight selection). MR is
# deferred, so mr_pval_cutoff and the remaining legacy params are kept declared
# for CLI stability but are not consumed by the S4 path.
parameter: coverage = "cs_coverage_0.95"
# Thresholds for rsq and CV p-value for imputability/best-model selection
parameter: rsq_cutoff = 0.01
parameter: rsq_pval_cutoff = 0.05
parameter: mr_pval_cutoff = 0.05
parameter: save_ctwas_data = True
parameter: save_mr_result = True
parameter: rsq_option = "rsq"
parameter: rsq_pval_option = ["adj_rsq_pval", "pval"]
# load by batches if memory resource is limited, default to load all at once
parameter: batch_load_memory = 500
parameter: event_filter_rules = path()
parameter: comment_string = "NULL"
parameter: rename_column = False
jobs = twas_manifest_rows
stop_if(len(jobs) == 0,
"No regions with overlapping xQTL weights found; skipping TWAS step.")
_weight_groups = [j['weight_files'].split(',') for j in jobs]
input: [f for grp in _weight_groups for f in grp], group_by = lambda x: _weight_groups, group_with = "jobs"
output: [f'{cwd:a}/{step_name}/{name}.{_jobs["region_id"]}.{gene}.twas.rds' for gene in sorted(set(_jobs["genes"].split(',')))]
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
bash: expand = '${ }', stdout = f"{_output[0]:n}.stdout", stderr = f"{_output[0]:n}.stderr", container = container
set -e
# SoS expands ${...}; bare $shellvars (no braces) are left for bash.
outdir=${cwd:a}/${step_name}
mkdir -p "$outdir"
region_id=${_jobs["region_id"]}
ld_block=${_jobs["chrom"]}:${_jobs["start"]}-${_jobs["stop"]}
# GWAS studies covering this region's chromosome (paths already resolved by
# get_analysis_regions). studies and gwas_tsvs iterate the same dict order.
studies="${_jobs["gwas_studies"]}"
gwas_tsvs="${_jobs["gwas_files"]}"
# (1) One multi-study GwasSumStats over the whole region.
gss="$outdir/${name}.$region_id.gwas_sumstats.rds"
Rscript ${modular_script_dir}/pecotmr_integration/gwas_sumstats_construct.R \
--study "$studies" \
--gwas-tsv "$gwas_tsvs" \
--ld-block "$ld_block" \
--ld-meta "${ld_meta_data}" \
--output "$gss"
# (2)+(3) Per unique gene in this region (first weight file per gene):
# run TWAS-Z + MR on the per-gene S4 weights. The expansion
# below emits a flat "gene file gene file ..." token list (first file per
# gene); the loop consumes it two tokens at a time.
n_files=${len(_input)}
n_genes=${len(set(_jobs["genes"].split(",")))}
if [ "$n_files" -ne "$n_genes" ]; then
echo "NOTE: $n_files weight file(s) across $n_genes gene(s) in region $region_id; using the first file per gene." >&2
fi
set -- ${" ".join([tok for i, (g, f) in enumerate(zip(_jobs["genes"].split(","), _input)) if g not in _jobs["genes"].split(",")[:i] for tok in (g, str(f))])}
while [ "$#" -ge 2 ]; do
gene="$1"
weight="$2"
shift 2
Rscript ${modular_script_dir}/pecotmr_integration/twas.R \
--twas-weights "$weight" \
--gwas-sumstats "$gss" \
--rsq-cutoff ${rsq_cutoff} \
--rsq-pval-cutoff ${rsq_pval_cutoff} \
--rsq-option ${rsq_option} \
--rsq-pval-option "${",".join(rsq_pval_option)}" \
--output "$outdir/${name}.$region_id.$gene.twas.rds"
done
[ctwas_1]
# cTWAS step 1 (assemble): build the per-LD-block inputs for one gene-bearing
# chromosome and assemble them via the pecotmr S4 wrappers (no inline analysis
# R; no Python helpers). Per chromosome: (1) ctwas_manifest.R enumerates the
# chromosome's LD blocks from ld_meta_data (just the block grid; pecotmr places
# each gene into its home block from the weight `region` provenance);
# (2) gwas_sumstats_construct.R
# builds one GwasSumStats per block (SNP background); (3) ctwas_assemble.R runs
# assembleCtwasInputs -> {name}.ctwas_inputs.rds. cTWAS runs over the whole
# chromosome's LD blocks, NOT the coarse twas analysis region.
# Fan-in consumes the get_analysis_regions manifest via its shared
# twas_manifest_rows, which auto-triggers that step -- so this runs exactly as
# before, e.g. sos run twas_ctwas.ipynb ctwas.
depends: sos_variable('twas_manifest_rows')
# chromosome to assemble (integer or 'chrN'); default: the lone gene-bearing chrom
parameter: chrom = ""
# weight pre-filters passed to assembleCtwasInputs
parameter: twas_weight_cutoff = 0.0
parameter: max_num_variants = "Inf"
parameter: cs_min_cor = 0.0
parameter: min_pip_cutoff = 0.0
# declared for CLI stability; consumed downstream / not by assembleCtwasInputs
parameter: thin = 1.0
parameter: maxSNP = 20000
parameter: multi_group = True
parameter: numThreads = 4
parameter: skip_assembly = False
# TwasWeights source (option c): empty -> the upstream twas step's per-gene
# weights ({cwd}/twas/{name}.{region}.{gene}.twas_weights.rds); override with a
# list of prebuilt TwasWeights RDS (e.g. the blessed ctwas weights extracted by
# the blessed S4 ctwas weights) when the toy twas weights carry no signal.
parameter: twas_weights = []
skip_if(skip_assembly == True, "Skip [ctwas_1] assemble.")
jobs = twas_manifest_rows
gene_chroms = sorted(set(j["chrom"] for j in jobs if j["genes"]))
ctwas_chrom = (f"chr{int(chrom)}" if str(chrom).isdigit() else str(chrom)) if str(chrom) else (gene_chroms[0] if len(gene_chroms) == 1 else "")
stop_if(not ctwas_chrom, f"Specify --chrom: expected one gene-bearing chromosome, found {gene_chroms}.")
ctwas_weights = list(twas_weights) if twas_weights else [f"{cwd:a}/twas/{name}.{j['region_id']}.{g}.twas_weights.rds" for j in jobs if j["chrom"] == ctwas_chrom for g in sorted(set(j["genes"].split(',')))]
_chrom_jobs = [j for j in jobs if j["chrom"] == ctwas_chrom]
ctwas_studies = _chrom_jobs[0]["gwas_studies"].split(',') if _chrom_jobs and _chrom_jobs[0]["gwas_studies"] else []
ctwas_gwas_files = _chrom_jobs[0]["gwas_files"] if _chrom_jobs else ""
stop_if(len(ctwas_studies) == 0, f"No GWAS study covers {ctwas_chrom}.")
output: f"{cwd:a}/ctwas/{name}.ctwas_inputs.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
outdir=${cwd:a}/ctwas
mkdir -p "$outdir"
manifest="$outdir/${name}.ctwas_manifest.${ctwas_chrom}.tsv"
# (1) enumerate the chromosome's LD-block grid (block enumeration only).
Rscript ${modular_script_dir}/pecotmr_integration/ctwas_manifest.R \
--ld-meta "${ld_meta_data}" \
--chrom "${ctwas_chrom}" \
--gwas-sumstats-dir "$outdir" \
--output "$manifest"
# (2) one GwasSumStats per LD block (all studies on this chromosome).
studies="${",".join(ctwas_studies)}"
gwas_tsvs="${ctwas_gwas_files}"
tail -n +2 "$manifest" | while IFS=$'\t' read -r region_id region gwas_rds; do
Rscript ${modular_script_dir}/pecotmr_integration/gwas_sumstats_construct.R \
--study "$studies" \
--gwas-tsv "$gwas_tsvs" \
--ld-block "$region" \
--ld-meta "${ld_meta_data}" \
--output "$gwas_rds"
done
# (3) assemble cTWAS inputs (pecotmr places each gene into its home block).
Rscript ${modular_script_dir}/pecotmr_integration/ctwas_assemble.R \
--manifest "$manifest" \
--twas-weights "${",".join([str(w) for w in ctwas_weights])}" \
--twas-weight-cutoff ${twas_weight_cutoff} \
--cs-min-cor ${cs_min_cor} \
--min-pip-cutoff ${min_pip_cutoff} \
--max-num-variants ${max_num_variants} \
--output ${_output}
[ctwas_2]
# cTWAS step 2 (estimate priors): estCtwasParam via ctwas_est.R. Reads the
# assembled inputs from [ctwas_1] on disk so it can be run independently.
# --fallback-to-prefit recovers the prefit estimate when the accurate EM
# cannot be estimated (e.g. the single-gene MWE), mirroring the legacy ctwas_2
# workaround.
parameter: run_param_est = False
parameter: skip_assembly = False
parameter: thin = 1.0
parameter: prior_var_structure = "shared_all"
parameter: niter = 50
# declared for CLI stability; not consumed by the S4 path
parameter: multi_group = True
parameter: numThreads = 8
skip_if(run_param_est == False, "Skip [ctwas_2] parameter estimation.")
input: f"{cwd:a}/ctwas/{name}.ctwas_inputs.rds"
output: f"{cwd:a}/ctwas/{name}.ctwas_est.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
Rscript ${modular_script_dir}/pecotmr_integration/ctwas_est.R \
--inputs ${_input} \
--thin ${thin} \
--niter ${niter} \
--group-prior-var-structure ${prior_var_structure} \
--min-group-size 1 \
--fallback-to-prefit \
--ncore ${numThreads} \
--output ${_output}
[ctwas_3]
# cTWAS step 3 (screen + finemap): screenCtwasRegions + finemapCtwasRegions via
# ctwas_finemap.R. Reads the estimated params from [ctwas_2] on disk so it can
# be run independently. Emits the per-gene cTWAS fine-mapping result.
parameter: run_finemapping = False
parameter: L = 5
parameter: min_nonSNP_PIP = 0.5
# Boundary-gene region merging (legacy ctwas_3 merge_regions, default-off): when
# True, after fine-mapping, each high-PIP boundary gene's adjacent LD blocks are
# merged and re-fine-mapped via mergeCtwasBoundaryRegions (require-in-CS, like
# the legacy susie_pip>0.5 & !is.na(cs) selection); maxSNP caps SNPs per merged
# region.
parameter: merge_regions = False
parameter: maxSNP = 20000
# retain the SNP background as a dedicated study=context="SNP" CtwasResult row;
# default True matches the legacy pipeline, which kept SNP rows in the finemap output
parameter: keep_snps = True
# declared for CLI stability; not consumed by the S4 path
parameter: thin = 1.0
parameter: max_iter = 0
parameter: prior_var_structure = "shared_all"
parameter: region_name = []
parameter: subset_context = []
parameter: multi_group = True
parameter: numThreads = 4
parameter: p_diff_thresh = 5e-8
parameter: alias = "NULL"
skip_if(run_finemapping == False, "Skip [ctwas_3] fine-mapping.")
input: f"{cwd:a}/ctwas/{name}.ctwas_est.rds"
output: f"{cwd:a}/ctwas/{name}.ctwas_finemap.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
Rscript ${modular_script_dir}/pecotmr_integration/ctwas_finemap.R \
--est ${_input} \
--L ${L} \
--min-nonsnp-pip ${min_nonSNP_PIP} \
${('--merge-regions --merge-filter-cs --max-snp ' + str(maxSNP)) if merge_regions else ''} \
${'--keep-snps' if keep_snps else ''} \
--ncore ${numThreads} \
--output ${_output}
[quantile_twas]
depends: sos_variable("filtered_regional_xqtl_files")
parameter: save_ctwas_data = True
parameter: save_mr_result = False
input: filtered_regional_xqtl_files, group_by = lambda x: group_by_region(x, filtered_regional_xqtl_files), group_with = "filtered_region_info"
output_files = [f'{cwd:a}/{step_name}/{name}.{_filtered_region_info[3]}.quantile_twas.tsv.gz']
if save_ctwas_data:
output_files.append(f'{cwd:a}/{step_name}/{name}.{_filtered_region_info[3]}.quantile_twas_data.rds')
if save_mr_result:
output_files.append(f'{cwd:a}/{step_name}/{name}.{_filtered_region_info[3]}.quantile_mr_result.tsv.gz')
output: output_files
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output[0]:bn}'
R: expand = '${ }', stdout = f"{_output[0]:n}.stdout", stderr = f"{_output[0]:n}.stderr", container = container
library(dplyr)
library(data.table)
library(pecotmr)
library(readr)
# get xQTL weight information
xqtl_meta_df <- fread("${xqtl_meta_data}") # Get related gene information from the xqtl_meta data table
xqtl_type_table <- if (isTRUE(file.exists("${xqtl_type_table}"))) fread("${xqtl_type_table}") else NULL
gwas_studies = c(${paths(regional_data["GWAS"].keys()):r,})
gwas_files = c(${paths([v[_filtered_region_info[0]] for k, v in regional_data["GWAS"].items()]):r,})
gene_list <- c(${', '.join([f"'{gene}'" for gene in _filtered_region_info[4]])})
# Initialize export_twas_weights_db
export_twas_weights_db <- list()
export_twas_weights_db[["${_filtered_region_info[3]}"]] <- list()
# Initialize twas_weights_results
twas_weights_results <- list()
# Create initial weight_db_list
weight_db_list <- c(${_input:r,})
names(weight_db_list) <- gene_list
# Split weight_db_list
weight_db_list <- split(weight_db_list, names(weight_db_list))
# Update weight_db_list
weight_db_list_update <- lapply(weight_db_list, function(file_list) {
valid_files <- Filter(function(file) {
if (file.size(file) > 200) {
content <- tryCatch({
data <- readRDS(file)
}, error = function(e) {
warning(paste("Failed to read RDS file:", file))
NULL
})
return(!is.null(content) && any(sapply(content, function(gene_data) {
is.list(gene_data) && any(sapply(gene_data, function(context_data) {
is.list(context_data) && "twas_variant_names" %in% names(context_data)
}))
})))
}
return(FALSE)
}, file_list)
if (length(valid_files) == 0) return(NULL)
return(valid_files)
})
weight_db_list_update <- Filter(Negate(is.null), weight_db_list_update)
if (length(weight_db_list_update) == 0 || all(sapply(weight_db_list_update, length) == 0)) {
message("No valid twas weight files found after filtering. Exiting the script.")
for (of in c(${_output:r,})) { if (grepl("[.]rds$", of)) saveRDS(list(), of) else readr::write_tsv(data.frame(message = character(0)), of) }
quit(save = "no", status = 0)
}
# Check if weight_db_list_update is empty
# Define tau_values
tau_values <- seq(0.01, 0.99, 0.01)
# Main processing loop
for (gene_db in names(weight_db_list_update)) {
weight_dbs <- weight_db_list_update[[gene_db]]
twas_weights_results[[gene_db]] <- load_quantile_twas_weights(
weight_db_files = weight_dbs,
tau_values = tau_values,
between_cluster = 0.8,
num_intervals = 3
)
if (!is.null(twas_weights_results[[gene_db]]) && !is.null(twas_weights_results[[gene_db]]$weights)) {
twas_weights_results[[gene_db]]$data_type <- setNames(
lapply(names(twas_weights_results[[gene_db]]$weights), function(context) {
xqtl_type_table$type[sapply(xqtl_type_table$context, function(x) grepl(x, context))]
}),
names(twas_weights_results[[gene_db]]$weights)
)
} else {
print(paste("Warning: No valid weights found for gene:", gene_db))
}
}
if (length(twas_weights_results) == 0 || all(sapply(twas_weights_results, is.null))) {
stop("twas_weights_results is empty or invalid. Exiting.")
}
saveRDS(twas_weights_results, "${_output[0]:nnn}.grouped_quantile_twas_weight.rds", compress='xz')
#twas_weights_results
#Step 2: twas analysis for imputable genes across contexts
twas_results_db <- twas_pipeline(twas_weights_data = twas_weights_results,
ld_meta_file_path = "${ld_meta_data}",
gwas_meta_file = "${gwas_meta_data}",
region_block = "${_filtered_region_info[3]}",
ld_reference_sample_size = ${ld_reference_sample_size},
quantile_twas = TRUE,
output_twas_data = ${"TRUE" if save_ctwas_data else "FALSE"} )
# Merging with xQTL meta-data
if (is.null(twas_results_db$twas_result) || nrow(twas_results_db$twas_result) == 0) {
message("twas_results_db$twas_result is NULL. Exiting script normally.")
for (of in c(${_output:r,})) { if (grepl("[.]rds$", of)) saveRDS(list(), of) else readr::write_tsv(data.frame(message = character(0)), of) }
quit(save = "no", status = 0)
}
message("Merging twas_result with xqtl_meta_df...")
# Check twas_results_db before merging
common_ids <- intersect(twas_results_db$twas_result$molecular_id, xqtl_meta_df$region_id)
if (length(common_ids) > 0) {
twas_results_db$twas_result <- merge(
twas_results_db$twas_result,
xqtl_meta_df[, c("region_id", "TSS", "start", "end")],
by.x = "molecular_id",
by.y = "region_id"
)
twas_results_db$twas_result <- unique(twas_results_db$twas_result)
} else {
warning("No common molecular_id and region_id. Skipping merge.")
}
fwrite(twas_results_db$twas_result[, c(2, 1, (ncol(twas_results_db$twas_result)-2):ncol(twas_results_db$twas_result), 3:(ncol(twas_results_db$twas_result)-3))], file = ${_output[0]:r}, sep = "\t", compress = "gzip")
# Step 3: reformat for follow up cTWAS analysis
if (${"TRUE" if save_ctwas_data else "FALSE"}) {
saveRDS(twas_results_db$twas_data, "${_output[0]:nnn}.quantile_twas_data.rds", compress='xz')
}
if (${"TRUE" if save_mr_result else "FALSE"}) {
fwrite(twas_results_db$mr_result, file = "${_output[0]:nnn}.quantile_mr_result.tsv.gz", sep = "\t", compress = "gzip")
}
message("quantile twas analysis is completed in this block.")