From 7dfadaf070a40eb0acb06c42cf84a55a58e5ee7a Mon Sep 17 00:00:00 2001 From: Brian Yee Date: Mon, 5 Oct 2026 05:37:31 -0700 Subject: [PATCH] Handle alternate GTF biotype keys --- README.md | 18 +++++---- annotator/annotate_bed.py | 24 +++++++++-- changelog.md | 10 +++++ unit_tests/test_biotype_keys.py | 71 +++++++++++++++++++++++++++++++++ 4 files changed, 112 insertions(+), 11 deletions(-) create mode 100644 unit_tests/test_biotype_keys.py diff --git a/README.md b/README.md index 07e8a4d..81bc26f 100644 --- a/README.md +++ b/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()`. diff --git a/annotator/annotate_bed.py b/annotator/annotate_bed.py index 559a17e..dc0a15b 100755 --- a/annotator/annotate_bed.py +++ b/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: diff --git a/changelog.md b/changelog.md index 868e0b2..dd86322 100644 --- a/changelog.md +++ b/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 diff --git a/unit_tests/test_biotype_keys.py b/unit_tests/test_biotype_keys.py new file mode 100644 index 0000000..8d1f0b3 --- /dev/null +++ b/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"