Fix add_sample_from_csv crash on unmatched transcript/gene ids (2.0.7)

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.0.7]
+
+* fixed: `add_sample_from_csv` crashed with a confusing `AttributeError` when one transcript_id/gene_id from the coverage csv wasn't found in the transcripts file; the warning now also suggests `infer_genes=True` when the file has no gene annotations (#25)
+
 ## [2.0.6]
 
 * fixed: `has_overlap`/`get_overlap` gave wrong results for reverse-strand features, where genomic coordinates are sometimes passed as `(end, start)` instead of `(start, end)`; this caused `add_domains_to_table` to silently miss overlapping domains on the reverse strand (#23)

VERSION.txt

@@ -1 +1 @@
-2.0.6
+2.0.7

src/isotools/_transcriptome_io.py

@@ -275,18 +275,31 @@ def add_sample_from_csv(
 
     if "gene_id" not in cov_tab:
         gene_id_dict = {tid: gid for gid, tids in transcripts.items() for tid in tids}
-        try:
-            cov_tab["gene_id"] = [gene_id_dict[tid] for tid in cov_tab.transcript_id]
-        except KeyError as e:
+        cov_tab["gene_id"] = [gene_id_dict.get(tid) for tid in cov_tab.transcript_id]
+        missing = cov_tab["gene_id"].isna()
+        if missing.any():
             logger.warning(
-                "transcript_id %s from csv file not found in gtf." % e.args[0]
+                "%d transcript_id(s) from %s not found in %s, e.g. %s",
+                missing.sum(),
+                coverage_csv_file,
+                transcripts_file,
+                cov_tab.loc[missing, "transcript_id"].iloc[0],
             )
+            cov_tab = cov_tab[~missing]
     if "chr" not in cov_tab:
         chrom_dict = {gid: chrom for chrom, gids in gene_infos.items() for gid in gids}
-        try:
-            cov_tab["chr"] = [chrom_dict[gid] for gid in cov_tab.gene_id]
-        except KeyError as e:
-            logger.warning("gene_id %s from csv file not found in gtf.", e.args[0])
+        cov_tab["chr"] = [chrom_dict.get(gid) for gid in cov_tab.gene_id]
+        missing = cov_tab["chr"].isna()
+        if missing.any():
+            logger.warning(
+                "%d gene_id(s) from %s not found in %s, e.g. %s -- if the file has no "
+                "gene annotations (only transcript/exon lines), try infer_genes=True",
+                missing.sum(),
+                coverage_csv_file,
+                transcripts_file,
+                cov_tab.loc[missing, "gene_id"].iloc[0],
+            )
+            cov_tab = cov_tab[~missing]
 
     used_transcripts = set()
     for _, row in cov_tab.iterrows():

tests/data/no_gene_lines_example.gtf

@@ -0,0 +1,4 @@
+chr2_part	PacBio	transcript	3250396	3250740	.	+	.	transcript_id "PB.7.1"; gene_id "PB.7"
+chr2_part	PacBio	exon	3250396	3250740	.	+	.	transcript_id "PB.7.1"; gene_id "PB.7";
+chr2_part	PacBio	transcript	3905348	3905724	.	-	.	transcript_id "PB.1.1"; gene_id "PB.1"
+chr2_part	PacBio	exon	3905348	3905724	.	-	.	transcript_id "PB.1.1"; gene_id "PB.1";

tests/data/no_gene_lines_example_coverage.csv

@@ -0,0 +1,4 @@
+transcript_id,sample1_coverage
+PB.7.1,5
+PB.1.1,3
+PB.999.1,1

tests/data_import_test.py

@@ -16,6 +16,25 @@ def test_import_gff():
     assert True
 
 
+def test_add_sample_from_csv_missing_gene_info():
+    # regression test for #25: a coverage csv row referencing a transcript_id
+    # not found in the transcripts file previously broke gene_id/chr column
+    # construction for *every* row, not just the unmatched one (an all-or-
+    # nothing list comprehension inside a try/except). PB.999.1 below is not
+    # in the gtf; the two real transcripts must still import successfully.
+    isoseq = Transcriptome.from_reference("tests/data/example.gff.gz")
+    id_map = isoseq.add_sample_from_csv(
+        "tests/data/no_gene_lines_example_coverage.csv",
+        transcripts_file="tests/data/no_gene_lines_example.gtf",
+        transcript_id_col="transcript_id",
+        reconstruct_genes=False,
+        infer_genes=True,
+        sep=",",
+    )
+    assert isoseq.n_transcripts == 2, "the two known transcripts should still import"
+    assert set(id_map) == {"PB.7", "PB.1"}
+
+
 def test_read_gff_progress_bar_plain_and_gzip():
     # regression test for #37: the progress bar's byte-position tracking
     # must work for both plain and gzip files. It previously crashed for