fixes for new pecotmr release

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",

code/script/mnm_analysis/mnm_methods/fsusie.R → 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")
 }
 

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

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

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

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

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)

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

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)

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:",

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")
 }