From 1d7c51963a7d29b88e33bac2821766f76af589da Mon Sep 17 00:00:00 2001 From: jfallmann Date: Wed, 21 Oct 2026 20:11:43 +0200 Subject: [PATCH 1/4] workound for missing gene_id in gtf in case of e.g. ERCC spikes --- scripts/lib/_lib.R | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/scripts/lib/_lib.R b/scripts/lib/_lib.R index 18c6de36..f5db3426 100644 --- a/scripts/lib/_lib.R +++ b/scripts/lib/_lib.R @@ -8,13 +8,16 @@ get_gene_name <- function(id, df) { message("WARNING: gene_name not found as colname, will be replaced by gene column, please make sure the gtf file is in the correct format") df$gene_name <- df$gene } - name_list <- df$gene_name[df["type"] == "gene" & df["gene_id"] == id] - if (length(unique(name_list)) == 1) { - return(name_list[1]) - } else { - message(paste("WARNING: ambigous gene id: ", id)) - return(paste(unique(name_list), sep = "|")) + name_list <- unique(df$gene_name[df["type"] == "gene" & df["gene_id"] == id]) + if (length(name_list) == 0) { + message(paste("WARNING: gene id not found in annotation, using id as name: ", id)) + return(id) + } + if (length(name_list) == 1) { + return(name_list) } + message(paste("WARNING: ambigous gene id: ", id)) + return(paste(name_list, collapse = "|")) } @@ -27,13 +30,16 @@ get_exon_name <- function(id, df) { message("WARNING: gene_name not found as colname, will be replaced by gene column, please make sure the gtf file is in the correct format") df$gene_name <- df$gene } - name_list <- df$gene_name[df["type"] == "exon" & df["gene_id"] == id] - if (length(unique(name_list)) == 1) { - return(name_list[1]) - } else { - message(paste("WARNING: ambigous gene id: ", id)) - return(paste(unique(name_list), sep = "|")) + name_list <- unique(df$gene_name[df["type"] == "exon" & df["gene_id"] == id]) + if (length(name_list) == 0) { + message(paste("WARNING: gene id not found in annotation, using id as name: ", id)) + return(id) + } + if (length(name_list) == 1) { + return(name_list) } + message(paste("WARNING: ambigous gene id: ", id)) + return(paste(name_list, collapse = "|")) } get_gene_coords <- function(id, df) { From 706ce0d8c5ee18b93626e32ce5f2d68d4e9a3cbd Mon Sep 17 00:00:00 2001 From: jfallmann Date: Wed, 21 Oct 2026 22:53:05 +0200 Subject: [PATCH 2/4] length fix for missing gene ids, tpms no longer na --- scripts/lib/_lib.R | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/scripts/lib/_lib.R b/scripts/lib/_lib.R index f5db3426..eae96f09 100644 --- a/scripts/lib/_lib.R +++ b/scripts/lib/_lib.R @@ -67,6 +67,10 @@ add_gene_coordinates <- function(df, gene_ids, gtf_df, after = NULL) { colnames(gtf_df)[grepl("id$", names(gtf_df), ignore.case = TRUE)][1] <- "gene_id" } df <- as.data.frame(df) + if (length(gene_ids) != nrow(df)) { + message(paste0("WARNING: gene_ids length (", length(gene_ids), ") does not match data rows (", nrow(df), "); falling back to rownames(df)")) + gene_ids <- rownames(df) + } g <- gtf_df[gtf_df["type"] == "gene", ] coords <- paste(g$seqnames, g$start, g$end, g$strand, sep = ":") names(coords) <- as.character(g$gene_id) @@ -102,10 +106,18 @@ calc_tpm <- function(counts, gtf) { gene_lengths <- gtf$width names(gene_lengths) <- gtf$gene_id matched_lengths <- gene_lengths[rownames(counts)] - gene_lengths_kb <- matched_lengths / 1000 - rpk <- counts / gene_lengths_kb + # ids without a length (e.g. ERCC spikes that only have exon rows in the GTF) + # must not enter the per-sample scaling sum, otherwise one NA poisons the whole column + has_len <- !is.na(matched_lengths) + if (any(!has_len)) { + message(paste0("WARNING: no gene length for ", sum(!has_len), " features, excluded from TPM (set to NA): ", paste(head(rownames(counts)[!has_len], 5), collapse = ", "), if (sum(!has_len) > 5) ", ..." else "")) + } + tpm <- matrix(NA_real_, nrow = nrow(counts), ncol = ncol(counts), dimnames = dimnames(counts)) + counts_l <- counts[has_len, , drop = FALSE] + gene_lengths_kb <- matched_lengths[has_len] / 1000 + rpk <- counts_l / gene_lengths_kb scaling_factors <- colSums(rpk) - tpm <- t(t(rpk) / scaling_factors * 1e6) + tpm[has_len, ] <- t(t(rpk) / scaling_factors * 1e6) return(tpm) } @@ -220,7 +232,7 @@ format_deseq2_results <- function(res, gtf_gene, shrink = TRUE, center = FALSE) message(paste0("Global shift (median log2FoldChange) subtracted in log2FoldChange_centered: ", shift)) res$log2FoldChange_centered <- res$log2FoldChange - shift } - add_gene_coordinates(res, res$Gene_ID, gtf_gene, after = "Gene_ID") + add_gene_coordinates(res, rownames(res), gtf_gene, after = "Gene_ID") } ## Format an edgeR result table (qlf$table or a topTags table) for export. From 03e11db34cd6eef2e5724234bf587188007486e7 Mon Sep 17 00:00:00 2001 From: jfallmann Date: Thu, 22 Oct 2026 13:22:17 +0200 Subject: [PATCH 3/4] sequential write to postde.rds to prevent OOM issues --- scripts/Analysis/DE/DESEQ2.R | 4 ++++ scripts/Analysis/DE/EDGER.R | 4 ++++ scripts/Analysis/PostDE/export.R | 36 ++++++++++++++++++++++++++++---- 3 files changed, 40 insertions(+), 4 deletions(-) diff --git a/scripts/Analysis/DE/DESEQ2.R b/scripts/Analysis/DE/DESEQ2.R index d4c15471..d8a694ea 100755 --- a/scripts/Analysis/DE/DESEQ2.R +++ b/scripts/Analysis/DE/DESEQ2.R @@ -55,6 +55,10 @@ if (postde_enabled) { } source(postde_lib) postde_bundles <- list() + postde_entry_ids <- character(0) + postde_cache_dir <- file.path(outdir, ".postde_cache") + dir.create(postde_cache_dir, recursive = TRUE, showWarnings = FALSE) + postde_cache_dir <- normalizePath(postde_cache_dir) } ## set thread-usage diff --git a/scripts/Analysis/DE/EDGER.R b/scripts/Analysis/DE/EDGER.R index 6beecf7d..52bc69f9 100755 --- a/scripts/Analysis/DE/EDGER.R +++ b/scripts/Analysis/DE/EDGER.R @@ -50,6 +50,10 @@ if (postde_enabled) { } source(postde_lib) postde_bundles <- list() + postde_entry_ids <- character(0) + postde_cache_dir <- file.path(outdir, ".postde_cache") + dir.create(postde_cache_dir, recursive = TRUE, showWarnings = FALSE) + postde_cache_dir <- normalizePath(postde_cache_dir) } ## plotMDS needs at least 3 samples/columns, fall back to a placeholder plot otherwise diff --git a/scripts/Analysis/PostDE/export.R b/scripts/Analysis/PostDE/export.R index 26ea063f..368a56bf 100644 --- a/scripts/Analysis/PostDE/export.R +++ b/scripts/Analysis/PostDE/export.R @@ -29,6 +29,10 @@ if (!is.null(postde_common) && file.exists(postde_common)) { } postde_bundles <- list() +postde_entry_ids <- character(0) +if (!exists("postde_cache_dir")) { + postde_cache_dir <- NULL +} postde_results_deseq2 <- function(res) { out <- data.frame( @@ -110,25 +114,49 @@ postde_capture <- function(engine, id, A, B, normalized, metadata, counts, expre results = results, mean_scale = mean_scale ) - postde_bundles[[length(postde_bundles) + 1L]] <<- entry - invisible(entry) + # Spill each entry to its own file immediately rather than keeping every + # contrast's counts/expression matrices resident for the whole run: with + # dozens of contrasts (and a _norm variant for each) the accumulated list + # has been large enough to exhaust job memory before the final save. Only + # the small id and file path are kept in memory; postde_write() streams + # the files back together one at a time for the single combined saveRDS. + if (!is.null(postde_cache_dir)) { + san_id <- postde_sanitize_id(id) + entry_path <- file.path(postde_cache_dir, paste0(san_id, ".rds")) + saveRDS(entry, entry_path) + postde_entry_ids <<- c(postde_entry_ids, id) + postde_bundles[[length(postde_bundles) + 1L]] <<- entry_path + rm(entry) + } else { + postde_entry_ids <<- c(postde_entry_ids, id) + postde_bundles[[length(postde_bundles) + 1L]] <<- entry + } + invisible(NULL) } postde_write <- function(outdir, combi, engine) { if (length(postde_bundles) == 0L) { stop("MONSDA_POSTDE=1 but no contrast entries were captured; refusing to write an empty bundle") } - ids <- vapply(postde_bundles, function(e) e$id, character(1)) + ids <- postde_entry_ids if (anyDuplicated(ids)) { stop(paste0("postde: duplicate entry ids: ", paste(ids[duplicated(ids)], collapse = ", "))) } postde_check_sanitized_ids(ids) + contrasts <- setNames(vector("list", length(postde_bundles)), ids) + for (i in seq_along(postde_bundles)) { + item <- postde_bundles[[i]] + contrasts[[i]] <- if (is.character(item)) readRDS(item) else item + } bundle <- list( schema_version = 1L, engine = engine, - contrasts = setNames(postde_bundles, ids) + contrasts = contrasts ) path <- file.path(outdir, paste0("DE_", engine, "_", combi, "_postde.rds")) saveRDS(bundle, path) + if (!is.null(postde_cache_dir)) { + unlink(postde_cache_dir, recursive = TRUE) + } invisible(path) } From c1ce9a2770193cd46c6e29240a96debc66ad9424 Mon Sep 17 00:00:00 2001 From: jfallmann Date: Thu, 22 Oct 2026 13:25:11 +0200 Subject: [PATCH 4/4] add extra option to build_all_local to enable single container builds --- containers/apptainer/build_all_local.sh | 33 ++++++++++++++++++++++++- 1 file changed, 32 insertions(+), 1 deletion(-) diff --git a/containers/apptainer/build_all_local.sh b/containers/apptainer/build_all_local.sh index 6375390f..aad72704 100755 --- a/containers/apptainer/build_all_local.sh +++ b/containers/apptainer/build_all_local.sh @@ -1,12 +1,43 @@ #!/usr/bin/env bash set -uo pipefail +usage() { + printf 'Usage: %s [-e environment[,environment...]] [repo] [out]\n' "$(basename "${BASH_SOURCE[0]}")" >&2 + printf ' -e Build only the given environment(s) instead of the full matrix.\n' >&2 + printf ' May be repeated or given as a comma-separated list, e.g.\n' >&2 + printf ' -e deseq2_DE -e edger_DE or -e deseq2_DE,edger_DE\n' >&2 +} + +only_environments=() +while getopts ":e:h" opt; do + case "$opt" in + e) IFS=',' read -ra envs <<< "$OPTARG"; only_environments+=("${envs[@]}") ;; + h) usage; exit 0 ;; + \?) printf 'Unknown option: -%s\n' "$OPTARG" >&2; usage; exit 1 ;; + :) printf 'Option -%s requires an argument\n' "$OPTARG" >&2; usage; exit 1 ;; + esac +done +shift $((OPTIND - 1)) + repo=${1:-${HOME}/Gits/Monsda} out=${2:-${HOME}/Work2/Container/MONSDA} script_dir=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd) version=$(cd "$repo" && python scripts/container_matrix.py --version | cut -d= -f2) revision=$(cd "$repo" && git rev-parse HEAD) -mapfile -t environments < <(cd "$repo" && python scripts/container_matrix.py | python -c 'import json,sys; print("\n".join(json.load(sys.stdin)["environment"]))') +if (( ${#only_environments[@]} > 0 )); then + mapfile -t all_environments < <(cd "$repo" && python scripts/container_matrix.py | python -c 'import json,sys; print("\n".join(json.load(sys.stdin)["environment"]))') + environments=() + for environment in "${only_environments[@]}"; do + [[ -z "$environment" ]] && continue + if [[ ! " ${all_environments[*]} " == *" $environment "* ]]; then + printf 'Unknown environment: %s (available: %s)\n' "$environment" "${all_environments[*]}" >&2 + exit 1 + fi + environments+=("$environment") + done +else + mapfile -t environments < <(cd "$repo" && python scripts/container_matrix.py | python -c 'import json,sys; print("\n".join(json.load(sys.stdin)["environment"]))') +fi mkdir -p "$out/logs" "$out/tmp" printf 'environment\tstatus\timage\n' > "$out/results.tsv"