ID version handling and GENCODE compatibility

[ccwang002]

Hi team, thanks again for maintaining this library.

We (together with @fkgruber) are trying to use pyensembl with a custom genome build based on GENCODE v48, but we’re running into issues where variant effects are not inferred correctly. For example,

from pyensembl import Genome
from varcode import Variant

g = Genome(
    reference_name="GRCh38",
    annotation_name="my_genome_features",
    gtf_path_or_url="gencode.v48.primary_assembly.annotation.gtf",
    protein_fasta_paths_or_urls="gencode.v48.pc_translations.fa.gz",
    transcript_fasta_paths_or_urls="gencode.v48.transcripts.fa.gz",
    cache_directory_path=".cache",
)

Variant("chr1", 12110124, "G", "A", ensembl=g).effects()
# <EffectCollection with 3 elements>
#  -- NoncodingTranscript(...)
#  -- NoncodingTranscript(...)
#  -- NoncodingTranscript(...)

# Expected output say using Ensembl release 114
# <EffectCollection with 3 elements>
#  -- Substitution(variant=..., transcript_name=TNFRSF8-201, transcript_id=ENST00000263932, effect_description=p.R199H)
#  -- Substitution(variant=..., transcript_name=TNFRSF8-203, transcript_id=ENST00000417814, effect_description=p.R88H)
#  -- NoncodingTranscript(variant=..., transcript_name=TNFRSF8-205, transcript_id=ENST00000514649)

So far we've identified some potential issues:

1. Protein/transcript ID version handling.
   It looks like pyensembl trims version suffixes from FASTA IDs (e.g. .1 in ENSPxxx.1), which seems to break mapping between GTF and FASTA entries.

   • GENCODE GTF includes versioned IDs directly (e.g. ENSTxxx.7, ENSPxxx.1)
   • Ensembl GTF separates version into a different field (transcript_version, protein_version)

   The version trimming seems to occur at these places:

   • https://github.com/openvax/pyensembl/blob/v2.6.9/pyensembl/fasta.py#L50-L61
   • https://github.com/openvax/pyensembl/blob/v2.6.9/pyensembl/transcript.py#L454-L463

2. Biotype field naming differences.
   GENCODE uses: gene_type transcript_type, while pyensembl expects: gene_biotype, transcript_biotype

Modifications to GENCODE GTF to make it work for pyensembl

zcat gencode.v48.primary_assembly.annotation.gtf.gz \
| sed -E '
  s/gene_type "/gene_biotype "/g;
  s/transcript_type "/transcript_biotype "/g;
  s/protein_id "ENSP([0-9]+)\.[0-9]+"/protein_id "ENSP\1"/g
' \
| gzip > gencode.v48.primary_assembly.annotation.biotype_protein_noversion.gtf.gz

Questions

1. Is there a recommended way to preserve transcript/protein versioning within pyensembl?
2. GENCODE compatibility (see also Compatability with Gencode GTF #168)
   • Are there official or recommended preprocessing steps for using GENCODE GTFs with pyensembl?
   • Beyond the fixes above, are there other known incompatibilities we should account for?

Happy to help test or contribute if useful. Thanks again for your work on this project.

[iskandr]

Part 2 shipped in v2.9.6 (PR #350).

The protein/transcript-ID half of this issue is now fixed:

• pyensembl's FASTA parser strips ENS `.N` suffixes from headers, so the dict key is unversioned.
• GENCODE GTFs embed the version directly in `protein_id` / `transcript_id`, so the previous literal lookup missed.
• A new internal helper sequence_lookup_with_ens_fallback is wired into Transcript.sequence, Transcript.protein_sequence, Genome.transcript_sequence(id), and Genome.protein_sequence(id). It tries the literal ID first (handles both styles) and falls back to the version-stripped form for ENS-prefix IDs on miss.

After updating you should no longer need the `sed` workaround in your repro for protein_id "ENSPxxx.1". `pip install -U pyensembl` to pick it up and please re-run the test variant — would be great to confirm the effects come back as `Substitution(...)` instead of `NoncodingTranscript(...)`.

The biotype-naming half (`gene_type` / `transcript_type` vs `gene_biotype` / `transcript_biotype`) is still open as part 1 — small follow-up PR.

[ccwang002]

Wow, this is amazing! Thanks @iskandr for your speedy work. Could you point us to the part 1 PR regarding the GTF field naming convention?

[iskandr]

Part 1 is stuck with some other refactoring going on but I'll let you know when it's done

[iskandr]

Both halves of this issue are now shipped:

Half	Version	PR
Versioned protein-ID FASTA matching	v2.9.6 (refined in pending #355)	#350
GENCODE biotype column aliases (gene_type → gene_biotype, transcript_type → transcript_biotype)	v2.9.8	#356

For the original repro (GENCODE GTF + protein FASTA), the pyensembl-side observables to confirm after pip install -U pyensembl:

from pyensembl import Genome

g = Genome(
    reference_name="GRCh38",
    annotation_name="my_genome_features",
    gtf_path_or_url="gencode.v48.primary_assembly.annotation.gtf",
    protein_fasta_paths_or_urls="gencode.v48.pc_translations.fa.gz",
    transcript_fasta_paths_or_urls="gencode.v48.transcripts.fa.gz",
    cache_directory_path=".cache",
)
g.index()

# Biotype now lands on the Ensembl canonical column (was empty before #356)
t = g.transcript_by_id("ENST00000263932.10")  # versioned IDs are preserved
assert t.biotype == "protein_coding"
assert t.is_protein_coding is True

# Protein sequence resolves against a GENCODE-style versioned protein_id
# (returned None before #350)
assert t.protein_sequence is not None
assert len(t.protein_sequence) > 0

# Coding sequence + translation pipeline works. The start codon is usually
# ATG but Ensembl also annotates alternative initiators (CTG, GTG, TTG,
# ACG, etc.) for transcripts where that is the biology, so don't assert
# on the literal — just that the call returns a non-empty string whose
# length is divisible by 3.
assert t.coding_sequence is not None
assert len(t.coding_sequence) % 3 == 0

You no longer need either of your sed workarounds (the gene_type → gene_biotype rewrite or the protein_id "ENSPxxx.\d+" version strip). The GTF and FASTA can be loaded as-is.

Downstream tools that consume these pyensembl accessors (varcode, isovar, vaxrank, etc.) will see the corrected data on the same upgrade. Closing — happy to reopen if you hit any residual rough edges on a re-run.