Постройте самый длинный транскрипт в GenomicRanges с помощью ggbio

Question

Я пытаюсь построить конкретную область c, используя ggbio. Я использую приведенный ниже код, который произвел мое желание, за исключением того, что он содержит несколько стенограмм. Можно ли построить только самую длинную стенограмму? Мне не удалось получить доступ к объекту genomi c range в пределах Homo.sapiens, который, как я полагаю, содержит эту информацию.

library(ggbio)
library(Homo.sapiens)

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

Answer 1

Вот решение, которое включает фильтрацию TxDb.Hsapiens.UCSC.hg19.knownGene в самой длинной транскрипции с помощью gene_id (что удаляет гены без gene_id):

suppressPackageStartupMessages({
    invisible(lapply(c("ggbio", "biovizBase", "data.table", 
        "TxDb.Hsapiens.UCSC.hg19.knownGene", 
        "org.Hs.eg.db"),
        require, character.only = TRUE))})
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

# retrieve transcript lengths
txlen <- transcriptLengths(txdb, with.utr5_len=TRUE, with.utr3_len=TRUE)
setDT(txlen)
txlen$len <- rowSums(as.matrix(txlen[, .(tx_len, utr5_len, utr3_len)]))
setkey(txlen, gene_id, len, tx_id)

# filter longesttranscript by gene_id
ltx <- txlen[!is.na(gene_id)][, tail(.SD,1), by=gene_id]$tx_id

# filter txdb object
txb <- as.list(txdb)
txb$transcripts <- txb$transcripts[txb$transcripts$tx_id %in% ltx, ]
txb$splicings <- txb$splicings[txb$splicings$tx_id %in% ltx,]
txb$genes <- txb$genes[txb$genes$tx_id %in% ltx,]
txb <- do.call(makeTxDb, txb)

# plot according to vignette, chapter 2.2.5
range <- GRanges("chr10", IRanges(start = 78000000 , end = 79000000))
gr.txdb <- crunch(txb, which = range)
#> Parsing transcripts...
#> Parsing exons...
#> Parsing cds...
#> Parsing utrs...
#> ------exons...
#> ------cdss...
#> ------introns...
#> ------utr...
#> aggregating...
#> Done
colnames(values(gr.txdb))[4] <- "model"
grl <- split(gr.txdb, gr.txdb$gene_id)
symbols <- select(org.Hs.eg.db, keys=names(grl), columns="SYMBOL", keytype="ENTREZID")
#> 'select()' returned 1:1 mapping between keys and columns
names(grl) <- symbols[match(symbols$ENTREZID, names(grl), nomatch=0),"SYMBOL"]
autoplot(grl, aes(type = "model"), gap.geom="chevron")
#> Constructing graphics...

^{Создано 29.05.2020 пакетом REPEX (v0.3.0)}

Изменить: Чтобы получить символы гена вместо гена ( или транскрипт) идентификаторов, просто замените имена grl соответствующими символами генов, например через org.Hs.eg.db, или любым другим ресурсом, который им соответствует.