Skip to contents

Causal TWAS Pipeline (cTWAS, multi LD block)

Source: R/ctwasPipeline.R
ctwasPipeline.Rd

Pipeline that hands a per-block set of GwasSumStats of GWAS Z-scores together with the matching per-block per-gene TWAS weights and LD sketches to ctwas::ctwas_sumstats, producing per-gene posterior inclusion probabilities for causal genes. Optionally accepts a precomputed TWAS-Z GRanges from causalInferencePipeline as the z_gene input so the per-gene Z is not recomputed inside ctwas.

Usage

ctwasPipeline(
  gwasSumStats,
  twasWeights,
  twasZ = NULL,
  fineMappingResult = NULL,
  method = NULL,
  thin = 0.1,
  niterPrefit = 3L,
  niter = 30L,
  L = 5L,
  groupPriorVarStructure = c("shared_type", "shared_context", "shared_nonSNP",
    "shared_all", "independent"),
  ncore = 1L,
  twasWeightCutoff = 0,
  csMinCor = 0.8,
  minPipCutoff = 0,
  maxNumVariants = Inf,
  fallbackToPrefit = FALSE,
  keepSnps = FALSE,
  mergeBoundary = FALSE,
  mergePipThresh = 0.5,
  mergeFilterCs = FALSE,
  mergeMaxSNP = Inf,
  ...
)

Arguments

gwasSumStats

A GwasSumStats whose elements are LD blocks (at least two), keyed by its blockId column, with getQcInfo() non-empty. Build it with loadGwasSumStatsFromManifest(..., ldBlocks = <blocks>) and pass it through summaryStatsQc.

twasWeights

The per-gene weight source. Either (a) a FLAT TwasWeights / QtlFineMappingResult (or a homogeneous list of them) carrying region provenance – each gene is placed into its home LD block internally by start(region) (matching cTWAS's p0 assignment rule); or (b) a pre-bucketed NAMED LIST keyed by region_id (keys a SUBSET of gwasSumStats's), used as-is. Blocks without any TWAS weights still contribute their SNP-level signal to ctwas's joint group prior estimate (the legacy whole-chromosome pattern where only a few of many LD blocks carry gene weights). A gene whose cis span straddles a block boundary is homed by its single anchor; the cross-block signal is cTWAS's boundary-gene concern (mergeCtwasBoundaryRegions), not placement.

twasZ

Optional GRanges of TWAS Z-scores (output of causalInferencePipeline). When supplied, the per-(trait, context) Z is used as the z_gene input to ctwas_sumstats so it is not recomputed.

fineMappingResult

Optional QtlFineMappingResult or GwasFineMappingResult carrying the per-variant PIP and credible-set membership data used by the CS / PIP rescue filters (csMinCor and minPipCutoff). When NULL (default) the smart filters are no-ops; only the magnitude filter (twasWeightCutoff) and the per-gene cap (maxNumVariants, ordered by |weight|) apply.

method

Optional character (length 1). Picks which TWAS method's weights to feed into ctwas for each (study, context, trait) gene. When NULL (default): use "ensemble" if that method is present across the weight sources; otherwise use the sole method when only one is present; otherwise run every method as an independent cTWAS run (one CtwasResult row-set per method). Passing the name explicitly (e.g. "mrash") restricts the run to that single method.

thin, niterPrefit, niter, L

Pass-throughs to ctwas::ctwas_sumstats.

groupPriorVarStructure

Pass-through (defaults "shared_type").

ncore

Number of cores. Default 1.

twasWeightCutoff

Numeric (length 1). Drop variants with |weight| < twasWeightCutoff from each gene's weight matrix before ctwas sees it. Default 0 (no filter).

csMinCor

Numeric (length 1). When fineMappingResult is provided, variants belonging to any 95% credible set with purity (min_abs_corr) >= csMinCor are marked as must-keep and survive the per-gene cap. Default 0.8. Ignored without a fineMappingResult.

minPipCutoff

Numeric (length 1). When fineMappingResult is provided, variants with PIP greater than minPipCutoff are marked as must-keep and survive the per-gene cap. Default 0 (no PIP rescue). Ignored without a fineMappingResult.

maxNumVariants

Numeric (length 1). Cap on per-gene variant count. When the gene has more variants than this, keep all must-keep variants and fill remaining slots by descending PIP (when available) or descending |weight|. Default Inf (no cap).

fallbackToPrefit

Logical (length 1). Forwarded to estCtwasParam. When TRUE, ctwas's accurate-EM NaN failure is recovered by falling back to the prefit estimates (mirrors the legacy ctwas_2 workaround on underpowered data). Default FALSE.

keepSnps

Logical (length 1). When TRUE, retain the context-agnostic SNP background of each run as one extra CtwasResult row (study = context = "SNP", mirroring cTWAS's own "SNP" group) so the full ctwas output is reconstructable from getFinemap / getSusieAlpha. Default FALSE – the SNP rows are the null background and are dropped from the structured gene-level result.

mergeBoundary

Logical (length 1). When TRUE, run mergeCtwasBoundaryRegions after fine-mapping each run: a high-PIP gene whose cis window straddles an LD-block boundary has its adjacent regions merged and re-fine-mapped (the legacy default-off ctwas_3 post-processing). Default FALSE.

mergePipThresh

Numeric (length 1). PIP threshold for selecting which boundary genes to merge (mergeCtwasBoundaryRegions pipThresh). Default 0.5. Ignored unless mergeBoundary = TRUE.

mergeFilterCs

Logical (length 1). Require the boundary gene to be in a credible set to be selected. Default FALSE. Ignored unless mergeBoundary = TRUE.

mergeMaxSNP

Numeric (length 1). Per-merged-region SNP cap. Default Inf. Ignored unless mergeBoundary = TRUE.

...

Additional arguments forwarded to ctwas::ctwas_sumstats.

Value

A CtwasResult collection: one row per (gwasStudy, study, context, method). A single-context run is one row per method; a multi-context (joint) run emits per-context rows sharing the same jointContexts set and the jointly-estimated group priors. Each row's CtwasResultEntry payload carries that context's per-gene fine-mapping posteriors (finemap), the run's param, and its regionInfo. For the raw ctwas::finemap_regions list (e.g. to feed mergeCtwasBoundaryRegions), call the granular assembleCtwasInputs \(\to\) estCtwasParam \(\to\) screenCtwasRegions \(\to\) finemapCtwasRegions path instead.

LD block convention

gwasSumStats is ONE GwasSumStats whose elements are LD blocks, keyed by its blockId column – build it with loadGwasSumStatsFromManifest(..., ldBlocks = <blocks>). Per-block region_info, LD_map, and snp_map entries are built automatically from the LD sketch and concatenated before the call to ctwas::ctwas_sumstats. A single-block input is rejected: cTWAS's EM cannot converge on a single region, so callers must supply at least two blocks.

LD-sketch identity check

Per block: getLdSketch(twasWeights) (when non-NULL) must match getLdSketch(gwasSumStats). Mismatch is a hard error.

Examples

data(ctwasWeightsExample)
ldDir <- system.file("extdata", "ld_reference", "chr22",
  package = "pecotmr")
ldStem <- file.path(ldDir, "protocol_example.LD.chr22")
gwasTsv <- system.file("extdata", "manifests",
  "protocol_example.twas.gwas_sumstats.chr22.tsv.gz", package = "pecotmr")
mani <- data.frame(study = "gwas1", sumStatsPath = gwasTsv)
blocks <- GenomicRanges::GRanges("chr22",
  IRanges::IRanges(c(10000000, 15000001), c(15000000, 19000000)),
  blockId = c("chr22_1", "chr22_2"))
gss <- loadGwasSumStatsFromManifest(manifest = mani, genome = "hg38",
  ldSketch = ldStem, region = "chr22:10000000-19000000", ldBlocks = blocks)
#> Warning: GwasSumStats[study=gwas1]: no effect-allele frequency declared (map `af: <col>` to export a directional af); top_loci$af will be NA. A directionless `maf`/`FRQ` is used for QC only, never as af.
gwasByRegion <- summaryStatsQc(gss, mafCutoff = 0.0025)
#> [gwas1] QC track: harmonization kept 20 of 20 variant(s) (corrected: sign-flipped 0, strand-flipped 0; dropped 0).
#> [gwas1] QC summary: 20 in -> 20 out | corrected: sign-flip 0, strand-flip 0
#> [gwas1] QC track: harmonization kept 155 of 155 variant(s) (corrected: sign-flipped 0, strand-flipped 0; dropped 0).
#> [gwas1] QC summary: 155 in -> 155 out | corrected: sign-flip 0, strand-flip 0
ctwasPipeline(gwasSumStats = gwasByRegion,
  twasWeights = list(ctwasWeightsExample), thin = 1, niterPrefit = 3,
  niter = 10, min_group_size = 1, min_p_single_effect = 0,
  fallbackToPrefit = TRUE)
#> 2026-09-04 23:05:28.661646 INFO::Computing gene z-scores...
#> 2026-09-04 23:05:28.666125 INFO::Assembling region_data...
#> 2026-09-04 23:05:28.669425 INFO::Number of regions in total: 2
#> 2026-09-04 23:05:28.669643 INFO::thin = 1
#> 2026-09-04 23:05:28.670388 INFO::Number of SNPs in snp_map: 175
#> 2026-09-04 23:05:28.6706 INFO::Number of SNPs in GWAS (z_snp$id): 175
#> 2026-09-04 23:05:28.6708 INFO::Average number of SNPs in weights per molecular trait: 18.00
#> 2026-09-04 23:05:28.670997 INFO::Total number of SNPs in weights: 18
#> 2026-09-04 23:05:28.671168 INFO::10.29% GWAS SNPs (z_snp$id) are in weights.
#> 2026-09-04 23:05:28.672325 INFO::Number of regions in chr22: 2
#> 2026-09-04 23:05:28.672659 INFO::Map gene regions...
#> 2026-09-04 23:05:28.673101 INFO::Number of boundary genes: 1
#> 2026-09-04 23:05:28.673373 INFO::Adjust region assignment for boundary genes
#> 2026-09-04 23:05:28.673716 INFO::Updating region z-scores...
#> 2026-09-04 23:05:28.674798 INFO::Get boundary genes...
#> 2026-09-04 23:05:28.675302 INFO::Map gene regions...
#> 2026-09-04 23:05:28.676119 INFO::Estimating parameters...
#> 2026-09-04 23:05:28.676368 INFO::thin = 1
#> 2026-09-04 23:05:28.677226 INFO::Skip 1 regions with < 1 genes.
#> 2026-09-04 23:05:28.677451 INFO::Run EM prefit iterations, getting rough estimates...
#> 2026-09-04 23:05:28.677622 INFO::group_prior_var_structure = 'shared_type'
#> 2026-09-04 23:05:28.677783 INFO::Using data in 1 regions
#> 2026-09-04 23:05:28.681031 INFO::Roughly estimated group_prior {bulk_rnaseq|bulk_rnaseq,SNP}: {0.01569,0.02411}
#> 2026-09-04 23:05:28.681302 INFO::Roughly estimated group_prior_var {bulk_rnaseq|bulk_rnaseq,SNP}: {2.618,4.831}
#> 2026-09-04 23:05:28.681765 INFO::Selected 1 regions with p(single effect) > 0 to run EM
#> 2026-09-04 23:05:28.682003 INFO::Run EM iterations, getting accurate estimates...
#> 2026-09-04 23:05:28.682186 INFO::Using data in 1 regions
#> 2026-09-04 23:05:28.684472 INFO::EM converged after 3 iterations.
#> 2026-09-04 23:05:28.685253 INFO::Estimated group_prior {bulk_rnaseq|bulk_rnaseq,SNP}: {0.01327,0.02401}
#> 2026-09-04 23:05:28.685495 INFO::Estimated group_prior_var {bulk_rnaseq|bulk_rnaseq,SNP}: {2.368,4.549}
#> 2026-09-04 23:05:28.685682 INFO::group_size {bulk_rnaseq|bulk_rnaseq,SNP}: {1,173}
#> 2026-09-04 23:05:28.686431 INFO::Screening regions...
#> 2026-09-04 23:05:28.687037 INFO::Skip 1 regions with < 1 genes.
#> 2026-09-04 23:05:28.68725 INFO::Selected 0 regions with GWAS min. pval < 5e-08.
#> 2026-09-04 23:05:28.687454 INFO::Screening 1 regions...
#> 2026-09-04 23:05:28.688623 INFO::Selected 0 regions with non-SNP PIP > 0.5.
#> 2026-09-04 23:05:28.688847 INFO::No regions selected after screening.
#> CtwasResult: 1 run(s)
#>   1 GWAS stud(y/ies), 1 method(s): susie