Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

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, or another molecular trait, 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 xQTLs 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). This module follows the multi-group cTWAS formulation of Qian et al. (2024+), in which each molecular context forms its own group. Quantile TWAS extends the scan by testing genetic effects at different quantiles of the trait distribution. Mendelian randomization uses fine-mapped xQTL variants as instrumental variables (Zhang et al., 2020) to estimate a causal effect for candidate genes: it runs inside the same pecotmr pipeline when an xQTL fine-mapping result is supplied, aggregates multiple instruments by fixed-effect meta-analysis, and excludes results with severe exclusion-restriction violations. The unit of analysis is a single gene-trait pair.

TWAS and cTWAS read prediction weights stored as pecotmr TwasWeights objects. Weights from the legacy FunGen-xQTL release on Synapse use an older list layout and are converted first by the optional twas_weights_conversion workflow; weights already stored as TwasWeights skip this step. The conversion copies the weights and cross-validation metrics unchanged and writes the xQTL meta and data-type tables that the other workflows read. It is optional: skip it when the weights are already TwasWeights. MR in twas also needs the QTL fine-mapping results of the same genes; the optional fine_mapping_conversion workflow converts legacy fine-mapping files the same way.

The module offers five workflows, in the order they run: twas_weights_conversion and fine_mapping_conversion (both optional) convert legacy weights and fine-mapping results, twas tests each gene and molecular context for association using the models that pass the cross-validation cutoffs, ctwas fine-maps genes and variants jointly in three chained steps, and quantile_twas is an alternative scan across trait quantiles.

When to run it. After xQTL fine-mapping and TWAS weight training (mnm_regression), with GWAS summary statistics and an LD reference in hand. Weights from the legacy FunGen-xQTL release need twas_weights_conversion first; weights already in pecotmr TwasWeights format do not. Cross-validated weights are strongly recommended. 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 output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv (one row per gene with its TwasWeights file, as written by twas_weights_conversion; for weights already in that format, a table of the same layout pointing at them. Columns #chr, start, end, region_id, TSS, original_data, contexts.)

#chr	start	end	region_id	TSS	original_data	contexts
chr22	0	18960000	ENSG00000130538	10685239	<path>/output/twas_weights_conversion/protocol_example/protocol_example.ENSG00000130538.twas_weights.s4.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 output/twas_weights_conversion/protocol_example.xqtl_type_table.tsv (maps each xQTL context to its modality, used when grouping contexts for cTWAS; written by twas_weights_conversion, or supplied directly)

context	type
bulk_rnaseq	eQTL

  • --legacy_weights tests/fixtures/twas/protocol_example.twas.legacy_weights.tsv (twas_weights_conversion only) (legacy FunGen-xQTL weight files to convert, one per row. Columns study, data_type, path; rows that share a study and a gene are merged into one weight file.)

  • --legacy_fine_mapping legacy_fine_mapping.tsv (fine_mapping_conversion only) (legacy FunGen-xQTL fine-mapping files to convert, one per gene. Columns study, path; study must match the study of the converted weights.)

  • --fine_mapping_dir output/fine_mapping_conversion (twas, optional) (folder written by fine_mapping_conversion. When set, twas also runs MR for each gene that has a fine-mapping file there; --mr_method (ivwPerVariant or csAware) and --mr_pip_cutoff (default 0.5) control it. Without it, twas reports TWAS only.)

study	data_type	path
protocol_example	eQTL	tests/fixtures/twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.univariate_twas_weights.rds

  • --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)

  • --gwas_study <study_id> (ctwas only: the one GWAS study to fine-map, required when --gwas_meta_data lists several; run each study with its own --cwd)

  • --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_weights_conversion/<study>/<study>.<gene>.twas_weights.s4.rds (twas_weights_conversion) (one pecotmr TwasWeights object per study and gene; one row per context and method, weights copied unchanged from the legacy file)

  • output/twas_weights_conversion/<name>.<study>.xqtl_meta.tsv (twas_weights_conversion) (--xqtl_meta_data input of twas and ctwas, one table per study)

#chr	start	end	region_id	TSS	original_data	contexts
chr22	0	18960000	ENSG00000130538	10685239	<path>/output/twas_weights_conversion/protocol_example/protocol_example.ENSG00000130538.twas_weights.s4.rds	bulk_rnaseq

  • output/twas_weights_conversion/<name>.xqtl_type_table.tsv (twas_weights_conversion) (--xqtl_type_table input: the data type of every converted context)

context	type
bulk_rnaseq	eQTL

  • output/fine_mapping_conversion/<study>/<study>.<gene>.fine_mapping.s4.rds (fine_mapping_conversion) (one pecotmr QtlFineMappingResult per study and gene, one row per context; the SuSiE fit and marginal QTL statistics are copied unchanged from the legacy file. Read by twas through --fine_mapping_dir.)

  • 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 and .gwas_sumstats.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.

twas_weights_conversion, below, is only needed for weights in the legacy FunGen-xQTL format, as the example weights here are; its tables are the --xqtl_meta_data and --xqtl_type_table inputs of the three methods. Weights already stored as pecotmr TwasWeights skip it and go straight to the three methods with a table of the same layout pointing at them.

Convert legacy TWAS weights

twas_weights_conversion converts legacy FunGen-xQTL weight files into pecotmr TwasWeights, one per study and gene, and writes the xQTL meta and data-type tables that the other workflows read.

Timing: ~20 sec (on toy dataset)

sos run pipeline/twas_ctwas.ipynb twas_weights_conversion \
    --cwd output \
    --name protocol_example \
    --legacy_weights tests/fixtures/twas/protocol_example.twas.legacy_weights.tsv

Convert legacy fine-mapping results

fine_mapping_conversion converts legacy FunGen-xQTL fine-mapping files (*.univariate_bvsr.rds) into pecotmr QtlFineMappingResult, one per study and gene, for MR in twas. No toy legacy fine-mapping file ships with this repository, so point --legacy_fine_mapping at a table of your own files (columns study, path).

Timing: ~7 sec per gene

sos run pipeline/twas_ctwas.ipynb fine_mapping_conversion \
    --cwd output \
    --name protocol_example \
    --legacy_fine_mapping legacy_fine_mapping.tsv

Convert GWAS summary statistics

gwas_sumstats_to_tsv is optional. It turns per-LD-block pecotmr GwasSumStats files of one study (for example the QC’d output of the RSS GWAS pipeline, <study>.<chr>_<start>_<end>.gwas_sumstats.rds) into one TSV per chromosome (<study>.chr<N>.tsv.gz) and writes the matching <study>.gwas_meta.tsv to pass as --gwas_meta_data in twas. Point --gwas_sumstats at a text file that lists the RDS files, one per line. No toy GwasSumStats files ship with this repository. A variant present in two blocks is kept once: from a block listed in the optional --prefer_blocks file first, then the measured copy over a RAISS-imputed one, then the block that starts first.

Timing: ~2 min for 1,358 LD blocks (one study)

sos run pipeline/twas_ctwas.ipynb gwas_sumstats_to_tsv \
    --cwd output \
    --gwas_sumstats AD_Bellenguez_main_2026.gwas_sumstats.list \
    --study AD_Bellenguez_main_2026

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. Add --fine_mapping_dir output/fine_mapping_conversion to also run MR, and --rsq_cutoff 0 --rsq_pval_cutoff inf to keep every method per gene instead of only the best one.

For a genome-wide scan, split the work across the cluster: either submit through SoS with -q and a host configuration (job_size, walltime, mem), or split --regions into one file per chromosome and run one job per chromosome, each with its own --cwd.

Timing: ~25 sec (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 output/twas_weights_conversion/protocol_example.protocol_example.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 output/twas_weights_conversion/protocol_example.xqtl_type_table.tsv \
    --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 output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv \
    --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.ctwas.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 output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv \
    --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.ctwas.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 output/twas_weights_conversion/protocol_example.protocol_example.xqtl_meta.tsv \
    --ld_meta_data tests/fixtures/ld_reference/ld_meta_file.ctwas.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 output/twas_weights_conversion/protocol_example.protocol_example.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 output/twas_weights_conversion/protocol_example.xqtl_type_table.tsv \
    --region-name chr22_10000000_19000000

Command Interface

sos run pipeline/twas_ctwas.ipynb -h
usage: sos run 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:
  twas_weights_conversion
  fine_mapping_conversion
  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
  twas_weights_conversion_1:
    Workflow Options:
      --legacy-weights . (as path)
                        Convert legacy FunGen-xQTL TWAS weight RDS files (the
                        list layout released on Synapse: gene ->
                        "<context>_<gene>" -> twas_weights / twas_cv_result /
                        region_info) into pecotmr TwasWeights, one RDS per
                        (study, gene), with twas_weights_conversion.R. Weights
                        and CV metrics are copied unchanged; no LD sketch is
                        attached. --legacy-weights is a TSV with columns study,
                        data_type, path (one legacy RDS per row). Rows sharing a
                        study and a gene (e.g. the univariate and multicontext
                        files of one gene) are merged into one TwasWeights,
                        because the twas step reads one weight file per gene.
                        The gene is the ENSG ID in the legacy file name.
  twas_weights_conversion_2:
  fine_mapping_conversion:
    Workflow Options:
      --legacy-fine-mapping . (as path)
                        Optional. Convert legacy FunGen-xQTL QTL fine-mapping
                        RDS files (*.univariate_bvsr.rds: gene ->
                        "<context>_<gene>" -> susie_result_trimmed / sumstats /
                        variant_names / region_info) into pecotmr
                        QtlFineMappingResult, one RDS per (study, gene), with
                        fine_mapping_conversion.R. The SuSiE fit, marginal QTL
                        effects and PIPs are copied unchanged; no LD sketch is
                        attached. Study and context labels match
                        twas_weights_conversion, so the twas step can pair each
                        fine-mapping row with its weights for MR: pass this
                        step's output folder as --fine_mapping_dir.
                        --legacy_fine_mapping is a TSV with columns study, path
                        (one legacy RDS per gene; the gene is the ENSG ID in the
                        name).
  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; set --rsq_cutoff 0
                        to keep every method) and the MR knobs mr_pval_cutoff,
                        mr_method, mr_pip_cutoff. 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)
      --fine-mapping-dir ''
                        Folder of converted QTL fine-mapping results
                        (fine_mapping_conversion output); empty = TWAS only, no
                        MR
      --mr-method ivwPerVariant
                        MR method (ivwPerVariant or csAware) and the PIP cutoff
                        for MR instruments
      --mr-pip-cutoff 0.5 (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)
      --ctwas-method ''
                        TWAS method fed to cTWAS (assembleCtwasInputs
                        weightSelectionParam=). Empty: each gene's best method
                        by cross-validated r-squared; a gene whose methods
                        carry no r-squared needs one, e.g. --ctwas-method susie.
      --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  paths()

                        TwasWeights source (option c): empty -> the upstream
                        twas step's per-gene weights (the per-gene TwasWeights
                        files in the TWAS manifest); override with prebuilt
                        TwasWeights RDS paths (e.g. the blessed S4 ctwas
                        weights) when the toy twas weights carry no signal. A
                        quoted wildcard is expanded.
      --twas-weights-list . (as path)
                        Or a file listing those RDS paths, one per line: a
                        genome-wide list is too long for the command line.
  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
# For thousands of regions use --regions instead: the names reach twas_manifest.R
# as one comma-joined argument, which Linux caps at 128 KiB.
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=[]


[twas_weights_conversion_1]
# Convert legacy FunGen-xQTL TWAS weight RDS files (the list layout released on
# Synapse: gene -> "<context>_<gene>" -> twas_weights / twas_cv_result /
# region_info) into pecotmr TwasWeights, one RDS per (study, gene), with
# twas_weights_conversion.R. Weights and CV metrics are copied unchanged; no LD
# sketch is attached. --legacy-weights is a TSV with columns study, data_type,
# path (one legacy RDS per row). Rows sharing a study and a gene (e.g. the
# univariate and multicontext files of one gene) are merged into one
# TwasWeights, because the twas step reads one weight file per gene. The gene
# is the ENSG ID in the legacy file name.
parameter: legacy_weights = path()
import csv, re
_legacy_rows = list(csv.DictReader(open(legacy_weights), delimiter = '\t'))
stop_if(len(_legacy_rows) == 0, f"No rows in {legacy_weights}.")
_conv_units = {}
for _r in _legacy_rows:
    _g = re.search(r'ENSG[0-9]+', path(_r['path']).name)
    stop_if(_g is None, f"No ENSG gene ID in legacy file name: {_r['path']}")
    _conv_units.setdefault((_r['study'], _g.group(0)), []).append(_r['path'])
_conv_keys = sorted(_conv_units)
_conv_types = {_r['study']: _r['data_type'] for _r in _legacy_rows}
input: [f for k in _conv_keys for f in _conv_units[k]], group_by = lambda x: [_conv_units[k] for k in _conv_keys]
output: f'{cwd:a}/twas_weights_conversion/{_conv_keys[_index][0]}/{_conv_keys[_index][0]}.{_conv_keys[_index][1]}.twas_weights.s4.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
    Rscript ${modular_script_dir}/pecotmr_integration/twas_weights_conversion.R \
        --legacy "${",".join([str(x) for x in _input])}" \
        --study "${_conv_keys[_index][0]}" \
        --data-type "${_conv_types[_conv_keys[_index][0]]}" \
        --output "${_output}" \
        --meta-row "${_output:n}.meta_row.tsv" \
        --type-rows "${_output:n}.type_rows.tsv"
[twas_weights_conversion_2]
# Collect the per-gene rows written by step 1 into one xqtl_meta table per study
# ({name}.<study>.xqtl_meta.tsv, the --xqtl-meta-data input of twas and ctwas)
# and one context -> data type table ({name}.xqtl_type_table.tsv, the
# --xqtl-type-table input). The meta tables are per study because twas uses
# the first weight file it finds for each gene: run twas once per study.
input: group_by = 'all'
_conv_studies = sorted(set(path(x).parent.name for x in _input))
output: [f'{cwd:a}/twas_weights_conversion/{name}.{s}.xqtl_meta.tsv' for s in _conv_studies] + [f'{cwd:a}/twas_weights_conversion/{name}.xqtl_type_table.tsv']
_ctx_type = {}
for _s, _out in zip(_conv_studies, _output[:-1]):
    _rows = []
    for _x in _input:
        if path(_x).parent.name != _s:
            continue
        _rows += [l.rstrip('\n').split('\t') for l in open(f'{_x:n}.meta_row.tsv') if l.strip()]
        for l in open(f'{_x:n}.type_rows.tsv'):
            if l.strip():
                _c, _t = l.rstrip('\n').split('\t')
                _ctx_type[_c] = _t
    _rows.sort(key = lambda r: (r[0], int(r[1]), r[3]))
    with open(_out, 'w') as _f:
        _f.write('#chr\tstart\tend\tregion_id\tTSS\toriginal_data\tcontexts\n')
        _f.writelines('\t'.join(r) + '\n' for r in _rows)
with open(_output[-1], 'w') as _f:
    _f.write('context\ttype\n')
    _f.writelines(f'{c}\t{t}\n' for c, t in sorted(_ctx_type.items()))
[fine_mapping_conversion]
# Optional. Convert legacy FunGen-xQTL QTL fine-mapping RDS files
# (*.univariate_bvsr.rds: gene -> "<context>_<gene>" -> susie_result_trimmed /
# sumstats / variant_names / region_info) into pecotmr QtlFineMappingResult,
# one RDS per (study, gene), with fine_mapping_conversion.R. The SuSiE fit,
# marginal QTL effects and PIPs are copied unchanged; no LD sketch is attached.
# Study and context labels match twas_weights_conversion, so the twas step can
# pair each fine-mapping row with its weights for MR: pass this step's output
# folder as --fine_mapping_dir. --legacy_fine_mapping is a TSV with columns
# study, path (one legacy RDS per gene; the gene is the ENSG ID in the name).
parameter: legacy_fine_mapping = path()
import csv, re
_fm_rows = list(csv.DictReader(open(legacy_fine_mapping), delimiter = '\t'))
stop_if(len(_fm_rows) == 0, f"No rows in {legacy_fine_mapping}.")
_fm_keys = []
for _r in _fm_rows:
    _g = re.search(r'ENSG[0-9]+', path(_r['path']).name)
    stop_if(_g is None, f"No ENSG gene ID in legacy file name: {_r['path']}")
    _fm_keys.append((_r['study'], _g.group(0)))
fail_if(len(set(_fm_keys)) < len(_fm_keys), f"Each (study, gene) must appear once in {legacy_fine_mapping}.")
input: [_r['path'] for _r in _fm_rows], group_by = 1
output: f'{cwd:a}/fine_mapping_conversion/{_fm_keys[_index][0]}/{_fm_keys[_index][0]}.{_fm_keys[_index][1]}.fine_mapping.s4.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
    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_conversion.R \
        --legacy "${_input}" \
        --study "${_fm_keys[_index][0]}" \
        --output "${_output}"

[gwas_sumstats_to_tsv]
# Optional. Export per-LD-block pecotmr GwasSumStats RDS files of one study (e.g. the
# QC'd output of the RSS GWAS pipeline) into one bgzipped, tabix-indexed TSV per
# chromosome plus a gwas_meta TSV that the twas step reads through --gwas_meta_data,
# with gwas_sumstats_to_tsv.R. --gwas_sumstats is a text file listing the RDS files,
# one per line; --study is the study_id written to the gwas_meta file.
# Optional --prefer_blocks: text file of block names (chr<N>_<start>_<end>) whose copy
# of a variant shared with another block is kept first.
parameter: gwas_sumstats = path()
parameter: study = str
parameter: prefer_blocks = ""
input: gwas_sumstats
output: f'{cwd:a}/gwas_sumstats_to_tsv/{study}.gwas_meta.tsv'
task: trunk_workers = 1, trunk_size = 1, 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
    Rscript ${modular_script_dir}/pecotmr_integration/gwas_sumstats_to_tsv.R \
        --input "${_input}" \
        --study "${study}" \
        --output-dir "${_output:d}" \
        --prefer-blocks "${prefer_blocks}" \
        --meta "${_output}"
[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 also runs when --fine_mapping_dir holds a
# converted fine-mapping file for the gene's study
# (<dir>/<study>/<study>.<gene>.fine_mapping.s4.rds, see
# fine_mapping_conversion); otherwise 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; set
# --rsq_cutoff 0 --rsq_pval_cutoff inf to keep every method) and the MR knobs mr_pval_cutoff,
# mr_method, mr_pip_cutoff. 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
# Folder of converted QTL fine-mapping results (fine_mapping_conversion output);
# empty = TWAS only, no MR
parameter: fine_mapping_dir = ""
# MR method (ivwPerVariant or csAware) and the PIP cutoff for MR instruments
parameter: mr_method = "ivwPerVariant"
parameter: mr_pip_cutoff = 0.5
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))])}
    # The per-gene paths go in a table, one row per gene: a chromosome-wide
    # region's lists are too long for the command line.
    gene_table="$outdir/${name}.$region_id.twas_genes.tsv"
    printf 'twas_weights\tfine_mapping_result\toutput\n' > "$gene_table"
    while [ "$#" -ge 2 ]; do
        gene="$1"
        weight="$2"
        shift 2
        study=$(basename "$weight" | cut -d. -f1)
        printf '%s\t%s\t%s\n' "$weight" \
            "${fine_mapping_dir}/$study/$study.$gene.fine_mapping.s4.rds" \
            "$outdir/${name}.$region_id.$gene.twas.rds" >> "$gene_table"
    done
    # One R call per region, so the region GwasSumStats (and its LD reference) is read once for every gene.
    Rscript ${modular_script_dir}/pecotmr_integration/twas.R \
        --gene-table "$gene_table" \
        --gwas-sumstats "$gss" \
        --rsq-cutoff ${rsq_cutoff} \
        --rsq-pval-cutoff ${"Inf" if rsq_pval_cutoff == float("inf") else rsq_pval_cutoff} \
        --rsq-option ${rsq_option} \
        --rsq-pval-option "${",".join(rsq_pval_option)}" \
        --mr-pval-cutoff ${mr_pval_cutoff} \
        --mr-method ${mr_method} \
        --mr-pip-cutoff ${mr_pip_cutoff}
[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
# TWAS method fed to cTWAS (assembleCtwasInputs weightSelectionParam=). Empty
# uses each gene's best method by cross-validated r-squared; a gene whose
# methods carry no r-squared needs this set explicitly, e.g. --ctwas-method susie.
parameter: ctwas_method = ""
# 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 (the per-gene TwasWeights files in the TWAS manifest); override with
# prebuilt TwasWeights RDS paths (e.g. the blessed S4 ctwas weights) when the
# toy twas weights carry no signal. A quoted wildcard is expanded.
parameter: twas_weights = paths()
# Or a file listing those RDS paths, one per line: a genome-wide list is too
# long for the command line.
parameter: twas_weights_list = path()
skip_if(skip_assembly == True, "Skip [ctwas_1] assemble.")
fail_if(len(twas_weights) > 0 and twas_weights_list != path(), "Pass --twas-weights or --twas-weights-list, not both.")
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}.")
# paths() expands wildcards in strings but not in path objects, so only
# --twas-weights globs; the list file's and the manifest's paths are taken as is.
if twas_weights_list != path():
    ctwas_weights = paths([path(x.strip()) for x in open(twas_weights_list) if x.strip()])
elif twas_weights:
    ctwas_weights = paths([str(w) for w in twas_weights])
else:
    ctwas_weights = paths([path(w) for w in dict.fromkeys(f for j in jobs if j["chrom"] == ctwas_chrom for f in j["weight_files"].split(',') if f)])
_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 ""
# assembleCtwasInputs() takes one GWAS per LD block, so cTWAS runs one study at
# a time: --gwas-study picks it when --gwas-meta-data lists several. The
# per-block GWAS files are named by region only, so give each study its own --cwd.
if gwas_study:
    _keep = [i for i, s in enumerate(ctwas_studies) if s in gwas_study]
    ctwas_gwas_files = ",".join([ctwas_gwas_files.split(",")[i] for i in _keep])
    ctwas_studies = [ctwas_studies[i] for i in _keep]
stop_if(len(ctwas_studies) == 0, f"No GWAS study covers {ctwas_chrom}.")
fail_if(len(ctwas_studies) > 1, f"cTWAS takes one GWAS study at a time; {ctwas_chrom} has {ctwas_studies}. Pass --gwas-study.")
# ctwasPipeline needs the LD panel as a real genotype file: a block built from the
# chromosome->path mapping of --ld-meta-data keeps only a "<chrom-meta>" placeholder.
# Resolve this chromosome's genotype prefix and pass it as --ld-sketch when unique.
import csv as _csv, os as _os
_ld_rows = list(_csv.DictReader(open(ld_meta_data), delimiter = "\t"))
_ld_chr = [x for x in ("#chr", "#chrom", "chr", "chrom") if _ld_rows and x in _ld_rows[0]]
_ld_dir = _os.path.dirname(_os.path.abspath(ld_meta_data))
_ld_paths = sorted(set(_os.path.join(_ld_dir, r["path"].split(",")[0]) for r in _ld_rows if _ld_chr and r[_ld_chr[0]].replace("chr", "") == ctwas_chrom.replace("chr", "")))
ctwas_ld_sketch = _ld_paths[0] if len(_ld_paths) == 1 and any(_os.path.exists(_ld_paths[0] + e) for e in (".pgen", ".bed")) else ""
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}" \
            ${f'--ld-sketch "{ctwas_ld_sketch}" ' if ctwas_ld_sketch else ""}--output "$gwas_rds"
    done

    # Drop LD blocks with no GWAS variants: gwas_sumstats_construct.R writes no
    # GwasSumStats for them, and cTWAS cannot use an empty block.
    kept="$outdir/${name}.ctwas_manifest.${ctwas_chrom}.nonempty.tsv"
    head -n 1 "$manifest" > "$kept"
    tail -n +2 "$manifest" | while IFS=$'\t' read -r region_id region gwas_rds; do
      if [ -s "$gwas_rds" ]; then
        printf '%s\t%s\t%s\n' "$region_id" "$region" "$gwas_rds" >> "$kept"
      else
        echo "Skipping LD block $region_id: no GWAS variants" >&2
      fi
    done

    # (3) assemble cTWAS inputs (pecotmr places each gene into its home block).
    # The weight paths go in a file: a genome-scale list overflows the 128 KB per-argument limit.
    printf '%s\n' ${ctwas_weights:r} > "${_output:n}.weights.txt"
    Rscript ${modular_script_dir}/pecotmr_integration/ctwas_assemble.R \
        --manifest "$kept" \
        --twas-weights-list "${_output:n}.weights.txt" \
        ${f"--method {ctwas_method} " if ctwas_method else ""}\
        --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
  1. 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
GWAS Integration
xQTL-GWAS pairwise enrichment and colocalization
GWAS Integration
TWAS and COLOC integration (INTACT)