This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Get the longest transcript in GenomicRanges

I am trying to plot an specific region using ggbio. I am using the below code that produces my desire output, except that it contains several transcript. Is it possible to only plot the longest transcript for each gene? I've not been able to access the genomic ranges object within Homo.sapiens that I assume contains this information.

library(ggbio)
library(Homo.sapiens)
data(genesymbol, package = "biovizBase")

range <- GRanges("chr10"  , IRanges(start = 78000000 , end = 79000000))
p.txdb <- autoplot(Homo.sapiens, which = range)
p.txdb
ggbio r genomicranges

Thanks for the help. I've checked which data can be retrieved from the Homo.sapiens object using columns(Homo.sapiens), and there is actually information on the transcript start and transcript end. I think I will need to get the width first, and then modify the Homo.sapiens object somehow?

I checked that, and a TxDb object does not contain gene names... I would like to show the name of the genes not the transcripts names.

That definitely works but there must be an easier (in-built way, I hope). Anyway, getting the longest transcript is simple:

library(ggbio)
library(Homo.sapiens)

data(genesymbol, package = "biovizBase")

GetLongest <- function(Ranges){

  ## list transcripts
  tx <- transcripts(Homo.sapiens)

  ## list genes
  gns <- genes(Homo.sapiens)

  ## get transcripts overlapping the ranges
  subs <- subsetByOverlaps(tx, Ranges)

  ## overlap with genes to get gene names
  olap <- findOverlaps(subs, gns)

  ## some transcripts might not have a matched gene in the database
  subs <- subs[olap@from]

  ## add gene ID to the metadata:
  elementMetadata(subs) <- data.frame(elementMetadata(subs), 
                                      GENEID = as.character(gns[olap@to]$GENEID))

  ## sort by length and take the unique genes.
  ## that gives you the longest transcript per annotated gene
  subs <- subs[order(width(subs), decreasing = TRUE)]
  longest <- subs[!duplicated(subs$GENEID)]
  return(longest)
}

GetLongest(GRanges("chr10"  , IRanges(start = 78000000 , end = 79000000)))

GRanges object with 2 ranges and 3 metadata columns:
      seqnames            ranges strand |   TXID     TXNAME      GENEID
         <Rle>         <IRanges>  <Rle> | <list>     <list> <character>
  [1]    chr10 77542519-78317126      + |  38363 uc001jxi.3       83938
  [2]    chr10 78629359-79397577      - |  40147 uc001jxj.2        3778
  -------
  seqinfo: 93 sequences (1 circular) from hg19 genome

Great, thanks! So, then, within the autoplot function, how do I select only those transcripts? I am trying to do that, but it keeps printing all transcripts.

I just found out that I can reduce the transcripts... which is not exactly what I want, as I think this is collapsing all the transcripts... i.e. p.txdb <- autoplot(Homo.sapiens, which = range, stat = "reduce")

I am trying to do that, but it keeps printing all transcripts.

Yeah, that is unfortunately what I did not find out yet.

0 answers

No answers yet.

Log in to answer this question.