hpages
commented
2 years ago Hi Janet,
The error message is a little bit misleading: it should say that makeTxDbFromGFF() does not support coding-transcripts that encode more than one protein.
In the Canonical Gene example shown in the GFF3 specs, transcript EDEN.3 encodes 2 proteins (edenprotein.3 and edenprotein.4), via 2 CDSs (cds00003 and cds00004). makeTxDbFromGFF(), and the underlying db schema used by TxDb objects, cannot handle this. Note that AFAIK the 2 major annotation providers UCSC and Ensembl don't handle this either. For example transcript tables at UCSC like knownGene or ncbiRefSeq use the cdsStart and cdsEnd fields to store a single CDS start and end value per transcript. IIRC the MySQL db used at Ensembl does something similar.
When we started working on TxDb objects back in 2010, our primary goal was to support UCSC and Ensembl annotations, so we designed a db schema that assumes a one-to-one relationship between coding transcripts and CDSs/proteins. So when you import the Canonical Gene example, only the first CDS on transcript EDEN.3 is imported:
library(txdbmaker)
EDEN <- system.file("extdata", "GFF3_files", "TheCanonicalGene_v1.gff3", package="txdbmaker")
txdb <- makeTxDbFromGFF(EDEN)
# Import genomic features from the file as a GRanges object ... OK
# Prepare the 'metadata' data frame ... OK
# Make the TxDb object ... OK
# Warning message:
# In .find_exon_cds(exons, cds) :
# The following transcripts have exons that contain more than one CDS
# (only the first CDS was kept for each exon): mRNA00003
cdsBy(txdb, by="tx", use.names=TRUE)
# GRangesList object of length 3:
# $EDEN.1
# GRanges object with 4 ranges and 3 metadata columns:
# seqnames ranges strand | cds_id cds_name exon_rank
# <Rle> <IRanges> <Rle> | <integer> <character> <integer>
# [1] ctg123 1201-1500 + | 1 edenprotein.1 1
# [2] ctg123 3000-3902 + | 3 edenprotein.1 2
# [3] ctg123 5000-5500 + | 5 edenprotein.1 3
# [4] ctg123 7000-7600 + | 8 edenprotein.1 4
# -------
# seqinfo: 1 sequence from an unspecified genome; no seqlengths
#
# $EDEN.2
# GRanges object with 3 ranges and 3 metadata columns:
# seqnames ranges strand | cds_id cds_name exon_rank
# <Rle> <IRanges> <Rle> | <integer> <character> <integer>
# [1] ctg123 1201-1500 + | 2 edenprotein.2 1
# [2] ctg123 5000-5500 + | 6 edenprotein.2 2
# [3] ctg123 7000-7600 + | 9 edenprotein.2 3
# -------
# seqinfo: 1 sequence from an unspecified genome; no seqlengths
#
# $EDEN.3
# GRanges object with 3 ranges and 3 metadata columns:
# seqnames ranges strand | cds_id cds_name exon_rank
# <Rle> <IRanges> <Rle> | <integer> <character> <integer>
# [1] ctg123 3301-3902 + | 4 edenprotein.3 2
# [2] ctg123 5000-5500 + | 7 edenprotein.3 3
# [3] ctg123 7000-7600 + | 10 edenprotein.3 4
# -------
# seqinfo: 1 sequence from an unspecified genome; no seqlengthsThe one-to-one relationship between coding transcripts and CDSs/proteins is a very early assumption that is at the root of the design/behavior of many functions in GenomicFeatures/txdbmaker, so is something that would be very difficult to change, and also very disruptive.
I knew that the GFF3 specs did support a one-to-many relationship between coding transcripts and CDSs/proteins, and that's a feature that I found surprising when I ran into it, because I was not aware of any real-world situation where this would actually happen, and because AFAIK none of the hundreds of thousands of GFF3 files provided by UCSC and Ensembl seemed to use this capability. Also makeTxDbFromGFF() has been used on hundreds or thousands of GFF3 files from FlyBase, WormBase, and other annotations providers, for many years, and AFAIK nobody has reported seeing that. So congratulations for finding and reporting the first case of such file :wink:! Out of curiosity, may I know where is this file coming from, if that's something you're willing to share?
I don't know how UCSC or Ensembl handle this situation, or how frequently it is observed in biology. Maybe they just use 2 distinct transcript ids to represent such event? In your case it means that the GFF3 file would need to be modified but I don't know how hard that would be. If this is a viral genome, then there's no exon, and each CDS can be represented by a single line (no CDS parts), so the file should be much simpler than in the general situation where the gene/transcript/exon/cds hierarchy can be complex and hard to figure out. Maybe makeTxDbFromGFF() should be able to do something this like e.g. by importing 2 transcripts instead of one and assigning the corresponding protein IDs to them, or by adding the .1 and .2 suffixes to the original transcript ID? I'm open to suggestions.
H.
jayoung
hi there,
I am working on a viral genome, where the annotations (based on ribosome profiling) fairly often include >1 CDS per exon. Using
makeTxDbFromGFFmeans that some of these CDSs get dropped with a warningThe following transcripts have exons that contain more than one CDS (only the first CDS was kept for each exon). I'd like to be able to keep all the CDSs.I'm not sure if this website is what determines 'official' gff3 specs, but it suggests that >1 CDS per exon should be allowed: https://github.com/The-Sequence-Ontology/Specifications/blob/master/gff3.md
Any change that makeTxDbFromGFF could support >1 CDS per exon? Not sure how tricky that would be. Example data below, from the website linked above.
Thanks!
Janet
Here's the example gff file from that page:
And here's the warning - we lose the edenprotein.4 annotation