Dev hotfix

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"
 

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

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

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)
 }

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) {
@@ -61,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)
@@ -96,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)
 }
 
@@ -214,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.