diff --git a/CHANGELOG.md b/CHANGELOG.md index 4f1d383..f9e2b2e 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/VERSION.txt b/VERSION.txt index 7ec1d6d..3e3c2f1 100644 --- a/VERSION.txt +++ b/VERSION.txt @@ -1 +1 @@ -2.1.0 +2.1.1 diff --git a/src/isotools/transcriptome.py b/src/isotools/transcriptome.py index 8330fd6..0948c1e 100644 --- a/src/isotools/transcriptome.py +++ b/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 diff --git a/tests/transcriptome_test.py b/tests/transcriptome_test.py new file mode 100644 index 0000000..8cf8bbc --- /dev/null +++ b/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"