Genome to protein does not work

[tprodanov]

I have R 4.5.1, emsembldb v2.32.0 and EnsDb.Hsapiens.v86 v2.99.0.
Currently, the code from this tutorial does not work for me:

library(ensembldb)
library(EnsDb.Hsapiens.v86)
edbx <- filter(EnsDb.Hsapiens.v86, filter = ~ seq_name == "X")
gnm <- GRanges("X:107716399-107716401")
genomeToProtein(gnm, edbx)

produces

IRanges object with 2 ranges and 8 metadata columns:
          start       end     width |           tx_id    cds_ok         exon_id exon_rank
      <integer> <integer> <integer> |     <character> <logical>     <character> <integer>
  [1]        -1        -1         1 | ENST00000372390      <NA> ENSE00001457675         1
  [2]        -1        -1         1 | ENST00000486554      <NA> ENSE00001927337         1
      seq_start   seq_end    seq_name  seq_strand
      <integer> <integer> <character> <character>
  [1] 107716399 107716401           X           *
  [2] 107716399 107716401           X           *
Warning message:
Provided coordinates for 'ENST00000372390', 'ENST00000486554' are not within the coding region

It also does not work for genes I am interested in, so I guess the problem is in general for all genes.

[jorainer]

According to the current genome build/transcript definitions from Ensembl (release 114) your selected region is within introns of the listed transcripts above (see here). That's why ensembldb can not map it to protein coordinates. ensembldb basically first checks whether the genomic region is within exons of transcripts, then it evaluates whether it is within the coding region of these transcripts and only then maps the coordinates to the associated proteins (i.e., their sequences).

Also, I would suggest to use annotations from a more recent Ensembl release (unless you work in fact with Ensembl release 86 annotations). For example (to use Ensembl release 113):

library(AnnotationHub)
ah <- AnnotationHub()
edb <- ah [["AH119325"]]

Then, selecting a genomic region that is within coding exons of transcripts, e.g., region X:102716420-102716422 you would get:

gnm <- GRanges("X:102716420-102716422")
genomeToProtein(gnm, edb)
IRangesList object of length 1:
[[1]]
IRanges object with 5 ranges and 8 metadata columns:
                      start       end     width |           tx_id    cds_ok
                  <integer> <integer> <integer> |     <character> <logical>
  ENSP00000339057       517       518         2 | ENST00000332262      TRUE
  ENSP00000507692       517       518         2 | ENST00000483720      TRUE
  ENSP00000508139       517       518         2 | ENST00000486814      TRUE
  ENSP00000437394       517       518         2 | ENST00000535209      TRUE
  ENSP00000437872       517       518         2 | ENST00000543253      TRUE
                          exon_id exon_rank seq_start   seq_end    seq_name
                      <character> <integer> <integer> <integer> <character>
  ENSP00000339057 ENSE00001923525         4 102716420 102716422           X
  ENSP00000507692 ENSE00001923525         5 102716420 102716422           X
  ENSP00000508139 ENSE00001890623         5 102716420 102716422           X
  ENSP00000437394 ENSE00002284792         4 102716420 102716422           X
  ENSP00000437872 ENSE00002284792         5 102716420 102716422           X
                   seq_strand
                  <character>
  ENSP00000339057           *
  ENSP00000507692           *
  ENSP00000508139           *
  ENSP00000437394           *
  ENSP00000437872           *

Note also that pre-filtering by chromosome (as you have done in your example) would not really be required - it could have a minimal performance advantage, but neglectible (I would say).