Test Ensembl transcripts without gene features

.gitignore

@@ -107,5 +107,6 @@ ENV/
 
 # big folders to ignore
 datasets/
-examples/
+examples/*
+!examples/ensembl_gtf/
 tests/

README.md

@@ -32,6 +32,9 @@ pip install .
 
 Use at your own risk! Not all GTF files are of the same structure.
 
+For a small Ensembl-format database with a transcript lacking a `gene` row,
+see [the example GTF, database, and BED file](examples/ensembl_gtf/README.md).
+
 # Annotator Example Usage:
 
 ```
@@ -191,4 +194,4 @@ differently formatted "attributes" terminologies (see the difference between a w
 Currently setting this flag to either 'ce10' or 'ce11' will assuming it is formatted the wormbase way. Otherwise
 the format will default to the 'gencode' nomenclature. To use another format, you will have to look at annotate_bed.get_keys(),
 which tells this package what keys to expect from the 'attributes' section of the GTF file.
--
+-

changelog.md

@@ -5,6 +5,11 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/)
 
 ## [Unreleased]
 
+### Added
+- A small Ensembl-format GTF, gffutils database, and BED example for
+  transcripts without a separate gene feature, including the transcript
+  named in issue #3. Tests verify the database and annotations.
+
 
 ## [0.1.0] - 2022-11-21
 ### Changed
@@ -76,4 +81,3 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/)
 - README now contains examples and default priorities/params
 
 [Unreleased]: https://github.com/byee4/annotator...HEAD
-

examples/ensembl_gtf/README.md

@@ -0,0 +1,43 @@
+# Small Ensembl GTF example
+
+This seven-feature GTF uses Ensembl's `gene_biotype` and
+`transcript_biotype` attributes. Neither transcript has a separate `gene`
+feature, which exercises [issue #3](https://github.com/byee4/annotator/issues/3).
+The coding transcript is illustrative. The second transcript uses the ID,
+coordinates, and exon structure returned by the
+[Ensembl GRCh37 lookup for ENST00000606034](https://grch37.rest.ensembl.org/lookup/id/ENST00000606034?expand=1),
+the transcript named in the issue; this is a curated subset, not the complete
+release 87 GTF.
+The committed `mini_ensembl.gtf.db` is a gffutils SQLite database built from
+`mini_ensembl.gtf`.
+
+From the repository root, rebuild the database with:
+
+```sh
+build_gffutils_db \
+  --annotation_file examples/ensembl_gtf/mini_ensembl.gtf \
+  --db_file examples/ensembl_gtf/mini_ensembl.gtf.db \
+  --force --disable_infer_genes --disable_infer_transcripts
+```
+
+Disabling inference keeps the database free of a synthetic `gene` feature.
+Run the example with:
+
+```sh
+annotator \
+  --input examples/ensembl_gtf/queries.bed \
+  --output /tmp/mini_ensembl.annotated.bed \
+  --gtfdb examples/ensembl_gtf/mini_ensembl.gtf.db \
+  --format ensembl
+```
+
+The `Genic Region` column (column 9) should read `CDS`, `five_prime_utr`,
+`intergenic`, and `noncoding_exon` for the four queries. The BED file uses
+zero-based starts; the GTF uses one-based starts. The source GTF and database
+are test fixtures, not a biological reference.
+
+To verify the committed database and annotations:
+
+```sh
+python -m pytest -q unit_tests/test_ensembl_example.py
+```

examples/ensembl_gtf/mini_ensembl.gtf

@@ -0,0 +1,8 @@
+## Synthetic Ensembl-style GTF for testing annotator; coordinates and IDs are illustrative.
+chr1	example	transcript	101	180	.	+	.	gene_id "ENSG_EXAMPLE_1"; transcript_id "ENST_EXAMPLE_1"; gene_name "EXAMPLE1"; gene_biotype "protein_coding"; transcript_biotype "protein_coding";
+chr1	example	exon	101	180	.	+	.	gene_id "ENSG_EXAMPLE_1"; transcript_id "ENST_EXAMPLE_1"; gene_name "EXAMPLE1"; gene_biotype "protein_coding"; transcript_biotype "protein_coding"; exon_number "1";
+chr1	example	five_prime_utr	101	110	.	+	.	gene_id "ENSG_EXAMPLE_1"; transcript_id "ENST_EXAMPLE_1"; gene_name "EXAMPLE1"; gene_biotype "protein_coding"; transcript_biotype "protein_coding";
+chr1	example	CDS	111	170	.	+	0	gene_id "ENSG_EXAMPLE_1"; transcript_id "ENST_EXAMPLE_1"; gene_name "EXAMPLE1"; gene_biotype "protein_coding"; transcript_biotype "protein_coding";
+chr1	example	three_prime_utr	171	180	.	+	.	gene_id "ENSG_EXAMPLE_1"; transcript_id "ENST_EXAMPLE_1"; gene_name "EXAMPLE1"; gene_biotype "protein_coding"; transcript_biotype "protein_coding";
+chr1	ensembl	transcript	931346	933431	.	-	.	gene_id "ENSG00000272512"; transcript_id "ENST00000606034"; gene_name "RP11-54O7.17"; gene_biotype "lincRNA"; transcript_biotype "lincRNA";
+chr1	ensembl	exon	931346	933431	.	-	.	gene_id "ENSG00000272512"; transcript_id "ENST00000606034"; gene_name "RP11-54O7.17"; gene_biotype "lincRNA"; transcript_biotype "lincRNA"; exon_number "1";

examples/ensembl_gtf/queries.bed

@@ -0,0 +1,4 @@
+chr1	119	120	example_cds	0	+
+chr1	104	105	example_5utr	0	+
+chr1	189	190	example_intergenic	0	+
+chr1	931399	931400	issue_3_transcript	0	-

unit_tests/test_ensembl_example.py

@@ -0,0 +1,58 @@
+"""Keep the committed Ensembl example database in sync with its source GTF."""
+
+from pathlib import Path
+
+import gffutils
+
+from annotator import annotate_bed
+from annotator.annotator import GENE_PRIORITY, TRANSCRIPT_PRIORITY
+from annotator.build_gffutils_db import build_db
+
+
+EXAMPLE = Path(__file__).resolve().parents[1] / "examples" / "ensembl_gtf"
+GTF = EXAMPLE / "mini_ensembl.gtf"
+DB = EXAMPLE / "mini_ensembl.gtf.db"
+BED = EXAMPLE / "queries.bed"
+
+
+def feature_records(db):
+    return sorted(
+        (
+            feature.featuretype,
+            feature.seqid,
+            feature.start,
+            feature.end,
+            feature.strand,
+            tuple(sorted((key, tuple(value)) for key, value in feature.attributes.items())),
+        )
+        for feature in db.all_features()
+    )
+
+
+def test_committed_ensembl_database_matches_gtf(tmp_path):
+    rebuilt_file = tmp_path / "rebuilt.gtf.db"
+    build_db(str(GTF), str(rebuilt_file))
+
+    committed = gffutils.FeatureDB(str(DB))
+    rebuilt = gffutils.FeatureDB(str(rebuilt_file))
+    assert feature_records(committed) == feature_records(rebuilt)
+    assert len(feature_records(committed)) == 7
+    assert list(committed.features_of_type("gene")) == []
+
+
+def test_ensembl_example_annotates_coding_utr_and_intergenic_queries(tmp_path):
+    output = tmp_path / "annotated.bed"
+    annotate_bed.annotate_bed(
+        str(DB), [str(BED)], [str(output)], True, [],
+        [row[:] for row in TRANSCRIPT_PRIORITY],
+        [row[:] for row in GENE_PRIORITY],
+        "ensembl", False, 0, 1,
+    )
+
+    rows = [line.rstrip("\n").split("\t") for line in output.read_text().splitlines()]
+    assert [(row[3], row[6], row[8]) for row in rows] == [
+        ("example_cds", "ENSG_EXAMPLE_1", "CDS"),
+        ("example_5utr", "ENSG_EXAMPLE_1", "five_prime_utr"),
+        ("example_intergenic", "intergenic", "intergenic"),
+        ("issue_3_transcript", "ENSG00000272512", "noncoding_exon"),
+    ]