From 2e1c3e15f350fcb423db1800050da21fceed4898 Mon Sep 17 00:00:00 2001
From: Jenny Anak Empawi <jennyempawi95@gmail.com>
Date: Wed, 21 Oct 2026 09:06:10 -0400
Subject: [PATCH 1/5] Handle weights stored for the QC variant subset and genes
 with no usable weights in twas_weights_conversion

---
 .../twas_weights_conversion.R                  | 18 ++++++++++++++++--
 1 file changed, 16 insertions(+), 2 deletions(-)

diff --git a/code/script/pecotmr_integration/twas_weights_conversion.R b/code/script/pecotmr_integration/twas_weights_conversion.R
index adaaea91c..d21da3f36 100644
--- a/code/script/pecotmr_integration/twas_weights_conversion.R
+++ b/code/script/pecotmr_integration/twas_weights_conversion.R
@@ -53,7 +53,11 @@ for (f in files) {
     for (wn in names(co$twas_weights)) {
       w <- co$twas_weights[[wn]]
       vids <- if (is.matrix(w)) rownames(w) else names(w)
-      if (is.null(vids)) vids <- co$variant_names
+      # Some releases store weights for the QC'd variant subset only, without row
+    # names; their IDs are then in preset_variants_result$variant_names
+    if (is.null(vids)) vids <- Find(function(v) length(v) == NROW(w),
+      list(co$variant_names, co$preset_variants_result$variant_names))
+    if (is.null(vids)) stop("No variant IDs of length ", NROW(w), " for ", cn, " ", wn, " in ", f)
       wv <- if (is.matrix(w)) as.numeric(w[, 1]) else as.numeric(w)
       if (all(is.na(wv)) || all(wv == 0, na.rm = TRUE)) next
       m <- sub("_weights$", "", wn)
@@ -67,7 +71,17 @@ for (f in files) {
     ctxs <- union(ctxs, ctx)
   }
 }
-if (!length(E)) stop("No usable (non-zero) weights in: ", a$legacy)
+if (!length(E)) {
+  # e.g. a gene stored as NULL or with all-zero weights: write an empty result
+  # and empty meta rows, so step 2 leaves the gene out instead of failing
+  msg <- paste("No usable (non-zero) weights in:", a$legacy)
+  warning(msg)
+  dir.create(dirname(a$output), showWarnings = FALSE, recursive = TRUE)
+  saveRDS(structure(list(), noResult = msg), a$output)
+  file.create(a$meta_row)
+  if (nzchar(a$type_rows)) file.create(a$type_rows)
+  quit(save = "no", status = 0)
+}
 
 tw <- TwasWeights(studyName = S, context = C, trait = Tr, method = M, entry = E,
                   traitPos = GenomicRanges::GRanges(TPc, IRanges::IRanges(TPp, width = 1L)))

From bd4b5a875623432ed76796121dda6c664145d7a7 Mon Sep 17 00:00:00 2001
From: Jenny Anak Empawi <jennyempawi95@gmail.com>
Date: Wed, 28 Oct 2026 10:36:24 -0400
Subject: [PATCH 2/5] Record noResult for genes in LD blocks with no GWAS
 variants

gwas_sumstats_construct.R writes no GwasSumStats file for a block without GWAS variants, so twas.R failed on readRDS and the whole region errored. Such genes now get an empty result with a noResult attribute.
---
 code/SoS/pecotmr_integration/twas_ctwas.ipynb | 15 +++-
 code/script/pecotmr_integration/twas.R        |  7 +-
 .../pecotmr_integration/twas_manifest.R       | 69 ++++++++++++++++---
 3 files changed, 79 insertions(+), 12 deletions(-)

diff --git a/code/SoS/pecotmr_integration/twas_ctwas.ipynb b/code/SoS/pecotmr_integration/twas_ctwas.ipynb
index 41f235e4f..2b1601a7e 100644
--- a/code/SoS/pecotmr_integration/twas_ctwas.ipynb
+++ b/code/SoS/pecotmr_integration/twas_ctwas.ipynb
@@ -976,6 +976,17 @@
     "# (e.g. sos run twas_ctwas.ipynb twas). Gene -> home-LD-block placement is NOT\n",
     "# done here: pecotmr's assembleCtwasInputs places each gene into its block from\n",
     "# the weight `region` provenance.\n",
+    "# TWAS analysis unit: gene (default) gives each gene its own window, the\n",
+    "# span of its weights from the xQTL meta, so GWAS and LD are read per gene\n",
+    "# and no LD blocks are involved; --regions / --region-name then only filter\n",
+    "# genes by TSS. region keeps the legacy grouping of genes by --regions.\n",
+    "# LD blocks enter only at the ctwas step.\n",
+    "parameter: twas_unit = \"gene\"\n",
+    "# Gene unit only: skip genes with no context/method passing the CV cutoffs\n",
+    "# (same values as the twas step), so their GWAS input is never built.\n",
+    "parameter: twas_prefilter = True\n",
+    "parameter: rsq_cutoff = 0.01\n",
+    "parameter: rsq_pval_cutoff = 0.05\n",
     "input: None\n",
     "twas_manifest_file = f\"{cwd:a}/get_analysis_regions/{name}.twas_manifest.tsv\"\n",
     "bash: expand = \"${ }\", container = container\n",
@@ -988,7 +999,9 @@
     "        ${(\"--gwas-name \" + \" \".join([str(x) for x in gwas_name])) if len(gwas_name) > 0 else \"\"} \\\n",
     "        ${(\"--gwas-data \" + \" \".join([str(x) for x in gwas_data])) if len(gwas_data) > 0 else \"\"} \\\n",
     "        ${(\"--column-mapping \" + \" \".join([str(x) for x in column_mapping])) if len(column_mapping) > 0 else \"\"} \\\n",
-    "        --output ${twas_manifest_file}"
+    "        --unit ${twas_unit} \\\n",
+    "    ${(\"--min-cv-rsq %s --max-cv-pval %s --threads %s\" % (rsq_cutoff, \"Inf\" if rsq_pval_cutoff == float(\"inf\") else rsq_pval_cutoff, numThreads)) if (twas_prefilter and twas_unit == \"gene\") else \"\"} \\\n",
+    "    --output ${twas_manifest_file}"
    ]
   },
   {
diff --git a/code/script/pecotmr_integration/twas.R b/code/script/pecotmr_integration/twas.R
index 42489f2b7..5b9a418ba 100644
--- a/code/script/pecotmr_integration/twas.R
+++ b/code/script/pecotmr_integration/twas.R
@@ -112,7 +112,9 @@ if (length(fmPaths) != length(twPaths)) stop("--fine-mapping-result needs one pa
 # this object's ldReference, per gene and scoped to that gene's matched
 # variants (and for a genotype panel it uses the LD factor, so the p x p matrix
 # is never formed). Nothing here computes or caches LD.
-gss <- readRDS(argv$gwas_sumstats)
+# gwas_sumstats_construct.R writes no GwasSumStats file for an LD block with
+# no GWAS variants; every gene of such a block gets an empty result.
+gss <- if (file.exists(argv$gwas_sumstats)) readRDS(argv$gwas_sumstats) else NULL
 
 readFmr <- function(p) if (nzchar(p) && p != "." && file.exists(p)) readRDS(p) else NULL
 
@@ -122,7 +124,8 @@ readFmr <- function(p) if (nzchar(p) && p != "." && file.exists(p)) readRDS(p) e
 # A gene with no scorable (QTL, GWAS) pair gets an empty list with a `noResult`
 # attribute, so the other genes of the block are still scored.
 for (i in seq_along(twPaths)) {
-  res <- tryCatch(causalInferencePipeline(
+  res <- if (is.null(gss)) structure(list(), noResult = "no GWAS variants in this LD block") else
+    tryCatch(causalInferencePipeline(
     gwasSumStats = gss,
     twasWeights = readRDS(twPaths[i]),
     fineMappingResult = readFmr(fmPaths[i]),
diff --git a/code/script/pecotmr_integration/twas_manifest.R b/code/script/pecotmr_integration/twas_manifest.R
index 9e5514ada..1c5ce0ff1 100644
--- a/code/script/pecotmr_integration/twas_manifest.R
+++ b/code/script/pecotmr_integration/twas_manifest.R
@@ -62,6 +62,14 @@ p <- add_argument(p, "--gwas-data", type = "character", nargs = Inf,
                   default = NA_character_, help = "Optional inline GWAS sumstats paths")
 p <- add_argument(p, "--column-mapping", type = "character", nargs = Inf,
                   default = NA_character_, help = "Optional inline GWAS column-mapping paths")
+p <- add_argument(p, "--unit", type = "character", default = "gene",
+  help = "Analysis unit: 'gene' (default; one unit per gene spanning its weight window, no LD blocks) or 'region' (legacy: genes grouped by --regions)")
+p <- add_argument(p, "--min-cv-rsq", type = "numeric", default = NA,
+  help = "gene unit only: skip genes where no context/method reaches this CV r2 (with --max-cv-pval); NA = no prefilter")
+p <- add_argument(p, "--max-cv-pval", type = "numeric", default = NA,
+  help = "gene unit only: CV p-value that must also be passed (Inf = no p-value check)")
+p <- add_argument(p, "--threads", type = "integer", default = 1L,
+  help = "Cores for reading weight files in the prefilter")
 p <- add_argument(p, "--output", type = "character", help = "Output manifest TSV path")
 argv <- parse_args(p)
 
@@ -139,8 +147,10 @@ if (length(rn) > 0L) {
     stringsAsFactors = FALSE)
 }
 regions <- unique(regions)
-if (nrow(regions) == 0L)
-  stop("No analysis regions: provide --regions and/or --region-name.")
+if (!argv$unit %in% c("gene", "region"))
+  stop("--unit must be 'gene' or 'region', got: ", argv$unit)
+if (argv$unit == "region" && nrow(regions) == 0L)
+  stop("--unit region needs --regions and/or --region-name.")
 
 # ---- xQTL weight metadata --------------------------------------------------
 if (is.null(argv$xqtl_meta) || !nzchar(argv$xqtl_meta) || !file.exists(argv$xqtl_meta))
@@ -152,12 +162,53 @@ xqtlChr <- chromAdd(xqtlDf[["#chr"]])
 xqtlTss <- suppressWarnings(as.integer(xqtlDf$TSS))
 
 # ---- per-region: overlapping genes/weights + covering GWAS -----------------
+# Analysis units. gene (default): one unit per xQTL meta row, spanning that
+# gene's weight window (meta start..end), so every weight variant is scored
+# with GWAS and LD read from the per-chromosome files; --regions /
+# --region-name, if given, only filter genes by TSS. region (legacy): genes
+# binned into --regions by TSS, one unit per region.
+inRegions <- function(j) nrow(regions) == 0L ||
+  any(regions$chr == xqtlChr[[j]] & !is.na(xqtlTss[[j]]) &
+      regions$start <= xqtlTss[[j]] & regions$stop >= xqtlTss[[j]])
+units <- if (argv$unit == "gene") {
+  keep <- Filter(inRegions, seq_len(nrow(xqtlDf)))
+  # Optional prefilter: drop genes that can have no imputable model, so the
+  # twas step does not build their GWAS input. Lenient on purpose: a gene is
+  # kept if ANY context/method row passes, a superset of what pecotmr's own
+  # best-method selection can keep; pecotmr still applies the exact rule.
+  if (!is.na(argv$min_cv_rsq)) {
+    maxP <- if (is.na(argv$max_cv_pval)) Inf else argv$max_cv_pval
+    passes <- function(j) {
+      f <- trimws(strsplit(xqtlDf$original_data[[j]], ",")[[1L]])
+      f <- vapply(f[nzchar(f)], adaptFilePath, character(1), referenceFile = argv$xqtl_meta, USE.NAMES = FALSE)
+      for (x in f) {
+        cv <- as.list(S4Vectors::mcols(readRDS(x))$cvResult)
+        for (r in cv) {
+          v <- r$metrics
+          if (is.null(v) || is.na(v["rsq"])) return(TRUE)  # no CV info: let pecotmr decide
+          pv <- if (!is.na(v["adj_rsq_pval"])) v["adj_rsq_pval"] else v["pval"]
+          if (v["rsq"] >= argv$min_cv_rsq && (is.infinite(maxP) || (!is.na(pv) && pv < maxP))) return(TRUE)
+        }
+      }
+      FALSE
+    }
+    ok <- unlist(parallel::mclapply(keep, passes, mc.cores = max(1L, argv$threads)))
+    message(sprintf("CV prefilter (r2 >= %g, p < %g): kept %d of %d genes.", argv$min_cv_rsq, maxP, sum(ok), length(keep)))
+    keep <- keep[ok]
+  }
+  lapply(keep, function(j) list(chrom = xqtlChr[[j]],
+    start = as.integer(xqtlDf$start[[j]]), stop = as.integer(xqtlDf$end[[j]]), hit = j))
+} else {
+  lapply(seq_len(nrow(regions)), function(i) {
+    chrom <- regions$chr[[i]]; start <- regions$start[[i]]; stop <- regions$stop[[i]]
+    list(chrom = chrom, start = start, stop = stop,
+         hit = which(xqtlChr == chrom & !is.na(xqtlTss) & xqtlTss >= start & xqtlTss <= stop))
+  })
+}
 rows <- list()
 skipped <- 0L
-for (i in seq_len(nrow(regions))) {
-  chrom <- regions$chr[[i]]; start <- regions$start[[i]]; stop <- regions$stop[[i]]
-  hit <- which(xqtlChr == chrom & !is.na(xqtlTss) &
-               xqtlTss >= start & xqtlTss <= stop)
+for (u in units) {
+  chrom <- u$chrom; start <- u$start; stop <- u$stop; hit <- u$hit
   genes <- character(0); weightFiles <- character(0)
   for (j in hit) {
     files <- trimws(strsplit(xqtlDf$original_data[[j]], ",")[[1L]])
@@ -185,7 +236,7 @@ for (i in seq_len(nrow(regions))) {
   }
 
   rows[[length(rows) + 1L]] <- data.frame(
-    region_id     = sprintf("%s_%d_%d", chrom, start, stop),
+    region_id     = if (argv$unit == "gene") genes[[1L]] else sprintf("%s_%d_%d", chrom, start, stop),
     chrom         = chrom, start = start, stop = stop,
     genes         = paste(genes, collapse = ","),
     weight_files  = paste(weightFiles, collapse = ","),
@@ -195,8 +246,8 @@ for (i in seq_len(nrow(regions))) {
     stringsAsFactors = FALSE)
 }
 if (skipped > 0L)
-  message(sprintf("Skipping %d of %d region(s): no overlapping xQTL weights.",
-                  skipped, nrow(regions)))
+  message(sprintf("Skipping %d of %d unit(s): no overlapping xQTL weights.",
+                  skipped, length(units)))
 if (length(rows) == 0L)
   stop("No regions with overlapping xQTL weights; nothing to analyse.")
 out <- do.call(rbind, rows)

From 34e6fd9eafb0cc91dbb7c022c22dba707c54c945 Mon Sep 17 00:00:00 2001
From: Jenny Anak Empawi <jennyempawi95@gmail.com>
Date: Wed, 28 Oct 2026 10:36:24 -0400
Subject: [PATCH 3/5] Add per-gene TWAS units and a CV prefilter

twas_manifest.R gains --unit gene (one unit per gene spanning its weight window, so GWAS and LD are read per gene) next to the existing region grouping, which stays the script default. An optional CV prefilter (--min-cv-rsq / --max-cv-pval) skips genes where no context/method passes, so their GWAS input is never built; weights without per-row CV results are kept for pecotmr to decide. twas_ctwas.ipynb uses the gene unit and the prefilter by default (--twas_unit, --twas_prefilter).
---
 code/script/pecotmr_integration/twas_manifest.R | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/code/script/pecotmr_integration/twas_manifest.R b/code/script/pecotmr_integration/twas_manifest.R
index 1c5ce0ff1..2526cd62c 100644
--- a/code/script/pecotmr_integration/twas_manifest.R
+++ b/code/script/pecotmr_integration/twas_manifest.R
@@ -62,8 +62,8 @@ p <- add_argument(p, "--gwas-data", type = "character", nargs = Inf,
                   default = NA_character_, help = "Optional inline GWAS sumstats paths")
 p <- add_argument(p, "--column-mapping", type = "character", nargs = Inf,
                   default = NA_character_, help = "Optional inline GWAS column-mapping paths")
-p <- add_argument(p, "--unit", type = "character", default = "gene",
-  help = "Analysis unit: 'gene' (default; one unit per gene spanning its weight window, no LD blocks) or 'region' (legacy: genes grouped by --regions)")
+p <- add_argument(p, "--unit", type = "character", default = "region",
+  help = "Analysis unit: 'region' (default; genes grouped by --regions) or 'gene' (one unit per gene spanning its weight window, no LD blocks; twas_ctwas.ipynb passes this by default)")
 p <- add_argument(p, "--min-cv-rsq", type = "numeric", default = NA,
   help = "gene unit only: skip genes where no context/method reaches this CV r2 (with --max-cv-pval); NA = no prefilter")
 p <- add_argument(p, "--max-cv-pval", type = "numeric", default = NA,
@@ -162,10 +162,10 @@ xqtlChr <- chromAdd(xqtlDf[["#chr"]])
 xqtlTss <- suppressWarnings(as.integer(xqtlDf$TSS))
 
 # ---- per-region: overlapping genes/weights + covering GWAS -----------------
-# Analysis units. gene (default): one unit per xQTL meta row, spanning that
+# Analysis units. gene: one unit per xQTL meta row, spanning that
 # gene's weight window (meta start..end), so every weight variant is scored
 # with GWAS and LD read from the per-chromosome files; --regions /
-# --region-name, if given, only filter genes by TSS. region (legacy): genes
+# --region-name, if given, only filter genes by TSS. region (default): genes
 # binned into --regions by TSS, one unit per region.
 inRegions <- function(j) nrow(regions) == 0L ||
   any(regions$chr == xqtlChr[[j]] & !is.na(xqtlTss[[j]]) &
@@ -182,7 +182,11 @@ units <- if (argv$unit == "gene") {
       f <- trimws(strsplit(xqtlDf$original_data[[j]], ",")[[1L]])
       f <- vapply(f[nzchar(f)], adaptFilePath, character(1), referenceFile = argv$xqtl_meta, USE.NAMES = FALSE)
       for (x in f) {
-        cv <- as.list(S4Vectors::mcols(readRDS(x))$cvResult)
+        w <- readRDS(x)
+      # Legacy list-layout weights carry no per-row cvResult column: keep the
+      # gene and let pecotmr apply the CV rule.
+      if (!methods::is(w, "TwasWeights")) return(TRUE)
+      cv <- as.list(S4Vectors::mcols(w)$cvResult)
         for (r in cv) {
           v <- r$metrics
           if (is.null(v) || is.na(v["rsq"])) return(TRUE)  # no CV info: let pecotmr decide

From 63acf05934ea1f28bccb92de102724868c3109fe Mon Sep 17 00:00:00 2001
From: Jenny Anak Empawi <jennyempawi95@gmail.com>
Date: Tue, 10 Nov 2026 11:19:37 -0500
Subject: [PATCH 4/5] Harden the gene-unit TWAS manifest and add regression
 tests

- twas_ctwas.ipynb keeps region as the default TWAS unit, so existing commands give the same outputs; the gene unit is opt-in with --twas_unit gene.
- CV prefilter: off by default in the notebook (--twas_prefilter turns it on for twas runs), since the ctwas steps read the same manifest and cTWAS uses every gene with weights. It now reads the per-method CV results of the pecotmr 0.9.5+ TwasWeights layout (it previously found no metrics there and kept every gene), ignores methods without a CV r2 such as ensemble, uses the twas step's --rsq-option and --rsq-pval-option, keeps genes with no CV metrics, no cvResult column or an unreadable weight file (reported on stderr), treats --min-cv-rsq <= 0 as off, and loads pecotmr before reading the weights.
- Gene unit: make region_id file-name safe and stop when the xQTL meta repeats a gene, since both would collide on output paths.
- twas_weights_conversion.R: skip all-zero/NA weights before resolving variant IDs, and prefer preset_variants_result IDs when their length fits.
- Tests: the notebook twas test pins the region unit; new tests cover the gene unit (same result as the region unit on the toy gene), gene-unit manifest rows, repeated genes, and the prefilter keep/drop rules.
---
 code/SoS/pecotmr_integration/twas_ctwas.ipynb | 27 ++++---
 .../pecotmr_integration/twas_manifest.R       | 76 ++++++++++++++++---
 .../twas_weights_conversion.R                 | 14 ++--
 .../pecotmr_integration/test_twas_ctwas.py    | 33 ++++++++
 .../pecotmr_integration/test_twas_manifest.py | 59 ++++++++++++++
 5 files changed, 180 insertions(+), 29 deletions(-)

diff --git a/code/SoS/pecotmr_integration/twas_ctwas.ipynb b/code/SoS/pecotmr_integration/twas_ctwas.ipynb
index 2b1601a7e..fb609c2de 100644
--- a/code/SoS/pecotmr_integration/twas_ctwas.ipynb
+++ b/code/SoS/pecotmr_integration/twas_ctwas.ipynb
@@ -976,17 +976,22 @@
     "# (e.g. sos run twas_ctwas.ipynb twas). Gene -> home-LD-block placement is NOT\n",
     "# done here: pecotmr's assembleCtwasInputs places each gene into its block from\n",
     "# the weight `region` provenance.\n",
-    "# TWAS analysis unit: gene (default) gives each gene its own window, the\n",
-    "# span of its weights from the xQTL meta, so GWAS and LD are read per gene\n",
-    "# and no LD blocks are involved; --regions / --region-name then only filter\n",
-    "# genes by TSS. region keeps the legacy grouping of genes by --regions.\n",
-    "# LD blocks enter only at the ctwas step.\n",
-    "parameter: twas_unit = \"gene\"\n",
-    "# Gene unit only: skip genes with no context/method passing the CV cutoffs\n",
-    "# (same values as the twas step), so their GWAS input is never built.\n",
-    "parameter: twas_prefilter = True\n",
+    "# TWAS analysis unit: region (default) groups genes by --regions as before.\n",
+    "# gene gives each gene its own window, the span of its weights from the\n",
+    "# xQTL meta, so GWAS and LD are read per gene and no LD blocks are involved;\n",
+    "# --regions / --region-name then only filter genes by TSS. LD blocks enter\n",
+    "# only at the ctwas step.\n",
+    "parameter: twas_unit = \"region\"\n",
+    "# Gene unit only, off by default: skip genes with no context/method passing\n",
+    "# the CV cutoffs (same cutoffs and metrics as the twas step; keep them in\n",
+    "# sync), so their GWAS input is never built. Turn it on (--twas_prefilter)\n",
+    "# for twas runs only: the ctwas steps read the same manifest, and cTWAS uses\n",
+    "# every gene with weights.\n",
+    "parameter: twas_prefilter = False\n",
     "parameter: rsq_cutoff = 0.01\n",
     "parameter: rsq_pval_cutoff = 0.05\n",
+    "parameter: rsq_option = \"rsq\"\n",
+    "parameter: rsq_pval_option = [\"adj_rsq_pval\", \"pval\"]\n",
     "input: None\n",
     "twas_manifest_file = f\"{cwd:a}/get_analysis_regions/{name}.twas_manifest.tsv\"\n",
     "bash: expand = \"${ }\", container = container\n",
@@ -1000,8 +1005,8 @@
     "        ${(\"--gwas-data \" + \" \".join([str(x) for x in gwas_data])) if len(gwas_data) > 0 else \"\"} \\\n",
     "        ${(\"--column-mapping \" + \" \".join([str(x) for x in column_mapping])) if len(column_mapping) > 0 else \"\"} \\\n",
     "        --unit ${twas_unit} \\\n",
-    "    ${(\"--min-cv-rsq %s --max-cv-pval %s --threads %s\" % (rsq_cutoff, \"Inf\" if rsq_pval_cutoff == float(\"inf\") else rsq_pval_cutoff, numThreads)) if (twas_prefilter and twas_unit == \"gene\") else \"\"} \\\n",
-    "    --output ${twas_manifest_file}"
+    "        ${(\"--min-cv-rsq %s --max-cv-pval %s --rsq-option %s --rsq-pval-option %s --threads %s\" % (rsq_cutoff, \"Inf\" if rsq_pval_cutoff == float(\"inf\") else rsq_pval_cutoff, rsq_option, \",\".join(rsq_pval_option), numThreads)) if (twas_prefilter and twas_unit == \"gene\") else \"\"} \\\n",
+    "        --output ${twas_manifest_file}"
    ]
   },
   {
diff --git a/code/script/pecotmr_integration/twas_manifest.R b/code/script/pecotmr_integration/twas_manifest.R
index 2526cd62c..e3ad1d4dc 100644
--- a/code/script/pecotmr_integration/twas_manifest.R
+++ b/code/script/pecotmr_integration/twas_manifest.R
@@ -68,6 +68,10 @@ p <- add_argument(p, "--min-cv-rsq", type = "numeric", default = NA,
   help = "gene unit only: skip genes where no context/method reaches this CV r2 (with --max-cv-pval); NA = no prefilter")
 p <- add_argument(p, "--max-cv-pval", type = "numeric", default = NA,
   help = "gene unit only: CV p-value that must also be passed (Inf = no p-value check)")
+p <- add_argument(p, "--rsq-option", type = "character", default = "rsq",
+  help = "gene unit prefilter: CV metric compared with --min-cv-rsq (same as twas.R --rsq-option)")
+p <- add_argument(p, "--rsq-pval-option", type = "character", default = "adj_rsq_pval,pval",
+  help = "gene unit prefilter: comma-separated CV p-value metrics, first present is used (same as twas.R --rsq-pval-option)")
 p <- add_argument(p, "--threads", type = "integer", default = 1L,
   help = "Cores for reading weight files in the prefilter")
 p <- add_argument(p, "--output", type = "character", help = "Output manifest TSV path")
@@ -176,28 +180,67 @@ units <- if (argv$unit == "gene") {
   # twas step does not build their GWAS input. Lenient on purpose: a gene is
   # kept if ANY context/method row passes, a superset of what pecotmr's own
   # best-method selection can keep; pecotmr still applies the exact rule.
-  if (!is.na(argv$min_cv_rsq)) {
+  # --min-cv-rsq <= 0 (e.g. --rsq_cutoff 0) disables the prefilter, as it
+  # disables the cutoff in twas.R.
+  if (!is.na(argv$min_cv_rsq) && argv$min_cv_rsq > 0) {
+    # The weights are pecotmr S4 objects: load it so their class and mcols()
+    # resolve here, instead of every gene falling back to "keep".
+    suppressPackageStartupMessages(loadNamespace("pecotmr"))
     maxP <- if (is.na(argv$max_cv_pval)) Inf else argv$max_cv_pval
+    pvalOpt <- trimws(strsplit(argv$rsq_pval_option, ",", fixed = TRUE)[[1L]])
+    # TRUE when one context/method row passes, or when the CV information
+    # needed to judge it is missing (pecotmr then applies the exact rule).
+    rowPasses <- function(v) {
+      r2 <- v[argv$rsq_option]
+      if (r2 < argv$min_cv_rsq) return(FALSE)
+      if (is.infinite(maxP)) return(TRUE)
+      pv <- v[intersect(pvalOpt, names(v))]
+      pv <- pv[!is.na(pv)]
+      if (length(pv) == 0L) return(TRUE)
+      pv[[1L]] < maxP
+    }
     passes <- function(j) {
       f <- trimws(strsplit(xqtlDf$original_data[[j]], ",")[[1L]])
       f <- vapply(f[nzchar(f)], adaptFilePath, character(1), referenceFile = argv$xqtl_meta, USE.NAMES = FALSE)
+      if (length(f) == 0L) return(TRUE)
       for (x in f) {
         w <- readRDS(x)
-      # Legacy list-layout weights carry no per-row cvResult column: keep the
-      # gene and let pecotmr apply the CV rule.
-      if (!methods::is(w, "TwasWeights")) return(TRUE)
-      cv <- as.list(S4Vectors::mcols(w)$cvResult)
+        # Legacy list-layout weights, or TwasWeights without a cvResult
+        # column, carry no per-row CV: keep the gene.
+        if (!methods::is(w, "TwasWeights")) return(TRUE)
+        cv <- as.list(S4Vectors::mcols(w)$cvResult)
+        if (length(cv) == 0L) return(TRUE)
+        # One cvResult per context row; it holds either the metrics directly
+        # or one list(metrics = ...) per method (pecotmr >= 0.9.5 layout).
+        # Methods without a CV r2 (e.g. ensemble) cannot be picked by the
+        # best-r2 rule and are ignored; a gene where no method has one is kept.
+        metrics <- list()
         for (r in cv) {
-          v <- r$metrics
-          if (is.null(v) || is.na(v["rsq"])) return(TRUE)  # no CV info: let pecotmr decide
-          pv <- if (!is.na(v["adj_rsq_pval"])) v["adj_rsq_pval"] else v["pval"]
-          if (v["rsq"] >= argv$min_cv_rsq && (is.infinite(maxP) || (!is.na(pv) && pv < maxP))) return(TRUE)
+          perMethod <- if (!is.null(r$metrics)) list(r) else as.list(r)
+          for (mm in perMethod) {
+            v <- mm$metrics
+            if (!is.null(v) && !is.na(v[argv$rsq_option])) metrics[[length(metrics) + 1L]] <- v
+          }
         }
+        if (length(metrics) == 0L) return(TRUE)
+        if (any(vapply(metrics, rowPasses, logical(1)))) return(TRUE)
+        rm(w)
       }
       FALSE
     }
-    ok <- unlist(parallel::mclapply(keep, passes, mc.cores = max(1L, argv$threads)))
-    message(sprintf("CV prefilter (r2 >= %g, p < %g): kept %d of %d genes.", argv$min_cv_rsq, maxP, sum(ok), length(keep)))
+    # A gene whose weight file cannot be read is kept (with a warning) so the
+    # twas step reports it, rather than failing the whole manifest.
+    # message() rather than warning(): warnings raised in mclapply's forked
+    # workers are lost, messages reach stderr.
+    safePasses <- function(j) tryCatch(passes(j), error = function(e) {
+      message(sprintf("Warning: CV prefilter could not read weights of %s (%s); keeping it.",
+                      xqtlDf$region_id[[j]], conditionMessage(e)))
+      TRUE
+    })
+    res <- parallel::mclapply(keep, safePasses, mc.cores = max(1L, argv$threads))
+    ok <- vapply(res, function(z) !isFALSE(z), logical(1))
+    message(sprintf("CV prefilter (%s >= %g, p < %g): kept %d of %d genes.",
+                    argv$rsq_option, argv$min_cv_rsq, maxP, sum(ok), length(keep)))
     keep <- keep[ok]
   }
   lapply(keep, function(j) list(chrom = xqtlChr[[j]],
@@ -255,6 +298,15 @@ if (skipped > 0L)
 if (length(rows) == 0L)
   stop("No regions with overlapping xQTL weights; nothing to analyse.")
 out <- do.call(rbind, rows)
+if (argv$unit == "gene") {
+  # region_id names the per-gene output files: make it file-name safe (sQTL
+  # IDs carry ':') and require one meta row per gene.
+  out$region_id <- gsub("[^A-Za-z0-9_.-]", "_", out$region_id)
+  dup <- unique(out$region_id[duplicated(out$region_id)])
+  if (length(dup) > 0L)
+    stop("--unit gene needs one xQTL meta row per gene; repeated: ",
+         paste(head(dup, 5L), collapse = ", "), if (length(dup) > 5L) ", ...")
+}
 
 writeManifest(out, argv$output)
-cat(sprintf("Wrote TWAS manifest with %d region(s) to %s\n", nrow(out), argv$output))
+cat(sprintf("Wrote TWAS manifest with %d %s unit(s) to %s\n", nrow(out), argv$unit, argv$output))
diff --git a/code/script/pecotmr_integration/twas_weights_conversion.R b/code/script/pecotmr_integration/twas_weights_conversion.R
index d21da3f36..acc03fa07 100644
--- a/code/script/pecotmr_integration/twas_weights_conversion.R
+++ b/code/script/pecotmr_integration/twas_weights_conversion.R
@@ -52,14 +52,16 @@ for (f in files) {
     names(perf) <- sub("_performance$", "", names(perf))
     for (wn in names(co$twas_weights)) {
       w <- co$twas_weights[[wn]]
-      vids <- if (is.matrix(w)) rownames(w) else names(w)
-      # Some releases store weights for the QC'd variant subset only, without row
-    # names; their IDs are then in preset_variants_result$variant_names
-    if (is.null(vids)) vids <- Find(function(v) length(v) == NROW(w),
-      list(co$variant_names, co$preset_variants_result$variant_names))
-    if (is.null(vids)) stop("No variant IDs of length ", NROW(w), " for ", cn, " ", wn, " in ", f)
       wv <- if (is.matrix(w)) as.numeric(w[, 1]) else as.numeric(w)
       if (all(is.na(wv)) || all(wv == 0, na.rm = TRUE)) next
+      vids <- if (is.matrix(w)) rownames(w) else names(w)
+      # Some releases store weights for the QC'd variant subset only, without
+      # row names; their IDs are then in preset_variants_result$variant_names.
+      # Prefer that subset when its length fits, since variant_names can have
+      # the same length but list the full set.
+      if (is.null(vids)) vids <- Find(function(v) length(v) == NROW(w),
+        list(co$preset_variants_result$variant_names, co$variant_names))
+      if (is.null(vids)) stop("No variant IDs of length ", NROW(w), " for ", cn, " ", wn, " in ", f)
       m <- sub("_weights$", "", wn)
       fits <- if (m %in% c("susie", "susie_inf", "susie_ash")) co$susie_weights_intermediate else NULL
       E[[length(E) + 1L]] <- TwasWeightsRow(variantIds = vids, weights = wv,
diff --git a/tests/notebooks/pecotmr_integration/test_twas_ctwas.py b/tests/notebooks/pecotmr_integration/test_twas_ctwas.py
index a68bcaeef..812fddc74 100644
--- a/tests/notebooks/pecotmr_integration/test_twas_ctwas.py
+++ b/tests/notebooks/pecotmr_integration/test_twas_ctwas.py
@@ -36,6 +36,9 @@ def test_twas(run_sos, repo_root, tmp_path):
             "xqtl_meta_data": meta,
             "ld_meta_data": fx / "ld_reference/ld_meta_file.tsv",
             "region-name": "chr22_10000000_19000000",
+            # The notebook defaults to the gene unit; this test covers the
+            # region unit (test_twas_gene_unit covers the default).
+            "twas_unit": "region",
             # twas.R's default, under which its expected output was made (the
             # notebook defaults to 0.05)
             "rsq_pval_cutoff": "inf",
@@ -55,6 +58,36 @@ def test_twas(run_sos, repo_root, tmp_path):
                             mode="tolerant", rtol=1e-6, atol=1e-8)
 
 
+
+def test_twas_gene_unit(run_sos, repo_root, tmp_path):
+    """The gene unit (--twas_unit gene, CV prefilter off by default) scores the toy
+    gene over its own weight window. The window equals the toy region here, so the result
+    must match the region-unit expected output."""
+    fx = repo_root / "tests" / "fixtures"
+    tw = fx / "twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.twas_weights.s4.rds"
+    meta = tmp_path / "xqtl_meta.tsv"
+    meta.write_text("#chr\tstart\tend\tregion_id\tTSS\toriginal_data\tcontexts\n"
+                    f"chr22\t10000000\t19000000\tENSG00000130538\t15528191\t{tw}\tbulk_rnaseq\n")
+    cwd = tmp_path / "out"
+    p = run_sos(
+        repo_root / "pipeline/twas_ctwas.ipynb", "twas",
+        {
+            "cwd": cwd,
+            "name": "protocol_example",
+            "gwas_meta_data": fx / "twas/protocol_example.twas.gwas_meta.tsv",
+            "xqtl_meta_data": meta,
+            "ld_meta_data": fx / "ld_reference/ld_meta_file.tsv",
+            "region-name": "chr22_10000000_19000000",
+            "rsq_pval_cutoff": "inf",
+            "twas_unit": "gene",
+            "modular_script_dir": repo_root / "code/script",
+        },
+        cwd=repo_root, timeout=600)
+    assert p.returncode == 0, p.stdout + p.stderr
+    out = cwd / "twas" / "protocol_example.ENSG00000130538.ENSG00000130538.twas.rds"
+    assert_matches_expected(out, fx / "twas/expected/twas.chr22.rds",
+                            mode="tolerant", rtol=1e-6, atol=1e-8)
+
 # The blessed cTWAS weights as committed S4 TwasWeights (region-provenance
 # stamped -- chrom/p0/p1 -- which assembleCtwasInputs needs for LD-block
 # placement; the legacy_ctwas_weights_to_s4 conversion step is retired).
diff --git a/tests/scripts/pecotmr_integration/test_twas_manifest.py b/tests/scripts/pecotmr_integration/test_twas_manifest.py
index a152b5faa..b1f4047c3 100644
--- a/tests/scripts/pecotmr_integration/test_twas_manifest.py
+++ b/tests/scripts/pecotmr_integration/test_twas_manifest.py
@@ -23,3 +23,62 @@ def test_twas_manifest(run_r, repo_root, tmp_path):
     assert rows[1][0] == "chr22_10000000_19000000" and "ENSG00000130538" in rows[1][4]
     assert_matches_expected(out, repo_root / EXPECTED, mode="tolerant",
                             rtol=1e-6, atol=1e-8, normalize_paths=True)
+
+
+TOY_WEIGHTS = "tests/fixtures/twas/protocol_example.twas.reshaped_toy.chr22_ENSG00000130538.twas_weights.s4.rds"
+
+
+def _gene_meta(path, rows):
+    path.write_text("#chr\tstart\tend\tregion_id\tTSS\toriginal_data\tcontexts\n"
+                    + "".join(f"chr22\t10000000\t19000000\t{g}\t15528191\t{w}\tbulk_rnaseq\n"
+                              for g, w in rows))
+    return path
+
+
+def _gene_manifest(run_r, repo_root, meta, out, *extra):
+    return run_r(repo_root / "code/script/pecotmr_integration/twas_manifest.R",
+                 ["--gwas-meta", repo_root / "tests/fixtures/twas/protocol_example.twas.gwas_meta.tsv",
+                  "--xqtl-meta", meta, "--unit", "gene", *extra, "--output", out], timeout=120)
+
+
+def test_twas_manifest_gene_unit(run_r, repo_root, tmp_path):
+    """--unit gene: one row per gene, named by the gene and spanning its meta window."""
+    meta = _gene_meta(tmp_path / "meta.tsv", [("ENSG00000130538", repo_root / TOY_WEIGHTS)])
+    out = tmp_path / "manifest.tsv"
+    p = _gene_manifest(run_r, repo_root, meta, out, "--min-cv-rsq", "0")
+    assert p.returncode == 0, p.stdout + p.stderr
+    rows = [l.split("\t") for l in out.read_text().splitlines()]
+    assert len(rows) == 2
+    assert rows[1][:4] == ["ENSG00000130538", "chr22", "10000000", "19000000"]
+
+
+def test_twas_manifest_gene_unit_rejects_repeated_gene(run_r, repo_root, tmp_path):
+    w = repo_root / TOY_WEIGHTS
+    meta = _gene_meta(tmp_path / "meta.tsv", [("ENSG00000130538", w), ("ENSG00000130538", w)])
+    p = _gene_manifest(run_r, repo_root, meta, tmp_path / "manifest.tsv")
+    assert p.returncode != 0
+    assert "one xQTL meta row per gene" in p.stdout + p.stderr
+
+
+def test_twas_manifest_prefilter_keeps_unreadable_weights(run_r, repo_root, tmp_path):
+    """A weight file the CV prefilter cannot read is kept (with a warning) for the
+    twas step to report, instead of failing the manifest."""
+    bad = tmp_path / "corrupt.rds"
+    bad.write_text("not an RDS file")
+    meta = _gene_meta(tmp_path / "meta.tsv", [("ENSG00000130538", repo_root / TOY_WEIGHTS),
+                                              ("ENSG_CORRUPT", bad)])
+    out = tmp_path / "manifest.tsv"
+    p = _gene_manifest(run_r, repo_root, meta, out, "--min-cv-rsq", "0.01", "--max-cv-pval", "Inf")
+    assert p.returncode == 0, p.stdout + p.stderr
+    assert "could not read weights of ENSG_CORRUPT" in p.stdout + p.stderr
+    ids = [l.split("\t")[0] for l in out.read_text().splitlines()[1:]]
+    assert "ENSG_CORRUPT" in ids
+
+
+def test_twas_manifest_prefilter_drops_failing_genes(run_r, repo_root, tmp_path):
+    """A CV r2 cutoff no model can reach drops the gene, so nothing is left."""
+    meta = _gene_meta(tmp_path / "meta.tsv", [("ENSG00000130538", repo_root / TOY_WEIGHTS)])
+    p = _gene_manifest(run_r, repo_root, meta, tmp_path / "manifest.tsv",
+                       "--min-cv-rsq", "2", "--max-cv-pval", "Inf")
+    assert p.returncode != 0
+    assert "kept 0 of 1 genes" in p.stdout + p.stderr

From 2afd16e0d1eea62816db5cf0d3977d40d9d6392d Mon Sep 17 00:00:00 2001
From: Jenny Anak Empawi <jennyempawi95@gmail.com>
Date: Tue, 10 Nov 2026 13:30:21 -0500
Subject: [PATCH 5/5] Tighten TWAS edge cases and document the new options

- twas_manifest.R: --unit gene stops early when a gene has no start/end in the xQTL meta, checks for repeated genes before making IDs file-name safe, and stops when two genes would share a file name. The CV prefilter skips cvResult entries that are not lists instead of reporting the file as unreadable.
- twas.R: a missing GwasSumStats file still means an empty window, but a missing folder now stops as a wrong path; the noResult text no longer assumes an LD block.
- twas_ctwas.ipynb: the twas step removes an old GwasSumStats file before rebuilding it, makes the gene part of output names file-name safe, and the Input and Command Interface sections list --twas_unit, --twas_prefilter, --rsq_option and --rsq_pval_option.
- twas_weights_conversion.R: report genes without usable weights with message(), since the warning was lost on quit().
- Tests: prefilter p-value cutoff (keep and drop, with the pval fall-back and ensemble skip), NA gene window, file-name safe IDs and clashes, missing GwasSumStats file and folder, and conversion of weights without names, with unmatched IDs and with all-zero weights.
---
 code/SoS/pecotmr_integration/twas_ctwas.ipynb | 24 ++++++-
 code/script/pecotmr_integration/twas.R        | 15 +++--
 .../pecotmr_integration/twas_manifest.R       | 27 ++++++--
 .../twas_weights_conversion.R                 |  7 ++-
 .../scripts/pecotmr_integration/test_twas.py  | 19 ++++++
 .../pecotmr_integration/test_twas_manifest.py | 46 ++++++++++++++
 .../test_twas_weights_conversion.py           | 63 +++++++++++++++++++
 7 files changed, 187 insertions(+), 14 deletions(-)
 create mode 100644 tests/scripts/pecotmr_integration/test_twas_weights_conversion.py

diff --git a/code/SoS/pecotmr_integration/twas_ctwas.ipynb b/code/SoS/pecotmr_integration/twas_ctwas.ipynb
index fb609c2de..688f0c16f 100644
--- a/code/SoS/pecotmr_integration/twas_ctwas.ipynb
+++ b/code/SoS/pecotmr_integration/twas_ctwas.ipynb
@@ -103,6 +103,9 @@
     "- `--ld_reference_sample_size` (sample size of the LD reference panel, used to scale the LD matrices)\n",
     "- `--prior_var_structure` (how the cTWAS prior variance is shared across groups)\n",
     "- `--gwas_study <study_id>` (`ctwas` only: the one GWAS study to fine-map, required when `--gwas_meta_data` lists several; run each study with its own `--cwd`)\n",
+    "- `--twas_unit` (`region`, the default, groups genes by `--regions`; `gene` gives each gene its own window, the span of its weights in the xQTL meta, and `--regions` / `--region-name` then only select genes by TSS; outputs are then named `<name>.<gene>.<gene>.twas.rds`)\n",
+    "- `--twas_prefilter` (gene unit only, off by default: skip genes where no context/method passes `--rsq_cutoff` / `--rsq_pval_cutoff`, so their GWAS input is never built; use it for `twas` runs only, since the `ctwas` steps read the same gene list)\n",
+    "- `--rsq_option`, `--rsq_pval_option` (which cross-validation metrics `--rsq_cutoff` and `--rsq_pval_cutoff` apply to; default `rsq` and `adj_rsq_pval pval`, the first one present is used)\n",
     "- `--rsq_cutoff` (minimum cross-validation r-squared for a weight model to be used)\n",
     "- `--rsq_pval_cutoff` (maximum cross-validation p-value for a weight model to be used)\n",
     "- `--twas_weight_cutoff` (drop weights below this magnitude before assembling cTWAS inputs)\n",
@@ -634,6 +637,18 @@
     "                        (one legacy RDS per gene; the gene is the ENSG ID in the\n",
     "                        name).\n",
     "  get_analysis_regions:\n",
+    "    Workflow Options:\n",
+    "      --twas-unit region\n",
+    "                        TWAS analysis unit: region (default) groups genes by\n",
+    "                        --regions as before. gene gives each gene its own\n",
+    "                        window, the span of its weights from the xQTL meta.\n",
+    "      --[no-]twas-prefilter (default to False)\n",
+    "                        Gene unit only: skip genes with no context/method\n",
+    "                        passing the CV cutoffs. Use for twas runs only.\n",
+    "      --rsq-cutoff 0.01 (as float)\n",
+    "      --rsq-pval-cutoff 0.05 (as float)\n",
+    "      --rsq-option rsq\n",
+    "      --rsq-pval-option adj_rsq_pval pval (as list)\n",
     "  twas:\n",
     "    Workflow Options:\n",
     "      --coverage 'cs_coverage_0.95'\n",
@@ -1066,7 +1081,7 @@
     "        \"No regions with overlapping xQTL weights found; skipping TWAS step.\")\n",
     "_weight_groups = [j['weight_files'].split(',') for j in jobs]\n",
     "input: [f for grp in _weight_groups for f in grp], group_by = lambda x: _weight_groups, group_with = \"jobs\"\n",
-    "output: [f'{cwd:a}/{step_name}/{name}.{_jobs[\"region_id\"]}.{gene}.twas.rds' for gene in sorted(set(_jobs[\"genes\"].split(',')))]\n",
+    "output: [f'{cwd:a}/{step_name}/{name}.{_jobs[\"region_id\"]}.{__import__(\"re\").sub(r\"[^A-Za-z0-9_.-]\", \"_\", gene)}.twas.rds' for gene in sorted(set(_jobs[\"genes\"].split(',')))]\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 = '${ }', stdout = f\"{_output[0]:n}.stdout\", stderr = f\"{_output[0]:n}.stderr\", container = container\n",
     "    set -e\n",
@@ -1083,6 +1098,9 @@
     "\n",
     "    # (1) One multi-study GwasSumStats over the whole region.\n",
     "    gss=\"$outdir/${name}.$region_id.gwas_sumstats.rds\"\n",
+    "    # Remove a file left by an earlier run: twas.R reads a missing file as\n",
+    "    # \"no GWAS variants here\", so a stale one must not stand in for it.\n",
+    "    rm -f \"$gss\"\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/gwas_sumstats_construct.R \\\n",
     "        --study \"$studies\" \\\n",
     "        --gwas-tsv \"$gwas_tsvs\" \\\n",
@@ -1109,9 +1127,11 @@
     "        weight=\"$2\"\n",
     "        shift 2\n",
     "        study=$(basename \"$weight\" | cut -d. -f1)\n",
+    "        # File-name safe gene ID (sQTL IDs carry ':'), matching the output: names.\n",
+    "        safe_gene=$(printf '%s' \"$gene\" | tr -c 'A-Za-z0-9_.-' '_')\n",
     "        printf '%s\\t%s\\t%s\\n' \"$weight\" \\\n",
     "            \"${fine_mapping_dir}/$study/$study.$gene.fine_mapping.s4.rds\" \\\n",
-    "            \"$outdir/${name}.$region_id.$gene.twas.rds\" >> \"$gene_table\"\n",
+    "            \"$outdir/${name}.$region_id.$safe_gene.twas.rds\" >> \"$gene_table\"\n",
     "    done\n",
     "    # One R call per region, so the region GwasSumStats (and its LD reference) is read once for every gene.\n",
     "    Rscript ${modular_script_dir}/pecotmr_integration/twas.R \\\n",
diff --git a/code/script/pecotmr_integration/twas.R b/code/script/pecotmr_integration/twas.R
index 5b9a418ba..738d5883a 100644
--- a/code/script/pecotmr_integration/twas.R
+++ b/code/script/pecotmr_integration/twas.R
@@ -112,9 +112,16 @@ if (length(fmPaths) != length(twPaths)) stop("--fine-mapping-result needs one pa
 # this object's ldReference, per gene and scoped to that gene's matched
 # variants (and for a genotype panel it uses the LD factor, so the p x p matrix
 # is never formed). Nothing here computes or caches LD.
-# gwas_sumstats_construct.R writes no GwasSumStats file for an LD block with
-# no GWAS variants; every gene of such a block gets an empty result.
-gss <- if (file.exists(argv$gwas_sumstats)) readRDS(argv$gwas_sumstats) else NULL
+# gwas_sumstats_construct.R writes no GwasSumStats file for a window with no
+# GWAS variants; every gene of such a window gets an empty result. A missing
+# folder is a wrong path, not an empty window, and stops here.
+gss <- if (file.exists(argv$gwas_sumstats)) readRDS(argv$gwas_sumstats) else {
+  if (!dir.exists(dirname(argv$gwas_sumstats)))
+    stop("--gwas-sumstats ", argv$gwas_sumstats, ": its folder does not exist.")
+  message("No GwasSumStats file at ", argv$gwas_sumstats,
+          " (none is written for a window without GWAS variants); every gene gets an empty result.")
+  NULL
+}
 
 readFmr <- function(p) if (nzchar(p) && p != "." && file.exists(p)) readRDS(p) else NULL
 
@@ -124,7 +131,7 @@ readFmr <- function(p) if (nzchar(p) && p != "." && file.exists(p)) readRDS(p) e
 # A gene with no scorable (QTL, GWAS) pair gets an empty list with a `noResult`
 # attribute, so the other genes of the block are still scored.
 for (i in seq_along(twPaths)) {
-  res <- if (is.null(gss)) structure(list(), noResult = "no GWAS variants in this LD block") else
+  res <- if (is.null(gss)) structure(list(), noResult = "no GWAS variants in this window") else
     tryCatch(causalInferencePipeline(
     gwasSumStats = gss,
     twasWeights = readRDS(twPaths[i]),
diff --git a/code/script/pecotmr_integration/twas_manifest.R b/code/script/pecotmr_integration/twas_manifest.R
index e3ad1d4dc..968aa4726 100644
--- a/code/script/pecotmr_integration/twas_manifest.R
+++ b/code/script/pecotmr_integration/twas_manifest.R
@@ -176,6 +176,14 @@ inRegions <- function(j) nrow(regions) == 0L ||
       regions$start <= xqtlTss[[j]] & regions$stop >= xqtlTss[[j]])
 units <- if (argv$unit == "gene") {
   keep <- Filter(inRegions, seq_len(nrow(xqtlDf)))
+  # The gene window comes from the meta start/end: stop now rather than pass
+  # an NA window to the twas step.
+  noWin <- keep[is.na(suppressWarnings(as.integer(xqtlDf$start[keep]))) |
+                is.na(suppressWarnings(as.integer(xqtlDf$end[keep])))]
+  if (length(noWin) > 0L)
+    stop("--unit gene needs a start and end for every gene in --xqtl-meta; missing for: ",
+         paste(head(xqtlDf$region_id[noWin], 5L), collapse = ", "),
+         if (length(noWin) > 5L) ", ...")
   # Optional prefilter: drop genes that can have no imputable model, so the
   # twas step does not build their GWAS input. Lenient on purpose: a gene is
   # kept if ANY context/method row passes, a superset of what pecotmr's own
@@ -216,9 +224,11 @@ units <- if (argv$unit == "gene") {
         # best-r2 rule and are ignored; a gene where no method has one is kept.
         metrics <- list()
         for (r in cv) {
-          perMethod <- if (!is.null(r$metrics)) list(r) else as.list(r)
+          if (!is.list(r)) next
+          perMethod <- if (!is.null(r[["metrics"]])) list(r) else as.list(r)
           for (mm in perMethod) {
-            v <- mm$metrics
+            if (!is.list(mm)) next
+            v <- mm[["metrics"]]
             if (!is.null(v) && !is.na(v[argv$rsq_option])) metrics[[length(metrics) + 1L]] <- v
           }
         }
@@ -299,13 +309,20 @@ if (length(rows) == 0L)
   stop("No regions with overlapping xQTL weights; nothing to analyse.")
 out <- do.call(rbind, rows)
 if (argv$unit == "gene") {
-  # region_id names the per-gene output files: make it file-name safe (sQTL
-  # IDs carry ':') and require one meta row per gene.
-  out$region_id <- gsub("[^A-Za-z0-9_.-]", "_", out$region_id)
+  # region_id names the per-gene output files: require one meta row per gene,
+  # then make it file-name safe (sQTL IDs carry ':') and check that no two
+  # genes end up with the same name.
   dup <- unique(out$region_id[duplicated(out$region_id)])
   if (length(dup) > 0L)
     stop("--unit gene needs one xQTL meta row per gene; repeated: ",
          paste(head(dup, 5L), collapse = ", "), if (length(dup) > 5L) ", ...")
+  safe <- gsub("[^A-Za-z0-9_.-]", "_", out$region_id)
+  clash <- out$region_id[duplicated(safe) | duplicated(safe, fromLast = TRUE)]
+  if (length(clash) > 0L)
+    stop("--unit gene: these gene IDs give the same file name once characters ",
+         "other than letters, digits, '_', '.' and '-' become '_': ",
+         paste(head(clash, 6L), collapse = ", "), if (length(clash) > 6L) ", ...")
+  out$region_id <- safe
 }
 
 writeManifest(out, argv$output)
diff --git a/code/script/pecotmr_integration/twas_weights_conversion.R b/code/script/pecotmr_integration/twas_weights_conversion.R
index acc03fa07..e03d1bf03 100644
--- a/code/script/pecotmr_integration/twas_weights_conversion.R
+++ b/code/script/pecotmr_integration/twas_weights_conversion.R
@@ -57,8 +57,8 @@ for (f in files) {
       vids <- if (is.matrix(w)) rownames(w) else names(w)
       # Some releases store weights for the QC'd variant subset only, without
       # row names; their IDs are then in preset_variants_result$variant_names.
-      # Prefer that subset when its length fits, since variant_names can have
-      # the same length but list the full set.
+      # That subset is tried first because it is what the weights were fitted
+      # on; variant_names is the fallback, used when its length fits.
       if (is.null(vids)) vids <- Find(function(v) length(v) == NROW(w),
         list(co$preset_variants_result$variant_names, co$variant_names))
       if (is.null(vids)) stop("No variant IDs of length ", NROW(w), " for ", cn, " ", wn, " in ", f)
@@ -77,7 +77,8 @@ if (!length(E)) {
   # e.g. a gene stored as NULL or with all-zero weights: write an empty result
   # and empty meta rows, so step 2 leaves the gene out instead of failing
   msg <- paste("No usable (non-zero) weights in:", a$legacy)
-  warning(msg)
+  # message(), not warning(): a top-level warning is lost when quit() runs first
+  message("Warning: ", msg)
   dir.create(dirname(a$output), showWarnings = FALSE, recursive = TRUE)
   saveRDS(structure(list(), noResult = msg), a$output)
   file.create(a$meta_row)
diff --git a/tests/scripts/pecotmr_integration/test_twas.py b/tests/scripts/pecotmr_integration/test_twas.py
index 54d74c18a..09837c70f 100644
--- a/tests/scripts/pecotmr_integration/test_twas.py
+++ b/tests/scripts/pecotmr_integration/test_twas.py
@@ -34,3 +34,22 @@ def test_twas(via_table, run_r, read_rds, repo_root, tmp_path):
     assert p.returncode == 0, p.stdout + p.stderr
     assert read_rds(out)["class"] == "GRanges"
     assert_matches_expected(out, repo_root / EXPECTED, mode="tolerant", rtol=1e-6, atol=1e-8)
+
+
+def test_twas_missing_gwas_sumstats(run_r, repo_root, tmp_path):
+    """No GwasSumStats file (gwas_sumstats_construct.R writes none for a window
+    without GWAS variants): every gene gets an empty result. A missing folder
+    is a wrong path and must fail instead."""
+    r = repo_root
+    out = tmp_path / "twas.rds"
+    p = run_r(r / "code/script/pecotmr_integration/twas.R",
+              ["--twas-weights", r / TW, "--output", out,
+               "--gwas-sumstats", tmp_path / "none.gwas_sumstats.rds"], timeout=300)
+    assert p.returncode == 0, p.stdout + p.stderr
+    assert out.exists()
+    assert "No GwasSumStats file" in p.stdout + p.stderr
+    p = run_r(r / "code/script/pecotmr_integration/twas.R",
+              ["--twas-weights", r / TW, "--output", tmp_path / "twas2.rds",
+               "--gwas-sumstats", tmp_path / "no_such_dir" / "gss.rds"], timeout=300)
+    assert p.returncode != 0
+    assert "folder does not exist" in p.stdout + p.stderr
diff --git a/tests/scripts/pecotmr_integration/test_twas_manifest.py b/tests/scripts/pecotmr_integration/test_twas_manifest.py
index b1f4047c3..bdd56aa5b 100644
--- a/tests/scripts/pecotmr_integration/test_twas_manifest.py
+++ b/tests/scripts/pecotmr_integration/test_twas_manifest.py
@@ -82,3 +82,49 @@ def test_twas_manifest_prefilter_drops_failing_genes(run_r, repo_root, tmp_path)
                        "--min-cv-rsq", "2", "--max-cv-pval", "Inf")
     assert p.returncode != 0
     assert "kept 0 of 1 genes" in p.stdout + p.stderr
+
+
+@pytest.mark.parametrize("max_p, kept", [("0.05", True), ("0.02", False)])
+def test_twas_manifest_prefilter_pvalue_cutoff(run_r, repo_root, tmp_path, max_p, kept):
+    """The CV p-value cutoff: the toy gene's best model (mrash, rsq 0.089, CV
+    p 0.022) passes p < 0.05 and no model passes p < 0.02. The toy weights have
+    no adj_rsq_pval, so this also covers the fall-back to pval, and they carry
+    an ensemble method without CV metrics, which must be skipped."""
+    meta = _gene_meta(tmp_path / "meta.tsv", [("ENSG00000130538", repo_root / TOY_WEIGHTS)])
+    p = _gene_manifest(run_r, repo_root, meta, tmp_path / "manifest.tsv",
+                       "--min-cv-rsq", "0.01", "--max-cv-pval", max_p)
+    out = p.stdout + p.stderr
+    if kept:
+        assert p.returncode == 0, out
+        assert "kept 1 of 1 genes" in out
+    else:
+        assert p.returncode != 0
+        assert "kept 0 of 1 genes" in out
+
+
+def test_twas_manifest_gene_unit_needs_window(run_r, repo_root, tmp_path):
+    meta = tmp_path / "meta.tsv"
+    meta.write_text("#chr\tstart\tend\tregion_id\tTSS\toriginal_data\tcontexts\n"
+                    f"chr22\tNA\t19000000\tENSG00000130538\t15528191\t{repo_root / TOY_WEIGHTS}\tbulk_rnaseq\n")
+    p = _gene_manifest(run_r, repo_root, meta, tmp_path / "manifest.tsv")
+    assert p.returncode != 0
+    assert "needs a start and end for every gene" in p.stdout + p.stderr
+
+
+def test_twas_manifest_gene_unit_file_name_safe(run_r, repo_root, tmp_path):
+    """sQTL-style IDs with ':' become file-name safe unit IDs."""
+    meta = _gene_meta(tmp_path / "meta.tsv", [("ENSG00000130538:clu_1", repo_root / TOY_WEIGHTS)])
+    out = tmp_path / "manifest.tsv"
+    p = _gene_manifest(run_r, repo_root, meta, out)
+    assert p.returncode == 0, p.stdout + p.stderr
+    row = out.read_text().splitlines()[1].split("\t")
+    assert row[0] == "ENSG00000130538_clu_1"
+    assert row[4] == "ENSG00000130538:clu_1"
+
+
+def test_twas_manifest_gene_unit_rejects_name_clash(run_r, repo_root, tmp_path):
+    w = repo_root / TOY_WEIGHTS
+    meta = _gene_meta(tmp_path / "meta.tsv", [("ENSG1:a", w), ("ENSG1_a", w)])
+    p = _gene_manifest(run_r, repo_root, meta, tmp_path / "manifest.tsv")
+    assert p.returncode != 0
+    assert "give the same file name" in p.stdout + p.stderr
diff --git a/tests/scripts/pecotmr_integration/test_twas_weights_conversion.py b/tests/scripts/pecotmr_integration/test_twas_weights_conversion.py
new file mode 100644
index 000000000..3d2e8597e
--- /dev/null
+++ b/tests/scripts/pecotmr_integration/test_twas_weights_conversion.py
@@ -0,0 +1,63 @@
+"""Tier B: twas_weights_conversion.R edge cases on small legacy weight files
+built here: weights without variant names, genes with no usable weights."""
+from __future__ import annotations
+
+import pytest
+
+BUILD = r'''
+args <- commandArgs(TRUE); out <- args[1]; mode <- args[2]
+ids <- c("chr22:15500100:A:G", "chr22:15500200:C:T", "chr22:15500300:G:A")
+w <- switch(mode,
+  zero    = list(enet_weights = matrix(0, 2, 1, dimnames = list(ids[1:2], "y"))),
+  preset  = list(enet_weights = c(0.1, -0.2)),
+  nomatch = list(enet_weights = c(0.1, -0.2, 0.3, 0.4)))
+perf <- matrix(c(0.3, 0.09, 0.08, 0.01, 1, 1), 1,
+               dimnames = list(NULL, c("corr", "rsq", "adj_rsq", "pval", "RMSE", "MAE")))
+co <- list(twas_weights = w,
+           twas_cv_result = list(performance = list(enet_performance = perf)),
+           region_info = list(grange = list(chrom = "22", start = 15000000L, end = 16000000L),
+                              region_coord = list(chrom = "22", start = 15500000L, end = 15500001L)),
+           variant_names = ids,
+           preset_variants_result = list(variant_names = ids[2:3]))
+saveRDS(list(ENSG00000000001 = list(bulk_rnaseq_ENSG00000000001 = co)), out)
+'''
+
+
+def _convert(run_r, repo_root, tmp_path, mode):
+    build = tmp_path / "build.R"
+    build.write_text(BUILD)
+    legacy = tmp_path / f"{mode}.ENSG00000000001.univariate_twas_weights.rds"
+    assert run_r(build, [legacy, mode], timeout=120).returncode == 0
+    out = tmp_path / "out" / "toy.ENSG00000000001.twas_weights.s4.rds"
+    meta = tmp_path / "meta_row.tsv"
+    p = run_r(repo_root / "code/script/pecotmr_integration/twas_weights_conversion.R",
+              ["--legacy", legacy, "--study", "toy", "--output", out, "--meta-row", meta],
+              timeout=300)
+    return p, out, meta
+
+
+def test_conversion_ids_from_preset_variants(run_r, read_rds, repo_root, tmp_path):
+    """Weights without names take their IDs from the vector of matching length
+    (here preset_variants_result, 2 of the 3 variant_names)."""
+    p, out, meta = _convert(run_r, repo_root, tmp_path, "preset")
+    assert p.returncode == 0, p.stdout + p.stderr
+    info = read_rds(out)
+    assert info["class"] == "TwasWeights"
+    assert info["nrow"] == 1
+    assert meta.read_text().split("\t")[3] == "ENSG00000000001"
+
+
+def test_conversion_stops_without_matching_ids(run_r, repo_root, tmp_path):
+    p, _, _ = _convert(run_r, repo_root, tmp_path, "nomatch")
+    assert p.returncode != 0
+    assert "No variant IDs of length 4" in p.stdout + p.stderr
+
+
+def test_conversion_gene_without_usable_weights(run_r, repo_root, tmp_path):
+    """All-zero weights: an empty result and an empty meta row, so the gene is
+    left out of the meta table instead of failing the conversion."""
+    p, out, meta = _convert(run_r, repo_root, tmp_path, "zero")
+    assert p.returncode == 0, p.stdout + p.stderr
+    assert "No usable (non-zero) weights" in p.stdout + p.stderr
+    assert out.exists()
+    assert meta.read_text() == ""