Bugs in 1.23.0/1.26.0: getMCI_inner complement silently equals all genes (changes MCI p-values); BioTIP.wrap crashes by default; unreproducible parallel defaults

[Phoenix12580]

Hi, thank you for BioTIP — the tipping-point framework has been genuinely useful in our work. While systematically profiling and validating the full workflow, we ran into several bugs in the current code (verified on GitHub master @ fb8b2f2, version 1.23.0, and confirmed present in the 1.26.0 release tarball from Bioconductor; R 4.4.3 on Linux). They compound: a default BioTIP.wrap() call cannot run at all, and even when fixed, results are not reproducible run-to-run. Reporting together with minimal reproducers and suggested fixes.

1. getMCI_inner: the "complement" silently equals ALL genes (affects MCI permutation p-values)

Inside getMCI_inner (R/BioTIP_update_06232025.R, ~line 2213), randomL is built with lapply(names(countsL), ...), so at that point the list has no names:

randomL = lapply(names(countsL), function(x) countsL[[x]][random_id, ])  # unnamed list!
comple = lapply(names(countsL), function(x)
    subset(countsL[[x]], !row.names(countsL[[x]]) %in% row.names(randomL[[x]])))
names(randomL) = names(comple) = names(countsL)   # <- assigned two lines too late

In the comple expression, randomL[[x]] does exact name lookup on an unnamed list, returns NULL, so rownames %in% NULL is all FALSE and the negation keeps every gene. Consequently:

PCCo_avg[i] <- mean(abs(M[rownames(comple[[i]]), rownames(randomL[[i]])]))

averages over all genes vs. the random module (including the module's own diagonal 1s) instead of the complement — this changes the MCI null distribution and hence the empirical p-values used to filter CTS candidates.

Minimal reproducer:

countsL <- local({
  set.seed(7); G <- 60; cells <- 90
  mat <- matrix(rpois(G * cells, 2), nrow = G)
  rownames(mat) <- paste0("g", 1:G); colnames(mat) <- paste0("c", 1:cells)
  lapply(split(colnames(mat), rep(paste0("S", 1:3), each = 30)), function(x) mat[, x])
})
rid <- c(58, 34, 56, 25, 16)

randomL <- lapply(names(countsL), function(x) countsL[[x]][rid, ])  # no names yet
randomL[["S1"]]                      # NULL

comple <- lapply(names(countsL), function(x)
  subset(countsL[[x]], !row.names(countsL[[x]]) %in% row.names(randomL[[x]])))
nrow(comple[[1]])                    # 60  — expected 55

# fix: name the list first
names(randomL) <- names(countsL)
comple <- lapply(names(countsL), function(x)
  countsL[[x]][!rownames(countsL[[x]]) %in% rownames(randomL[[x]]), , drop = FALSE])
nrow(comple[[1]])                    # 55  — correct

Suggested fix: set names(randomL) immediately after creating it (or build comple positionally). Note this also affects simulationMCI, which calls getMCI_inner per permutation — published p-values obtained with the current behavior would change after the fix, so it may deserve a NEWS entry.

2. BioTIP.wrap is not exported, and the export list is stale

After the 2025-06-23 restructure, NAMESPACE was not regenerated: BioTIP.wrap is missing from the exports (users must resort to BioTIP:::BioTIP.wrap), while export(getReadthrough) still refers to a function that no longer exists in the package. Regenerating NAMESPACE with roxygen2 fixes both.

3. BioTIP.wrap crashes with its default arguments

The signature default is permutation.method = c('gene','both','sample') — a length-3 vector. The body then does if (!permutation.method %in% c("gene","both","sample")), which errors with "the condition has length > 1" before any computation. The default should presumably be 'gene' (and/or match.arg(permutation.method)).

4. Parallel defaults are not reproducible (no worker seeding)

optimize.sd_selection and simulation_Ic default to doParallel = TRUE and call sample() inside %dopar% blocks without clusterSetRNGStream() — outcomes depend on which iteration lands on which worker, and all workers share the same default RNG state. On our benchmark (Schiebinger et al. 2019 GSE115943 subset, 38 states), the number of CTS candidates varied between identical-looking runs. Two robust fixes: (a) draw all random numbers up front in the main process before distributing (then any parallel plan reproduces exactly), or (b) seed workers via clusterSetRNGStream().

5. Two smaller ones

• In BioTIP.wrap (~line 3375), getTopMCI.gene.maxsiz is a typo for getTopMCI.gene.maxsize — an undefined-object error raised at the very end of an otherwise successful run (the run prints "BioTIP is done" and then fails before returning).
• In getMCI (df-given branch, ~line 2078), PCC_avg[j] <- mean(abs(U)) averages a logical matrix (the fraction of upper-triangular cells, a constant) instead of mean(abs(tmp[U])). This branch is not reached by BioTIP.wrap defaults, but it looks unintended.

Environment

• BioTIP 1.23.0 (GitHub master @ fb8b2f2) and 1.26.0 release tarball; R 4.4.3, Ubuntu 22.04 (WSL2)
• Found during a bit-exact re-validation: our re-implementation of the documented behavior disagreed with the package output, which is how bug 1 surfaced

Happy to share our benchmark artifacts (locked fixture, seeds, before/after outputs) if that would help. Thanks again for the package!

---

Attribution: this report was drafted by ZCode (an AI coding agent, powered by the GLM model) as part of an automated refactoring-and-verification task on BioTIP, and reviewed & submitted by @Phoenix12580. All findings and reproducers were empirically verified against the package as installed.