Export fSuSiE curves and bands in the canonical fine-mapping workflow

code/SoS/mnm_analysis/mnm_methods/mnm_regression.ipynb

@@ -1286,7 +1286,6 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "62f20d88",
    "metadata": {
     "kernel": "SoS"
    },
@@ -1310,12 +1309,17 @@
     "        --output ${_output}\n",
     "\n",
     "[fsusie_2]\n",
+    "# Empty preserves standalone behavior; Snakemake passes its chromosome scope.\n",
+    "parameter: chromosomes = []\n",
     "# Functional SuSiE per locus (fsusieR::susiF) over the pre-built QtlDataset.\n",
     "# fsusie-specific knobs ride on --method-args; --susie-top-pc>0 additionally\n",
     "# fine-maps each context's top principal components with univariate SuSiE\n",
     "# (usePCA; ports the legacy fsusie.R susie_on_top_pc). Replaces fsusie.R.\n",
     "parameter: prior = \"mixture_normal\"\n",
     "parameter: max_scale = 10\n",
+    "# Credible-set purity. pecotmr supplies susiF's `min_purity` from\n",
+    "# CredibleSetParam, the same knob every other method uses, so this rides\n",
+    "# on --min-abs-corr rather than the fsusie kwargs.\n",
     "parameter: min_purity = 0.5\n",
     "parameter: post_processing = \"TI\"\n",
     "parameter: max_SNP_EM = 100\n",
@@ -1329,13 +1333,16 @@
     "import csv, json\n",
     "manifest = path(f\"{cwd:a}/fsusie/{name}.fsusie_region_manifest.tsv\")\n",
     "jobs = list(csv.DictReader(open(manifest), delimiter='\\t'))\n",
+    "if chromosomes:\n",
+    "    allowed_chromosomes = {str(chrom).removeprefix(\"chr\") for chrom in chromosomes}\n",
+    "    jobs = [job for job in jobs if str(job[\"chr\"]).removeprefix(\"chr\") in allowed_chromosomes]\n",
     "stop_if(len(jobs) == 0, \"fsusie: empty region manifest; check --region-name / --customized-association-windows.\")\n",
     "qtl_dataset = path(f\"{cwd:a}/qtl_dataset/{name}.qtl_dataset.rds\")\n",
-    "fsusie_args = fsusie_method_args if fsusie_method_args else json.dumps({\"fsusie\": {\"prior\": prior, \"max_scale\": max_scale, \"min_purity\": min_purity, \"post_processing\": post_processing, \"max_SNP_EM\": max_SNP_EM, \"cor_small\": bool(small_sample_correction)}})\n",
+    "fsusie_args = fsusie_method_args if fsusie_method_args else json.dumps({\"fsusie\": {\"prior\": prior, \"max_scale\": max_scale, \"post_processing\": post_processing, \"max_SNP_EM\": max_SNP_EM, \"cor_small\": bool(small_sample_correction)}})\n",
     "input: qtl_dataset, for_each = \"jobs\"\n",
-    "output: f\"{cwd:a}/fsusie/{name}.{_jobs['region_id']}.fsusie.rds\"\n",
-    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output:bn}\"\n",
-    "bash: expand = \"${ }\", stderr = f\"{_output:n}.stderr\", stdout = f\"{_output:n}.stdout\", container = container, entrypoint = entrypoint\n",
+    "output: f\"{cwd:a}/fsusie/{name}.{_jobs['region_id']}.fsusie.rds\", f\"{cwd:a}/fsusie/{name}.{_jobs['region_id']}.fsusie.exported.bed.gz\"\n",
+    "task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f\"{step_name}_{_output[0]:bn}\"\n",
+    "bash: expand = \"${ }\", stderr = f\"{_output[0]:n}.stderr\", stdout = f\"{_output[0]:n}.stdout\", container = container, entrypoint = entrypoint\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping.R \\\n",
     "        --qtl-dataset ${_input} \\\n",
     "        --region ${_jobs['ld_block']} \\\n",
@@ -1344,12 +1351,36 @@
     "        ${(\"--secondary-coverage \" + \",\".join([str(x) for x in coverage[1:]])) if len(coverage) > 1 else \"\"} \\\n",
     "        --pip-cutoff ${pip_cutoff} \\\n",
     "        --pip-cutoff-to-skip ${pip_cutoff_to_skip} \\\n",
+    "        --L ${L} \\\n",
+    "        --min-abs-corr ${min_purity} \\\n",
     "        --seed ${seed} \\\n",
     "        ${(\"--contexts \" + contexts) if contexts else \"\"} \\\n",
     "        --method-args '${fsusie_args}' \\\n",
     "        ${(\"--use-pca --n-pcs \" + str(susie_top_pc)) if susie_top_pc > 0 else \"\"} \\\n",
-    "        --output ${_output}\n"
-   ]
+    "        --output ${_output[0]}\n"
+   ],
+   "id": "62f20d88"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "kernel": "SoS"
+   },
+   "outputs": [],
+   "source": [
+    "[fsusie_3]\n",
+    "# Assemble the same regional tables into one indexed context table.\n",
+    "input: output_from(\"fsusie_2\")\n",
+    "bed_inputs = paths([p for p in _input if str(p).endswith(\".exported.bed.gz\")])\n",
+    "output: f\"{cwd:a}/fsusie/{name}.exported.bed.gz\", f\"{cwd:a}/fsusie/{name}.exported.bed.gz.tbi\"\n",
+    "task: trunk_workers=1, trunk_size=1, mem=\"40G\", walltime=\"1h\", cores=1\n",
+    "bash: expand=\"${ }\", stderr=f\"{_output[0]}.stderr\", stdout=f\"{_output[0]}.stdout\", container=container, entrypoint=entrypoint\n",
+    "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_bed_merge.R \\\n",
+    "        --input ${bed_inputs:r} \\\n",
+    "        --output ${_output[0]:r}"
+   ],
+   "id": "fsusie-export"
   },
   {
    "cell_type": "markdown",

code/SoS/mnm_analysis/mnm_postprocessing.ipynb

@@ -1708,24 +1708,9 @@
     "output: f\"{cwd}/summary/{name}.{qtl_type}{_tag}.top_loci.bed.gz\"\n",
     "task: trunk_workers = 1, walltime = '2h', trunk_size = 1, mem = '8G', cores = 1, tags = f'{_output:bn}'\n",
     "bash: expand = \"${ }\", container = container, stderr = f'{_output:n}.stderr', stdout = f'{_output:n}.stdout', entrypoint=entrypoint\n",
-    "    set -euo pipefail\n",
-    "    mkdir -p \"$(dirname \"${_output}\")\"\n",
-    "    GZ=$(command -v bgzip || echo gzip)\n",
-    "    {\n",
-    "        # header from first non-empty file\n",
-    "        for f in ${_input}; do\n",
-    "            if [ \"$(zcat -f \"$f\" 2>/dev/null | wc -l)\" -ge 2 ]; then\n",
-    "                zcat -f \"$f\" | head -1\n",
-    "                break\n",
-    "            fi\n",
-    "        done\n",
-    "        # data rows from every non-empty file\n",
-    "        for f in ${_input}; do\n",
-    "            if [ \"$(zcat -f \"$f\" 2>/dev/null | wc -l)\" -ge 2 ]; then\n",
-    "                zcat -f \"$f\" | tail -n +2\n",
-    "            fi\n",
-    "        done\n",
-    "    } | $GZ > ${_output}\n"
+    "    Rscript ${modular_script_dir}/pecotmr_integration/fine_mapping_bed_merge.R \\\n",
+    "        --input ${_input:r} \\\n",
+    "        --output ${_output:r}"
    ]
   },
   {

code/script/data_preprocessing/phenotype/phenotype_formatting.R

@@ -11,6 +11,8 @@
 #   phenotype_by_region_2     — build the per-region file manifest
 #   gct_extract_samples       — filter samples from a GCT by a keep-list
 #   phenotype_annotate_by_tad — build a phenotype-per-TAD region list
+#   qtl_dataset_manifest      — single-context QtlDataset phenotype manifest
+#                               (+ region-id and sample-id side files)
 #   bam_subsetting            — subset a BAM to the requested regions
 #
 # Conventions: read with vroom, manipulate with dplyr, bgzip/tabix/bam via
@@ -59,6 +61,24 @@ parser <- add_argument(parser, "--TAD_list",
 parser <- add_argument(parser, "--phenotype-per-tad",
                        help = "[phenotype_annotate_by_tad] Min phenotypes per TAD to keep",
                        type = "numeric", default = 2)
+parser <- add_argument(parser, "--covariate-file",
+                       help = "[qtl_dataset_manifest] Covariate TSV (samples as columns)",
+                       type = "character", default = "")
+parser <- add_argument(parser, "--context",
+                       help = "[qtl_dataset_manifest] Context name for the `cond` column",
+                       type = "character", default = "")
+parser <- add_argument(parser, "--chromosomes", nargs = Inf,
+                       help = "[qtl_dataset_manifest] Chromosomes to keep, in output order",
+                       type = "character")
+parser <- add_argument(parser, "--phenotype-id-column",
+                       help = "[qtl_dataset_manifest] 4th BED column name", type = "character",
+                       default = "ID")
+parser <- add_argument(parser, "--output-region-ids",
+                       help = "[qtl_dataset_manifest] One feature id per line",
+                       type = "character", default = "")
+parser <- add_argument(parser, "--output-sample-ids",
+                       help = "[qtl_dataset_manifest] One shared sample id per line",
+                       type = "character", default = "")
 parser <- add_argument(parser, "--numThreads", help = "Number of threads",
                        type = "numeric", default = 1)
 argv <- parse_args(parser)
@@ -224,6 +244,57 @@ phenotype_annotate_by_tad <- function(argv) {
   message(sprintf("Written: %s (%d of %d TADs kept)", out_file, nrow(out), nrow(tad)))
 }
 
+# --- steps: QtlDataset manifest --------------------------------------------
+# One context -> one phenotype matrix, which is what
+# pecotmr::loadQtlDatasetFromManifest() resolves: every feature row points at
+# the same BED and the same covariate file. Emits the manifest plus the
+# feature-id and shared-sample-id lists the fan-out steps consume.
+qtl_dataset_manifest <- function(argv) {
+  need(argv, c("phenoFile", "covariate_file", "context", "chromosomes", "output",
+               "output_region_ids", "output_sample_ids"))
+  id_col <- argv$phenotype_id_column
+  # Absolute paths: the manifest is read from other working directories.
+  bed <- normalizePath(argv$phenoFile, mustWork = TRUE)
+  cov <- normalizePath(argv$covariate_file, mustWork = TRUE)
+
+  canonical <- function(x) tolower(sub("^chr", "", x, ignore.case = TRUE))
+  wanted <- canonical(argv$chromosomes)
+
+  header <- strsplit(read_header_line(bed), "\t", fixed = TRUE)[[1L]]
+  expected <- c("#chr", "start", "end", id_col)
+  if (!identical(header[1:4], expected))
+    stop(sprintf("%s must begin with %s", bed, paste(expected, collapse = ", ")))
+
+  pheno <- vroom::vroom(bed, delim = "\t", col_select = 1:4,
+                        col_types = cols(.default = "c"), show_col_types = FALSE)
+  names(pheno) <- c("chr", "start", "end", "ID")
+  rows <- pheno %>%
+    mutate(start = as.integer(.data$start), end = as.integer(.data$end),
+           .ord = match(canonical(.data$chr), wanted)) %>%
+    filter(!is.na(.data$.ord)) %>%
+    arrange(.data$.ord, .data$start, .data$end, .data$ID) %>%
+    transmute(chr, start, end, ID, path = bed, cond = argv$context, cov_path = cov)
+  if (nrow(rows) == 0L)
+    stop(sprintf("No features on the configured chromosomes: %s", bed))
+
+  # Samples carried by BOTH files; the covariate file names samples from its
+  # second column on (the first holds the covariate name).
+  covHeader <- strsplit(read_header_line(cov), "\t", fixed = TRUE)[[1L]]
+  samples <- sort(intersect(header[-(1:4)], covHeader[-1L]))
+  if (length(samples) == 0L)
+    stop("Phenotype BED and covariate file share no sample IDs")
+
+  out <- rows
+  names(out)[names(out) == "chr"] <- "#chr"
+  for (f in c(argv$output, argv$output_region_ids, argv$output_sample_ids))
+    dir.create(dirname(f), showWarnings = FALSE, recursive = TRUE)
+  vroom::vroom_write(out, argv$output, delim = "\t")
+  writeLines(rows$ID, argv$output_region_ids)
+  writeLines(samples, argv$output_sample_ids)
+  message(sprintf("Written: %s (%d features, %d shared samples)",
+                  argv$output, nrow(rows), length(samples)))
+}
+
 # --- steps: BAM subsetting -------------------------------------------------
 bam_subsetting <- function(argv) {
   need(argv, c("phenoFile", "region", "output"))
@@ -252,6 +323,7 @@ switch(argv$step,
   phenotype_by_region_2     = phenotype_by_region_2(argv),
   gct_extract_samples       = gct_extract_samples(argv),
   phenotype_annotate_by_tad = phenotype_annotate_by_tad(argv),
+  qtl_dataset_manifest      = qtl_dataset_manifest(argv),
   bam_subsetting            = bam_subsetting(argv),
   stop(sprintf("Unknown step: '%s'", argv$step))
 )

code/script/pecotmr_integration/fine_mapping.R

@@ -454,7 +454,7 @@ if (has_gwas || has_qss) {
   label <- if (has_region) paste0("region '", argv$region, "'")
            else paste0("gene '", argv$gene_id, "'")
   run_fm <- function() if (has_region) {
-    do.call(fineMappingPipeline, c(qtl_args, list(region = argv$region)))
+    do.call(fineMappingPipeline, c(qtl_args, list(region = GenomicRanges::GRanges(argv$region))))
   } else {
     do.call(fineMappingPipeline, c(qtl_args,
                                    list(traitId = argv$gene_id,
@@ -500,11 +500,12 @@ if (has_gwas && !is.null(res)) {
       # 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
+      # credibleSetSummary() IS the per-credible-set table, one row each, so
+      # the count is its nrow. Counting distinct cs_95 labels by hand got this
+      # wrong: the "not captured" sentinel is "<method>_0", not a bare 0, so
+      # setdiff(..., 0) never removed it and every row reported one credible
+      # set too many (1 where there were none).
+      ncs <- nrow(as.data.frame(credibleSetSummary(row_i, coverage = 0.95)))
       fit <- susieFit(row_i)
       ser <- if (!is.null(fit$serFallbackUsed))
                isTRUE(fit$serFallbackUsed) else NA
@@ -520,3 +521,11 @@ if (has_gwas && !is.null(res)) {
      warning("fine_mapping.R: GWAS summary skipped (", conditionMessage(e), ")",
              call. = FALSE)))
 }
+
+if (has_qtl && has_region && "fsusie" %in% methods) {
+    # The published fSuSiE table (schema, packing, sort, bgzip + tabix) is
+    # pecotmr's exportToBed(); region_ID / region_* come off the result's
+    # own fitRegion() provenance, so --region is not passed again here.
+    exportToBed(res, sub("\\.rds$", ".exported.bed.gz", argv$output),
+                schema = "fsusie")
+}

code/script/pecotmr_integration/fine_mapping_bed_merge.R

@@ -0,0 +1,80 @@
+#!/usr/bin/env Rscript
+# fine_mapping_bed_merge.R
+#
+# Concatenate the per-region published tables written by fine_mapping.R /
+# fine_mapping_top_loci_bed.R (pecotmr::exportToBed) into one sorted,
+# compressed table for a study or context.
+#
+# Reads the per-region tables rather than the per-region FineMappingResult
+# RDSes on purpose: exportToBed() can fold many results into one call,
+# but that holds every fit in memory at once, and a genome-wide run is
+# thousands of regions. The tables are small and already carry the published
+# schema, so merging them keeps the peak footprint flat.
+#
+# Every table pecotmr publishes is a half-open BED now -- #chr / start / end
+# for the top-loci schema as much as for the fSuSiE one, which is the point of
+# the two sharing one vocabulary -- so both are tabix-indexable. The header is
+# still sniffed, because a hand-assembled position table (chr / pos / ...) can
+# reach this step from outside the pipeline, and indexing one with the tabix
+# "bed" preset would read its third column, an allele, as an end coordinate.
+#
+# Inputs:
+#   --input <TABLE> [...]  Per-region tables (all sharing one header).
+#   --output <PATH>        Merged output; bgzipped when it ends in .gz/.bgz.
+
+suppressPackageStartupMessages({
+  library(argparser)
+  library(vroom)
+  library(dplyr)
+  library(Rsamtools)
+  library(pecotmr)
+})
+
+parser <- arg_parser("Merge per-region fine-mapping tables into one sorted table")
+parser <- add_argument(parser, "--input", help = "Per-region table paths",
+                       type = "character", nargs = Inf)
+parser <- add_argument(parser, "--output", help = "Merged output path",
+                       type = "character")
+argv <- parse_args(parser)
+
+inputs <- as.character(argv$input)
+if (length(inputs) == 0L) stop("--input requires at least one table path.")
+
+# col_types = "c" throughout: the packed ";"-joined columns must survive
+# untouched and the table is re-emitted verbatim. vroom reads the shared
+# header once and skips it in the remaining files.
+merged <- vroom(inputs, delim = "\t", col_types = vroom::cols(.default = "c"),
+                progress = FALSE)
+
+cols <- names(merged)
+# A half-open BED starts #chr/start/end; anything else came from elsewhere and
+# is treated as a position table.
+is_bed <- length(cols) >= 3L &&
+  grepl("^#?chr", cols[[1L]], ignore.case = TRUE) &&
+  identical(tolower(cols[[2L]]), "start") &&
+  identical(tolower(cols[[3L]]), "end")
+
+if (nrow(merged) > 0L) {
+  # genomicOrder() (pecotmr) is the single definition of genomic order, shared
+  # with the writers: chromosomes naturally (2 before 10), positions
+  # numerically. `end` only breaks ties, and is ignored when it is not a
+  # position -- which is what makes one call serve both schemas.
+  merged <- merged[genomicOrder(merged[[1L]], merged[[2L]],
+                                if (is_bed) merged[[3L]] else NULL), ,
+                   drop = FALSE]
+}
+
+dir.create(dirname(argv$output), showWarnings = FALSE, recursive = TRUE)
+compress <- grepl("\\.(gz|bgz)$", argv$output)
+plain <- if (compress) sub("\\.(gz|bgz)$", "", argv$output) else argv$output
+vroom_write(merged, plain, delim = "\t", quote = "none", na = "NA",
+            progress = FALSE)
+if (compress) {
+  bgzip(plain, dest = argv$output, overwrite = TRUE)
+  file.remove(plain)
+  # Only the real BED is tabix-indexable, and an empty table has no records.
+  if (is_bed && nrow(merged) > 0L) indexTabix(argv$output, format = "bed")
+}
+cat(sprintf("Merged %d table(s), %d rows, into %s%s\n",
+            length(inputs), nrow(merged), argv$output,
+            if (is_bed && compress && nrow(merged) > 0L) " (tabix-indexed)" else ""))

code/script/pecotmr_integration/fine_mapping_cis_db_export.R

@@ -7,7 +7,7 @@
 # purity), so per region this wrapper only:
 #   1. combines the per-study/context FineMappingResult RDS(es) into the
 #      cis_results_db (itself an FMR) via combineFineMappingResults, and saves it;
-#   2. saves the marginal sumstats table (getMarginalEffects) as the
+#   2. saves the marginal sumstats table (marginalEffects) as the
 #      combined_data_sumstats db;
 #   3. emits the enrichment meta row SuSiE_enloc reads and the pip_sum table.
 #