genome browser | ggplot | 基因组可视化 | R | transcript | isoform

2021年08月24日 更新

给每个exon编排序号,方便管理

cat gencode.v37.annotation.gtf | grep ENSG00000011304 | cut -f1-5 | grep exon | sort -k 4 | uniq | less -S

 


更新:又用了一段时间,发现是真的好用,可惜就是服务器上装不上,只能在Mac上用了。

 

核心代码很少:

library(ggbio)

# hg19
library(EnsDb.Hsapiens.v75)
ensdb <- EnsDb.Hsapiens.v75

options(repr.plot.width=15, repr.plot.height=9)
autoplot(ensdb, GeneNameFilter("PKM")) +
    # ME of exon 9 and 10
    geom_vline(xintercept = c(72492996,72494795), linetype="dashed", color = "red", size=0.5) +
    theme_bw()

# hg38
library(EnsDb.Hsapiens.v86)
ensdb <- EnsDb.Hsapiens.v86

options(repr.plot.width=15, repr.plot.height=9)
autoplot(ensdb, GeneNameFilter("PKM")) +
    # ME of exon 9 and 10
    geom_vline(xintercept = c(72200655,72202454), linetype="dashed", color = "red", size=0.5) +
    theme_bw()

  

版本问题

ensembl版本列表:https://asia.ensembl.org/info/website/archives/assembly.html

NCBI上的基因都有hg19和hg38两个版本的坐标,可以辅助使用。

 


只搞基因表达还没啥,现在搞到了isoform,即exon级别,想要随时查看底层的数据就很困难了。

比如下面就是AS分析的结果,想要手动去检查数据,完全不知所措。

一个AS event的命名:

isoform1=junction:chr11:237145-244160:+|isoform2=junction:chr11:237145-238997:+@exon:chr11:238998-239076:+@junction:chr11:239077-244160:+

  

里面都是具体的基因组位置了,UCSC Genome Browser和IGV都无法胜任,因为无法标明单个碱基的位置。

 

这时候就发现了一个神器,ggbio: An R implementation for extending the Grammar of Graphics for Genomic Data

tutorial

http://bioconductor.org/packages/release/bioc/html/ggbio.html

https://github.com/tengfei/ggbio

 

另一个发现:spliceclust

 

下面描述其具体用法:

library(ggbio)
library(Homo.sapiens)
class(Homo.sapiens)

data(genesymbol, package = "biovizBase")
wh <- genesymbol[c("PTBP1")]
wh <- range(wh, ignore.strand = TRUE)
p.txdb <- autoplot(Homo.sapiens, which = wh)
p.txdb

options(repr.plot.width=8, repr.plot.height=4)
p.txdb

autoplot(Homo.sapiens, which = wh, label.color = "black", color = "brown", fill = "brown")

autoplot(Homo.sapiens, which = wh, gap.geom = "chevron")

autoplot(Homo.sapiens, which = wh, stat = "reduce")
--------------------------------------------------
library(EnsDb.Hsapiens.v75)
ensdb <- EnsDb.Hsapiens.v75
autoplot(ensdb, GeneNameFilter("PTBP1"))

options(repr.plot.width=10, repr.plot.height=6)
autoplot(ensdb, GeneNameFilter("PTBP1"))

autoplot(ensdb, ~ symbol == "PTBP1", names.expr="gene_name")

autoplot(ensdb, ~ symbol == "PSMD13", names.expr="gene_name")
--------------------------------------------------
## We specify "*" as strand, thus we query for genes encoded on both strands
gr <- GRanges(seqnames = 11, IRanges(238998, 239076), strand = "+")
autoplot(ensdb, GRangesFilter(gr), names.expr = "gene_name")

options(repr.plot.width=15, repr.plot.height=9)
autoplot(ensdb, GRangesFilter(gr), names.expr = "gene_name") +
    #geom_vline(xintercept = c(237145,244160), linetype="dashed", color = "blue", size=0.5) +
    #geom_vline(xintercept = c(237145,238997), linetype="dashed", color = "red", size=0.5) +
    geom_vline(xintercept = c(238998,239076), linetype="dashed", color = "green", size=0.5) +
    #geom_vline(xintercept = c(239077,244160), linetype="dashed", color = "black", size=0.5) +
    theme_bw()

  

最后的图片:

 

稍微搞清楚了的概念:

junction:从一个exon的结束到另一个exon的起点,这就是一个junction,俗称跳跃点、接合点,因为这两个点注定要连在一起。

AS event:一个可变剪切的事件,如何定义呢?一个exon,以及三个junction即可定义,要满足等式相等原则【一个完整的junction = 一个sub junction + exon + 一个sub junction】,即一个exon的inclusion/exclusion。

 

可视化的威力真是无穷,直观明了的揭示了数据之间的关系,这是肉眼直接看数据所做不到的!!!

posted @ 2021-03-05 17:57  Life·Intelligence  阅读(905)  评论(0编辑  收藏  举报
TOP