Warn on ambiguous gene name lookup (2.1.1)

CHANGELOG.md

@@ -6,6 +6,10 @@
 * planned new feature: during import of long reads, (optionally) correct for short exon alignment issues.
 * separate new read import and classification of isoforms.
 
+## [2.1.1]
+
+* fixed: `Transcriptome[...]` gene lookup by name silently returned an arbitrary gene when the name was shared by multiple genes (common for duplicated gene symbols); now warns clearly, separately from the existing gene id ambiguity check (which itself had a bug: it checked against the combined id+name index instead of ids alone) (#27)
+
 ## [2.1.0]
 
 * **breaking**: renamed `tpm` to `cpm` throughout the API (`transcript_table(cpm=..., cpm_pseudocount=...)`, `Gene.cpm()`, `estimate_cpm_threshold`, `_cpm`/`_sum_cpm` columns). The values were always counts per million (no transcript-length normalization) despite the "tpm" name -- correct for full-length long reads, where read count is already a direct proxy for molecule count, but mislabeled. No deprecated alias; update `tpm=` to `cpm=` in existing code.

VERSION.txt

@@ -1 +1 @@
-2.1.0
+2.1.1

src/isotools/transcriptome.py

@@ -201,14 +201,24 @@ def _extract_reference(self):
     def make_index(self):
         """Updates the index of gene names and ids (e.g. used by the the [] operator)."""
         idx = dict()
+        seen_ids = set()
+        seen_names = set()
         for gene in self:
-            if gene.id in idx:  # at least id should be unique - maybe raise exception?
+            if gene.id in seen_ids:
                 logger.warning(
-                    "%s seems to be ambigous: %s vs %s",
+                    "gene id %r is ambiguous: shared by multiple genes, e.g. %s vs %s",
                     gene.id,
                     str(idx[gene.id]),
                     str(gene),
                 )
+            seen_ids.add(gene.id)
+            if gene.name in seen_names:
+                logger.warning(
+                    "gene name %r is ambiguous: shared by multiple genes; "
+                    "looking it up by name will return an arbitrary one -- use the gene id instead",
+                    gene.name,
+                )
+            seen_names.add(gene.name)
             idx[gene.name] = gene
             idx[gene.id] = gene
         self._idx = idx

tests/transcriptome_test.py

@@ -0,0 +1,43 @@
+import logging
+from intervaltree import IntervalTree
+from isotools import Gene, Transcriptome
+
+
+def test_make_index_warns_on_ambiguous_name(caplog):
+    # regression test for #27: only gene *id* collisions were checked in
+    # make_index(); gene *name* collisions (common in real annotation --
+    # paralogs, duplicated symbols) were completely silent, with the
+    # second gene silently overwriting the first in the name->gene index.
+    gene1 = Gene(
+        0, 100, {"chr": "chr1", "strand": "+", "ID": "GENE1", "name": "DUP"}, None
+    )
+    gene2 = Gene(
+        200, 300, {"chr": "chr1", "strand": "+", "ID": "GENE2", "name": "DUP"}, None
+    )
+
+    with caplog.at_level(logging.WARNING, logger="isotools"):
+        Transcriptome(
+            data={"chr1": IntervalTree([gene1, gene2])},
+            infos={"reference_file": "test"},
+        )
+
+    messages = [r.getMessage() for r in caplog.records]
+    assert any(
+        "DUP" in m and "ambiguous" in m for m in messages
+    ), "expected a warning about the ambiguous gene name"
+
+    # id collisions must still be reported, and must not be confused with
+    # name collisions (they are tracked separately)
+    gene3 = Gene(
+        400, 500, {"chr": "chr1", "strand": "+", "ID": "GENE1", "name": "OTHER"}, None
+    )
+    caplog.clear()
+    with caplog.at_level(logging.WARNING, logger="isotools"):
+        Transcriptome(
+            data={"chr1": IntervalTree([gene1, gene3])},
+            infos={"reference_file": "test"},
+        )
+    messages = [r.getMessage() for r in caplog.records]
+    assert any(
+        "GENE1" in m and "ambiguous" in m for m in messages
+    ), "expected a warning about the ambiguous gene id"