使用 ggbio 绘制 GenomicRanges 中最长的转录本

Question

我正在尝试使用 ggbio 绘制特定区域。我正在使用下面的代码来产生我想要的输出，除了它包含几个抄本。是否可以只绘制最长的成绩单？我无法访问 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

这是一个解决方案，涉及通过 gene_id 对最长的转录本进行过滤 TxDb.Hsapiens.UCSC.hg19.knownGene（确实会删除没有 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...

^{由 reprex package (v0.3.0)}

于 2020-05-29 创建

编辑： 要获取基因符号而不是基因（或转录本）id，只需将 grl 的名称替换为相关的基因符号，例如通过 org.Hs.eg.db，或任何其他匹配它们的资源。

使用 ggbio 绘制 GenomicRanges 中最长的转录本

Plot the longest transcript in GenomicRanges with ggbio

r

genomicranges