diff --git a/CHANGELOG.md b/CHANGELOG.md index d425ddf..494d5a4 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.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) diff --git a/VERSION.txt b/VERSION.txt index 157e54f..f1547e6 100644 --- a/VERSION.txt +++ b/VERSION.txt @@ -1 +1 @@ -2.0.6 +2.0.7 diff --git a/src/isotools/_transcriptome_io.py b/src/isotools/_transcriptome_io.py index 1c95d62..6ed4c34 100644 --- a/src/isotools/_transcriptome_io.py +++ b/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(): diff --git a/tests/data/no_gene_lines_example.gtf b/tests/data/no_gene_lines_example.gtf new file mode 100644 index 0000000..976c05d --- /dev/null +++ b/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"; diff --git a/tests/data/no_gene_lines_example_coverage.csv b/tests/data/no_gene_lines_example_coverage.csv new file mode 100644 index 0000000..b92d1f9 --- /dev/null +++ b/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 diff --git a/tests/data_import_test.py b/tests/data_import_test.py index f100b9b..9f8d0b2 100644 --- a/tests/data_import_test.py +++ b/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