Hi,

I am having a problem with AnnotationHub when trying to convert transcript coordinates to genome coordinates.

Some information first on the context of where the data comes from: I started with a GTF file assembled from RNA-Seq data using the Ensembl reference transcriptome (v. 99) as reference. I then translated these transcripts in 3 frames to generate a theoretical translatome. This protein database was then used to search mass spec data to identify proteins in my samples. I next mapped the identified proteins back to their transcript IDs and sequences and was able to get transcript coordinates for each protein sequence. Now, I would like to convert the transcript positions to genomic positions.

I used the Annotation Hub package with the following code:

    ah = AnnotationHub()
    ahDb <- query(ah, pattern = c("Homo sapiens", "EnsDb", 99))
    ahEdb <- ahDb[[1]]

    df_protein_coding['genome.start'] <- NA
    df_protein_coding['genome.end'] <- NA
    row.names(df_protein_coding) <- NULL

    for (i in 1:nrow(df_protein_coding)){

      start = df_protein_coding['tx_start'][[1]][i]
      end = df_protein_coding['tx_end'][[1]][i]
      tx_id=df_protein_coding['transcript.id'][[1]][i]
      chr = df_protein_coding['Chromosome.scaffold.name'][[1]][i]

      rng_tx <- IRanges(start = start, end = end, names = tx_id)
      edbx <- filter(ahEdb, filter = ~ seq_name == chr)

      rng_gnm <- transcriptToGenome(rng_tx, edbx)
      a <- GRangesList(rng_gnm)
      df_protein_coding['genome.start'][[1]][i] = start(a)[[1]]
      df_protein_coding['genome.end'][[1]][i] = end(a)[[1]]
    }

However, I run into this error:

Error in df_protein_coding["genome.start"][[1]][i] <- start(a)[[1]] : 
  replacement has length zero
In addition: There were 24 warnings (use warnings() to see them)

The code seems to have worked for the first few entires. Here is the truncated version for the entries that worked:

The table has been transposed for easier formatting.

The entry that AnnotationHub has an issue with is a sequence from position 8698 to 8760 of ENST00000515408. Not sure where I have gone wrong here as it seems that the transcript coordinates are out of bounds from the one that AnnotationHub has. Does anyone have any suggestions on how to solve this? Would appreciate any help with this issue.

The length of transcript ENST00000515408 is only 8029 nucleotides,

transcriptLengths(filter(ahEdb, filter = ~ tx_id == "ENST00000515408"))
                          tx_id         tx_name         gene_id nexon tx_len
ENST00000515408 ENST00000515408 ENST00000515408 ENSG00000172795    10   8029

but you try to map nucleotides 8698 to 8760 to the genome. And this results in an empty GRanges:

rng_tx <- IRanges(start = 8698, end = 8760, names = "ENST00000515408")
transcriptToGenome(rng_tx, ahEdb)
GRangesList object of length 1:
$ENST00000515408
GRanges object with 0 ranges and 0 metadata columns:
   seqnames    ranges strand
      <Rle> <IRanges>  <Rle>
  -------
  seqinfo: no sequences

Warning messages:
1: The within transcript/CDS coordinates are outside the region defined by the provided exons 
2: In transcriptToGenome(rng_tx, ahEdb) :
  1 transcript(s) could either not be found in the database or the specified range is outside the transcript's sequence

I would suggest you to re-check how you calculate the within-transcript coordinates. Also, always ensure that you're using the same Ensembl version for all your steps - could be that the 5' and 3' UTRs for some transcripts change over Ensembl releases (although I guess that's happening much less frequent than some years ago).

best, jo