diff --git a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb index e9dc2f16e..e435ffd48 100644 --- a/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb +++ b/code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb @@ -1073,6 +1073,12 @@ "# weights pass then reuses its SuSiE fits via --fine-mapping-result\n", "# (pecotmr_integration/twas_weights.R -> twasWeightsPipeline). Replaces the\n", "# legacy regional_data + susie_twas.R path.\n", + "#\n", + "# Both passes get the same --cv-folds: a cross-validated fine-mapping-derived\n", + "# weight method (susie) needs each fold's OWN fine-mapping fit, so the\n", + "# fine-mapping pass must be cross-validated for the TWAS pass to reuse it.\n", + "# pecotmr then reuses that fold partition for every other weight method, so\n", + "# all of them are scored on the same folds.\n", "parameter: cis_window = 1000000\n", "parameter: fine_mapping_methods = \"susie\"\n", "parameter: twas_methods = \"susie,mrash,enet,lasso,mcp,scad,l0learn,bayes_r,bayes_c\"\n", @@ -1140,6 +1146,8 @@ " --methods ${fine_mapping_methods} \\\n", " --coverage ${fine_mapping_coverage} \\\n", " --pip-cutoff-to-skip ${pip_cutoff_to_skip} \\\n", + " --cv-folds ${cv_folds} \\\n", + " --cv-threads ${cv_threads} \\\n", " ${seed_arg} ${filter_overrides} ${fm_maf_xvar} \\\n", " ${(\"--contexts \" + contexts) if contexts else \"\"} \\\n", " ${(\"--method-args '\" + fine_mapping_method_args + \"'\") if fine_mapping_method_args else \"\"} \\\n", diff --git a/code/script/mnm_analysis/mnm_methods/fsusie.R b/code/script/graveyard/fsusie.R similarity index 93% rename from code/script/mnm_analysis/mnm_methods/fsusie.R rename to code/script/graveyard/fsusie.R index 0bfd7a9ee..490f7d293 100644 --- a/code/script/mnm_analysis/mnm_methods/fsusie.R +++ b/code/script/graveyard/fsusie.R @@ -1,4 +1,9 @@ #!/usr/bin/env Rscript +# RETIRED: superseded by pecotmr_integration/fine_mapping.R (fineMappingPipeline's +# fsusie / mvsusie methods, and --use-pca for susie_on_top_pc). Nothing invokes this +# file. It targets the PRE-camelCase pecotmr API -- load_regional_functional_data(), +# parse_region(), susie_wrapper(), susie_post_processor(), twas_weights_pipeline(), +# fsusie_wrapper() -- none of which exist any more, so it cannot run as written. # fsusie.R — Functional SuSiE fine-mapping for epigenomic QTL # Called directly by current SoS notebook task blocks. # @@ -131,7 +136,10 @@ if (length(fdat$Y_coordinates) == 0) { } if (opt[["save-data"]]) { - saveRDS(list(opt$region = fdat), + # One element NAMED BY THE VALUE of opt$region -- list(opt$region = ...) is a + # syntax error, since the tag has to be a literal name (the same intent the + # result[[opt$region]] <- idiom above expresses). + saveRDS(setNames(list(fdat), opt$region), paste0(prefix, ".", mark_thresh, "_marks.dataset.rds"), compress="xz") } diff --git a/code/script/molecular_phenotypes/h3k9ac_calling/macs2Peaks/scripts/lib/encodeMetrics.R b/code/script/molecular_phenotypes/h3k9ac_calling/macs2Peaks/scripts/lib/encodeMetrics.R index 56be0f606..ef47126ff 100644 --- a/code/script/molecular_phenotypes/h3k9ac_calling/macs2Peaks/scripts/lib/encodeMetrics.R +++ b/code/script/molecular_phenotypes/h3k9ac_calling/macs2Peaks/scripts/lib/encodeMetrics.R @@ -19,6 +19,9 @@ encodeMetrics <- function(sample, bed, peaks, blacklist) { require(rtracklayer) + # seqlevelsStyle()/keepSeqlevels() are GenomeInfoDb's; GenomicRanges now + # Depends on Seqinfo instead, so rtracklayer no longer attaches them. + require(GenomeInfoDb) fragments <- import(bed, format="BED") peaks <- import(peaks, format="BED") blacklist <- import(blacklist, format="BED") @@ -82,6 +85,9 @@ encodeMetrics <- function(sample, bed, peaks, blacklist) { crossCorrelation <- function(bed, blacklist, rmDup=FALSE) { require(rtracklayer) + # seqlevelsStyle()/keepSeqlevels() are GenomeInfoDb's; GenomicRanges now + # Depends on Seqinfo instead, so rtracklayer no longer attaches them. + require(GenomeInfoDb) require(epigenomix) fragments <- import(bed, format="BED") blacklist <- import(blacklist, format="BED") @@ -129,6 +135,9 @@ crossCorrelation <- function(bed, blacklist, rmDup=FALSE) { chromatinStateOverlap <- function(sample, bed, chromState, blacklist, rmDup=FALSE) { require(rtracklayer) + # seqlevelsStyle()/keepSeqlevels() are GenomeInfoDb's; GenomicRanges now + # Depends on Seqinfo instead, so rtracklayer no longer attaches them. + require(GenomeInfoDb) epigenome <- strsplit(basename(chromState), "_")[[1]][1] fragments <- import(bed, format="BED") chromState <- import(chromState, format="BED") diff --git a/code/script/molecular_phenotypes/h3k9ac_calling/values/counts/scripts/01_h3k9acDomains.R b/code/script/molecular_phenotypes/h3k9ac_calling/values/counts/scripts/01_h3k9acDomains.R index 5e11cc005..ccfcf0394 100644 --- a/code/script/molecular_phenotypes/h3k9ac_calling/values/counts/scripts/01_h3k9acDomains.R +++ b/code/script/molecular_phenotypes/h3k9ac_calling/values/counts/scripts/01_h3k9acDomains.R @@ -2,6 +2,7 @@ suppressPackageStartupMessages({ library(argparser) library(vroom) library(GenomicRanges) + library(GenomeInfoDb) # seqlevelsStyle(), sortSeqlevels() library(rtracklayer) }) diff --git a/code/script/molecular_phenotypes/h3k9ac_calling/values/domains/scripts/annotateDomains.R b/code/script/molecular_phenotypes/h3k9ac_calling/values/domains/scripts/annotateDomains.R index 765987c85..81b858c6b 100644 --- a/code/script/molecular_phenotypes/h3k9ac_calling/values/domains/scripts/annotateDomains.R +++ b/code/script/molecular_phenotypes/h3k9ac_calling/values/domains/scripts/annotateDomains.R @@ -3,6 +3,8 @@ suppressPackageStartupMessages({ library(rtracklayer) library(MatrixGenerics) library(GenomicRanges) + library(GenomeInfoDb) # seqlevelsStyle() + library(SummarizedExperiment) # assay() }) p <- arg_parser("Annotate H3K9ac domains with count stats + Roadmap chromatin-state overlap and cluster them") diff --git a/code/script/molecular_phenotypes/h3k9ac_calling/values/residuals/scripts/calculateResiduals.R b/code/script/molecular_phenotypes/h3k9ac_calling/values/residuals/scripts/calculateResiduals.R index ebbb82d2b..1e65ffbd7 100644 --- a/code/script/molecular_phenotypes/h3k9ac_calling/values/residuals/scripts/calculateResiduals.R +++ b/code/script/molecular_phenotypes/h3k9ac_calling/values/residuals/scripts/calculateResiduals.R @@ -7,6 +7,7 @@ suppressPackageStartupMessages({ library(argparser) + library(SummarizedExperiment) # assay(), colData(), rowData() library(limma) library(edgeR) }) diff --git a/code/script/pecotmr_integration/coloc.R b/code/script/pecotmr_integration/coloc.R index 79f4a717a..c474cd7c5 100644 --- a/code/script/pecotmr_integration/coloc.R +++ b/code/script/pecotmr_integration/coloc.R @@ -6,8 +6,8 @@ # inline fine-mapping when given sumstats), and calls # `pecotmr::colocPipeline()`. When `--enrichment` points at the output # of `qtl_enrichment.R`, the per-pair `p12` prior is scaled by -# `(1 + enrichment)` (capped at `--p12-max`), the "enloc" variant of -# coloc. +# `(1 + enrichment)`, the "enloc" variant of coloc. (pecotmr 0.9.0 removed +# the cap on that adjusted prior, so the former --p12-max flag is gone.) # # Inputs: # --qtl-fine-mapping Path to S4 QtlFineMappingResult RDS @@ -20,7 +20,6 @@ # --p1 Per-variant QTL prior (default 1e-4) # --p2 Per-variant GWAS prior (default 1e-4) # --p12 Per-variant shared prior (default 5e-6) -# --p12-max Cap on enrichment-adjusted p12 (default 1e-3) # --no-adjust-pips Disable adjustPips (default: PIPs renormalized to overlap) # --finemapping-methods Comma-separated list for inline GWAS FM when --gwas-input is GwasSumStats (default 'susie') # --output Output RDS path (ColocResult collection) @@ -56,9 +55,6 @@ parser <- add_argument(parser, "--p2", parser <- add_argument(parser, "--p12", help = "Per-variant shared prior", type = "numeric", default = 5e-6) -parser <- add_argument(parser, "--p12-max", - help = "Cap on enrichment-adjusted p12", - type = "numeric", default = 1e-3) parser <- add_argument(parser, "--no-adjust-pips", help = "Disable adjustPips (use FMRs as supplied)", flag = TRUE) @@ -90,19 +86,37 @@ enrich <- if (nzchar(argv$enrichment) && argv$enrichment != "." && methods <- trimws(strsplit(argv$finemapping_methods, ",", fixed = TRUE)[[1L]]) -res <- colocPipeline( - qtlFineMappingResult = qtlFmr, - gwasInput = gwasIn, - filterLbfCs = argv$filter_lbf_cs, - filterLbfCsSecondary = if (is.na(argv$filter_lbf_cs_secondary)) NULL - else argv$filter_lbf_cs_secondary, - p1 = argv$p1, - p2 = argv$p2, - p12 = argv$p12, - finemappingMethods = methods, - enrichment = enrich, - p12Max = argv$p12_max, - adjustPips = !argv$no_adjust_pips) +# pecotmr 0.9.0 groups these into ColocPriorParam / ColocLbfFilterParam / +# GwasFineMappingParam. `priors` is passed only WITHOUT an enrichment table: +# 0.9.0 refuses both at once (with an enrichment table the priors are enloc's, +# derived from the genome-wide QTL results, so supplying them here too would +# mean scoring with one set and enriching with another), and it detects +# "supplied" as missing(priors) -- so even default values must be withheld. +# Priors the caller actually chose are not silently dropped: that combination +# is rejected here, with pecotmr's own reasoning. +prior_defaults <- c(p1 = 1e-4, p2 = 1e-4, p12 = 5e-6) +prior_values <- c(p1 = argv$p1, p2 = argv$p2, p12 = argv$p12) +priors_set <- !isTRUE(all.equal(prior_values, prior_defaults)) +if (!is.null(enrich) && priors_set) + stop("--p1/--p2/--p12 cannot be combined with --enrichment: with an ", + "enrichment table the colocalization priors are enloc's, derived from ", + "the genome-wide QTL results by qtl_enrichment.R. Drop the prior flags, ", + "or drop --enrichment to use your own.") + +coloc_args <- list( + qtlFineMappingResult = qtlFmr, + gwasInput = gwasIn, + lbfFilterArgs = ColocLbfFilterParam( + filterLbfCs = argv$filter_lbf_cs, + secondary = if (is.na(argv$filter_lbf_cs_secondary)) NULL + else argv$filter_lbf_cs_secondary), + gwasFineMappingArgs = GwasFineMappingParam(methods = methods), + enrichment = enrich, + adjustPips = !argv$no_adjust_pips) +if (is.null(enrich)) + coloc_args$priors <- ColocPriorParam(p1 = argv$p1, p2 = argv$p2, p12 = argv$p12) + +res <- do.call(colocPipeline, coloc_args) dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) saveRDS(res, argv$output) diff --git a/code/script/pecotmr_integration/colocboost.R b/code/script/pecotmr_integration/colocboost.R index aac283e99..b5562c520 100644 --- a/code/script/pecotmr_integration/colocboost.R +++ b/code/script/pecotmr_integration/colocboost.R @@ -32,7 +32,7 @@ # --separate-gwas Flag: run per-GWAS-study separate-focal coloc # --pip-cutoff-to-skip # Optional per-context single-effect skip cutoff passed to -# colocboostPipeline(pipCutoffToSkip = ). Either a single +# colocboostPipeline's SER pre-screen. Either a single # number applied to every context, or comma-separated # context=value pairs (e.g. "context1=0.1,context2=0"). # A negative value selects the data-driven 3/ncol(X) @@ -80,7 +80,10 @@ parser <- add_argument(parser, "--output", argv <- parse_args(parser) # Parse --pip-cutoff-to-skip into the scalar / context-named-vector shape -# colocboostPipeline(pipCutoffToSkip = ) accepts. Mirrors susie_twas.R. +# the SER pre-screen accepts. Mirrors susie_twas.R. (pecotmr 0.9.0 takes it as +# SignalScreenParam(pip = ); a field set to 0 there is off, which is this +# flag's documented 0 = off, and a negative value still selects the adaptive +# 3/nVariants threshold.) parse_pip_cutoff <- function(s) { if (is.null(s) || !nzchar(s)) return(0) pairs <- strsplit(s, ",")[[1L]] @@ -129,7 +132,7 @@ res <- if (has_region) { xqtlColoc = xqtl_coloc, jointGwas = joint_gwas, separateGwas = separate_gwas, - pipCutoffToSkip = pip_cutoff_to_skip) + signalScreenArgs = SignalScreenParam(pip = pip_cutoff_to_skip)) } else { colocboostPipeline( qd, @@ -140,7 +143,7 @@ res <- if (has_region) { xqtlColoc = xqtl_coloc, jointGwas = joint_gwas, separateGwas = separate_gwas, - pipCutoffToSkip = pip_cutoff_to_skip) + signalScreenArgs = SignalScreenParam(pip = pip_cutoff_to_skip)) } # Production output is emitted verbatim — including colocboostPipeline()'s wall-clock diff --git a/code/script/pecotmr_integration/ctwas_assemble.R b/code/script/pecotmr_integration/ctwas_assemble.R index 709408e39..fe54ace64 100644 --- a/code/script/pecotmr_integration/ctwas_assemble.R +++ b/code/script/pecotmr_integration/ctwas_assemble.R @@ -119,16 +119,20 @@ fmrPaths <- split_paths(argv$fine_mapping_results) fmr <- if (length(fmrPaths) > 0L) combineFineMappingResults(lapply(fmrPaths, readRDS)) else NULL +# pecotmr 0.9.0 grouped the four per-gene variant-pruning knobs into one +# bundle, with two renamed: minPipCutoff -> pipRescue (same meaning -- variants +# above it are must-keep) and maxNumVariants -> maxVariants. inputs <- assembleCtwasInputs( gwasSumStats = gwasSumStats, twasWeights = weightObjs, twasZ = tz, fineMappingResult = fmr, method = if (nzchar(argv$method)) argv$method else NULL, - twasWeightCutoff = argv$twas_weight_cutoff, - csMinCor = argv$cs_min_cor, - minPipCutoff = argv$min_pip_cutoff, - maxNumVariants = argv$max_num_variants) + variantPruningArgs = VariantPruningParam( + weightCutoff = argv$twas_weight_cutoff, + csMinCor = argv$cs_min_cor, + pipRescue = argv$min_pip_cutoff, + maxVariants = argv$max_num_variants)) dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) saveRDS(inputs, argv$output) diff --git a/code/script/pecotmr_integration/ctwas_est.R b/code/script/pecotmr_integration/ctwas_est.R index 2f21218e8..fdfa38562 100644 --- a/code/script/pecotmr_integration/ctwas_est.R +++ b/code/script/pecotmr_integration/ctwas_est.R @@ -3,7 +3,7 @@ # # cTWAS step 2: estimate group prior + prior variance. # Reads the assembled inputs from ctwas_assemble.R and calls -# pecotmr::estCtwasParam() (which wraps ctwas::assemble_region_data + +# pecotmr::estCtwasGroupPriors() (which wraps ctwas::assemble_region_data + # ctwas::est_param). Saves the augmented state — `inputs` + `region_data` # + `boundary_genes` + `z_gene` + `param` — for ctwas_finemap.R to # consume downstream. @@ -58,20 +58,34 @@ parser <- add_argument(parser, "--output", argv <- parse_args(parser) inputs <- readRDS(argv$inputs) -est <- estCtwasParam( +# pecotmr 0.9.0: estCtwasParam() -> estCtwasGroupPriors(), with the EM prior +# knobs in CtwasPriorParam() (groupPriorVarStructure -> varStructure), ncore -> +# numThreads, and the two ctwas::est_param pass-throughs in the est-param +# options bundle under ctwas's own names. +est <- estCtwasGroupPriors( inputs, - thin = argv$thin, - niterPrefit = argv$niter_prefit, - niter = argv$niter, - groupPriorVarStructure = argv$group_prior_var_structure, - ncore = argv$ncore, - fallbackToPrefit = argv$fallback_to_prefit, - min_group_size = argv$min_group_size, - min_p_single_effect = argv$min_p_single_effect) + ctwasPriorArgs = CtwasPriorParam( + thin = argv$thin, + niterPrefit = argv$niter_prefit, + niter = argv$niter, + varStructure = argv$group_prior_var_structure, + fallbackToPrefit = argv$fallback_to_prefit), + numThreads = argv$ncore, + # min_group_size is a formal of BOTH ctwas::assemble_region_data and + # ctwas::est_param. Pre-0.9.0 a single `...` was splatted into every step, so + # one value reached both; the per-step bundles must each carry it or + # assemble_region_data silently keeps ctwas's own default (100) and drops + # every group in a small run. min_p_single_effect is est_param's alone. + methodArgs = CtwasOptions( + regionData = CtwasRegionDataOptions( + min_group_size = argv$min_group_size), + estParam = CtwasEstParamOptions( + min_group_size = argv$min_group_size, + min_p_single_effect = argv$min_p_single_effect))) dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) saveRDS(est, argv$output) -cat(sprintf("Wrote estCtwasParam result to %s\n", argv$output)) +cat(sprintf("Wrote estCtwasGroupPriors result to %s\n", argv$output)) cat("group_prior: ", paste(format(est$param$group_prior, digits = 4), collapse = " | "), "\n") cat("group_prior_var:", diff --git a/code/script/pecotmr_integration/ctwas_finemap.R b/code/script/pecotmr_integration/ctwas_finemap.R index 74390bbda..c8078a7e5 100644 --- a/code/script/pecotmr_integration/ctwas_finemap.R +++ b/code/script/pecotmr_integration/ctwas_finemap.R @@ -15,7 +15,6 @@ # screen/finemap. Mirrors the legacy ctwas_3 # escape hatch for NaN-on-iter-2 EM divergence. # --L Pass-through (default 5) -# --no-filter-L Flag: disable ctwas's internal L >= 1 region screen # --min-nonsnp-pip Pass-through (default 0.5) # --merge-regions Flag: after fine-mapping, merge boundary genes' # adjacent LD blocks and re-fine-map the merged @@ -47,9 +46,6 @@ parser <- add_argument(parser, "--param-override", parser <- add_argument(parser, "--L", help = "Max number of credible sets per region", type = "integer", default = 5L) -parser <- add_argument(parser, "--no-filter-L", - help = "Disable ctwas's internal L >= 1 region screen (toy data)", - flag = TRUE) parser <- add_argument(parser, "--min-nonsnp-pip", help = "min_nonSNP_PIP threshold for screen_regions", type = "numeric", default = 0.5) @@ -89,27 +85,38 @@ if (nzchar(argv$param_override) && argv$param_override != "." && cat("Applied param override from ", argv$param_override, "\n", sep = "") } +# pecotmr 0.9.0: ncore -> numThreads everywhere, and L splits by step. The +# screening step has no L formal of its own -- its three knobs are +# ctwas::screen_regions() arguments, which travel under ctwas's own names in +# CtwasScreenOptions() -- while the fine-mapping and merge steps take L as the +# named `maxNumSingleEffects` formal. screened <- screenCtwasRegions( est, - L = argv$L, - ncore = argv$ncore, - filter_L = !argv$no_filter_L, - min_nonSNP_PIP = argv$min_nonsnp_pip) + numThreads = argv$ncore, + # Only min_nonSNP_PIP: the installed ctwas's screen_regions() has neither an + # `L` nor a `filter_L` argument. Pre-0.9.0 both were accepted here and then + # silently discarded (pecotmr filtered pass-throughs to the callee's + # formals); 0.9.0's options bundle validates against those formals instead, + # which is what surfaced it. L still matters for the fine-map / merge steps + # below, where it is a real argument. + methodArgs = CtwasOptions( + screen = CtwasScreenOptions( + min_nonSNP_PIP = argv$min_nonsnp_pip))) final <- finemapCtwasRegions( screened, - L = argv$L, - ncore = argv$ncore) + maxNumSingleEffects = argv$L, + numThreads = argv$ncore) # Optional step 4: boundary-gene region merging + re-fine-map. if (argv$merge_regions) { final <- mergeCtwasBoundaryRegions( final, - pipThresh = argv$merge_pip_cutoff, - filterCs = argv$merge_filter_cs, - maxSNP = argv$max_snp, - L = argv$L, - ncore = argv$ncore) + pipThresh = argv$merge_pip_cutoff, + filterCs = argv$merge_filter_cs, + maxSNP = argv$max_snp, + maxNumSingleEffects = argv$L, + numThreads = argv$ncore) cat("Applied boundary-region merging (merge_regions).\n") } @@ -121,11 +128,11 @@ result <- asCtwasResult(final, keepSnps = argv$keep_snps) dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) saveRDS(result, argv$output) cat(sprintf("Wrote CtwasResult (%d row(s)) to %s\n", nrow(result), argv$output)) -fm <- getFinemap(result) +fm <- ctwasPosteriors(result) if (!is.null(fm) && nrow(fm) > 0L) { g <- fm[fm$type != "SNP", , drop = FALSE] cat(sprintf(" fine-mapped rows: %d (%d gene-level, %d SNP-level)\n", nrow(fm), nrow(g), nrow(fm) - nrow(g))) } else { - cat(" no fine-mapped genes (no regions surviving filter_L >= 1)\n") + cat(" no fine-mapped genes (no regions survived screening)\n") } diff --git a/code/script/pecotmr_integration/fine_mapping.R b/code/script/pecotmr_integration/fine_mapping.R index d1ef895e3..58d50c884 100644 --- a/code/script/pecotmr_integration/fine_mapping.R +++ b/code/script/pecotmr_integration/fine_mapping.R @@ -89,23 +89,23 @@ parser <- add_argument(parser, "--coverage", help = "SuSiE primary credible-set coverage", type = "numeric", default = 0.95) parser <- add_argument(parser, "--secondary-coverage", - help = "Comma-separated secondary credible-set coverages (fineMappingPipeline secondaryCoverage)", + help = "Comma-separated secondary credible-set coverages (CredibleSetParam secondaryCoverage)", type = "character", default = "0.7,0.5") parser <- add_argument(parser, "--min-abs-corr", - help = "Credible-set purity threshold (fineMappingPipeline minAbsCorr)", + help = "Credible-set purity threshold (CredibleSetParam minAbsCorr)", type = "numeric", default = 0.8) parser <- add_argument(parser, "--median-abs-corr", - help = "Optional median-abs-corr purity, OR-logic with --min-abs-corr (fineMappingPipeline medianAbsCorr; requires pecotmr step-1). Omit to leave off.", + help = "Optional median-abs-corr purity, OR-logic with --min-abs-corr (CredibleSetParam medianAbsCorr). Omit to leave off.", type = "numeric", default = NA) parser <- add_argument(parser, "--pip-cutoff", - help = "PIP signal cutoff (fineMappingPipeline signalCutoff)", + help = "PIP signal cutoff (CredibleSetParam signalCutoff)", type = "numeric", default = 0.025) parser <- add_argument(parser, "--L", - help = "SuSiE number of single effects (susie L); pipeline default 10", + help = "SuSiE number of single effects (CredibleSetParam maxNumSingleEffects); pipeline default 10", type = "integer", default = 10L) parser <- add_argument(parser, "--L-greedy", help = paste("SuSiE greedy init count (susie L_greedy): a positive integer to run", - "the greedy-L loop, or 'none'/'off' to disable it; pipeline default", + "the greedy-L loop (CredibleSetParam maxNumSingleEffectsGreedy), or 'none'/'off' to disable it; pipeline default", "none (greedy off)"), type = "character", default = "none") parser <- add_argument(parser, "--method-args", @@ -115,13 +115,13 @@ parser <- add_argument(parser, "--seed", help = "Integer RNG seed set before fitting (reproducibility); unset = no seeding", type = "integer", default = NA) parser <- add_argument(parser, "--cv-folds", - help = "Cross-validation folds for the mash-prior / predictive refits (fineMappingPipeline cvFolds); 0 = off", + help = "Cross-validation folds for the mash-prior / predictive refits (CrossValidationParam folds); 0 = off", type = "integer", default = 0L) parser <- add_argument(parser, "--cv-threads", - help = "Parallel workers for the cross-validation fold refits (fineMappingPipeline cvThreads); -1 = all cores", + help = "Parallel workers for the cross-validation fold refits (CrossValidationParam numThreads); -1 = all cores", type = "integer", default = 1L) parser <- add_argument(parser, "--pip-cutoff-to-skip", - help = "Single-effect (SER) pre-screen cutoff (fineMappingPipeline pipCutoffToSkip), QTL mode; 0 = off, <0 = adaptive 3/nVariants", + help = "Single-effect (SER) pre-screen cutoff (SignalScreenParam pip), QTL mode; 0 = off, <0 = adaptive 3/nVariants", type = "numeric", default = 0) # --- GWAS-mode SuSiE-RSS knobs (fineMappingPipeline, GwasSumStats path). SER # fallback is ON by default (reproduces the single-panel notebook): fit a @@ -130,19 +130,19 @@ parser <- add_argument(parser, "--pip-cutoff-to-skip", # mode only); rFinite forwards only when set. susie_rss_control() settings # (susieR >= 0.16.6) are passed as a JSON object via --rss-control. parser <- add_argument(parser, "--ser-fallback", - help = "GWAS mode: TRUE/FALSE. Fall back to the single-effect (SER) result when susieR flags R unreliable (fineMappingPipeline serFallback). Default TRUE.", + help = "GWAS mode: TRUE/FALSE. Fall back to the single-effect (SER) result when susieR flags R unreliable (SusieRssParam serFallback). Default TRUE.", type = "character", default = "TRUE") parser <- add_argument(parser, "--r-finite", - help = "GWAS mode: finite-sample R correction size (fineMappingPipeline rFinite). Empty = auto (LD-panel sample size).", + help = "GWAS mode: finite-sample R correction size (SusieRssParam rFinite). Empty = auto (LD-panel sample size).", type = "character", default = "") parser <- add_argument(parser, "--r-mismatch", - help = "GWAS mode: LD-mismatch correction (fineMappingPipeline rMismatch): 'none', 'eb', or 'eb_mix' (residual-mixture EB, susieR >= 0.16.6). Default 'eb_mix'.", + help = "GWAS mode: LD-mismatch correction (SusieRssParam rMismatch): 'none', 'eb', or 'eb_mix' (residual-mixture EB, susieR >= 0.16.6). Default 'eb_mix'.", type = "character", default = "eb_mix") parser <- add_argument(parser, "--rss-control", - help = "GWAS mode: JSON object of susieR::susie_rss_control() settings (e.g. '{\"check_prior\":true,\"mismatch_estimator\":\"map\"}'), forwarded as fineMappingPipeline rssControl. Empty = susie_rss_control() defaults.", + help = "GWAS mode: JSON object of susieR::susie_rss_control() settings (e.g. '{\"check_prior\":true,\"mismatch_estimator\":\"map\"}'), forwarded as SusieRssParam control. Empty = susie_rss_control() defaults.", type = "character", default = "") parser <- add_argument(parser, "--keep-original-fit", - help = "GWAS mode: TRUE/FALSE. Retain the pre-fallback multi-effect fit (fineMappingPipeline keepFullFit) on SER-fallback regions, where it genuinely differs from the reported SER result. Default FALSE (keep none).", + help = "GWAS mode: TRUE/FALSE. Retain the pre-fallback multi-effect fit (SusieRssParam keepFullFit) on SER-fallback regions, where it genuinely differs from the reported SER result. Default FALSE (keep none).", type = "character", default = "FALSE") # --- Multivariate / joint-fit knobs (QTL mode; mvsusie / fsusie). Each is # opt-in and omitted from the pipeline call when left at its default, so this @@ -151,7 +151,7 @@ parser <- add_argument(parser, "--joint-specification", help = "Comma-separated joint-fit axes (context,trait,study) for mvsusie/fsusie; empty = implicit per-(context,trait) branch", type = "character", default = "") parser <- add_argument(parser, "--twas-weights", - help = "Optional TwasWeights RDS (preceding mr.mash run) supplying the mvSuSiE data-driven prior; omit for the canonical prior", + help = "Optional TwasWeights RDS (preceding mr.mash run) supplying the mvSuSiE data-driven prior (fineMappingPipeline mrmashPrior); omit for the canonical prior", type = "character", default = "") parser <- add_argument(parser, "--fine-mapping-result", help = "Optional existing FineMappingResult RDS used as a resume cache (fineMappingPipeline fineMappingResult); (study, context, trait, method) tuples already present are not refit", @@ -202,15 +202,20 @@ parser <- add_argument(parser, "--output", help = "Output RDS path", type = "character") argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + # Read a whitespace-delimited ID file into a unique character vector; NULL when # the path is empty / "." / missing. read_ids <- function(p) if (nzchar(p) && p != "." && file.exists(p)) unique(trimws(unlist(strsplit(readLines(p), "[[:space:]]+")))) else NULL -# Build the per-call genotype-filter overrides as pipeline ARGUMENTS (NULL = -# leave the QtlDataset's construct-time slot untouched). The pipeline applies -# these to a validated copy; the wrapper no longer mutates @slots directly. -# Mirrored in twas_weights.R. +# Build the per-call genotype-filter overrides as one GenotypeFilterParam() +# bundle (pecotmr 0.9.0 collapsed the seven loose arguments into it). An UNSET +# field leaves the QtlDataset's construct-time slot untouched, which is exactly +# the override role the bundle documents, so the NULLs are dropped rather than +# passed. The pipeline applies these to a validated copy; the wrapper never +# mutates @slots directly. Mirrored in twas_weights.R. qd_filter_overrides <- function(argv) { ov <- list( mafCutoff = if (!is.na(argv$maf_cutoff)) argv$maf_cutoff else NULL, @@ -220,7 +225,8 @@ qd_filter_overrides <- function(argv) { keepIndel = if (isTRUE(argv$drop_indel)) FALSE else NULL, keepSamples = read_ids(argv$keep_samples), keepVariants = read_ids(argv$keep_variants)) - ov[!vapply(ov, is.null, logical(1))] + ov <- ov[!vapply(ov, is.null, logical(1))] + list(genotypeFilterParam = do.call(GenotypeFilterParam, ov)) } # Parse --method-args into a nested named list of per-method kwargs. @@ -255,9 +261,8 @@ median_abs_corr <- if (length(argv$median_abs_corr) != 1L || is.na(argv$median_a # L_greedy: 'none'/'off'/'null'/empty (case-insensitive) -> NULL, i.e. greedy-L # off (susie's own L_greedy default is NULL). Otherwise a positive integer. -# The NULL is kept via the inline list() in cs_args below, which preserves a -# named NULL through c()/do.call to fineMappingPipeline; only an x$Lgreedy <- NULL -# assignment would drop it (re-enabling pecotmr's default 5). +# NULL is also CredibleSetParam(maxNumSingleEffectsGreedy = )'s own default, so +# dropping it from the bundle below leaves greedy-L off either way. l_greedy_raw <- trimws(argv[["L_greedy"]]) l_greedy_val <- if (tolower(l_greedy_raw) %in% c("none", "off", "null", "")) { NULL @@ -291,7 +296,11 @@ if (!has_qtl && !has_gwas && !has_qss) contexts_arg <- if (nzchar(argv$contexts) && argv$contexts != ".") trimws(strsplit(argv$contexts, ",", fixed = TRUE)[[1L]]) else NULL -methods <- trimws(strsplit(argv$methods, ",", fixed = TRUE)[[1L]]) +# Snake_case spellings stay the public interface on this side (they are the +# notebook's documented --methods values); pecotmr took its tokens camelCase, +# so translate. Without this, --methods susie_inf / susie_ash reached +# fineMappingPipeline verbatim and died on "unknown method token". +methods <- split_pecotmr_names(argv$methods, FM_METHOD_ALIASES, split = ",") # Joint-fit axes for mvsusie/fsusie: "" -> NULL (implicit per-tuple branch). joint_spec <- if (nzchar(argv$joint_specification) && argv$joint_specification != ".") @@ -307,9 +316,13 @@ fmr_obj <- if (nzchar(argv$fine_mapping_result) && argv$fine_mapping_result != " # Build the `methods` argument: the bare tokens, or the named-list {token: kwargs} # form when --method-args is given. SuSiE L / L_greedy defaults + susie-family # routing live in pecotmr (.fmNormalizeMethods); the wrapper only forwards the -# --L / --L-greedy values as top-level args (see cs_args). --method-args keys -# must be among the --methods tokens (no silent typos). +# --L / --L-greedy values through the credible-set bundle (see cs_param_args). +# --method-args keys must be among the --methods tokens (no silent typos). if (!is.null(parsed_method_args)) { + # Keyed by method token, so normalise the keys the same way or the subset + # check below would reject an alias the user spelled consistently. + names(parsed_method_args) <- apply_pecotmr_aliases(names(parsed_method_args), + FM_METHOD_ALIASES) unknown <- setdiff(names(parsed_method_args), methods) if (length(unknown) > 0L) stop("--method-args has keys not listed in --methods (got '", @@ -321,26 +334,39 @@ methods_arg <- if (is.null(parsed_method_args)) methods else if (tk %in% names(parsed_method_args)) parsed_method_args[[tk]] else list()), methods) -# Credible-set / coverage knobs common to both modes. medianAbsCorr is added -# only when set (NULL -> omitted), so the call also works against a pecotmr -# that predates that argument. L / L_greedy are forwarded top-level; pecotmr -# seeds the susie-family tokens and applies explicit --method-args overrides. -cs_args <- list(methods = methods_arg, - L = argv$L, - Lgreedy = l_greedy_val, - coverage = argv$coverage, - secondaryCoverage = secondary_cov, - signalCutoff = argv$pip_cutoff, - minAbsCorr = argv$min_abs_corr, - # Per-CS variant-level export: default emits within_cs_pip only; - # --full-fit widens per-CS matrices (alpha; +logbf/effect/effect-var - # with --full-fit-all; +filtered CS with --include-all-cs). - fullFit = isTRUE(argv[["full_fit"]]), - fullFitAlphaOnly = !isTRUE(argv[["full_fit_all"]]), - includeAllCs = isTRUE(argv[["include_all_cs"]])) -if (!is.null(median_abs_corr)) cs_args$medianAbsCorr <- median_abs_corr +# Per-CS variant-level export. pecotmr 0.9.0 replaced the fullFit / +# fullFitAlphaOnly pair with one perCsColumns choice: "none" emits only the +# within_cs_pip scalar (the default), "alpha" also widens the per-CS alpha +# matrices, "full" widens logbf / effect / effect-variance as well. +per_cs_columns <- if (!isTRUE(argv[["full_fit"]])) { + "none" +} else if (isTRUE(argv[["full_fit_all"]])) { + "full" +} else { + "alpha" +} + +# Credible-set / coverage knobs common to both modes, as the single +# CredibleSetParam() bundle 0.9.0 expects (L -> maxNumSingleEffects, +# L_greedy -> maxNumSingleEffectsGreedy). medianAbsCorr and the greedy count +# are added only when set, so an unset field keeps the bundle's own default. +# The susie-family token seeding and the explicit --method-args overrides still +# live in pecotmr (.fmNormalizeMethods). +cs_param_args <- list(coverage = argv$coverage, + secondaryCoverage = secondary_cov, + signalCutoff = argv$pip_cutoff, + minAbsCorr = argv$min_abs_corr, + perCsColumns = per_cs_columns, + includeAllCs = isTRUE(argv[["include_all_cs"]]), + maxNumSingleEffects = argv$L) +if (!is.null(l_greedy_val)) + cs_param_args$maxNumSingleEffectsGreedy <- l_greedy_val +if (!is.null(median_abs_corr)) cs_param_args$medianAbsCorr <- median_abs_corr + +cs_args <- list(methods = methods_arg, + credibleSetParam = do.call(CredibleSetParam, cs_param_args)) # Resume cache applies to both QTL and GWAS modes (both pipeline methods accept -# fineMappingResult). Added only when supplied, for pecotmr-version tolerance. +# fineMappingResult). Added only when supplied. if (!is.null(fmr_obj)) cs_args$fineMappingResult <- fmr_obj if (has_gwas || has_qss) { @@ -359,14 +385,13 @@ if (has_gwas || has_qss) { # (keepFullFit = "fallback"); FALSE (default) drops it (keepFullFit = "none"). # pecotmr's "all" mode is intentionally not exposed. keep_full_fit <- if (isTRUE(keep_original_fit)) "fallback" else "none" - # GWAS-only SuSiE-RSS knobs on top of the shared cs_args. rFinite forwards - # only when set; --rss-control (JSON) becomes the rssControl named list of - # susie_rss_control() settings. - gwas_args <- c(cs_args, - list(serFallback = ser_fallback, - rMismatch = argv$r_mismatch, - keepFullFit = keep_full_fit)) - if (nzchar(argv$r_finite)) gwas_args$rFinite <- as.numeric(argv$r_finite) + # SuSiE-RSS knobs, as the SusieRssParam() bundle 0.9.0 expects (rssControl -> + # control). rFinite and control forward only when set, so an unset field keeps + # the bundle's own default. + rss_param_args <- list(serFallback = ser_fallback, + rMismatch = argv$r_mismatch, + keepFullFit = keep_full_fit) + if (nzchar(argv$r_finite)) rss_param_args$rFinite <- as.numeric(argv$r_finite) if (nzchar(argv$rss_control) && argv$rss_control != "." && argv$rss_control != "{}") { rc <- tryCatch(jsonlite::fromJSON(argv$rss_control, simplifyVector = TRUE), @@ -376,9 +401,12 @@ if (has_gwas || has_qss) { if (!is.list(rc) || is.null(names(rc)) || any(!nzchar(names(rc)))) stop("--rss-control must be a JSON object with named fields, e.g. ", '\'{"check_prior":true,"mismatch_estimator":"map"}\'.') - gwas_args$rssControl <- rc + rss_param_args$control <- rc } - if (!is.null(seed_val)) gwas_args$seed <- seed_val + # No seed= here: neither RSS method takes one (pre-0.9.0 they ended in `...`, + # which silently swallowed it), so RSS mode is seeded by the set.seed() above. + gwas_args <- c(cs_args, + list(susieRssParam = do.call(SusieRssParam, rss_param_args))) res <- do.call(fineMappingPipeline, c(list(gss), gwas_args)) label <- if (has_gwas) paste0("GwasSumStats '", basename(argv$gwas_sumstats), "'") else paste0("QtlSumStats '", basename(argv$qtl_sumstats), "'") @@ -397,17 +425,25 @@ if (has_gwas || has_qss) { # cisWindow is a gene-mode (traitId) knob — it expands each trait's own # coordinates. In region mode the literal variant window comes from --region, # and passing both is an error, so cisWindow rides on the traitId branch only. + # pipCutoffToSkip is now the `pip` field of the SER pre-screen bundle (a + # field set to 0 is off there, which is this flag's documented "0 = off", and + # a negative value still selects the adaptive 3/nVariants threshold); + # cvFolds / cvThreads are now the cross-validation bundle's folds / + # numThreads. qtl_args <- c(list(qd), cs_args, qd_filter_overrides(argv), - list(contexts = contexts_arg, - pipCutoffToSkip = argv$pip_cutoff_to_skip, - cvFolds = argv$cv_folds, - cvThreads = argv$cv_threads)) + list(contexts = contexts_arg, + signalScreenArgs = SignalScreenParam( + pip = argv$pip_cutoff_to_skip), + crossValidationArgs = CrossValidationParam( + folds = argv$cv_folds, numThreads = argv$cv_threads))) # Opt-in multivariate / joint knobs. Each is added only when the user set it, # so the call stays compatible with a pecotmr that predates the argument # (mirrors the medianAbsCorr handling above). if (!is.null(joint_spec)) qtl_args$jointSpecification <- joint_spec if (!is.null(twas_weights_obj)) { - qtl_args$twasWeights <- twas_weights_obj + # 0.9.0 renamed this argument mrmashPrior; it still takes the TwasWeights + # object a preceding mr.mash run produced. + qtl_args$mrmashPrior <- twas_weights_obj qtl_args$dataDrivenPriorWeightsCutoff <- argv$data_driven_prior_weights_cutoff } if (isTRUE(argv$use_pca)) { @@ -461,11 +497,16 @@ if (has_gwas && !is.null(res)) { study_i <- as.character(res$study[i]) region_i <- as.character(res$blockId[i]) method_i <- as.character(res$method[i]) - tl <- tryCatch(as.data.frame(getTopLoci(row_i)), error = function(e) NULL) - ncs <- if (!is.null(tl) && "cs_95" %in% names(tl)) + # No per-row guards: the block-level one below already keeps a + # diagnostics failure from touching the run, and a second layer here + # would turn a real accessor break into silent "NA" columns. + # topLoci() is a GRanges from 0.9.0 on; as.data.frame() still exposes + # cs_95 (it becomes a column alongside the range columns). + tl <- as.data.frame(topLoci(row_i)) + ncs <- if ("cs_95" %in% names(tl)) length(setdiff(unique(tl$cs_95), 0)) else NA - fit <- tryCatch(getSusieFit(row_i), error = function(e) NULL) - ser <- if (!is.null(fit) && !is.null(fit$serFallbackUsed)) + fit <- susieFit(row_i) + ser <- if (!is.null(fit$serFallbackUsed)) isTRUE(fit$serFallbackUsed) else NA cat(sprintf("[FM %s %s] method=%s CS(cs_95)=%s SER-fallback=%s\n", study_i, region_i, method_i, @@ -473,5 +514,9 @@ if (has_gwas && !is.null(res)) { if (isTRUE(ser)) "YES" else if (identical(ser, FALSE)) "no" else "NA")) } }, error = function(e) - message("fine_mapping.R: GWAS summary skipped (", conditionMessage(e), ")"))) + # warning(), not message(): this runs AFTER saveRDS, so it must never fail + # the run -- but a skipped summary means an accessor broke, which should be + # loud enough to see in a CI log rather than reading as ordinary output. + warning("fine_mapping.R: GWAS summary skipped (", conditionMessage(e), ")", + call. = FALSE))) } diff --git a/code/script/pecotmr_integration/fine_mapping_cis_db_export.R b/code/script/pecotmr_integration/fine_mapping_cis_db_export.R index b7a45f96c..305d4c92f 100644 --- a/code/script/pecotmr_integration/fine_mapping_cis_db_export.R +++ b/code/script/pecotmr_integration/fine_mapping_cis_db_export.R @@ -14,10 +14,10 @@ # Meta schema (one tab-separated row), exactly as SuSiE_enloc consumes it: # chr start end region_id TSS original_data combined_data # combined_data_sumstats conditions conditions_top_loci -# TSS = start(getTraitPosition(db)) (NA when the trait position was never +# TSS = start(traitPosition(db)) (NA when the trait position was never # supplied, e.g. a QtlSumStats-derived FMR); conditions = getContexts; # conditions_top_loci = contexts with >=1 top locus; chr/start/end default to -# getRegion(db) when not supplied via --chr/--start/--end. +# genomicRegion(db) when not supplied via --chr/--start/--end. suppressPackageStartupMessages({ library(argparser) @@ -46,6 +46,9 @@ parser <- add_argument(parser, "--pip-sum-output", help = "Optional pip_sum TSV type = "character", default = NA) argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + inputs <- as.character(argv$input) if (length(inputs) == 0L) stop("--input requires at least one FineMappingResult RDS path.") minPurity <- if (is.na(argv$min_purity)) NULL else argv$min_purity @@ -67,27 +70,28 @@ saveRDS(db, argv$combined_data_output) # ---- 2. Marginal sumstats db (optional) ------------------------------------ if (!is.na(argv$combined_data_sumstats_output)) { - me <- tryCatch(as.data.frame(getMarginalEffects(db)), error = function(e) NULL) + me <- as.data.frame(marginalEffects(db)) dir.create(dirname(argv$combined_data_sumstats_output), showWarnings = FALSE, recursive = TRUE) saveRDS(me, argv$combined_data_sumstats_output) } # ---- 3. FMR-derived meta fields -------------------------------------------- -tp <- getTraitPosition(db) +tp <- traitPosition(db) TSS <- if (methods::is(tp, "GRanges") && length(tp) > 0L) GenomicRanges::start(tp)[1L] else NA -ctx <- getContexts(db) +ctx <- contexts(db) conditions <- if (is.null(ctx)) "" else paste(unique(as.character(ctx)), collapse = ",") -tl <- tryCatch(as.data.frame(getTopLoci(db, signalCutoff = argv$signal_cutoff, minPurity = minPurity)), - error = function(e) NULL) -conditions_top_loci <- if (!is.null(tl) && nrow(tl) > 0L && "context" %in% names(tl)) +# Unguarded: an empty top-loci set is a 0-row table, which the next line +# handles; an accessor error is a real break and must surface. +tl <- as_variant_df(topLoci(db, signalCutoff = argv$signal_cutoff, minPurity = minPurity)) +conditions_top_loci <- if (nrow(tl) > 0L && "context" %in% names(tl)) paste(unique(as.character(tl$context)), collapse = ",") else "" regionId <- if (!is.na(argv$region_id)) argv$region_id else { tr <- unique(as.character(db$trait)); if (length(tr) == 1L) tr else NA } -rg <- tryCatch(getRegion(db), error = function(e) NULL) +rg <- genomicRegion(db) gr1 <- function(fn, i) if (methods::is(rg, "GRanges") && length(rg) > 0L) fn(rg)[i] else NA chr <- if (!is.na(argv$chr)) argv$chr else as.character(gr1(GenomicRanges::seqnames, 1L)) start <- if (!is.na(argv$start)) argv$start else gr1(GenomicRanges::start, 1L) @@ -110,17 +114,27 @@ cat(sprintf("Wrote cis_results_db (%d FMR input(s) -> %d rows) + meta (region_id # ---- 4. pip_sum (optional) ------------------------------------------------- if (!is.na(argv$pip_sum_output)) { + # The identity columns live in mcols(): on these collection classes `[[` + # selects an ELEMENT, not a column, so db[["study"]] is NULL and every + # selector below used to come out NA -- which silently produced an EMPTY + # pip_sum (the error that NA selection raises was being swallowed). + idm <- S4Vectors::mcols(db) rows <- lapply(seq_len(nrow(db)), function(i) { # NULL-safe per-row identity read: a GwasFineMappingResult has no context / - # trait column, so db$context is NULL and as.character(db$context[i]) would be - # character(0) -- length 0 -- which makes the data.frame() below error on any - # region with signal. Substitute a length-1 NA when the column is absent. - .col1 <- function(nm) { v <- db[[nm]]; if (is.null(v)) NA_character_ else as.character(v[i]) } + # trait column, so that column is absent from mcols() entirely. Substitute + # a length-1 NA, and drop those fields from the selection -- selecting on + # NA matches no entry. The remaining ones come straight off this row, so + # together they always select it. + .col1 <- function(nm) { + v <- idm[[nm]] + if (is.null(v)) NA_character_ else as.character(v[i]) + } sel <- list(study = .col1("study"), context = .col1("context"), trait = .col1("trait"), method = .col1("method")) - pip <- tryCatch(do.call(getPip, c(list(db), sel)), error = function(e) NULL) - if (is.null(pip)) return(NULL) - data.frame(pip_sum = sum(pip[pip > 0], na.rm = TRUE), condition = sel$context, + context <- sel$context + sel <- sel[!vapply(sel, is.na, logical(1))] + pipVec <- do.call(pip, c(list(db), sel)) + data.frame(pip_sum = sum(pipVec[pipVec > 0], na.rm = TRUE), condition = context, stringsAsFactors = FALSE) }) pipSum <- do.call(rbind, Filter(Negate(is.null), rows)) diff --git a/code/script/pecotmr_integration/fine_mapping_export.R b/code/script/pecotmr_integration/fine_mapping_export.R index be74316cc..790d6ac48 100644 --- a/code/script/pecotmr_integration/fine_mapping_export.R +++ b/code/script/pecotmr_integration/fine_mapping_export.R @@ -38,6 +38,9 @@ parser <- add_argument(parser, "--output", help = "Output TSV path", type = "character") argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + view <- match.arg(argv$view, c("topLoci", "cs", "cs_summary", "pip", "marginals", "lbf", "credible_band", "affected_regions")) inputs <- as.character(argv$input) @@ -51,14 +54,14 @@ minPurity <- if (is.na(argv$min_purity)) NULL else argv$min_purity # wrapper only tags each row with its source RDS and concatenates across inputs. view_fn <- switch(view, topLoci = function(fmr) - as.data.frame(getTopLoci(fmr, signalCutoff = argv$signal_cutoff, - minPurity = minPurity)), - cs = function(fmr) as.data.frame(getCs(fmr, minPurity = minPurity)), + as_variant_df(topLoci(fmr, signalCutoff = argv$signal_cutoff, + minPurity = minPurity)), + cs = function(fmr) as.data.frame(credibleSets(fmr, minPurity = minPurity)), # cs_summary: one row per credible set (size / purity / V / logBF / lead). - cs_summary = function(fmr) as.data.frame(getCredibleSetSummary(fmr)), + cs_summary = function(fmr) as.data.frame(credibleSetSummary(fmr)), # lbf: wide variant x effect log Bayes factor matrix (lbf_L1..lbf_LL). - lbf = function(fmr) as.data.frame(getLbf(fmr)), - marginals = function(fmr) as.data.frame(getMarginalEffects(fmr)), + lbf = function(fmr) as.data.frame(lbf(fmr)), + marginals = function(fmr) as.data.frame(marginalEffects(fmr)), # fSuSiE functional views (require untrimmed fits; degrade to empty otherwise). # credible_band = fitted effect curve + band; affected_regions = GRanges of # the intervals where a CS's band excludes zero (coerced to a flat table). @@ -66,7 +69,7 @@ view_fn <- switch(view, affected_regions = function(fmr) as.data.frame(fsusieAffectedRegions(fmr)), # pip view: the (identity + variant_id + pip) projection of topLoci. pip = function(fmr) { - tl <- as.data.frame(getTopLoci(fmr, signalCutoff = argv$signal_cutoff)) + tl <- as_variant_df(topLoci(fmr, signalCutoff = argv$signal_cutoff)) tl[, intersect(c("study", "context", "trait", "blockId", "method", "variant_id", "pip"), names(tl)), drop = FALSE] }) @@ -88,12 +91,12 @@ for (path in inputs) { warning("Skipping non-FineMappingResult input: ", path) next } - df <- tryCatch(view_fn(fmr), - error = function(e) { - message(basename(path), ": ", conditionMessage(e)) - NULL - }) - if (is.null(df) || nrow(df) == 0L) next + # Unguarded: every view returns an EMPTY table rather than erroring when it + # does not apply (the fSuSiE views give a 0-row tibble / empty GRanges on a + # susie-only fit), and the next line skips those. An accessor error is a real + # break and must surface instead of silently dropping the input. + df <- view_fn(fmr) + if (nrow(df) == 0L) next df$source <- basename(path) df$block_id <- .blockIdFromName(basename(path)) pieces[[length(pieces) + 1L]] <- df diff --git a/code/script/pecotmr_integration/fine_mapping_overlap.R b/code/script/pecotmr_integration/fine_mapping_overlap.R index 3556279b1..decc9c2c2 100644 --- a/code/script/pecotmr_integration/fine_mapping_overlap.R +++ b/code/script/pecotmr_integration/fine_mapping_overlap.R @@ -44,6 +44,9 @@ parser <- add_argument(parser, "--meta-output", help = "Output path for the updated meta row", type = "character", default = NA) argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + # A QTL region may span several per-study FMRs (its meta original_data); combine # them so overlapTopLoci sees the region's full top-loci set. qtlPaths <- as.character(argv$qtl) @@ -63,13 +66,11 @@ for (gp in gwasPaths) { warning("Skipping non-GwasFineMappingResult input: ", gp) next } - ov <- tryCatch( - as.data.frame(overlapTopLoci(qtl, gwas, signalCutoff = argv$signal_cutoff)), - error = function(e) { - message(basename(gp), ": ", conditionMessage(e)) - NULL - }) - if (is.null(ov) || nrow(ov) == 0L) next + # Unguarded: "no shared variants" is a 0-row table, which the next line + # skips. The class check above is what rejects an unusable input; an error + # from the accessor itself is a real break and must surface. + ov <- as_variant_df(overlapTopLoci(qtl, gwas, signalCutoff = argv$signal_cutoff)) + if (nrow(ov) == 0L) next ov$gwas_source <- basename(gp) pieces[[length(pieces) + 1L]] <- ov blockTopLoci <- c(blockTopLoci, unique(as.character(gwas$blockId))) diff --git a/code/script/pecotmr_integration/fine_mapping_pip_plot.R b/code/script/pecotmr_integration/fine_mapping_pip_plot.R index 8c88336ed..e75901981 100644 --- a/code/script/pecotmr_integration/fine_mapping_pip_plot.R +++ b/code/script/pecotmr_integration/fine_mapping_pip_plot.R @@ -33,6 +33,9 @@ parser <- add_argument(parser, "--width", help = "Plot width in inches", type = parser <- add_argument(parser, "--height", help = "Per-PIP-panel height in inches", type = "numeric", default = 3) argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + fmr <- readRDS(argv$input) if (!methods::is(fmr, "FineMappingResultBase")) stop("--input must be a FineMappingResultBase subclass (got '", class(fmr)[[1L]], "').") @@ -53,10 +56,12 @@ writeStub <- function(msg) { if (nrow(fmr) == 0L) { writeStub("Empty FineMappingResult"); quit(save = "no") } -# getTopLoci(signalCutoff = 0): every fitted variant with pos, pip, CS membership +# topLoci(signalCutoff = 0): every fitted variant with pos, pip, CS membership # (cs_95 = "_"), and the row-identity columns -> one facet panel. -tl <- tryCatch(as.data.frame(getTopLoci(fmr, signalCutoff = 0)), error = function(e) NULL) -if (is.null(tl) || nrow(tl) == 0L) { writeStub("No usable PIP vectors"); quit(save = "no") } +# Unguarded on purpose: "nothing to plot" is a 0-row table, which the next line +# handles; an error here is a real break and must not become a blank stub. +tl <- as_variant_df(topLoci(fmr, signalCutoff = 0)) +if (nrow(tl) == 0L) { writeStub("No usable PIP vectors"); quit(save = "no") } isGwas <- methods::is(fmr, "GwasFineMappingResult") panel <- if (isGwas) { diff --git a/code/script/pecotmr_integration/fine_mapping_top_loci_bed.R b/code/script/pecotmr_integration/fine_mapping_top_loci_bed.R index 1cf2cd67e..58bf37e31 100644 --- a/code/script/pecotmr_integration/fine_mapping_top_loci_bed.R +++ b/code/script/pecotmr_integration/fine_mapping_top_loci_bed.R @@ -31,6 +31,9 @@ parser <- add_argument(parser, "--output", help = "Output BED path", type = "character") argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + inputs <- as.character(argv$input) if (length(inputs) == 0L) stop("--input requires at least one RDS path.") minPurity <- if (is.na(argv$min_purity)) NULL else argv$min_purity @@ -91,11 +94,13 @@ for (path in inputs) { if (!methods::is(fmr, "FineMappingResultBase")) { warning("Skipping non-FineMappingResult input: ", path); next } - tl <- tryCatch(as.data.frame(getTopLoci(fmr, signalCutoff = argv$signal_cutoff, - minPurity = minPurity)), - error = function(e) { message(basename(path), ": ", conditionMessage(e)); NULL }) - if (is.null(tl) || nrow(tl) == 0L) next - me <- as.data.frame(getMarginalEffects(fmr)) + # Unguarded: "no top loci past the cutoff" is a 0-row table, which the next + # line skips. The class check above rejects an unusable input; an accessor + # error is a real break and must surface rather than silently emptying the BED. + tl <- as_variant_df(topLoci(fmr, signalCutoff = argv$signal_cutoff, + minPurity = minPurity)) + if (nrow(tl) == 0L) next + me <- as.data.frame(marginalEffects(fmr)) piece <- bedFrom(tl, me) piece$block_id <- .blockIdFromName(basename(path)) pieces[[length(pieces) + 1L]] <- piece diff --git a/code/script/pecotmr_integration/fine_mapping_upset.R b/code/script/pecotmr_integration/fine_mapping_upset.R index c830d1d82..f73b75452 100644 --- a/code/script/pecotmr_integration/fine_mapping_upset.R +++ b/code/script/pecotmr_integration/fine_mapping_upset.R @@ -38,6 +38,9 @@ parser <- add_argument(parser, "--height", type = "numeric", default = 6) argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + if (!requireNamespace("UpSetR", quietly = TRUE)) stop("fine_mapping_upset.R requires the UpSetR package: ", "install.packages('UpSetR')") @@ -47,7 +50,7 @@ if (length(inputs) == 0L) stop("--input requires at least one RDS path.") # Build a named list of {set-label -> variant_ids} for UpSetR::fromList(). -# getCs(fmr) aggregates every entry's credible sets into one table already +# credibleSets(fmr) aggregates every entry's credible sets into one table already # tagged with the row identity (study/context/trait/blockId/method), so each # per-tuple set is a split() on the label -- no per-entry loop. sets <- list() @@ -57,8 +60,10 @@ for (path in inputs) { warning("Skipping non-FineMappingResult input: ", path) next } - cs <- tryCatch(as.data.frame(getCs(fmr)), error = function(e) NULL) - if (is.null(cs) || nrow(cs) == 0L || !"variant_id" %in% names(cs)) next + # Unguarded: an FMR with no credible sets yields a 0-row table, which the + # next line skips; an accessor error is a real break and must surface. + cs <- as.data.frame(credibleSets(fmr)) + if (nrow(cs) == 0L || !"variant_id" %in% names(cs)) next label <- if (methods::is(fmr, "GwasFineMappingResult")) { paste(cs$study, cs$method, cs$blockId, sep = "|") } else { diff --git a/code/script/pecotmr_integration/gwas_rss_plot.R b/code/script/pecotmr_integration/gwas_rss_plot.R index 9f8e4fc76..3985a900a 100644 --- a/code/script/pecotmr_integration/gwas_rss_plot.R +++ b/code/script/pecotmr_integration/gwas_rss_plot.R @@ -32,6 +32,9 @@ parser <- add_argument(parser, "--height", type = "numeric", default = 3) argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + write_empty <- function(msg) { message(msg, "; writing an empty plot to ", argv$output) dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) @@ -48,12 +51,13 @@ if (nrow(fmr) == 0L) { quit(save = "no") } -# One aggregated PIP table for the collection. getTopLoci(signalCutoff = 0) +# One aggregated PIP table for the collection. topLoci(signalCutoff = 0) # returns every variant with its decoded `pos` and the row-identity columns # (study / method / region_id), so the panel label and position come straight -# off the table -- one panel per (study, method, region_id). -tl <- getTopLoci(fmr, signalCutoff = 0) -if (is.null(tl) || nrow(tl) == 0L) { +# off the table -- one panel per (study, method, region_id). It is a GRanges +# from 0.9.0 on, so flatten it: nrow() on a GRanges is NULL, not a count. +tl <- as_variant_df(topLoci(fmr, signalCutoff = 0)) +if (nrow(tl) == 0L) { write_empty("No PIP vectors in input") quit(save = "no") } diff --git a/code/script/pecotmr_integration/gwas_sumstats_construct.R b/code/script/pecotmr_integration/gwas_sumstats_construct.R index 824ffb87f..708eec36a 100644 --- a/code/script/pecotmr_integration/gwas_sumstats_construct.R +++ b/code/script/pecotmr_integration/gwas_sumstats_construct.R @@ -137,6 +137,9 @@ parser <- add_argument(parser, "--output", help = "Output RDS path", type = "character") argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + splitCsv <- function(s) { if (is.null(s) || !nzchar(s)) return(character(0)) trimws(strsplit(s, ",", fixed = TRUE)[[1L]]) @@ -307,18 +310,20 @@ gss_out <- if (argv$skip_qc) { effN <- as.logical(argv$effective_n) if (is.na(effN)) stop("--effective-n must be TRUE or FALSE (got: ", argv$effective_n, ")") - qc_call_args <- c(list(gss, - zMismatchQc = qc_method, - impute = isTRUE(argv$impute), - pipCutoffToSkip = argv$pip_cutoff_to_skip, - mafCutoff = argv$maf, - effectiveN = effN), - if (isTRUE(argv$allele_flip_kriging)) - list(alleleFlipKriging = TRUE) else list(), - if (!is.null(skip_region_vec)) - list(skipRegion = skip_region_vec) else list(), - qc_extra) - do.call(summaryStatsQc, qc_call_args) + # Still assembled under the FLAT argument names -- the wrapper flags and + # --qc-args both speak them -- then routed into pecotmr 0.9.0's bundles + # (panelFilterParam / signalScreenArgs / ...) by summary_stats_qc_args(). + qc_flat <- c(list(zMismatchQc = qc_method, + impute = isTRUE(argv$impute), + pipCutoffToSkip = argv$pip_cutoff_to_skip, + mafCutoff = argv$maf, + effectiveN = effN), + if (isTRUE(argv$allele_flip_kriging)) + list(alleleFlipKriging = TRUE) else list(), + if (!is.null(skip_region_vec)) + list(skipRegion = skip_region_vec) else list(), + qc_extra) + do.call(summaryStatsQc, c(list(gss), summary_stats_qc_args(qc_flat))) } dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) @@ -331,7 +336,9 @@ cat(sprintf("Wrote %sGwasSumStats (%d stud%s: %s) over %s to %s\n", # ----- Per-entry QC diagnostics summary (to the step log) ------------------- # Surface the qcInfo audit so the SoS log shows what QC did, not just the path. if (!argv$skip_qc) { - audit <- tryCatch(getQcInfo(gss_out)$entryAudit, error = function(e) NULL) + # Unguarded: qcInfo() is a plain slot accessor, and an absent audit is a + # NULL element, which the loop below skips. An error here is a real break. + audit <- qcInfo(gss_out)$entryAudit for (i in seq_along(audit)) { a <- audit[[i]] if (is.null(a)) next diff --git a/code/script/pecotmr_integration/mash_covariance.R b/code/script/pecotmr_integration/mash_covariance.R index d9ff2d72f..bd1a02f6e 100644 --- a/code/script/pecotmr_integration/mash_covariance.R +++ b/code/script/pecotmr_integration/mash_covariance.R @@ -50,7 +50,7 @@ alpha <- if (toupper(argv$effect_model) == "EZ") 1 else 0 dat <- readRDS(argv$data) strong <- qtlSumStatsFromBetaMatrix( - as.matrix(dat$strong.b), as.matrix(dat$strong.s), study = "mash") + as.matrix(dat$strong.b), as.matrix(dat$strong.s), studyName = "mash") components <- split_pecotmr_names(argv$component, MASH_COMPONENT_ALIASES) nPcs <- if (is.na(argv$npc)) NULL else argv$npc diff --git a/code/script/pecotmr_integration/mash_fit.R b/code/script/pecotmr_integration/mash_fit.R index bd1d83226..5cd1f6c3f 100644 --- a/code/script/pecotmr_integration/mash_fit.R +++ b/code/script/pecotmr_integration/mash_fit.R @@ -49,7 +49,7 @@ U <- if (is.list(prior) && !is.null(prior$U)) prior$U else prior # Wrap the random (fit) subset as a beta-scale QtlSumStats. ldSketch stays NULL: # mash needs no LD reference. random <- qtlSumStatsFromBetaMatrix( - as.matrix(dat$random.b), as.matrix(dat$random.s), study = "mash") + as.matrix(dat$random.b), as.matrix(dat$random.s), studyName = "mash") model <- mashModelFit(list(random = random), alpha = alpha, priorCovariances = U, vhat = vhat, diff --git a/code/script/pecotmr_integration/mash_posterior.R b/code/script/pecotmr_integration/mash_posterior.R index 50230e102..22a96c1e2 100644 --- a/code/script/pecotmr_integration/mash_posterior.R +++ b/code/script/pecotmr_integration/mash_posterior.R @@ -53,7 +53,7 @@ model <- if (is.list(mm) && !is.null(mm$mash_model)) mm$mash_model else mm bhat <- as.matrix(dat[[argv$bhat_key]]) shat <- as.matrix(dat[[argv$shat_key]]) -target <- qtlSumStatsFromBetaMatrix(bhat, shat, study = "mash") +target <- qtlSumStatsFromBetaMatrix(bhat, shat, studyName = "mash") conds <- colnames(bhat) # Resolve --exclude-condition tokens: numeric tokens are 1-based column indices diff --git a/code/script/pecotmr_integration/mash_prior.R b/code/script/pecotmr_integration/mash_prior.R index 88753bdb2..eacf0a0d1 100644 --- a/code/script/pecotmr_integration/mash_prior.R +++ b/code/script/pecotmr_integration/mash_prior.R @@ -56,11 +56,15 @@ argv <- parse_args(p) .d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) source(file.path(.d, "pecotmr_aliases.R")) +# cov_ed / ud / ud_ted -> the 0.9.0 engine spelling (a name, or the udr +# options object that carries the ED-vs-TED choice). +engine_arg <- mash_prior_engine(argv$engine) + alpha <- if (toupper(argv$effect_model) == "EZ") 1 else 0 dat <- readRDS(argv$data) strong <- qtlSumStatsFromBetaMatrix( - as.matrix(dat$strong.b), as.matrix(dat$strong.s), study = "mash") + as.matrix(dat$strong.b), as.matrix(dat$strong.s), studyName = "mash") vhat <- if (nzchar(argv$vhat_data)) readRDS(argv$vhat_data) else NULL @@ -71,13 +75,13 @@ if (nzchar(argv$component_files)) { priorComponents <- do.call(c, lapply(files, readRDS)) prior <- mashPriorCovariances(list(strong = strong), alpha = alpha, vhat = vhat, priorComponents = priorComponents, - engine = argv$engine, setSeed = argv$seed) + engine = engine_arg, setSeed = argv$seed) } else { # Self-contained mode: build the components here. components <- split_pecotmr_names(argv$components, MASH_COMPONENT_ALIASES) nPcs <- if (is.na(argv$npc)) NULL else argv$npc prior <- mashPriorCovariances(list(strong = strong), alpha = alpha, vhat = vhat, - components = components, engine = argv$engine, + components = components, engine = engine_arg, nPcs = nPcs, setSeed = argv$seed) } diff --git a/code/script/pecotmr_integration/mash_sumstats_construct.R b/code/script/pecotmr_integration/mash_sumstats_construct.R index 4194d8381..1780367e9 100644 --- a/code/script/pecotmr_integration/mash_sumstats_construct.R +++ b/code/script/pecotmr_integration/mash_sumstats_construct.R @@ -100,7 +100,7 @@ if (nrow(zMat) == 0L) # Wrap the region's z-matrix as a multi-context QtlSumStats (one column per # condition). qtlSumStatsFromZMatrix stamps a non-empty qcInfo, so mashInput() # accepts it directly; ldSketch stays NULL (mash needs no LD reference). -qss <- qtlSumStatsFromZMatrix(zMat, study = argv$study, context = conditions, +qss <- qtlSumStatsFromZMatrix(zMat, studyName = argv$study, context = conditions, trait = "mash") dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) diff --git a/code/script/pecotmr_integration/mash_vhat.R b/code/script/pecotmr_integration/mash_vhat.R index f5fee7e57..b96de8f2d 100644 --- a/code/script/pecotmr_integration/mash_vhat.R +++ b/code/script/pecotmr_integration/mash_vhat.R @@ -55,7 +55,7 @@ dat <- readRDS(argv$data) mk <- function(bKey, sKey) { if (is.null(dat[[bKey]])) return(NULL) qtlSumStatsFromBetaMatrix(as.matrix(dat[[bKey]]), as.matrix(dat[[sKey]]), - study = "mash") + studyName = "mash") } ssl <- Filter(Negate(is.null), list( strong = mk("strong.b", "strong.s"), diff --git a/code/script/pecotmr_integration/pecotmr_aliases.R b/code/script/pecotmr_integration/pecotmr_aliases.R index f6f45f9a6..6f26d4be7 100644 --- a/code/script/pecotmr_integration/pecotmr_aliases.R +++ b/code/script/pecotmr_integration/pecotmr_aliases.R @@ -16,7 +16,8 @@ # here. Both spellings are accepted, so callers already using the pecotmr # names work unchanged. # -# This is argument marshalling only; no analysis logic (that lives in pecotmr). +# This is argument / return-value marshalling only; no analysis logic (that +# lives in pecotmr). # mashCovarianceComponents(components=) / mashPriorCovariances(components=) MASH_COMPONENT_ALIASES <- c(flash_nonneg = "flashNonneg") @@ -24,6 +25,37 @@ MASH_COMPONENT_ALIASES <- c(flash_nonneg = "flashNonneg") # mashResidualCorrelation(method=) MASH_VHAT_ALIASES <- c(simple_specific = "simpleSpecific") +# mashPriorCovariances(engine=). 0.9.0 renamed the engines camelCase AND folded +# udr's ED-vs-TED choice out of the engine name into the options object, so the +# three names this side exposes (and mixture_prior.ipynb's [ud] / +# [ud_unconstrained] / [ed_bovy] steps pass) no longer map one-to-one: +# +# cov_ed -> "covEd" mashr's exported cov_ed(), unchanged +# ud -> covUdr + unconstrained.update = "ed" +# ud_ted -> covUdr + unconstrained.update = "ted" +# +# Both udr engines pin the update explicitly. Left unset, 0.9.0 hands the +# choice to udr, which picks TED on z-scale data and ED on the beta scale -- +# correct on its own terms, but it would silently merge the `ud` and `ud_ted` +# steps into the same fit on z-scale input, and those are two different +# documented engines on this side. +mash_prior_engine <- function(value) { + switch(value, + cov_ed = "covEd", + ud = CovUdrOptions(fit = UdFitOptions(unconstrained.update = "ed")), + ud_ted = CovUdrOptions(fit = UdFitOptions(unconstrained.update = "ted")), + # Already a pecotmr spelling (or an unknown one): hand it over untouched + # and let pecotmr report it. + value) +} + +# fineMappingPipeline(methods=). Its token set is susie / susieInf / susieAsh / +# ser / mvsusie / fsusie -- only the two multiword ones need translating. +FM_METHOD_ALIASES <- c( + susie_ash = "susieAsh", + susie_inf = "susieInf" +) + # twasWeightsPipeline(methods=). Mirrors pecotmr's # .twasKnownMethodLookupNames(); tokens that are already single words # (susie, mrash, enet, lasso, scad, mcp, l0learn, mvsusie, mrmash) need no alias. @@ -55,3 +87,87 @@ split_pecotmr_names <- function(value, aliases, split = "[ ,]+") { parts <- parts[nzchar(parts)] apply_pecotmr_aliases(parts, aliases) } + +# pecotmr 0.9.0 returns the top-loci tables (topLoci(), overlapTopLoci()) as +# GRanges rather than data.frames. Every consumer here wants the flat +# per-variant table, and as.data.frame() on a GRanges PREPENDS five range +# columns (seqnames/start/end/width/strand) that are not part of the published +# schema -- chrom and pos are already mcols. mcols() takes the data columns +# only, so the frame keeps exactly the shape it had before 0.9.0. +# +# Only for those variant tables: a GRanges whose RANGES are the data +# (fsusieAffectedRegions()) must keep going through as.data.frame(). +as_variant_df <- function(x) { + if (methods::is(x, "GenomicRanges")) return(as.data.frame(S4Vectors::mcols(x))) + as.data.frame(x) +} + +# ---- summaryStatsQc() argument routing ------------------------------------- +# +# pecotmr 0.9.0 split summaryStatsQc()'s flat argument list into four bundles. +# This side still speaks the flat names: they are the wrapper flags, and +# --qc-args is a documented JSON passthrough whose own examples +# ({"mafCutoff":0.01,"nCutoff":0,"alleleFlipKriging":true}, rss_analysis.ipynb) +# span three different destinations. So route each flat name to whichever +# bundle now owns it, and leave everything else at top level -- where +# alleleFlipKriging / effectiveN / impute / skipRegion / keepVariants / +# matchMinProp still live, and where summaryStatsQc reports an unknown name +# itself rather than this table swallowing it. +SUMSTATS_QC_RENAMES <- c(zMismatchQc = "ldMismatchQcMethod") + +SUMSTATS_QC_BUNDLE_CTOR <- c( + sumstatsFilterArgs = "SumstatsFilterParam", + panelFilterParam = "PanelFilterParam", + signalScreenArgs = "SignalScreenParam", + sumstatsCleaningArgs = "SumstatsCleaningParam") + +# flat name -> "/" +SUMSTATS_QC_ROUTING <- c( + removeIndels = "sumstatsFilterArgs/removeIndels", + removeStrandAmbiguous = "sumstatsFilterArgs/removeStrandAmbiguous", + infoCutoff = "sumstatsFilterArgs/infoCutoff", + nCutoff = "sumstatsFilterArgs/nCutoff", + mafCutoff = "panelFilterParam/mafCutoff", + macCutoff = "panelFilterParam/macCutoff", + imissCutoff = "panelFilterParam/imissCutoff", + # the *CutoffToSkip screens became one bundle, one field per metric + pipCutoffToSkip = "signalScreenArgs/pip", + absZCutoffToSkip = "signalScreenArgs/absZ", + bfCutoffToSkip = "signalScreenArgs/bf", + logBfCutoffToSkip = "signalScreenArgs/logBf", + coerceNumeric = "sumstatsCleaningArgs/coerceNumeric", + normalizeChr = "sumstatsCleaningArgs/normalizeChr", + dropNonstandardChr = "sumstatsCleaningArgs/dropNonstandardChr", + dropMissData = "sumstatsCleaningArgs/dropMissData", + dropPOutOfRange = "sumstatsCleaningArgs/dropPOutOfRange", + clampSmallP = "sumstatsCleaningArgs/clampSmallP", + smallPFloor = "sumstatsCleaningArgs/smallPFloor", + dropZeroEffect = "sumstatsCleaningArgs/dropZeroEffect", + dropNonpositiveSe = "sumstatsCleaningArgs/dropNonpositiveSe") + +# Turn a flat named list of summaryStatsQc() arguments into the 0.9.0 call +# shape: the routed ones grouped into their bundle constructors, the rest +# passed straight through. +summary_stats_qc_args <- function(flat) { + nms <- names(flat) + hit <- match(nms, names(SUMSTATS_QC_RENAMES)) + names(flat) <- ifelse(is.na(hit), nms, unname(SUMSTATS_QC_RENAMES[hit])) + + out <- list() + bundled <- list() + for (nm in names(flat)) { + # [ ] not [[ ]]: on a named character vector a missing name is an + # out-of-bounds ERROR with [[ ]], but NA with [ ]. + route <- unname(SUMSTATS_QC_ROUTING[nm]) + if (is.na(route)) { + out[nm] <- list(flat[[nm]]) # out[nm] <- list(x) keeps a NULL value + next + } + part <- strsplit(route, "/", fixed = TRUE)[[1L]] + bundled[[part[[1L]]]][[part[[2L]]]] <- flat[[nm]] + } + for (b in names(bundled)) + out[[b]] <- do.call(SUMSTATS_QC_BUNDLE_CTOR[[b]], bundled[[b]]) + out +} + diff --git a/code/script/pecotmr_integration/qtl_association_postprocessing.R b/code/script/pecotmr_integration/qtl_association_postprocessing.R index 5386a3fdd..e7af510ff 100644 --- a/code/script/pecotmr_integration/qtl_association_postprocessing.R +++ b/code/script/pecotmr_integration/qtl_association_postprocessing.R @@ -118,14 +118,19 @@ entries <- lapply(genes, mkEntry) nVar <- as.numeric(regional$n_variants) nVarFcol <- if (!is.null(nvFilt)) as.numeric(nvFilt[genes]) else NULL -qssArgs <- list( - study = rep(argv$study, length(genes)), context = rep(argv$context, length(genes)), - trait = genes, entry = entries, genome = argv$genome, - n_variants = nVar, p_beta = as.numeric(regional$p_beta), +# pecotmr 0.9.0: `study` -> `studyName`, and the per-trait statistic columns +# travel in `extraCols` instead of `...` (which the constructor no longer has). +extraCols <- list( + n_variants = nVar, + p_beta = as.numeric(regional$p_beta), beta_shape1 = as.numeric(regional$beta_shape1), beta_shape2 = as.numeric(regional$beta_shape2)) -if (!is.null(nVarFcol)) qssArgs$n_variants_filtered <- nVarFcol -qss <- do.call(QtlSumStats, qssArgs) +if (!is.null(nVarFcol)) extraCols$n_variants_filtered <- nVarFcol +qss <- QtlSumStats( + studyName = rep(argv$study, length(genes)), + context = rep(argv$context, length(genes)), + trait = genes, entry = entries, genome = argv$genome, + extraCols = extraCols) filtering <- (argv$maf_cutoff > 0 || argv$cis_window > 0) && !is.null(nVarFcol) r <- qtlAssociationPostprocess( @@ -157,7 +162,7 @@ if (!is.na(argv$output_dir)) { # Per-method EVENT-significance column (which genes pass): permutation and the # q-value SNP method gate events on the Storey q of the beta-approx permutation # p (q_beta); Bonferroni gates on its BH-FDR of the per-gene min. This matches - # the variant-level rules getSignificantQtls() applies. + # the variant-level rules significantQtls() applies. flavours <- c(permutation = "q_beta", bonferroni_original = "fdr_bonferroni_min_original", bonferroni_filtered = "fdr_bonferroni_min_filtered", @@ -172,10 +177,13 @@ if (!is.na(argv$output_dir)) { if (nrow(events) > 0) readr::write_tsv(events, file.path(od, paste0(base, ".significant_events.", m, ".tsv.gz"))) # significant variant-level QTL (derived, never stored) - sig <- tryCatch(getSignificantQtls(r, m, threshold = argv$fdr_threshold), - error = function(e) NULL) + # Unguarded: significantQtls() aborts only when the object did not come + # from qtlAssociationPostprocess(), and `r` is that function's own output, + # so the one catchable condition cannot arise here. "No significant QTL" + # is a zero-length return, which the next line handles. + sig <- significantQtls(r, m, threshold = argv$fdr_threshold) nqtl <- 0L - if (!is.null(sig) && length(sig) > 0) { + if (length(sig) > 0) { sdf <- as.data.frame(sig) readr::write_tsv(sdf, file.path(od, paste0(base, ".significant_qtl.", m, ".tsv.gz"))) nqtl <- length(sig) diff --git a/code/script/pecotmr_integration/qtl_dataset_construct.R b/code/script/pecotmr_integration/qtl_dataset_construct.R index b3c60f9a6..397cc3cb1 100644 --- a/code/script/pecotmr_integration/qtl_dataset_construct.R +++ b/code/script/pecotmr_integration/qtl_dataset_construct.R @@ -117,25 +117,31 @@ genotype_prefix <- abs_path(argv$genotype_prefix) # The manifest's cond/path/cov_path columns are recognised via the loader's # snake_case aliases (cond -> context, path -> phenotypePath, # cov_path -> covariatePath). +# The seven genotype filters travel as one GenotypeFilterParam() bundle. Every +# field is passed explicitly: to QtlDataset() the bundle is a SPECIFICATION, so +# an unset field would take the constructor's own default rather than this +# script's CLI default -- which happen to coincide today, but pinning them keeps +# the wrapper's documented defaults authoritative. qd <- loadQtlDatasetFromManifest( manifest = argv$phenotype_manifest, - study = argv$study, + studyName = argv$study, genotypes = genotype_prefix, genotypeCovariates = genotype_covariates, scaleResiduals = !isTRUE(argv$no_scale_residuals), - mafCutoff = argv$maf_cutoff, - macCutoff = argv$mac_cutoff, - xvarCutoff = argv$xvar_cutoff, - imissCutoff = argv$imiss_cutoff, - keepSamples = keep_samples, - keepVariants = keep_variants, - keepIndel = !isTRUE(argv$drop_indel), + genotypeFilterParam = GenotypeFilterParam( + mafCutoff = argv$maf_cutoff, + macCutoff = argv$mac_cutoff, + xvarCutoff = argv$xvar_cutoff, + imissCutoff = argv$imiss_cutoff, + keepSamples = keep_samples, + keepVariants = keep_variants, + keepIndel = !isTRUE(argv$drop_indel)), transposeCovariates = isTRUE(argv$transpose_covariates)) dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) saveRDS(qd, argv$output) -contexts <- getContexts(qd) +ctx <- contexts(qd) cat(sprintf("Wrote QtlDataset for study '%s' (%d contexts: %s) to %s\n", - argv$study, length(contexts), - paste(contexts, collapse = ", "), + argv$study, length(ctx), + paste(ctx, collapse = ", "), argv$output)) diff --git a/code/script/pecotmr_integration/qtl_enrichment.R b/code/script/pecotmr_integration/qtl_enrichment.R index f2c4df755..668e7f6cb 100644 --- a/code/script/pecotmr_integration/qtl_enrichment.R +++ b/code/script/pecotmr_integration/qtl_enrichment.R @@ -78,13 +78,16 @@ gwasFmr <- .loadCombine(argv$gwas_fine_mapping) # compatibility with a pecotmr that predates the argument. seed_val <- if (length(argv$seed) == 1L && !is.na(argv$seed)) as.integer(argv$seed) else NULL +# pecotmr 0.9.0 groups the estimator's own knobs into QtlEnrichmentParam(); +# numThreads and seed stay top-level pipeline arguments. enr_args <- list( gwasFineMappingResult = gwasFmr, qtlFineMappingResult = qtlFmr, - numGwas = if (is.na(argv$num_gwas)) NULL else argv$num_gwas, - piQtl = if (is.na(argv$pi_qtl)) NULL else argv$pi_qtl, - lambda = argv$lambda, - impN = argv$imp_n, + methodArgs = QtlEnrichmentParam( + numGwas = if (is.na(argv$num_gwas)) NULL else argv$num_gwas, + piQtl = if (is.na(argv$pi_qtl)) NULL else argv$pi_qtl, + lambda = argv$lambda, + impN = argv$imp_n), numThreads = argv$ncore) if (!is.null(seed_val)) enr_args$seed <- seed_val res <- do.call(qtlEnrichmentPipeline, enr_args) diff --git a/code/script/pecotmr_integration/qtl_sumstats_construct.R b/code/script/pecotmr_integration/qtl_sumstats_construct.R index 4b8e4a5c8..aec1eb38f 100755 --- a/code/script/pecotmr_integration/qtl_sumstats_construct.R +++ b/code/script/pecotmr_integration/qtl_sumstats_construct.R @@ -73,6 +73,9 @@ parser <- add_argument(parser, "--output", type = "character", default = "", help = "Output RDS") argv <- parse_args(parser) +.d <- dirname(sub("^--file=", "", grep("^--file=", commandArgs(FALSE), value = TRUE)[1L])) +source(file.path(.d, "pecotmr_aliases.R")) + for (a in c("sumstats", "study", "context", "trait", "ld_sketch", "output")) if (!nzchar(argv[[a]])) stop("--", gsub("_", "-", a), " is required") @@ -100,18 +103,23 @@ out <- if (argv$skip_qc) { message("--skip-qc set; serialising the raw QtlSumStats without summaryStatsQc().") qss } else { - summaryStatsQc(qss, - mafCutoff = argv$maf, - macCutoff = argv$mac, - imissCutoff = argv$imiss, - zMismatchQc = argv$z_mismatch_qc, - pipCutoffToSkip = argv$pip_cutoff_to_skip) + # Flat flag names routed into pecotmr 0.9.0's bundles: maf/mac/imiss -> + # PanelFilterParam, zMismatchQc -> ldMismatchQcMethod, pipCutoffToSkip -> + # SignalScreenParam(pip =). See summary_stats_qc_args() in pecotmr_aliases.R. + do.call(summaryStatsQc, + c(list(qss), + summary_stats_qc_args(list( + mafCutoff = argv$maf, + macCutoff = argv$mac, + imissCutoff = argv$imiss, + zMismatchQc = argv$z_mismatch_qc, + pipCutoffToSkip = argv$pip_cutoff_to_skip)))) } dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) saveRDS(out, argv$output, compress = "xz") cat(sprintf("Wrote QtlSumStats for %s/%s/%s (%d variant(s), LD sketch %s) to %s\n", argv$study, argv$context, argv$trait, length(out[[1L]]), - if (is.null(getLdSketch(out))) "NULL" else - paste(dim(getLdSketch(out)), collapse = " x "), + if (is.null(ldSketch(out))) "NULL" else + paste(dim(ldSketch(out)), collapse = " x "), argv$output)) diff --git a/code/script/pecotmr_integration/sldsc_postprocess.R b/code/script/pecotmr_integration/sldsc_postprocess.R index 89e49f819..b537328ef 100644 --- a/code/script/pecotmr_integration/sldsc_postprocess.R +++ b/code/script/pecotmr_integration/sldsc_postprocess.R @@ -130,7 +130,7 @@ if (argv$maf_cutoff > 0) { frq <- readSldscFrq(argv$frqfile_dir, plinkName = argv$plink_name) } -sldscData <- SldscData(annot = annot, frq = frq, traits = traitList) +sldscData <- SldscData(annotData = annot, frqData = frq, traits = traitList) res <- sldscPostprocessingPipeline( sldscData, diff --git a/code/script/pecotmr_integration/twas.R b/code/script/pecotmr_integration/twas.R index 1fd20da4d..72b527075 100644 --- a/code/script/pecotmr_integration/twas.R +++ b/code/script/pecotmr_integration/twas.R @@ -83,18 +83,24 @@ fmr <- if (nzchar(fmr_path) && fmr_path != "." && file.exists(fmr_path)) { NULL } +# pecotmr 0.9.0 grouped the eight rsq*/mr* arguments into two bundles. Pure +# regrouping -- the CLI values carry over unchanged (--rsq-option feeds +# `metric`, --rsq-pval-option feeds `pvalMetric`, whose default c("adj_rsq_pval", +# "pval") is this flag's default too). res <- causalInferencePipeline( gwasSumStats = gss, twasWeights = tw, fineMappingResult = fmr, - rsqCutoff = argv$rsq_cutoff, - rsqPvalCutoff = argv$rsq_pval_cutoff, - rsqOption = argv$rsq_option, - rsqPvalOption = rsq_pval_option, - mrPipCutoff = argv$mr_pip_cutoff, - mrMethod = argv$mr_method, - mrCpipCutoff = argv$mr_cpip_cutoff, - mrPvalCutoff = argv$mr_pval_cutoff, + weightSelectionArgs = WeightSelectionParam( + cutoff = argv$rsq_cutoff, + pvalCutoff = argv$rsq_pval_cutoff, + metric = argv$rsq_option, + pvalMetric = rsq_pval_option), + mrArgs = MrParam( + method = argv$mr_method, + pipCutoff = argv$mr_pip_cutoff, + cpipCutoff = argv$mr_cpip_cutoff, + pvalCutoff = argv$mr_pval_cutoff), combineMethods = combine_methods) dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE) diff --git a/code/script/pecotmr_integration/twas_weights.R b/code/script/pecotmr_integration/twas_weights.R index 3c039dcee..567e73f21 100644 --- a/code/script/pecotmr_integration/twas_weights.R +++ b/code/script/pecotmr_integration/twas_weights.R @@ -148,7 +148,10 @@ qd_filter_overrides <- function(argv) { keepIndel = if (isTRUE(argv$drop_indel)) FALSE else NULL, keepSamples = read_ids(argv$keep_samples), keepVariants = read_ids(argv$keep_variants)) - ov[!vapply(ov, is.null, logical(1))] + # pecotmr 0.9.0: one GenotypeFilterParam() bundle, where an UNSET field + # leaves the dataset's construct-time value alone -- the override role. + ov <- ov[!vapply(ov, is.null, logical(1))] + list(genotypeFilterParam = do.call(GenotypeFilterParam, ov)) } # Seed up front for reproducible fits (mirrors the legacy susie_twas set.seed). @@ -264,13 +267,15 @@ fmr <- if (nzchar(fmr_path) && fmr_path != "." && file.exists(fmr_path)) { # --maf-cutoff/--xvar-cutoff/... args fine-mapping uses (variant QC is a data # property, applied identically); the cv* knobs control the cross-validated # predictive-performance refits. +# maxCvVariants / cvFolds / cvThreads are now one CrossValidationParam bundle. tw_args <- c(list(methods = methods_arg, mashPrior = mash_prior, cisWindow = argv$cis_window, contexts = contexts_arg, - maxCvVariants = argv$max_cv_variants, - cvFolds = argv$cv_folds, - cvThreads = argv$cv_threads, + crossValidationArgs = CrossValidationParam( + folds = argv$cv_folds, + numThreads = argv$cv_threads, + maxVariants = argv$max_cv_variants), fineMappingResult = fmr), qd_filter_overrides(argv)) # Opt-in joint axis / resume-from-checkpoint, added only when supplied (keeps @@ -282,10 +287,14 @@ res <- if (has_qss) { # RSS mode: no trait/region selection, and the individual-level-only knobs # (cisWindow + the genotype filters) do not apply to a summary-statistics # collection, so only the shared arguments are forwarded. + # No mashPrior / seed here: the QtlSumStats method takes neither (pre-0.9.0 + # its signature ended in `...`, which collected both into a dots list that + # the summary-statistics worker never read, so dropping them changes + # nothing). Seeding on this path is the set.seed() above. rss_args <- tw_args[intersect( names(tw_args), - c("methods", "mashPrior", "contexts", "fineMappingResult", - "jointSpecification", "twasWeights", "seed"))] + c("methods", "contexts", "fineMappingResult", + "jointSpecification", "twasWeights"))] do.call(twasWeightsPipeline, c(list(qd), rss_args)) } else if (has_region) { do.call(twasWeightsPipeline, diff --git a/tests/fixtures/ctwas/expected/ctwas_assemble.rds b/tests/fixtures/ctwas/expected/ctwas_assemble.rds index b908582d7..cf95c1727 100644 Binary files a/tests/fixtures/ctwas/expected/ctwas_assemble.rds and b/tests/fixtures/ctwas/expected/ctwas_assemble.rds differ diff --git a/tests/fixtures/ctwas/expected/ctwas_est.rds b/tests/fixtures/ctwas/expected/ctwas_est.rds index 626974c55..3f4bb5356 100644 Binary files a/tests/fixtures/ctwas/expected/ctwas_est.rds and b/tests/fixtures/ctwas/expected/ctwas_est.rds differ diff --git a/tests/fixtures/ctwas/expected/ctwas_finemap.rds b/tests/fixtures/ctwas/expected/ctwas_finemap.rds index 83f33e77f..9b6eb1796 100644 Binary files a/tests/fixtures/ctwas/expected/ctwas_finemap.rds and b/tests/fixtures/ctwas/expected/ctwas_finemap.rds differ diff --git a/tests/fixtures/ctwas/expected/expected_manifest.tsv b/tests/fixtures/ctwas/expected/expected_manifest.tsv index 7d8f8b3d2..368b09bdb 100644 --- a/tests/fixtures/ctwas/expected/expected_manifest.tsv +++ b/tests/fixtures/ctwas/expected/expected_manifest.tsv @@ -1,6 +1,6 @@ # expected-output regression manifest — ctwas # columns: output mode rtol atol notes ctwas_manifest.tsv tolerant 1e-6 1e-8 Per-LD-block cTWAS grid from ctwas_manifest.R (21 rows: header + 20 chr22 blocks); gwas_sumstats_rds is an absolute tmp path -> normalize_paths (basename compare). -ctwas_assemble.rds tolerant 1e-6 1e-8 Assembled cTWAS inputs list (z_snp/z_gene/weights/region_info/snp_map/LD_map/...) from ctwas_assemble.R over the blessed 20-block fixture; LD/snp maps embed absolute reference paths -> normalize_paths. -ctwas_est.rds tolerant 1e-6 1e-8 estCtwasParam result list (assemble inputs + region_data/boundary_genes/param) from ctwas_est.R with --fallback-to-prefit --niter 50 (deterministic prefit path); embeds absolute reference paths -> normalize_paths. -ctwas_finemap.rds tolerant 1e-6 1e-8 CtwasResult S4 (1 gene-level + 108 SNP-level rows; geneMaxPip=1, geneTopZ=3.168) from ctwas_finemap.R --L 5 --min-nonsnp-pip 0.5; deterministic single-region susie fine-map, no embedded paths. +ctwas_assemble.rds tolerant 1e-6 1e-8 Assembled cTWAS inputs list (z_snp/z_gene/weights/region_info/snp_map/LD_map/...) from ctwas_assemble.R over the blessed 20-block fixture; LD/snp maps embed absolute reference paths -> normalize_paths. Re-baselined for pecotmr 0.9.0: the ONLY delta is that each per-region snp_map data.frame now carries reset row names (1..n) instead of the parent table's indices; contents identical at 1e-6. +ctwas_est.rds tolerant 1e-6 1e-8 estCtwasParam result list (assemble inputs + region_data/boundary_genes/param) from ctwas_est.R with --fallback-to-prefit --niter 50 (deterministic prefit path); embeds absolute reference paths -> normalize_paths. Re-baselined for pecotmr 0.9.0 (estCtwasParam -> estCtwasGroupPriors): inherits the snp_map row-name reset from ctwas_assemble; everything else, group_prior and group_prior_var included, is identical at 1e-6. +ctwas_finemap.rds tolerant 1e-6 1e-8 CtwasResult S4 (1 gene-level + 108 SNP-level rows; geneMaxPip=1, geneTopZ=3.168) from ctwas_finemap.R --L 5 --min-nonsnp-pip 0.5; deterministic single-region susie fine-map, no embedded paths. Re-baselined for pecotmr 0.9.0: CtwasResultEntry's slots changed (posteriors/susieAlpha/groupPriors/regionInfo), so rds_compare cannot even traverse a pre-0.9.0 serialization. The science is unchanged -- still 1 gene-level + 108 SNP-level rows, geneMaxPip=1, geneTopZ=3.168. diff --git a/tests/fixtures/mash/expected/mash_input.fmr.rds b/tests/fixtures/mash/expected/mash_input.fmr.rds index 945e2326f..973930e1b 100644 Binary files a/tests/fixtures/mash/expected/mash_input.fmr.rds and b/tests/fixtures/mash/expected/mash_input.fmr.rds differ diff --git a/tests/fixtures/mash/expected/mash_input.indep.rds b/tests/fixtures/mash/expected/mash_input.indep.rds index e168961c0..c8245043c 100644 Binary files a/tests/fixtures/mash/expected/mash_input.indep.rds and b/tests/fixtures/mash/expected/mash_input.indep.rds differ diff --git a/tests/fixtures/mash/expected/mash_model.EE.rds b/tests/fixtures/mash/expected/mash_model.EE.rds index be3b8d090..22733edcd 100644 Binary files a/tests/fixtures/mash/expected/mash_model.EE.rds and b/tests/fixtures/mash/expected/mash_model.EE.rds differ diff --git a/tests/fixtures/mash/expected/mixture_prior.EE.prior.rds b/tests/fixtures/mash/expected/mixture_prior.EE.prior.rds index 6ccb39f73..ce21a9266 100644 Binary files a/tests/fixtures/mash/expected/mixture_prior.EE.prior.rds and b/tests/fixtures/mash/expected/mixture_prior.EE.prior.rds differ diff --git a/tests/fixtures/mash/expected/posterior.EE.rds b/tests/fixtures/mash/expected/posterior.EE.rds index d1deffaa7..f58e41fc0 100644 Binary files a/tests/fixtures/mash/expected/posterior.EE.rds and b/tests/fixtures/mash/expected/posterior.EE.rds differ diff --git a/tests/fixtures/mash/expected/prior.cov_ed.EE.rds b/tests/fixtures/mash/expected/prior.cov_ed.EE.rds index f50ec147a..aca7b8c9f 100644 Binary files a/tests/fixtures/mash/expected/prior.cov_ed.EE.rds and b/tests/fixtures/mash/expected/prior.cov_ed.EE.rds differ diff --git a/tests/fixtures/mash/expected/vhat.corshrink.EE.rds b/tests/fixtures/mash/expected/vhat.corshrink.EE.rds index 240c94ff5..decce9eea 100644 Binary files a/tests/fixtures/mash/expected/vhat.corshrink.EE.rds and b/tests/fixtures/mash/expected/vhat.corshrink.EE.rds differ diff --git a/tests/fixtures/mash/expected_manifest.tsv b/tests/fixtures/mash/expected_manifest.tsv index dbacff8ea..14e98315c 100644 --- a/tests/fixtures/mash/expected_manifest.tsv +++ b/tests/fixtures/mash/expected_manifest.tsv @@ -1,20 +1,20 @@ # expected-output regression manifest — mash # columns: output mode rtol atol notes -mixture_prior.EE.prior.rds tolerant 1e-6 1e-8 MASH prior list (U/w/loglik) from mixture_prior.ipynb ed_bovy (cov_ed + canonical/pca, EE, vhat=simple); deterministic, notebook-tier (sos run). x86/ARM CI may need a looser rtol. +mixture_prior.EE.prior.rds tolerant 1e-6 1e-8 MASH prior list (U/w/loglik) from mixture_prior.ipynb ed_bovy (cov_ed + canonical/pca, EE, vhat=simple); deterministic, notebook-tier (sos run). x86/ARM CI may need a looser rtol. Re-baselined for pecotmr 0.9.0: 32 -> 16 components, the duplicated UNREFINED copies are gone (this path supplies components via --component-files, which pecotmr treats as all-data-driven, so every one is ED-refined). The 16 refined matrices are IDENTICAL to the old fixture at 1e-6. vhat.identity.EE.rds tolerant 1e-6 1e-8 mash_vhat.R residual-correlation matrix (identity, EE); deterministic. vhat.simple.EE.rds tolerant 1e-6 1e-8 mash_vhat.R residual-correlation matrix (simple, EE); deterministic. -vhat.corshrink.EE.rds tolerant 1e-6 1e-8 mash_vhat.R residual-correlation matrix (corshrink, EE); deterministic. +vhat.corshrink.EE.rds tolerant 1e-6 1e-8 mash_vhat.R residual-correlation matrix (corshrink, EE); deterministic. Re-baselined for pecotmr 0.9.0: it stopped hardcoding ash.control=list(mixcompdist='halfuniform') into CorShrinkData and now forwards CorShrinkOptions() (empty by default), so this records CorShrink's own default. Verified the whole delta is that one setting -- re-running with method=CorShrinkOptions(ash.control=list(mixcompdist='halfuniform')) reproduces the old fixture exactly. vhat.simple_specific.EE.rds tolerant 1e-6 1e-8 mash_vhat.R residual-correlation matrix (simple_specific, EE); deterministic. -prior.cov_ed.EE.rds tolerant 1e-6 1e-8 mash_prior.R prior list (cov_ed engine, canonical+pca, npc 3, EE); deterministic. +prior.cov_ed.EE.rds tolerant 1e-6 1e-8 mash_prior.R prior list (cov_ed engine, canonical+pca, npc 3, EE); deterministic. Re-baselined for pecotmr 0.9.0: the prior no longer carries BOTH the raw and the ED-refined copy of every component -- canonical components are fixed structural hypotheses and go to the fit unrefined, only the data-driven set is ED-refined (the split mashr's own eQTL vignette makes). 34 -> 17 components; the 13 canonical ones are identical to the old fixture at 1e-6, the 4 ED_PCA/tPCA ones move because ED now fits 4 matrices instead of 17. cov.canonical.EE.rds tolerant 1e-6 1e-8 mash_covariance.R canonical component list (13 matrices, EE); deterministic. cov.pca.EE.rds tolerant 1e-6 1e-8 mash_covariance.R pca component list (4 matrices, npc 3, EE); deterministic. cov.flash.EE.rds tolerant 1e-6 1e-8 mash_covariance.R flash component list (EE); flashier self-seeds — empty on the weak toy input (0 factors), a valid regression baseline. cov.flash_nonneg.EE.rds tolerant 1e-6 1e-8 mash_covariance.R flash_nonneg component list (EE); empty on weak toy input (0 factors). mash_sumstats.region1.rds tolerant 1e-6 1e-8 QtlSumStats S4 (40 variants x 2 conditions) from mash_sumstats_construct.R over tensorqtl z-score TSVs; deterministic, no embedded paths. Re-baselined for pecotmr 0.7.8 (#591 class cleanup): the genome build moved off the SumStatsBase `genome` slot into seqinfo(), so this fixture now records genome=GRCh38; values are unchanged (old-vs-new all.equal 1e-6 with the build neutralized). -mash_model.EE.rds tolerant 1e-6 1e-8 mash_fit.R fitted model list(mash_model, vhat_file, prior_file) over the mash_model_chain (vhat=simple -> prior cov_ed canonical+pca npc 3 -> fit), EE; reproducible. vhat_file/prior_file embed the chain tmp paths -> compared with normalize_paths=True (rds_compare normPaths basenames the character leaves). -posterior.EE.rds tolerant 1e-6 1e-8 mash_posterior.R posterior list (PosteriorMean/SD/lfdr/NegativeProb/lfsr/PosteriorCov; 17 variants x 8 conditions) on strong.b/strong.s given the fitted model, EE; deterministic, no embedded paths. +mash_model.EE.rds tolerant 1e-6 1e-8 Re-baselined for pecotmr 0.9.0, downstream of the prior change above: fitted_g$Ulist 34 -> 17 components and pi 749 -> 375, loglik essentially unchanged (459.37 -> 459.06). mash_fit.R fitted model list(mash_model, vhat_file, prior_file) over the mash_model_chain (vhat=simple -> prior cov_ed canonical+pca npc 3 -> fit), EE; reproducible. vhat_file/prior_file embed the chain tmp paths -> compared with normalize_paths=True (rds_compare normPaths basenames the character leaves). +posterior.EE.rds tolerant 1e-6 1e-8 Re-baselined for pecotmr 0.9.0, downstream of the prior change above (the posterior moments move 17-20%). mash_posterior.R posterior list (PosteriorMean/SD/lfdr/NegativeProb/lfsr/PosteriorCov; 17 variants x 8 conditions) on strong.b/strong.s given the fitted model, EE; deterministic, no embedded paths. posterior_contrast.rds tolerant 1e-6 1e-8 mash_posterior_contrast.R contrast data.frame (17 variants x 108 contrasts) from the mash_posterior fixture, cells ALL..Oli; deterministic. row.names are stable variant ids (mash::mash::varN), no embedded paths. mash_input.qss.rds tolerant 1e-6 1e-8 mash_preprocessing.R mash_input from a QtlSumStats region (mash_sumstats_construct.R over Mic/Ast tensorqtl TSVs), n-random/n-null 15, sig-p-cutoff 0.1, seed 1; seeded-reproducible, no embedded paths. -mash_input.fmr.rds tolerant 1e-6 1e-8 mash_preprocessing.R mash_input from the committed QtlFineMappingResult, n-random/n-null 15, coverage 0.95, sig-p-cutoff 0.1, seed 1; seeded-reproducible. -mash_input.indep.rds tolerant 1e-6 1e-8 mash_preprocessing.R mash_input from the QtlFineMappingResult with --independent-variant-list, n-random/n-null 10, coverage 0.95, sig-p-cutoff 0.1, seed 1; seeded-reproducible. +mash_input.fmr.rds tolerant 1e-6 1e-8 mash_preprocessing.R mash_input from the committed QtlFineMappingResult, n-random/n-null 15, coverage 0.95, sig-p-cutoff 0.1, seed 1; seeded-reproducible. Re-baselined for pecotmr 0.9.0: strong.b/strong.s/strong.z now carry variant-id rownames (they were NULL); values identical at 1e-6 with dimnames stripped. +mash_input.indep.rds tolerant 1e-6 1e-8 mash_preprocessing.R mash_input from the QtlFineMappingResult with --independent-variant-list, n-random/n-null 10, coverage 0.95, sig-p-cutoff 0.1, seed 1; seeded-reproducible. Re-baselined for pecotmr 0.9.0: same new strong.* rownames as mash_input.fmr.rds; values identical at 1e-6 with dimnames stripped. feature_score_merge.nsig.csv tolerant 1e-6 1e-8 mash_feature_score_merge.R merged CSV of an nsig feature-score chunk (GENE1 x 8 conditions, nsig_ratio) over the mash_posterior contrast; deterministic. write.csv(row.names=FALSE) so no path-derived row names. (The per-method mash_feature_score.R RDS is NOT value-compared: its row.names embed the tmp contrast path, unfixable by the frozen rds_compare canon.) diff --git a/tests/fixtures/mnm_regression/expected/expected_manifest.tsv b/tests/fixtures/mnm_regression/expected/expected_manifest.tsv index 97b5b972b..02afcbb5a 100644 --- a/tests/fixtures/mnm_regression/expected/expected_manifest.tsv +++ b/tests/fixtures/mnm_regression/expected/expected_manifest.tsv @@ -1,5 +1,5 @@ # expected-output regression manifest — mnm_regression (notebook tier) # columns: output mode rtol atol notes -univariate_bvsr.rds tolerant 1e-6 1e-8 susie_twas cell -> fine_mapping.R (QtlFineMappingResult for test_study ENSG00000283047 over qtl_mini) run with --seed 1 (the notebook `seed` parameter is forwarded to fine_mapping.R). Bit-identical across two seeded reruns; no embedded absolute paths. The x86/ARM CI matrix may need a looser rtol (calibrate from the first CI run). Re-baselined for PR #1422 (SuSiE L default 10, greedy-L off): the fit now uses a FIXED 10 single effects where the greedy search previously settled on 5, so alpha goes 5xP -> 10xP and the diffuse PIP mass scales with L. The locus yields no credible sets either way, so no conclusion changes. Greedy-off is what makes this platform-stable: with greedy on the retained count was 5 on macOS and 10 on Linux, a shape difference no tolerance could absorb. +univariate_bvsr.rds tolerant 1e-5 1e-8 susie_twas cell -> fine_mapping.R (QtlFineMappingResult for test_study ENSG00000283047 over qtl_mini) run with --seed 1 (the notebook `seed` parameter is forwarded to fine_mapping.R). Bit-identical across two seeded reruns; no embedded absolute paths. The x86/ARM CI matrix may need a looser rtol (calibrate from the first CI run). Re-baselined for PR #1422 (SuSiE L default 10, greedy-L off): the fit now uses a FIXED 10 single effects where the greedy search previously settled on 5, so alpha goes 5xP -> 10xP and the diffuse PIP mass scales with L. The locus yields no credible sets either way, so no conclusion changes. Greedy-off is what makes this platform-stable: with greedy on the retained count was 5 on macOS and 10 on Linux, a shape difference no tolerance could absorb. Re-baselined again for pecotmr 0.9.0: the susie_twas cell now passes --cv-folds/--cv-threads to the fine-mapping pass as well, because a cross-validated fine-mapping-derived weight method (susie) needs each fold's own fine-mapping fit -- 0.9.0 aborts otherwise (.twasCheckFmCvPresent). The only delta is that cvResult is now populated (samplePartition/prediction/performance/foldFits) where it was NULL; everything else is identical at 1e-6. rtol LOOSENED 1e-6 -> 1e-5 after the first CI matrix run on this fixture (the note above anticipated it): with cvResult now populated, each CV fold's SuSiE prior variance V (~6e-4) differs by 1.02e-6 and 1.04e-6 relative on ubuntu-24.04-arm while matching exactly on ubuntu-latest and macos-latest (actions run 38183809313). univariate_twas_weights.rds tolerant 1e-6 1e-8 susie_twas cell -> twas_weights.R (TwasWeights, same seed=1) for the same gene. Bit-identical across two seeded reruns; no embedded absolute paths. multicontext_bvsr.rds tolerant 1e-6 1e-8 mnm cell -> fine_mapping.R --methods mvsusie --joint-specification context (QtlFineMappingResult for test_study ENSG00000283047 over qtl_mini, contexts context1;context2) run with --seed 1. Byte-identical across two seeded reruns; no embedded absolute paths. diff --git a/tests/fixtures/mnm_regression/expected/univariate_bvsr.rds b/tests/fixtures/mnm_regression/expected/univariate_bvsr.rds index deaa90a11..3d2314c6b 100644 Binary files a/tests/fixtures/mnm_regression/expected/univariate_bvsr.rds and b/tests/fixtures/mnm_regression/expected/univariate_bvsr.rds differ diff --git a/tests/fixtures/qtl_dataset/expected/qtl_dataset.rds b/tests/fixtures/qtl_dataset/expected/qtl_dataset.rds index c90dfeaf4..8da43c950 100644 Binary files a/tests/fixtures/qtl_dataset/expected/qtl_dataset.rds and b/tests/fixtures/qtl_dataset/expected/qtl_dataset.rds differ diff --git a/tests/fixtures/qtl_rss_analysis/expected/expected_manifest.tsv b/tests/fixtures/qtl_rss_analysis/expected/expected_manifest.tsv index fa91b6dee..1a195ed7d 100644 --- a/tests/fixtures/qtl_rss_analysis/expected/expected_manifest.tsv +++ b/tests/fixtures/qtl_rss_analysis/expected/expected_manifest.tsv @@ -1,5 +1,5 @@ # expected-output regression manifest — qtl_rss_analysis (notebook tier) # columns: output mode rtol atol notes -qtl_sumstats.rds tolerant 1e-6 1e-8 generate_qtl_sumstats -> qtl_sumstats_construct.R: QtlSumStats for test_study/context1/ENSG00000283047 (200 variants of the gene's cis window, filtered out of the 16-gene TensorQTL nominal table by molecular_trait_id) with the qtl_mini genotypes attached as the ldSketch and summaryStatsQc applied. Deterministic (no RNG); byte-identical across two runs. normalize_paths: the ldSketch embeds the genotype path, basename-normalized by rds_compare.R. +qtl_sumstats.rds tolerant 1e-6 1e-8 generate_qtl_sumstats -> qtl_sumstats_construct.R: QtlSumStats for test_study/context1/ENSG00000283047 (200 variants of the gene's cis window, filtered out of the 16-gene TensorQTL nominal table by molecular_trait_id) with the qtl_mini genotypes attached as the ldSketch and summaryStatsQc applied. Deterministic (no RNG); byte-identical across two runs. normalize_paths: the ldSketch embeds the genotype path, basename-normalized by rds_compare.R. Re-baselined for pecotmr 0.9.0: the only delta is the qcInfo provenance key zMismatchQc -> ldMismatchQcMethod (same value); every other recorded option and data slot is identical. qtl_rss_finemap.rds tolerant 1e-4 1e-8 qtl_rss_fine_mapping -> fine_mapping.R --qtl-sumstats: QtlFineMappingResult from susieR::susie_rss over the same window, --seed 999. Byte-identical across two seeded runs on one machine; the x86/ARM CI matrix may need a looser rtol (in-sample LD from 49 samples is rank-deficient, so the RSS fit is closer to the tolerance edge than the individual-level one). Compared at rtol=1e-4, not 1e-6: the LD is in-sample from 49 samples, so a 200-variant panel has rank <= 48 and the solve sits near-singular, where log-Bayes-factor sums amplify last-bit BLAS differences. Measured macOS-vs-Linux drift 3.5e-6 (logBF) / 5.3e-6 (susieFit lbf_variable) -- bounded rounding, ~20x headroom, still four orders tighter than a real regression (the SuSiE L default change moved this comparison by 0.49 relative). qtl_rss_twas_weights.rds tolerant 1e-4 1e-8 qtl_rss_twas_weights -> twas_weights.R --qtl-sumstats: TwasWeights from the RSS lasso implementation (lassosumRssWeights) over the same collection, --seed 999. Byte-identical across two seeded runs. Compared at rtol=1e-4, not 1e-6: the LD is in-sample from 49 samples, so a 200-variant panel has rank <= 48 and the solve sits near-singular, where log-Bayes-factor sums amplify last-bit BLAS differences. Measured macOS-vs-Linux drift 3.5e-6 (logBF) / 5.3e-6 (susieFit lbf_variable) -- bounded rounding, ~20x headroom, still four orders tighter than a real regression (the SuSiE L default change moved this comparison by 0.49 relative). diff --git a/tests/fixtures/qtl_rss_analysis/expected/qtl_sumstats.rds b/tests/fixtures/qtl_rss_analysis/expected/qtl_sumstats.rds index b62559d37..1302b7826 100644 Binary files a/tests/fixtures/qtl_rss_analysis/expected/qtl_sumstats.rds and b/tests/fixtures/qtl_rss_analysis/expected/qtl_sumstats.rds differ diff --git a/tests/fixtures/rss_analysis/expected/expected_manifest.tsv b/tests/fixtures/rss_analysis/expected/expected_manifest.tsv index 22f2c4217..ab071ef7b 100644 --- a/tests/fixtures/rss_analysis/expected/expected_manifest.tsv +++ b/tests/fixtures/rss_analysis/expected/expected_manifest.tsv @@ -1,4 +1,4 @@ # expected-output regression manifest — rss_analysis (notebook tier) # columns: output mode rtol atol notes -gwas_sumstats.rds tolerant 1e-6 1e-8 generate_gwas_sumstats cell -> gwas_sumstats_construct.R (GwasSumStats for AD_Bellenguez_2022 chr22:49355984-50799822 over the toy chr22 GWAS + ld_reference panel). Deterministic QC/reshaping (no RNG); bit-identical run-to-run. normalize_paths: the ldSketch slot embeds the absolute ld_reference LD-panel path (resolved from --ld-meta at the checkout location), basename-normalized by rds_compare.R. Re-baselined for pecotmr 0.7.8 (#591 class cleanup): the genome build moved off the SumStatsBase `genome` slot into seqinfo(), so this fixture now records genome=GRCh38; values are unchanged (old-vs-new all.equal 1e-6 with the build neutralized). +gwas_sumstats.rds tolerant 1e-6 1e-8 generate_gwas_sumstats cell -> gwas_sumstats_construct.R (GwasSumStats for AD_Bellenguez_2022 chr22:49355984-50799822 over the toy chr22 GWAS + ld_reference panel). Deterministic QC/reshaping (no RNG); bit-identical run-to-run. normalize_paths: the ldSketch slot embeds the absolute ld_reference LD-panel path (resolved from --ld-meta at the checkout location), basename-normalized by rds_compare.R. Re-baselined for pecotmr 0.7.8 (#591 class cleanup): the genome build moved off the SumStatsBase `genome` slot into seqinfo(), so this fixture now records genome=GRCh38; values are unchanged (old-vs-new all.equal 1e-6 with the build neutralized). Re-baselined for pecotmr 0.9.0: the only delta is the qcInfo provenance key zMismatchQc -> ldMismatchQcMethod (same value); the other 23 recorded options and every data slot are identical. gwas_finemap.rds tolerant 1e-6 1e-8 gwas_fine_mapping cell -> fine_mapping.R (GwasFineMappingResult, susie). The notebook cell forwards no --seed but susie_rss is deterministic for this fixed input+init: bit-identical across two full pipeline reruns on this machine. Same embedded ld_reference path -> normalize_paths. The x86/ARM CI matrix may need a looser rtol (calibrate from the first CI run), mirroring fine_mapping/expected_manifest.tsv. Re-baselined for PR #1422 (SuSiE L default 10, greedy-L off): the fit now uses a FIXED 10 single effects where the greedy search previously settled on 5, so alpha goes 5xP -> 10xP and the diffuse PIP mass scales with L. The locus yields no credible sets either way, so no conclusion changes. Greedy-off is what makes this platform-stable: with greedy on the retained count was 5 on macOS and 10 on Linux, a shape difference no tolerance could absorb. diff --git a/tests/fixtures/rss_analysis/expected/gwas_sumstats.rds b/tests/fixtures/rss_analysis/expected/gwas_sumstats.rds index c26b6c497..92b8a5c8b 100644 Binary files a/tests/fixtures/rss_analysis/expected/gwas_sumstats.rds and b/tests/fixtures/rss_analysis/expected/gwas_sumstats.rds differ diff --git a/tests/fixtures/twas/expected/gwas_sumstats.chr22.rds b/tests/fixtures/twas/expected/gwas_sumstats.chr22.rds index 7eaee65af..9b082966e 100644 Binary files a/tests/fixtures/twas/expected/gwas_sumstats.chr22.rds and b/tests/fixtures/twas/expected/gwas_sumstats.chr22.rds differ diff --git a/tests/fixtures/twas/expected_manifest.tsv b/tests/fixtures/twas/expected_manifest.tsv index c34a1c27c..f3371785e 100644 --- a/tests/fixtures/twas/expected_manifest.tsv +++ b/tests/fixtures/twas/expected_manifest.tsv @@ -1,4 +1,4 @@ # expected-output regression manifest — twas # columns: output mode rtol atol notes -gwas_sumstats.chr22.rds tolerant 1e-6 1e-8 GwasSumStats S4 (one chr22 LD block) from gwas_sumstats_construct.R; deterministic. Compared with normalize_paths=True — the S4 embeds the machine-specific LD-reference path. Re-baselined for pecotmr 0.7.8 (#591 class cleanup): the genome build moved off the SumStatsBase `genome` slot into seqinfo(), so this fixture now records genome=GRCh38; values are unchanged (old-vs-new all.equal 1e-6 with the build neutralized). +gwas_sumstats.chr22.rds tolerant 1e-6 1e-8 GwasSumStats S4 (one chr22 LD block) from gwas_sumstats_construct.R; deterministic. Compared with normalize_paths=True — the S4 embeds the machine-specific LD-reference path. Re-baselined for pecotmr 0.7.8 (#591 class cleanup): the genome build moved off the SumStatsBase `genome` slot into seqinfo(), so this fixture now records genome=GRCh38; values are unchanged (old-vs-new all.equal 1e-6 with the build neutralized). Re-baselined again for pecotmr 0.9.0: the ONLY delta is the qcInfo provenance key zMismatchQc -> ldMismatchQcMethod (same value); the other 23 recorded options and every data slot are identical. Note 0.9.0 still records the FLAT option names in qcInfo$options even though the call now takes bundles, which is what makes that one-key diff the whole story. twas.chr22.rds tolerant 1e-6 1e-8 TWAS-Z GRanges from twas.R (causalInferencePipeline) over the committed S4 TwasWeights + a built GwasSumStats; deterministic, no embedded paths. diff --git a/tests/helpers/rds_compare.R b/tests/helpers/rds_compare.R index 5a624e504..284f724eb 100644 --- a/tests/helpers/rds_compare.R +++ b/tests/helpers/rds_compare.R @@ -48,6 +48,9 @@ norm_paths <- function(v) { # identical see the contents. Ordinary lists recurse; data.frames and atomics are # left as-is (all.equal already compares those element-wise with tolerance), except # under normPaths where character leaves / columns are basename-normalized. +# Environments open on the current canon() descent (see the closure branch). +.open_envs <- list() + canon <- function(x) { # Strip run-to-run-varying TIME metadata to a fixed sentinel (length preserved so a # genuine structural diff still surfaces). Wall-clock timings and timestamps are @@ -68,6 +71,16 @@ canon <- function(x) { if (is.null(e) || environmentName(e) != "" || isNamespace(e)) { return("") } + # A closure can be SELF-REFERENTIAL: pecotmr 0.9.0's ctwas LD_loader_fun + # keeps a copy of itself (`fn`) in its own frame, so descending into that + # frame reaches the same environment again and recurses until the node + # stack overflows. Collapse an environment already open on this path -- + # both sides hit it at the same place, so the comparison stays fair. + if (any(vapply(.open_envs, identical, logical(1), e))) { + return("") + } + .open_envs[[length(.open_envs) + 1L]] <<- e + on.exit(.open_envs[[length(.open_envs)]] <<- NULL, add = TRUE) return(list(.closure = canon(as.list(e, all.names = TRUE)))) } out <- if (isS4(x)) { diff --git a/tests/helpers/rds_probe.R b/tests/helpers/rds_probe.R index 9b49e2718..531171b78 100644 --- a/tests/helpers/rds_probe.R +++ b/tests/helpers/rds_probe.R @@ -4,8 +4,8 @@ # Emit a compact jsonlite structural summary of an RDS on stdout so the Python # test side can assert on a wrapper's output without an R testing stack. Reports # the S4 class, row count (for DFrame-subclass collections), and whichever -# pecotmr accessors resolve (method names / contexts / study). Any accessor that -# does not apply to the object is simply omitted. +# pecotmr accessors resolve (method names / contexts / study name). Any accessor +# that does not apply to the object is simply omitted. suppressWarnings(suppressMessages({ library(pecotmr) library(jsonlite) @@ -23,13 +23,18 @@ if (!is.null(nr)) summ$nrow <- as.integer(nr) # Optional accessors: include each only when it resolves without error and # returns something non-NULL. Keeps the probe generic across QtlDataset / # FineMappingResult / TwasWeights / CtwasResult / ... -for (acc in c("getStudy", "getContexts", "getMethodNames")) { - if (!exists(acc)) next +# pecotmr 0.9.0 dropped the `get` prefix from every accessor, so the reported +# key can no longer be derived from the function name -- it is named here, which +# also keeps the JSON keys stable for the Python assertions. No exists() guard: +# an accessor missing from the installed pecotmr is a real break and must fail. +accessors <- c(studyName = "Study", contexts = "Contexts", + methodNames = "MethodNames") +for (acc in names(accessors)) { v <- tryCatch(do.call(acc, list(obj)), error = function(e) NULL) if (!is.null(v) && length(v) > 0L) # I() keeps these as JSON arrays even for length 1, so the Python side can # assert uniformly (["susie"], not the auto-unboxed scalar "susie"). - summ[[sub("^get", "", acc)]] <- I(as.character(v)) + summ[[accessors[[acc]]]] <- I(as.character(v)) } # For a plain named list (e.g. assembled ctwas inputs), report the top-level @@ -47,8 +52,8 @@ if (is.character(cn) && length(cn) > 0L) summ$colnames <- I(cn) # CtwasResult: surface the top gene-level signal (max susie_pip + its z), so a # smoke test can assert numeric parity independent of the SNP-background rows. -if (methods::is(obj, "CtwasResult") && exists("getFinemap")) { - fm <- tryCatch(as.data.frame(getFinemap(obj)), error = function(e) NULL) +if (methods::is(obj, "CtwasResult")) { + fm <- tryCatch(as.data.frame(ctwasPosteriors(obj)), error = function(e) NULL) if (!is.null(fm) && "susie_pip" %in% names(fm) && nrow(fm) > 0L) { genes <- if ("type" %in% names(fm)) fm[fm$type != "SNP", , drop = FALSE] else fm if (nrow(genes) > 0L) { diff --git a/tests/helpers/wrapper_api_check.R b/tests/helpers/wrapper_api_check.R new file mode 100644 index 000000000..594e1a97b --- /dev/null +++ b/tests/helpers/wrapper_api_check.R @@ -0,0 +1,285 @@ +#!/usr/bin/env Rscript +# wrapper_api_check.R [ ...] +# +# Static API-contract check for the R wrapper scripts. Parses every wrapper +# WITHOUT running it and reports two classes of break: +# +# UNRESOLVED a function is called by a name that resolves nowhere -- not in +# the file, not in a file it source()s, not in any package it +# loads, not in base R. This is what an upstream rename looks +# like from here (pecotmr 0.9.0 dropping the `get` prefix turned +# every getTopLoci() call into one of these). +# BAD-ARG a named argument that is not in the callee's formals, for a +# pecotmr function whose signature has no `...` to absorb it. +# +# Exists because a runtime test only covers a wrapper that a test actually +# runs: most of the pecotmr_integration wrappers have no runtime test, and a +# wrapper that wraps its accessor call in tryCatch() can swallow a rename and +# still exit 0. This check sees all of them. +# +# Known blind spots (deliberate, documented rather than silently missed): +# * do.call(fn, args) -- the argument names are a runtime value, not syntax. +# * source() -- the helper file is chosen at run time, so for +# those callers any definition in the scanned tree is accepted. +# * The `...`-dispatching S4 generics (fineMappingPipeline, twasWeightsPipeline, +# colocboostPipeline, qtlAssociationPostprocess): findMethods() rewrites a +# method's formals to the generic's `(data, ...)`, so their real argument +# lists are only reachable through pecotmr internals. All four have runtime +# tests, which is what covers them. +suppressWarnings(suppressMessages(library(pecotmr))) + +args <- commandArgs(trailingOnly = TRUE) +roots <- if (length(args)) args else "code/script" + +ns <- asNamespace("pecotmr") +pecotmrExports <- getNamespaceExports(ns) +# Generics whose methods' formals are rewritten to the generic's signature. +DISPATCH_BLIND <- c("fineMappingPipeline", "twasWeightsPipeline", + "colocboostPipeline", "qtlAssociationPostprocess") + +isEmptyArg <- function(x) identical(x, quote(expr = )) + +# ---- AST walking ----------------------------------------------------------- + +# Apply `fn` to every call in an expression tree. +walkCalls <- function(e, fn) { + if (is.call(e)) { + fn(e) + l <- as.list(e) + for (i in seq_along(l)) { + if (isEmptyArg(l[[i]])) next + walkCalls(l[[i]], fn) + } + } else if (is.pairlist(e) || is.list(e)) { + l <- as.list(e) + for (i in seq_along(l)) { + if (isEmptyArg(l[[i]])) next + walkCalls(l[[i]], fn) + } + } + invisible(NULL) +} + +# Every name bound by assignment anywhere in the tree (`x <- `, `x = `, `x <<- `, +# assign("x", ...)) plus for-loop variables and function formals: the set of +# names that are NOT free variables at the point of use. +boundNames <- function(exprs) { + out <- character(0) + collect <- function(e) { + if (!is.call(e)) return(invisible()) + op <- as.character(e[[1]])[1] + if (op %in% c("<-", "=", "<<-") && length(e) >= 2L && + (is.name(e[[2]]) || is.character(e[[2]]))) + out <<- c(out, as.character(e[[2]])) + if (op == "assign" && length(e) >= 2L && is.character(e[[2]])) + out <<- c(out, e[[2]]) + if (op == "for" && length(e) >= 2L && is.name(e[[2]])) + out <<- c(out, as.character(e[[2]])) + if (op == "function" && is.pairlist(e[[2]])) + out <<- c(out, names(e[[2]])) + } + for (e in exprs) walkCalls(e, collect) + unique(out) +} + +# Packages the file attaches, plus every pkg used as pkg::fn. +packagesUsed <- function(exprs) { + out <- character(0) + collect <- function(e) { + if (!is.call(e)) return(invisible()) + op <- as.character(e[[1]])[1] + if (op %in% c("library", "require", "requireNamespace", "loadNamespace") && + length(e) >= 2L) { + a <- e[[2]] + if (is.name(a) || is.character(a)) out <<- c(out, as.character(a)) + } + if (op %in% c("::", ":::") && length(e) >= 2L) + out <<- c(out, as.character(e[[2]])) + } + for (e in exprs) walkCalls(e, collect) + unique(out) +} + +# TRUE when the file source()s a path it computes at run time (e.g. +# `source(argv$lib)`), which no static pass can resolve to a file. +sourcesDynamically <- function(exprs) { + hit <- FALSE + collect <- function(e) { + if (!is.call(e) || as.character(e[[1]])[1] != "source") return(invisible()) + a <- as.list(e)[-1] + if (length(a) == 0L) return(invisible()) + lits <- FALSE + walkCalls(e, function(x) { + for (v in as.list(x)[-1]) { + if (isEmptyArg(v)) next + if (is.character(v)) lits <<- TRUE + } + }) + if (!lits) hit <<- TRUE + } + for (e in exprs) walkCalls(e, collect) + hit +} + +# Sibling files this one source()s, resolved by basename against the same tree. +sourcedFiles <- function(exprs, self, allFiles) { + out <- character(0) + collect <- function(e) { + if (!is.call(e) || as.character(e[[1]])[1] != "source") return(invisible()) + lits <- character(0) + walkCalls(e, function(x) { + for (a in as.list(x)[-1]) { + if (isEmptyArg(a)) next + if (is.character(a)) lits <<- c(lits, a) + } + }) + for (s in lits) { + hit <- allFiles[basename(allFiles) == basename(s)] + out <<- c(out, setdiff(hit, self)) + } + } + for (e in exprs) walkCalls(e, collect) + unique(out) +} + +# library(pkg) also ATTACHES everything in pkg's Depends:, and those exports +# are visible to the script too -- library(GenomicRanges) is what puts +# IRanges() and seqnames() on the search path. Walk that closure. +# A meta-package attaches its members from .onAttach, which is invisible to +# the DESCRIPTION walk below; list what they put on the search path. +META_PACKAGES <- list( + tidyverse = c("dplyr", "tidyr", "readr", "purrr", "tibble", "stringr", + "forcats", "ggplot2", "lubridate") +) + +attachClosure <- function(pkgs) { + seen <- character(0) + todo <- unique(c(pkgs, unlist(META_PACKAGES[intersect(pkgs, names(META_PACKAGES))], + use.names = FALSE))) + while (length(todo)) { + pk <- todo[[1L]] + todo <- todo[-1L] + if (pk %in% seen) next + seen <- c(seen, pk) + dep <- tryCatch(utils::packageDescription(pk)$Depends, error = function(...) NULL) + if (is.null(dep) || is.na(dep)) next + dep <- trimws(strsplit(dep, ",")[[1L]]) + dep <- sub("\\s*\\(.*$", "", dep) + dep <- dep[nzchar(dep) & dep != "R"] + todo <- c(todo, dep) + } + seen +} + +# ---- the check ------------------------------------------------------------- + +files <- unique(unlist(lapply(roots, function(r) + list.files(r, pattern = "\\.R$", recursive = TRUE, full.names = TRUE)))) +files <- files[!grepl("/graveyard/", files)] +files <- sort(files) + +findings <- list() +note <- function(file, kind, detail) + findings[[length(findings) + 1L]] <<- + data.frame(file = file, kind = kind, detail = detail, stringsAsFactors = FALSE) + +parsed <- list() +for (f in files) { + p <- tryCatch(as.list(parse(f)), error = function(e) { + note(f, "PARSE", gsub("[\r\n]+", " ", conditionMessage(e))) + NULL + }) + if (!is.null(p)) parsed[[f]] <- p +} + +# A file that is source()d runs inside the sourcing script's session, so it +# sees whatever that script attached. pecotmr_aliases.R attaches nothing of its +# own and would otherwise look like it calls pecotmr out of nowhere. +sourcedBy <- list() +for (f in names(parsed)) { + for (sf in sourcedFiles(parsed[[f]], f, files)) { + sourcedBy[[sf]] <- unique(c(sourcedBy[[sf]], f)) + } +} + +for (f in names(parsed)) { + exprs <- parsed[[f]] + pkgs <- packagesUsed(exprs) + for (parent in sourcedBy[[f]]) { + pp <- parsed[[parent]] + if (!is.null(pp)) pkgs <- unique(c(pkgs, packagesUsed(pp))) + } + # Base + attached packages make up the visible search space. all.names + # reaches the dot-prefixed base objects (.lm.fit, .Machine, ...). + base_pkgs <- c("base", "stats", "utils", "methods", "graphics", "grDevices", + "tools", "parallel", "compiler", "datasets") + visible <- unique(unlist(lapply(base_pkgs, function(p) + tryCatch(ls(getNamespace(p), all.names = TRUE), error = function(...) character(0))))) + for (pk in attachClosure(pkgs)) { + nm <- tryCatch(getNamespaceExports(pk), error = function(...) character(0)) + visible <- c(visible, nm) + } + local_defs <- boundNames(exprs) + for (sf in sourcedFiles(exprs, f, files)) { + sp <- parsed[[sf]] + if (!is.null(sp)) local_defs <- c(local_defs, boundNames(sp)) + } + # A run-time source() path (source(argv$lib)) names a file only the caller + # knows, so accept anything the scanned tree defines rather than reporting + # the helper it pulls in. Loose by necessity, and only for these files. + if (sourcesDynamically(exprs)) { + for (sp in parsed) local_defs <- c(local_defs, boundNames(sp)) + } + known <- unique(c(visible, local_defs)) + + seenUnresolved <- character(0) + check <- function(e) { + fn <- e[[1]] + if (!is.name(fn)) return(invisible()) + nm <- as.character(fn) + # Operators and control flow are calls too; they are all in base. + if (!(nm %in% known) && !(nm %in% seenUnresolved)) { + seenUnresolved <<- c(seenUnresolved, nm) + note(f, "UNRESOLVED", paste0(nm, "()")) + return(invisible()) + } + if (!(nm %in% pecotmrExports) || nm %in% DISPATCH_BLIND) return(invisible()) + # A local binding of the same name shadows the pecotmr one: do not + # arg-check against a signature the script is not calling. + if (nm %in% local_defs) return(invisible()) + target <- get(nm, envir = ns) + fml <- if (is(target, "genericFunction")) { + unique(unlist(c(lapply(findMethods(nm, where = ns), function(m) names(formals(m))), + list(names(formals(args(target))))))) + } else { + names(formals(args(target))) + } + if ("..." %in% fml) return(invisible()) + argn <- names(as.list(e)) + if (length(argn) < 2L) return(invisible()) + argn <- argn[-1] + argn <- argn[!is.na(argn) & nzchar(argn)] + bad <- setdiff(argn, fml) + if (length(bad)) + note(f, "BAD-ARG", paste0(nm, "(): ", paste(bad, collapse = ", "))) + } + for (e in exprs) walkCalls(e, check) +} + +if (length(findings) == 0L) { + cat("wrapper API check: OK -", length(parsed), "scripts, no findings\n") + quit(save = "no", status = 0L) +} + +df <- unique(do.call(rbind, findings)) +df <- df[order(df$file, df$kind, df$detail), ] +cat("wrapper API check: ", nrow(df), " finding(s) in ", + length(unique(df$file)), " of ", length(files), " scripts\n\n", sep = "") +for (f in unique(df$file)) { + d <- df[df$file == f, ] + cat(sub("^\\./", "", f), "\n") + for (i in seq_len(nrow(d))) cat(sprintf(" %-11s %s\n", d$kind[i], d$detail[i])) +} +cat("\nUNRESOLVED = called by a name nothing provides (an upstream rename, or a typo).\n") +cat("BAD-ARG = named argument absent from the callee's formals.\n") +quit(save = "no", status = 1L) diff --git a/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py b/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py index bb6a5b4c9..f357c56e5 100644 --- a/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py +++ b/tests/notebooks/mnm_analysis/mnm_methods/test_mnm_regression.py @@ -50,8 +50,15 @@ def test_susie_twas(run_sos, read_rds, repo_root, qtl_mini, tmp_path): # closes the gap is unresolved. Checked for existence + class (above) only. # See memory: cross-platform-numeric-divergence. exp = repo_root / "tests/fixtures/mnm_regression/expected" + # rtol 1e-5, not the 1e-6 used elsewhere: the FMR now carries the CV fold + # fits (the susie_twas cell cross-validates the fine-mapping pass so the + # TWAS pass can reuse per-fold fits), and each fold's SuSiE prior variance + # V -- values of order 6e-4 -- lands ~1.0e-6 relative away on ubuntu-24.04-arm + # while matching on x86 Linux and macOS (run 38183809313). Last-bit BLAS + # divergence amplified by the EM, not a result change; 1e-5 keeps ~10x + # headroom over the observed 1.04e-6 and still pins the fit. assert_matches_expected(fmr, exp / "univariate_bvsr.rds", mode="tolerant", - rtol=1e-6, atol=1e-8) + rtol=1e-5, atol=1e-8) def test_mnm(run_sos, read_rds, repo_root, qtl_mini, tmp_path): diff --git a/tests/scripts/pecotmr_integration/test_fine_mapping_fullfit.py b/tests/scripts/pecotmr_integration/test_fine_mapping_fullfit.py index ad7bbf8c2..9aa9a8372 100644 --- a/tests/scripts/pecotmr_integration/test_fine_mapping_fullfit.py +++ b/tests/scripts/pecotmr_integration/test_fine_mapping_fullfit.py @@ -5,8 +5,9 @@ genuine fit. The committed enloc MiGA_eQTL fit is a real chr22 SuSiE fit with a credible set AND the full per-effect posterior matrices (alpha / mu / mu2 / lbf_variable / X_column_scale_factors). We reprocess it through pecotmr's -buildTopLoci(fullFit=TRUE, fullFitAlphaOnly=FALSE) and assert every wide per-CS -column lands with real (unscaled) values — an integration guard that the +buildTopLoci(credibleSetParam = CredibleSetParam(perCsColumns = "full")) -- the +0.9.0 spelling of the old fullFit=TRUE, fullFitAlphaOnly=FALSE pair -- and +assert every wide per-CS column lands with real (unscaled) values — an integration guard that the installed pecotmr the wrappers run against actually carries the feature. The reprocessing R is written to tmp at test time (not a committed wrapper): it @@ -37,7 +38,8 @@ attr(csTables, "coverage") <- c(0.95, 0.7, 0.5) tl <- as.data.frame(pecotmr:::buildTopLoci(fit, csTables, variantNames = vn, method = "susie", signalCutoff = 0, - fullFit = TRUE, fullFitAlphaOnly = FALSE)) + credibleSetParam = CredibleSetParam(perCsColumns = "full", + includeAllCs = FALSE))) write.table(tl, a[2], sep = "\t", quote = FALSE, row.names = FALSE, na = "") """ diff --git a/tests/scripts/pecotmr_integration/test_fine_mapping_pip_plot.py b/tests/scripts/pecotmr_integration/test_fine_mapping_pip_plot.py index 00c7edd55..a2a62293f 100644 --- a/tests/scripts/pecotmr_integration/test_fine_mapping_pip_plot.py +++ b/tests/scripts/pecotmr_integration/test_fine_mapping_pip_plot.py @@ -1,6 +1,16 @@ -"""Tier B: fine_mapping_pip_plot.R renders a non-empty PIP PNG from an FMR.""" +"""Tier B: fine_mapping_pip_plot.R renders a REAL PIP PNG from an FMR. + +Asserting on the wrapper's own summary line, not just a non-zero file size: the +script's fallback path writes a valid PNG carrying only the words "No usable PIP +vectors", so a size check passes on a blank plot. That is exactly what happened +when pecotmr 0.9.0 renamed getTopLoci() -> topLoci() -- the accessor error was +swallowed, every panel vanished, and this test stayed green. The "Wrote PIP +landscape" line is printed only on the real path. +""" from __future__ import annotations +import re + import pytest @@ -10,3 +20,7 @@ def test_fine_mapping_pip_plot(fmr, run_r, repo_root, tmp_path): ["--input", fmr, "--output", out], timeout=120) assert p.returncode == 0, p.stdout + p.stderr assert out.exists() and out.stat().st_size > 0 + assert "Wrote PIP landscape" in p.stdout, ( + f"stub plot instead of a real one:\n{p.stdout}\n{p.stderr}") + panels = int(re.search(r"Wrote PIP landscape \((\d+) panel", p.stdout).group(1)) + assert panels > 0, f"rendered {panels} panels" diff --git a/tests/scripts/pecotmr_integration/test_fine_mapping_upset.py b/tests/scripts/pecotmr_integration/test_fine_mapping_upset.py index a0fc58b12..4cb4c1627 100644 --- a/tests/scripts/pecotmr_integration/test_fine_mapping_upset.py +++ b/tests/scripts/pecotmr_integration/test_fine_mapping_upset.py @@ -1,5 +1,12 @@ """Tier B: fine_mapping_upset.R renders a valid PNG even for an FMR with no -credible sets (the empty-upset path).""" +credible sets (the empty-upset path). + +The toy FMR genuinely has zero credible sets, so the empty plot IS the expected +result -- but the test asserts the wrapper's reason for it rather than just the +file's existence. Without that, "empty because the input has no CS" is +indistinguishable from "empty because the accessor blew up and was swallowed", +which is how the pecotmr 0.9.0 getCs() -> credibleSets() rename slipped through. +""" from __future__ import annotations import pytest @@ -11,3 +18,5 @@ def test_fine_mapping_upset(fmr, run_r, repo_root, tmp_path): ["--input", fmr, "--output", out], timeout=120) assert p.returncode == 0, p.stdout + p.stderr assert out.exists() and out.stat().st_size > 0 + assert "No credible sets found across any input" in p.stderr, ( + f"expected the no-CS path for this fixture:\n{p.stdout}\n{p.stderr}") diff --git a/tests/scripts/test_wrapper_api_contract.py b/tests/scripts/test_wrapper_api_contract.py new file mode 100644 index 000000000..696374e4a --- /dev/null +++ b/tests/scripts/test_wrapper_api_contract.py @@ -0,0 +1,48 @@ +"""Guard: every pecotmr call in a wrapper must still resolve against the INSTALLED pecotmr. + +A runtime test only covers a wrapper some test actually runs, and most of the +``pecotmr_integration`` wrappers have none (the suite prints the untested list after +every run). Worse, a wrapper that wraps its accessor call in ``tryCatch`` can swallow an +upstream rename and still exit 0 — pecotmr 0.9.0 dropping the ``get`` prefix from every +accessor turned ``fine_mapping_pip_plot.R`` into a blank "No usable PIP vectors" stub +while its test stayed green. + +This check parses every wrapper WITHOUT running it, so it sees all of them, and fails on: + +* ``UNRESOLVED`` — a function called by a name nothing provides: not defined in the file, + not in a file it ``source()``s, not exported by a package it attaches (following the + ``Depends`` closure that ``library()`` attaches too), not in base R. An upstream rename + looks exactly like this. +* ``BAD-ARG`` — a named argument that is not in the callee's formals, for a pecotmr + function whose signature has no ``...`` to absorb it. + +SCOPE is ``code/script/pecotmr_integration`` — the wrappers that hold the pecotmr API +contract. The checker itself takes any root, so ``Rscript tests/helpers/wrapper_api_check.R +code/script`` sweeps the whole tree on demand (it reports real findings in the vendored +h3k9ac scripts, which nothing runs and which are out of scope here). + +Blind spots, documented rather than silently missed: ``do.call(fn, args)`` argument names +are a runtime value, not syntax; and the four ``...``-dispatching S4 generics +(``fineMappingPipeline``, ``twasWeightsPipeline``, ``colocboostPipeline``, +``qtlAssociationPostprocess``) have their methods' formals rewritten to the generic's +``(data, ...)``, so only their runtime tests cover their argument lists. +""" +from __future__ import annotations + +import subprocess + +from helpers.r_runner import rscript_bin + +SCOPE = "code/script/pecotmr_integration" + + +def test_wrappers_match_installed_pecotmr_api(repo_root): + r = subprocess.run( + [rscript_bin(), str(repo_root / "tests/helpers/wrapper_api_check.R"), SCOPE], + capture_output=True, text=True, cwd=repo_root, timeout=600) + # 0 = clean, 1 = findings printed on stdout; anything else = the checker broke + assert r.returncode in (0, 1), ( + f"wrapper_api_check.R errored:\n{r.stdout}\n{r.stderr}") + assert r.returncode == 0, ( + "wrapper(s) call a pecotmr API that the installed pecotmr does not provide:\n" + + r.stdout)