Use ORFdata with varcode, isovar, topiary and vaxrank

Those tools read a PyEnsembl genome: a GTF of transcript models, plus transcript and protein FASTAs whose identifiers match the GTF. ORFdata writes exactly that for the two things a reference annotation leaves out — noncanonical ORFs and viral references — so a variant in either gets a coding consequence from the same tools that annotate reference transcripts.

Three layers, three jobs

ORFdata keeps the molecular layers apart and names them the same way in every command and API:

Layer	Molecule	What it answers	Who reads it
genome	DNA	Where a read or a variant call sits	Aligners; the contigs a VCF names
transcriptome	RNA	What a variant does to a reading frame	varcode, isovar
proteome	protein	Which peptide a patient could present	topiary, vaxrank, MS search

A DNA reference such as GDC's GRCh38.d1.vd1 is the genome layer: human chromosomes plus the viral sequences as extra contigs. It carries no transcript models and no proteins, which is what the other two layers supply.

Viral transcripts in a tumour sample

Install the collection, or one cancer subset of it, and open it as a genome:

orfdata download viral --species human --subset cancer-causing
orfdata viral summary --subset cancer-causing

from varcode import Variant

from orfdata import select_viral_references, viral_genome

references = select_viral_references(
    subsets=["cancer-causing"], species="human", download_if_missing=True
)
genome = viral_genome(references)  # A PyEnsembl genome over the viral contigs.

# EBV (AJ507799.2), codon 11 of BNRF1
variant = Variant(contig="AJ507799.2", start=1766, ref="C", alt="T", genome=genome)
print(variant.effects().top_priority_effect().short_description)

p.P11S

viral_genome() indexes the set's annotation.gtf, transcripts.fasta and proteins.fasta into a PyEnsembl database kept beside the cache, never inside the checksummed reference directory. The same object goes to isovar with the tumour RNA alignment, and to topiary or vaxrank with an MHC predictor, wherever those take a genome.

Coding references by exact sequence

Isovar asks ORFdata for coding references by exact genome sequence, not through a PyEnsembl genome. viral_references(), viral_orfs() and viral_panel_coverage() read every installed collection as well as the bundled HPV16 catalogue, so a collection installed with orfdata data fetch supplies Isovar's viral models for every contig it holds. A CDS is offered when it translates exactly to its source protein, wherever it lies: on either strand, across the origin of a circular genome, with a frameshift that skips bases (recorded as a subcodon_intron), or partial at either end (the ORF says which boundaries the source establishes). viral_references() says in each reference's scope how many were left out, and they remain in the transcriptome layer. In the GDC collection that is eight of 1,758: three with an ambiguous base, two whose source protein does not match their codons, two HTLV-1 ribosomal-slippage CDSs that read a base twice (#65), and an EBV gene spliced back across the terminal repeat.

An installed reference carries what its own GenBank record says: its group is the genus in the record's lineage, and host_taxon_ids and host_source name the host the sequence was isolated from, with the accession as the source. That is an isolation host, not a host range, and it is empty when the record names none, so viral_references(host="human") returns the references whose record says human and not every virus a human can carry. ictv_species stays empty until curated (#104).

Declaring the genome for vaxrank

A PyEnsembl genome over viral contigs has no species, so vaxrank cannot check it against the analysis host. viral_reference_context(references) states it, for the species the collection targets, and validates against viral_genome():

from orfdata import viral_genome, viral_reference_context

genome = viral_genome(references)
context = viral_reference_context(references)  # Pass as vaxrank's reference_context.

Four things are worth knowing before the first run:

• Contig names are the source accession.version, such as AY446894.2, because that is the record the sequences and the annotation come from. genomes.tsv keeps the inventory's source_genbank and source_refseq spellings beside it, so you can map the names your alignment used.
• A non-standard genetic code is carried, as transl_table on the transcript and CDS rows, so PyEnsembl translates with the code the source record declares.
• A partial or frameshifted source CDS gets no start or stop codon rows, so effect predictors will not translate it. annotation_complete says which transcripts are whole and partial says which end the source marks as running past the record: 69 of 1,758, mostly HPV E4 and E5 ORFs that are really expressed from spliced transcripts, plus HIV-1's pol, which is a Gag-Pol frameshift fusion. HPV16's E6, E7, E1, E2, L1 and L2 are complete, and HBV's polymerase is joined back together across the circular origin, so variants on either side of it get a coding consequence.
• A coding transcript here is the CDS span, not a UTR-containing mRNA. The source records annotate 1,678 CDS features but only 39 mRNA features, so the transcriptome layer is built from the CDS features and says so rather than inventing transcript bounds.
• Three source records annotate no CDS at all — HPV2, HPV8 and HPV51 — so their coding features come from PaVE instead, on the identical sequence, and every transcript records which source annotated it. Use --no-annotation-fill to keep strictly to the GenBank submissions; orfdata viral summary then names the genomes that are left without transcripts. For the opposite trade — one source's uniform, spliced gene models everywhere it publishes the same molecule, in place of each submission's own naming — use --alternate-annotations, which leaves the genome FASTA byte-identical and adds 235 coding transcripts across the GDC collection. orfdata viral completeness shows what either choice leaves on the table.
• A transcript's initiation says how the source starts it (atg, alternative, declared, non_canonical, partial, unknown) and recoding names any frameshift the source claims. Effect predictors handle a non-canonical start correctly, because the protein a source publishes is the one in proteins.fasta; they do not handle a variant at the base a ribosomal-frameshift join reads twice, which is two positions in HTLV-1.

Noncanonical ORFs

Build a proteome as usual and ask for the matching annotation:

orfdata proteome --species human --release 115 --nc-orf-sources ms \
  --output proteins.fasta --annotation ncorfs/

That writes ncorfs/annotation.gtf, ncorfs/transcripts.fasta and ncorfs/proteins.fasta: one coding transcript per ncORF, placed on its host gene, with the spliced coding sequence and the translation. In Python the same is proteome.write_annotation("ncorfs").

Catalogue imports preserve short translated ORFs for variant annotation. If a particular peptide-search workflow requires a length floor, add --min-protein-length 5; the generated proteome records the explicit policy and its omitted count, while the underlying catalogue remains complete.

from pyensembl import Genome

genome = Genome(
    reference_name="GRCh38",
    annotation_name="orfdata-ncorfs",
    annotation_version="1",
    gtf_path_or_url="ncorfs/annotation.gtf",
    transcript_fasta_paths_or_urls=["ncorfs/transcripts.fasta"],
    protein_fasta_paths_or_urls=["ncorfs/proteins.fasta"],
    cache_directory_path="ncorfs/index",
)
genome.index()

A variant inside an ncORF now has a protein consequence — a stop gained in a lncRNA ORF reads as p.K2* rather than disappearing into "non-coding transcript".

To keep the reference transcripts as well, concatenate the files instead of replacing them. PyEnsembl takes several FASTAs, and the ncORF GTF writes no gene rows precisely so it can be appended to an annotation that has them:

cat Homo_sapiens.GRCh38.115.gtf ncorfs/annotation.gtf > combined.gtf

genome = Genome(
    reference_name="GRCh38", annotation_name="ensembl-plus-ncorfs", annotation_version="115",
    gtf_path_or_url="combined.gtf",
    transcript_fasta_paths_or_urls=["Homo_sapiens.GRCh38.cdna.all.fa", "ncorfs/transcripts.fasta"],
    protein_fasta_paths_or_urls=["Homo_sapiens.GRCh38.pep.all.fa", "ncorfs/proteins.fasta"],
    cache_directory_path="combined/index",
)

For the models alone, without building a proteome, the query command writes the same GTF rows for any selection:

orfdata --evidence ms --format gtf > ncorfs.gtf
orfdata --gene LARP4B --format gtf

Two deliberate choices in that file:

• ncORF transcripts are marked transcript_biotype "protein_coding", because each has a complete CDS with start and stop codons and effect predictors only translate transcripts marked that way. The source's own category stays in orf_biotype, and the host transcript in host_transcript_id.
• The transcript is the ORF's coding path, so its identifier is the catalogue's orf_id. The protein is <orf_id>.p. Nothing is renamed to look like an Ensembl accession.

What this does not claim

An ORF model here is a catalogue ORF placed on an assembly, with the evidence the source reported; a predicted coding consequence is not evidence that a stable protein exists or that a peptide is presented. Viral collection and context membership record where a reference was selected, not that a virus drives a particular tumour. The layers stay separate so each of those questions can be asked of the right one.