TWAS, cTWAS and MR - FunGen-xQTL Consortium

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_data tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv (one row per GWAS study. Columns study_id, chrom, file_path, sample_size; file_path points 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_data tests/fixtures/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	tests/fixtures/twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.univariate_twas_weights.rds	bulk_rnaseq

--ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv (LD blocks and where the matrices live. Columns #chr, start, end, path. Produced by ld_reference_generation.)

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

--regions tests/fixtures/twas/protocol_example.twas.LD_blocks.chr22.bed (the LD blocks to analyze. Columns chr, start, stop.)

chr	start	stop
chr22	10000000	19000000

--xqtl_type_table tests/fixtures/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 existing ctwas_inputs.rds instead of rebuilding it)

File formats:

GWAS summary statistics are tab-delimited and tabix-indexed by chrom/pos; the first four columns are chrom, pos, A1, A2. For MR, effect_allele_frequency and a sample-size column are also required.
The GWAS meta-file may carry an optional column_mapping_file column pointing at a YAML that renames columns: required chrom, pos, A1, A2, z (or betahat/sebetahat); optional n, var_y. chrom 0 means genome-wide.
The LD reference meta-file gives the LD-matrix path as a comma-separated matrix,bim pair, 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. Columns gwas_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_method marks 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 tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \
    --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \
    --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \
    --ld_reference_sample_size 17000 \
    --regions tests/fixtures/twas/protocol_example.twas.LD_blocks.chr22.bed \
    --xqtl_type_table tests/fixtures/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 tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \
    --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \
    --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \
    --regions tests/fixtures/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 tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \
    --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \
    --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \
    --regions tests/fixtures/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 tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \
    --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \
    --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \
    --regions tests/fixtures/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 tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv \
    --xqtl_meta_data tests/fixtures/twas/protocol_example.twas.xqtl_meta.tsv \
    --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.tsv \
    --ld_reference_sample_size 17000 \
    --regions tests/fixtures/twas/protocol_example.twas.LD_blocks.chr22.bed \
    --xqtl_type_table tests/fixtures/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.")

References¶

Zhao, S., Crouse, W., Qian, S., Luo, K., Stephens, M., & He, X. (2024). Adjusting for genetic confounders in transcriptome-wide association studies improves discovery of risk genes of complex traits. Nature Genetics, 56(2), 336–347. 10.1038/s41588-023-01648-9