From 60e29b2b667d89537250bcda871e810b168c12f2 Mon Sep 17 00:00:00 2001 From: Yining97 Date: Tue, 11 Aug 2026 01:22:52 -0400 Subject: [PATCH] Fix GWAS fine-mapping postprocessing export - pip_sum crashed on a GwasFineMappingResult (no context/trait column) for any region with signal; read identity fields NULL-safely so absent columns are NA. - export_top_loci produced a silent empty bed under --gwas: export_suffix stayed cis_results_db while step 1 writes block_results_db. Derive it from the gwas flag. - export_top_loci_3 combine now uses zcat -f so it doesn't abort on macOS once a non-empty bed exists to concatenate (no-op on Linux). --- code/SoS/mnm_analysis/mnm_postprocessing.ipynb | 14 +++++++++----- .../fine_mapping_cis_db_export.R | 9 +++++++-- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/code/SoS/mnm_analysis/mnm_postprocessing.ipynb b/code/SoS/mnm_analysis/mnm_postprocessing.ipynb index 89ea6be1c..e495554d7 100644 --- a/code/SoS/mnm_analysis/mnm_postprocessing.ipynb +++ b/code/SoS/mnm_analysis/mnm_postprocessing.ipynb @@ -958,7 +958,10 @@ "# `export_suffix` deliberately differs from step 1's `--suffix` (input rds) to avoid CLI collision.\n", "parameter: export_path = cwd\n", "parameter: export_prefix = '' # falls back to step-1 `name` if empty\n", - "parameter: export_suffix = 'cis_results_db.rds'\n", + "parameter: export_suffix = '' # empty => derived from `gwas` below\n", + "# Re-declared so the suffix can default by mode: step 1 writes\n", + "# block_results_db.rds under --gwas True, else cis_results_db.rds.\n", + "parameter: gwas = False\n", "parameter: fsusie_prefix = ''\n", "parameter: preset_top_loci = False\n", "# PIP cutoff for getTopLoci; optional independent CS purity (min.abs.corr) cutoff (<0 => none).\n", @@ -969,8 +972,9 @@ "\n", "_pfx = export_prefix if export_prefix else name\n", "import glob as _glob\n", - "_rds_files = sorted(_glob.glob(f\"{export_path}/{_pfx}.{fsusie_prefix}*.{export_suffix}\"))\n", - "stop_if(len(_rds_files) == 0, f\"No RDS files matched: {export_path}/{_pfx}.{fsusie_prefix}*.{export_suffix}\")\n", + "_suffix = export_suffix if export_suffix else ('block_results_db.rds' if gwas else 'cis_results_db.rds')\n", + "_rds_files = sorted(_glob.glob(f\"{export_path}/{_pfx}.{fsusie_prefix}*.{_suffix}\"))\n", + "stop_if(len(_rds_files) == 0, f\"No RDS files matched: {export_path}/{_pfx}.{fsusie_prefix}*.{_suffix}\")\n", "_purity_arg = f\"--min-purity {min_purity}\" if float(min_purity) >= 0 else \"\"\n", "\n", "input: _rds_files, group_by = 1\n", @@ -1019,14 +1023,14 @@ " # 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\" | head -1\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\" | tail -n +2\n", + " zcat -f \"$f\" | tail -n +2\n", " fi\n", " done\n", " } | $GZ > ${_output}\n" diff --git a/code/script/pecotmr_integration/fine_mapping_cis_db_export.R b/code/script/pecotmr_integration/fine_mapping_cis_db_export.R index 6a7280d2a..b7a45f96c 100644 --- a/code/script/pecotmr_integration/fine_mapping_cis_db_export.R +++ b/code/script/pecotmr_integration/fine_mapping_cis_db_export.R @@ -111,8 +111,13 @@ cat(sprintf("Wrote cis_results_db (%d FMR input(s) -> %d rows) + meta (region_id # ---- 4. pip_sum (optional) ------------------------------------------------- if (!is.na(argv$pip_sum_output)) { rows <- lapply(seq_len(nrow(db)), function(i) { - sel <- list(study = as.character(db$study[i]), context = as.character(db$context[i]), - trait = as.character(db$trait[i]), method = as.character(db$method[i])) + # NULL-safe per-row identity read: a GwasFineMappingResult has no context / + # trait column, so db$context is NULL and as.character(db$context[i]) would be + # character(0) -- length 0 -- which makes the data.frame() below error on any + # region with signal. Substitute a length-1 NA when the column is absent. + .col1 <- function(nm) { v <- db[[nm]]; if (is.null(v)) NA_character_ else as.character(v[i]) } + sel <- list(study = .col1("study"), context = .col1("context"), + trait = .col1("trait"), method = .col1("method")) pip <- tryCatch(do.call(getPip, c(list(db), sel)), error = function(e) NULL) if (is.null(pip)) return(NULL) data.frame(pip_sum = sum(pip[pip > 0], na.rm = TRUE), condition = sel$context,