Causal TWAS Pipeline (cTWAS, multi LD block)

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 references 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,
  maxNumSingleEffects = 5L,
  numThreads = 1L,
  ctwasPriorParam = CtwasPriorParam(),
  variantPruningParam = VariantPruningParam(),
  boundaryMergeParam = BoundaryMergeParam(),
  keepSnps = FALSE,
  methodArgs = CtwasOptions()
)

Arguments

gwasSumStats: A GwasSumStats whose elements are LD blocks (at least two), keyed by its blockId column, with qcInfo() 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.
maxNumSingleEffects: Integer. Max number of single effects in cTWAS fine-mapping. Default 5. Pass-through to ctwas::ctwas_sumstats.
numThreads: Number of cores. Default 1.
ctwasPriorParam: The EM stage that estimates cTWAS's group priors, built with CtwasPriorParam: thin, niterPrefit, niter, varStructure and fallbackToPrefit. Handed whole to estCtwasGroupPriors.
variantPruningParam: Which variants survive into each gene's weight matrix, built with VariantPruningParam: weightCutoff, csMinCor, pipRescue and maxVariants. Handed whole to assembleCtwasInputs. The csMinCor / pipRescue rescues need a fineMappingResult.
boundaryMergeParam: Whether a high-PIP gene straddling an LD-block boundary has its regions merged and re-fine-mapped, built with BoundaryMergeParam: enabled (default FALSE), pipThreshold, filterCs and maxSnp.
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 ctwasPosteriors / susieAlpha. Default FALSE – the SNP rows are the null background and are dropped from the structured gene-level result.
methodArgs: Additional arguments forwarded to ctwas, built with CtwasOptions. Names are checked against what the ctwas steps accept between them.

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\) estCtwasGroupPriors \(\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 reference 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.

cTWAS's log

cTWAS narrates each of its steps through the logging package, which writes to the console with cat(), out of reach of suppressMessages(). The pipeline keeps that narration off the console, cTWAS's warnings excepted. To keep the full log, give a step a log file through its options, e.g. CtwasEstParamOptions(logfile = "est_param.log").

LD reference compatibility check

Per block: ldReference(twasWeights) (when non-NULL) must be the same LD reference as ldReference(gwasSumStats). The two must be the same kind of reference. Two genotype panels must share their samples and allele orientation; two LdMatrix objects must agree on the LD they share; two GenotypeMixture objects must average the same panels with the same weights. The two need NOT carry the same variants; a partial overlap, which is what QC-ing the two sides separately produces, only warns. A mismatch of any of those, or no shared variant at all, 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",
  ldReference = 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,
  panelFilterParam = PanelFilterParam(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),
  ctwasPriorParam = CtwasPriorParam(thin = 1, niterPrefit = 3, niter = 10,
    fallbackToPrefit = TRUE),
  methodArgs = CtwasOptions(
    regionData = CtwasRegionDataOptions(min_group_size = 1),
    estParam = CtwasEstParamOptions(
      min_group_size = 1, min_p_single_effect = 0
    )
  ))
#> CtwasResult: 1 run(s)
#>   1 GWAS stud(y/ies), 1 method(s): susie