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。

 

可視化的威力真是無窮，直觀明了的揭示了數據之間的關系，這是肉眼直接看數據所做不到的！！！