Add bridge ortholog mapping

[ivokwee]

Add bridge orthology search for poorly annotated species. Bridge orthology search uses an intermediate species (BRIDGE_SPECIES = dmelanogaster, drerio, celegans, scerevisiae, athalian) to map an organism first, then using the bridge species maps to human orthology next. This can improve orthology mapping for poorly annotated species (like Fall armyworm). Better orthology mapping is useful for geneset analysis, drug connectivity and experiment similarity in OPG.

[ESCRI11]

current state of PR does not even parse

Error in FUN(X[[i]], ...) : 
  /playbase/R/pgx-annot.R:1443:1: unexpected '}'
1442:   return(top_match)
1443: }
      ^

please fix before i can try to run it and review

[ivokwee]

@ESCRI11 ok. try again.

[ESCRI11]

Review: bridge ortholog mapping

Nice idea — bridging through a better-annotated intermediate is the right move for species g:Profiler barely knows. I ran the branch against real datasets rather than reviewing statically, so everything below is reproducible. All reprexes need a clone of bigomics/opg-exampledata (@ c453f3d); datasets used are tomatoe/, trout-uniprot/, hamster/ and human-symbol/.

Three bugs I'd like fixed before merge, plus one design question about whether the feature activates at all.

---

1. getHumanOrtholog() crashes on any single-gene input — R/pgx-annot.R:1067

apply(df, 2, df.concat) simplifies to a plain vector when the input collapses to one unique symbol, so df2$human is NULL, mean.mapped at :1073 is NaN, and the if at :1075 errors.

suppressMessages(devtools::load_all("path/to/playbase", quiet = TRUE))

## gene ID from opg-exampledata/tomatoe/counts.csv, row 1
getHumanOrtholog("Solanum lycopersicum", "Solyc02g050330", verbose = 0)
#> WARNING: bridge conversion with dmelanogaster failed
#> WARNING: bridge conversion with drerio failed
#> WARNING: bridge conversion with celegans failed
#> WARNING: bridge conversion with scerevisiae failed
#> WARNING: bridge conversion with athaliana failed
#> Error in `if (mean.mapped == 0) ...`:
#> ! missing value where TRUE/FALSE needed

## not tomato-specific -- any input collapsing to ONE unique symbol does it
## (protein ID from opg-exampledata/trout-uniprot/counts.csv)
getHumanOrtholog("Oncorhynchus mykiss", "A0A060PIE3", verbose = 0)
#> Error in `if (mean.mapped == 0) ...`:
#> ! missing value where TRUE/FALSE needed

apply(df, 2, df.concat, simplify = FALSE) fixes it.

---

2. human_ortholog becomes a ;-joined list — R/pgx-annot.R:1105 and :1272

gorth(mthreshold = Inf) is one-to-many, and the two paths resolve that ambiguity in two different ways — neither of which downstream code can consume:

• Direct path (res[jj,] <- out[,], :1105): every duplicated input matches the same row index, so the assignment silently keeps whichever candidate came last. A0A060VMM0 has 9 candidate orthologs and arrives as AKR1C8 for no reason other than row order.
• Bridged path (:1272): pastes them all with collapse = ";".

suppressMessages(devtools::load_all("path/to/playbase", quiet = TRUE))

## protein IDs from opg-exampledata/trout-uniprot/counts.csv
g <- c("A0A060VMM0", "A0A060VMQ9", "A0A060PIE3")

## gprofiler2 is one-to-many: 9 candidate orthologs for A0A060VMM0
playbase:::.query_orthologs(g, "omykiss", method = "gprofiler2", verbose = 0)
#>         input ortholog
#> 1  A0A060VMM0   AKR7A2
#> 2  A0A060VMM0   AKR1C2
#> 3  A0A060VMM0   AKR7A3
#> 4  A0A060VMM0   AKR1C1
#> 5  A0A060VMM0   AKR1C3
#> 6  A0A060VMM0  AKR1B10
#> 7  A0A060VMM0   AKR1C4
#> 8  A0A060VMM0  AKR1B15
#> 9  A0A060VMM0   AKR1C8
#> 10 A0A060VMQ9     ADH6
#> 11 A0A060VMQ9    ADH1A
#> ...
#> 17 A0A060PIE3     BTF3

## ...but convert_orthologs keeps only the LAST row of each duplicated input
playbase:::convert_orthologs(g, "omykiss", methods = "gprofiler2", bridge = FALSE, verbose = 0)
#>        input ortholog     method
#> 1 A0A060VMM0   AKR1C8 gprofiler2
#> 2 A0A060VMQ9    ADH1C gprofiler2
#> 3 A0A060PIE3     BTF3 gprofiler2

## the bridged path collapses the same ambiguity the other way, with ";"
b <- playbase:::.query_orthologs_bridged(g, "omykiss", bridge_species = "drerio", verbose = 0)
b
#>                 input      bridge                              ortholog
#> A0A060VMM0 A0A060VMM0        <NA>                                  <NA>
#> A0A060VMQ9 A0A060VMQ9 adh8a;adh8b ADH1A;ADH1B;ADH1C;ADH4;ADH5;ADH6;ADH7
#> A0A060PIE3 A0A060PIE3        <NA>                                  <NA>

## which is what lands in pgx$genes$human_ortholog -- and breaks exact-match consumers
playdata::GENE_SUMMARY[b$ortholog[2]]
#> <NA>
#>   NA

That ;-joined value propagates to annot$ORTHOGENE (pgx-annot.R:512), df$human_ortholog (:641) and pgx$genes$human_ortholog (pgx-init.R:210), and every downstream consumer matches on it exactly:

• playdata::GENE_SUMMARY[annot$human_ortholog] → NA (pgx-annot-info.R:39)
• data.table::chmatch(rownames(X), pgx$genes$human_ortholog) → no match (compute2-extra.R:500)
• rename_by2(..., new_id = "human_ortholog") writes "ADH1A;ADH1B;..." as a rowname into the tiledb / drug-connectivity matrices (pgx-tiledb.R:152)

Both paths need to agree. Either pick one symbol deterministically (and document the rule), or keep human_ortholog single-valued and add a separate human_orthologs column for the full set — a column that is sometimes one symbol and sometimes seven can't be consumed either way.

---

3. cluster_columns fix is only correct when splitx = NULL — R/gx-heatmap.r:429

The relaxed is.logical() guard now admits a precomputed hclust/dendrogram, but the per-group loop at :610 subsets the columns (gx0 <- gx[, jj, drop = FALSE]) and then hands the unmodified whole-matrix object straight to Heatmap() at :643.

suppressMessages(devtools::load_all("path/to/playbase", quiet = TRUE))
D <- "path/to/opg-exampledata"

## opg-exampledata/human-symbol : 18 samples, 6 groups
counts  <- as.matrix(read.csv(file.path(D, "human-symbol/counts.csv"), row.names = 1))
samples <- read.csv(file.path(D, "human-symbol/samples.csv"), row.names = 1)

gx <- log2(1 + counts)
gx <- gx[head(order(-apply(gx, 1, sd)), 50), ]   # 50 most variable genes
splitx <- factor(samples[colnames(gx), "group"])

dim(gx)
#> [1] 50 18
table(splitx)
#> splitx
#> act12h act24h act48h act72h act96h notact
#>      3      1      2      3      2      7

## a precomputed column dendrogram over ALL samples -- the thing the new
## is.logical() guard at R/gx-heatmap.r:429 was relaxed to allow
hc <- hclust(dist(t(gx)), method = "ward.D2")
length(hc$order)
#> [1] 18

## (A) splitx = NULL -> single Heatmap over all 18 columns: works, renders fine
r <- try(gx.splitmap(gx, splitx = NULL, cluster_columns = hc), silent = TRUE)
inherits(r, "try-error")
#> [1] FALSE

## (B) splitx set -> gx-heatmap.r:610 subsets columns per group,
##     gx-heatmap.r:643 still passes the FULL 18-leaf dendrogram to Heatmap(gx0)
try(gx.splitmap(gx, splitx = splitx, cluster_columns = hc))
#> Error : The length of the column clustering object is not the same as the
#> number of matrix columns

## the mismatch, stated directly:
sapply(split(colnames(gx), splitx), length)  # columns each per-group Heatmap() gets
#> act12h act24h act48h act72h act96h notact
#>      3      1      2      3      2      7
length(hc$order)                             # leaves in the dendrogram handed to it
#> [1] 18

Either subset/reorder the dendrogram per group, or keep rejecting non-logical cluster_columns when splitx is set.

---

Design question: the bridge never fires on the species it was built for

This is the one I'd most like your read on. The auto-trigger is bridge <- (orth.ratio < 0.10) (:1116), but the actual poorly-annotated datasets sit above that threshold — so they leave most features unmapped and never attempt a bridge:

suppressMessages(devtools::load_all("path/to/playbase", quiet = TRUE))
D <- "path/to/opg-exampledata"

## trout-uniprot: 5312 proteins, Oncorhynchus mykiss
g <- rownames(read.csv(file.path(D, "trout-uniprot/counts.csv"), row.names = 1))
res <- getHumanOrtholog("Oncorhynchus mykiss", g, verbose = 1)
#> [getHumanOrtholog] processing 5590 genes in 4 batches
#> [convert_orthologs] pass 1: orth.ratio = 0.183333333333333
#> [convert_orthologs] pass 1: orth.ratio = 0.164666666666667
#> [convert_orthologs] pass 1: orth.ratio = 0.109333333333333
#> [convert_orthologs] pass 1: orth.ratio = 0.132110091743119
#> [getHumanOrtholog] total ratio mapped  = 15.6627%
c(features = length(g), mapped = sum(!is.na(res$human)), bridged = sum(grepl("bridged", res$source)))
#> features   mapped  bridged
#>     5312      832        0

## hamster: Cricetulus griseus
gh <- rownames(read.csv(file.path(D, "hamster/counts.csv"), row.names = 1, nrows = 400))
rh <- getHumanOrtholog("Cricetulus griseus", gh, verbose = 1)
#> [convert_orthologs] pass 1: orth.ratio = 0.324189526184539
#> [getHumanOrtholog] total ratio mapped  = 32.5%
c(features = length(gh), mapped = sum(!is.na(rh$human)), bridged = sum(grepl("bridged", rh$source)))
#> features   mapped  bridged
#>      400      130        0

## forcing the bridge on the same trout genes DOES add mappings -- the gate is the problem
table(playbase:::convert_orthologs(g[1:300], "omykiss", bridge = TRUE, verbose = 0)$method)
#> WARNING: bridge conversion with celegans failed
#> WARNING: bridge conversion with scerevisiae failed
#> WARNING: bridge conversion with athaliana failed
#>
#> bridged:dmelanogaster        bridged:drerio             gprofiler
#>                     3                     1                    47
#>            gprofiler2
#>                     4

Across trout-uniprot, hamster and celegans-wormbase the bridge contributed orthologs in zero default runs. Trout finishes at 15.7% mapped — 4480 of 5312 features unmapped — without a single bridge attempt, because 15.7% > 10%. The machinery works when forced (4 extra mappings on 300 genes), so this is the gate, not the implementation.

Two related observations:

• The threshold is evaluated per batch. :1116 runs inside each convert_orthologs() call, and trout's four batches came back 0.183 / 0.165 / 0.109 / 0.132 — one batch is a hair from bridging while the others aren't. The same dataset can bridge inconsistently depending on feature order. Worth deciding bridge once, on a sample, above the batching loop.
• tomatoe triggers the bridge but can't use it. g:Profiler doesn't recognise Solyc* IDs in any namespace and returns N/A for every gene even to athaliana, so all five bridges fail and the result is 0 mapped. Might be worth detecting an unrecognised-namespace input rather than fanning out five failing round-trips.

Would it make more sense to trigger on absolute unmapped count, or to bridge whenever orth.ratio is below something much higher (~0.5)?

[ESCRI11]

@ivokwee commisioned claude review and steered it manually until reproducible valid bugs where found