Handle alternate GTF biotype keys

README.md

@@ -42,7 +42,7 @@ annotator \
 --input BED6_FILE \
 --output OUTPUT_FILE \
 --gtfdb gencode.v19.annotation.gtf.db \
---species hg19
+--format gencode
 ```
 
 ### Output file:
@@ -119,6 +119,12 @@ however if strand is specified in the BED file, we'll try to look for
 correctly stranded features first. If strand is not specified, we'll
 prioritize positive stranded features first.
 
+```--format``` selects the expected GTF attribute names. Use `gencode` for
+`gene_type`/`transcript_type` files and `ensembl` for
+`gene_biotype`/`transcript_biotype` files. If the selected format and the
+database use opposite names for these two biotypes, annotator also checks the
+alternate names.
+
 ```--limit-chroms-to``` will limit the dictionary build to only
 include these chromosomes for faster processing and less memory
 footprint. Leave blank to hash all chromosomes in the db file
@@ -189,9 +195,7 @@ where MATURE is an optional column linking the precursor accession to its proces
 - You can use the ```build_gffutils_db``` script to create a gffdb from a gff file downloaded from mirbase.
 
 # General notes:
-- The ```--species``` flag is only important for setting the GTF file nomenclature; different GTF/GFF files have
-differently formatted "attributes" terminologies (see the difference between a wormbase.org and a gencode GTF file).
-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.
--
+- `annotator --format` selects the GTF attribute convention. The `ce10` and
+  `ce11` values use WormBase keys; the default is GENCODE. Other commands may
+  still use `--species`. For additional formats, see
+  `annotation_functions.get_keys()`.

annotator/annotate_bed.py

@@ -651,14 +651,30 @@ def annotate(
                 to_append = to_append[:-1] + delim
             else:
                 to_append = to_append + '-' + delim
-            if keys['transcript_type'] in feature.attributes.keys():
-                for t in feature.attributes[keys['transcript_type']]:
+            transcript_types = feature.attributes.get(keys['transcript_type'])
+            if not transcript_types:
+                alternate_key = (
+                    'transcript_biotype'
+                    if keys['transcript_type'] == 'transcript_type'
+                    else 'transcript_type'
+                )
+                transcript_types = feature.attributes.get(alternate_key)
+            if transcript_types:
+                for t in transcript_types:
                     to_append += '{},'.format(t)
                 to_append = to_append[:-1] + delim
             else:
                 to_append = to_append + '-' + delim
-            if keys['gene_type'] in feature.attributes.keys():
-                for t in feature.attributes[keys['gene_type']]:
+            gene_types = feature.attributes.get(keys['gene_type'])
+            if not gene_types:
+                alternate_key = (
+                    'gene_biotype'
+                    if keys['gene_type'] == 'gene_type'
+                    else 'gene_type'
+                )
+                gene_types = feature.attributes.get(alternate_key)
+            if gene_types:
+                for t in gene_types:
                     to_append += '{},'.format(t)
                 to_append = to_append[:-1] + delim
             else:

changelog.md

@@ -10,6 +10,16 @@ The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/)
   transcripts without a separate gene feature, including the transcript
   named in issue #3. Tests verify the database and annotations.
 
+### Fixed
+- Accept either `transcript_type`/`gene_type` or
+  `transcript_biotype`/`gene_biotype` attributes when annotating GTF databases,
+  so coding regions are not misclassified when the selected format uses the
+  other pair of keys.
+
+### Tests
+- Verify coding classification under both biotype naming conventions and
+  preserve noncoding classification for a noncoding transcript of a coding gene.
+
 
 ## [0.1.0] - 2022-11-21
 ### Changed

unit_tests/test_biotype_keys.py

@@ -0,0 +1,71 @@
+"""Regression tests for coding and noncoding GTF biotype classification."""
+
+import gffutils
+import pytest
+
+from annotator import annotate_bed
+from annotator.annotator import GENE_PRIORITY, TRANSCRIPT_PRIORITY
+
+
+def annotate_features(tmp_path, gtf_format, gene_key, transcript_key, transcript_type,
+                      include_cds=True):
+    gene_attrs = (
+        f'gene_id "ENSG1"; gene_name "Example"; '
+        f'{gene_key} "protein_coding";'
+    )
+    transcript_attrs = (
+        f'gene_id "ENSG1"; transcript_id "ENST1"; gene_name "Example"; '
+        f'{gene_key} "protein_coding"; {transcript_key} "{transcript_type}";'
+    )
+    features = [("gene", gene_attrs), ("transcript", transcript_attrs),
+                ("exon", transcript_attrs)]
+    if include_cds:
+        features.append(("CDS", transcript_attrs))
+    gtf = tmp_path / "reference.gtf"
+    gtf.write_text("\n".join(
+        f"chr1\ttest\t{feature}\t100\t200\t.\t+\t.\t{attributes}"
+        for feature, attributes in features
+    ) + "\n")
+    db_file = tmp_path / "reference.gtf.db"
+    gffutils.create_db(
+        str(gtf), dbfn=str(db_file), force=True,
+        keep_order=True, merge_strategy="merge", sort_attribute_values=True,
+        disable_infer_genes=False, disable_infer_transcripts=False,
+    )
+    _, _, cds, indexed_features, keys = annotate_bed.create_definitions(
+        str(db_file), gtf_format=gtf_format,
+    )
+    return annotate_bed.annotate(
+        "chr1", 119, 120, "query", "0", "+", True,
+        [priority[:] for priority in TRANSCRIPT_PRIORITY],
+        [priority[:] for priority in GENE_PRIORITY],
+        indexed_features, cds, keys,
+    )
+
+
+@pytest.mark.parametrize(
+    "gtf_format,gene_key,transcript_key",
+    [
+        pytest.param("gencode", "gene_type", "transcript_type",
+                     id="gencode-standard-keys"),
+        pytest.param("gencode", "gene_biotype", "transcript_biotype",
+                     id="gencode-setting-with-biotype-keys"),
+        pytest.param("ensembl", "gene_type", "transcript_type",
+                     id="ensembl-setting-with-type-keys"),
+    ],
+)
+def test_protein_coding_cds_is_selected(
+    tmp_path, gtf_format, gene_key, transcript_key,
+):
+    result = annotate_features(
+        tmp_path, gtf_format, gene_key, transcript_key, "protein_coding",
+    )
+    assert result[6:9] == ("ENSG1", "Example", "CDS")
+
+
+def test_noncoding_transcript_in_coding_gene_stays_noncoding(tmp_path):
+    result = annotate_features(
+        tmp_path, "gencode", "gene_type", "transcript_type",
+        "retained_intron", include_cds=False,
+    )
+    assert result[8] == "noncoding_exon"