diff --git a/.gitignore b/.gitignore index 25854a1..56c9597 100644 --- a/.gitignore +++ b/.gitignore @@ -107,5 +107,6 @@ ENV/ # big folders to ignore datasets/ -examples/ +examples/* +!examples/ensembl_gtf/ tests/ diff --git a/README.md b/README.md index 69e44b4..07e8a4d 100644 --- a/README.md +++ b/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. -- \ No newline at end of file +- diff --git a/changelog.md b/changelog.md index 16e2836..868e0b2 100644 --- a/changelog.md +++ b/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 - diff --git a/examples/ensembl_gtf/README.md b/examples/ensembl_gtf/README.md new file mode 100644 index 0000000..7576b4e --- /dev/null +++ b/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 +``` diff --git a/examples/ensembl_gtf/mini_ensembl.gtf b/examples/ensembl_gtf/mini_ensembl.gtf new file mode 100644 index 0000000..c812a67 --- /dev/null +++ b/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"; diff --git a/examples/ensembl_gtf/mini_ensembl.gtf.db b/examples/ensembl_gtf/mini_ensembl.gtf.db new file mode 100644 index 0000000..9d46397 Binary files /dev/null and b/examples/ensembl_gtf/mini_ensembl.gtf.db differ diff --git a/examples/ensembl_gtf/queries.bed b/examples/ensembl_gtf/queries.bed new file mode 100644 index 0000000..f2ab79a --- /dev/null +++ b/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 - diff --git a/unit_tests/test_ensembl_example.py b/unit_tests/test_ensembl_example.py new file mode 100644 index 0000000..0b2bed3 --- /dev/null +++ b/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"), + ]