Fix die_test crash on isoforms unexpressed in both compared groups (2.1.3)

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.3]
+
+* fixed: `die_test` crashed with `ValueError: The internally computed table of expected frequencies has a zero element` for genes with isoforms only covered in samples outside the two compared groups (#29)
+
 ## [2.1.2]
 
 * fixed: gff3 import silently dropped genes when the file's chromosome/seqid naming didn't match the genome (e.g. RefSeq-style accessions like `NC_000001.11` vs a plain `1`/`chr1` genome FASTA); chromosome name aliasing via `region` feature lines is now resolved again, restoring behavior lost when the tabix-based reader was replaced to fix #28 (#36)

VERSION.txt

@@ -1 +1 @@
-2.1.2
+2.1.3

src/isotools/gene.py

@@ -1322,6 +1322,14 @@ def die_test(self, groups, min_cov=25, n_isoforms=10):
 
         if np.any(cov.sum(0) < min_cov):
             return np.nan, np.nan, []
+        # isoforms of the gene that are only covered in samples outside the two
+        # groups being compared have zero reads in both groups here; such an
+        # all-zero row makes chi2_contingency's expected-frequency table singular
+        # (raises ValueError "expected frequencies has a zero element") -- drop
+        # them first, keeping track of their original indices for the returned ids
+        expressed = cov.sum(1) > 0
+        orig_idx = np.flatnonzero(expressed)
+        cov = cov[expressed]
         # if there are more than 'numIsoforms' isoforms of the gene, all additional least expressed get summarized.
         if cov.shape[0] > n_isoforms:
             idx = np.argpartition(
@@ -1330,11 +1338,12 @@ def die_test(self, groups, min_cov=25, n_isoforms=10):
             additional = cov[idx[n_isoforms:]].sum(0)
             cov = cov[idx[:n_isoforms]]
             cov[n_isoforms - 1] += additional
+            idx = orig_idx[idx]
             idx[n_isoforms - 1] = -1  # this isoform gets all other - I give it index
         elif cov.shape[0] < 2:
             return np.nan, np.nan, []
         else:
-            idx = np.array(range(cov.shape[0]))
+            idx = orig_idx
         try:
             _, pval, _, _ = chi2_contingency(cov)
         except ValueError:

tests/die_test_test.py

@@ -0,0 +1,30 @@
+import numpy as np
+from isotools import Gene
+
+
+def test_die_test_ignores_isoforms_unexpressed_in_both_groups():
+    # regression test for #29: a gene's isoform list is the union of isoforms
+    # seen across *all* samples in the Transcriptome, not just the two groups
+    # being compared. An isoform only covered in a sample outside both groups
+    # ends up as an all-zero row here, which made chi2_contingency raise
+    # ValueError("expected frequencies has a zero element") since a zero row
+    # makes the expected-frequency table singular. min_cov (checked on group
+    # totals) does not catch this, since the other isoforms still carry
+    # enough coverage.
+    coverage = np.array(
+        [
+            [706, 3, 0, 2, 0, 0, 0, 0],  # group 0 samples
+            [1218, 0, 0, 0, 1, 0, 0, 0],  # group 1 samples
+        ]
+    )
+    gene = Gene(
+        100,
+        500,
+        {"chr": "chr1", "strand": "+", "ID": "GENE1", "coverage": coverage},
+        None,
+    )
+    pval, deltaPI, transcript_ids = gene.die_test([[0], [1]], min_cov=1)
+    assert np.isfinite(pval)
+    assert np.isfinite(deltaPI)
+    # the zero-coverage isoforms (original indices 2, 5, 6, 7) must not appear
+    assert set(transcript_ids).issubset({0, 1, 3, 4})