fix some tests after pecotmr update

code/SoS/pecotmr_integration/SuSiE_enloc.ipynb

@@ -389,6 +389,9 @@
     "  get_overlapped_analysis_regions: get overlapped regions from flatten table, to get the regions with overlapped variants and select those\n",
     "                        regions for coloc analysis\n",
     "  xqtl_gwas_enrichment:\n",
+    "    Workflow Options:\n",
+    "      --seed -1 (as int)\n",
+    "                        Reproducibility: integer RNG seed for the multiple-imputation sampler; -1 = unset (nondeterministic)\n",
     "  susie_coloc:\n",
     "    Workflow Options:\n",
     "      --[no-]skip-enrich (default to False)\n",
@@ -552,6 +555,11 @@
    "source": [
     "[xqtl_gwas_enrichment]\n",
     "depends: sos_variable(\"enrichment_rows\")\n",
+    "# Reproducibility: integer RNG seed for the enrichment multiple-imputation\n",
+    "# sampler, which draws in C++ and so is unreachable from R's set.seed; -1 = unset\n",
+    "# (nondeterministic).\n",
+    "parameter: seed = -1\n",
+    "seed_arg = f\"--seed {seed}\" if seed >= 0 else \"\"\n",
     "jobs = enrichment_rows\n",
     "stop_if(len(jobs) == 0, f'No files left for analysis')\n",
     "_qtl_groups = [j['qtl_files'].split(',') for j in jobs]\n",
@@ -567,6 +575,7 @@
     "        --qtl-fine-mapping ${\" \".join([str(x) for x in _input])} \\\n",
     "        --gwas-fine-mapping ${_jobs[\"gwas_files\"].replace(\",\", \" \")} \\\n",
     "        --ncore ${numThreads} \\\n",
+    "        ${seed_arg} \\\n",
     "        --output ${_output}\n"
    ]
   },

code/script/molecular_phenotypes/calling/RNA_calling.sh

@@ -661,7 +661,11 @@ _picard_qc_star_align_2() {
     local metrics_prefix ri_file
     local -a java_args=()
     metrics_prefix="${PICARD_METRICS%.*}"
-    ri_file="${INPUT_BAM}.RI"
+    # Derived intermediate (ribosomal interval list), consumed only by
+    # CollectRnaSeqMetrics below. It MUST live beside the picard OUTPUTS, not beside
+    # $INPUT_BAM -- writing next to the input mutates a read-only input tree (and, in
+    # the test suite, dirtied the committed BAM fixture directory on every run).
+    ri_file="${metrics_prefix}.RI"
     # portable mapfile: macOS ships bash 3.2 (no mapfile). picard_java_args emits
     # 0 or 1 lines (-Xmx...), so read them into the array with a while loop.
     while IFS= read -r _ja_line; do java_args+=("$_ja_line"); done < <(picard_java_args)

code/script/pecotmr_integration/fine_mapping.R

@@ -243,7 +243,14 @@ median_abs_corr <- if (length(argv$median_abs_corr) != 1L || is.na(argv$median_a
                      NULL else argv$median_abs_corr
 
 # Seed up front for reproducible fits (mirrors the legacy susie_twas set.seed).
-if (length(argv$seed) == 1L && !is.na(argv$seed)) set.seed(as.integer(argv$seed))
+# NULL when --seed is unset. In QTL mode the value is ALSO forwarded to the
+# pipeline as seed= (below): set.seed() only touches this process's RNG, so it
+# cannot reach the BiocParallel workers that fit the cross-validation folds --
+# only the forwarded seed makes a --cv-threads > 1 run reproducible. The
+# GwasSumStats method takes no seed= argument, so GWAS mode relies on set.seed().
+seed_val <- if (length(argv$seed) == 1L && !is.na(argv$seed))
+              as.integer(argv$seed) else NULL
+if (!is.null(seed_val)) set.seed(seed_val)
 
 has_qtl  <- nzchar(argv$qtl_dataset)
 has_gwas <- nzchar(argv$gwas_sumstats)
@@ -366,6 +373,7 @@ if (has_gwas) {
     qtl_args$usePCA <- TRUE
     qtl_args$nPCs   <- argv$n_pcs
   }
+  if (!is.null(seed_val)) qtl_args$seed <- seed_val
   label <- if (has_region) paste0("region '", argv$region, "'")
            else paste0("gene '", argv$gene_id, "'")
   run_fm <- function() if (has_region) {

code/script/pecotmr_integration/qtl_enrichment.R

@@ -4,7 +4,7 @@
 # xQTL-GWAS enrichment worker. Loads an S4 `QtlFineMappingResult` and
 # an S4 `GwasFineMappingResult`, calls
 # `pecotmr::qtlEnrichmentPipeline()`, and saves the resulting
-# per-(gwasStudy, qtlContext) enrichment data.frame. The output feeds
+# per-(gwasStudy, qtlContext) enrichment tibble. The output feeds
 # downstream into `coloc.R` (via colocPipeline's `enrichment` arg) as
 # the prior-adjustment factor for colocalization.
 #
@@ -17,7 +17,9 @@
 #   --lambda               Pass-through (default 1)
 #   --imp-n                Pass-through (default 25)
 #   --ncore                Number of threads (default 1)
-#   --output               Output RDS path (per-pair enrichment data.frame)
+#   --seed                 Integer RNG seed for the multiple-imputation sampler
+#                          (reproducibility); unset = nondeterministic
+#   --output               Output RDS path (per-pair enrichment tibble)
 
 suppressPackageStartupMessages({
   library(argparser)
@@ -46,6 +48,10 @@ parser <- add_argument(parser, "--imp-n",
 parser <- add_argument(parser, "--ncore",
                        help = "Pass-through ncore",
                        type = "integer", default = 1L)
+parser <- add_argument(parser, "--seed",
+                       help = paste("Integer RNG seed for the multiple-imputation sampler",
+                                    "(reproducibility); unset = nondeterministic"),
+                       type = "integer", default = NA)
 parser <- add_argument(parser, "--output",
                        help = "Output RDS path", type = "character")
 argv <- parse_args(parser)
@@ -64,14 +70,24 @@ argv <- parse_args(parser)
 qtlFmr  <- .loadCombine(argv$qtl_fine_mapping)
 gwasFmr <- .loadCombine(argv$gwas_fine_mapping)
 
-res <- qtlEnrichmentPipeline(
+# The impN imputation rounds are drawn in C++ (qtl_enrichment.h), so a
+# main-process set.seed() cannot reach them -- the base seed has to be forwarded,
+# and each round derives its own seed from it so the result is independent of
+# which OpenMP thread runs which round. Unset --seed keeps the historical
+# nondeterministic std::random_device behaviour. Added only when supplied, for
+# 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
+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,
   numThreads            = argv$ncore)
+if (!is.null(seed_val)) enr_args$seed <- seed_val
+res <- do.call(qtlEnrichmentPipeline, enr_args)
 
 dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE)
 saveRDS(res, argv$output)

code/script/pecotmr_integration/twas_weights.R

@@ -146,7 +146,13 @@ qd_filter_overrides <- function(argv) {
 }
 
 # Seed up front for reproducible fits (mirrors the legacy susie_twas set.seed).
-if (length(argv$seed) == 1L && !is.na(argv$seed)) set.seed(as.integer(argv$seed))
+# NULL when --seed is unset. The value is ALSO forwarded to the pipeline as
+# seed= (below): set.seed() only touches this process's RNG, so it cannot reach
+# the BiocParallel workers that fit the per-method weights and the CV folds --
+# only the forwarded seed makes a --cv-threads > 1 run reproducible.
+seed_val <- if (length(argv$seed) == 1L && !is.na(argv$seed))
+              as.integer(argv$seed) else NULL
+if (!is.null(seed_val)) set.seed(seed_val)
 
 # Parse --method-args into a nested named list of per-method kwargs.
 parsed_method_args <- if (nzchar(argv$method_args) && argv$method_args != "." &&
@@ -248,6 +254,7 @@ tw_args <- c(list(methods           = methods_arg,
 # the call compatible with a pecotmr that predates the argument).
 if (!is.null(joint_spec))       tw_args$jointSpecification <- joint_spec
 if (!is.null(twas_weights_obj)) tw_args$twasWeights        <- twas_weights_obj
+if (!is.null(seed_val))         tw_args$seed               <- seed_val
 res <- if (has_region) {
   do.call(twasWeightsPipeline,
           c(list(qd), tw_args, list(region = argv$region)))