align genomic features with phylogenetic tree
biostars
上有人问,如何画出下面这种图,把alignment和进化树对齐画在一起。
这个一点都不难,我可以写一个geom来画这个alignment,然后上facet_plot来对齐。
但我不想重新造轮子,因为ggbio包有很多geom画genomic data,我希望ggtree和ggbio可以互作,这也是R的可爱之处,软件包之间可以互补,让对方更加完整。
以下是基因BRCA1和NBR1的位置信息:
library(ggbio)library(biovizBase)
library(Homo.sapiens)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
data(genesymbol, package = "biovizBase")
wh <- genesymbol[c("BRCA1", "NBR1")]
wh <- range(wh, ignore.strand = TRUE)
gr.txdb <- crunch(txdb, which = wh)
## change column to model
colnames(values(gr.txdb))[4] <- "model"
grl <- split(gr.txdb, gr.txdb$tx_id)
set.seed(2016-10-25)
names(grl) <- sample(LETTERS, size = length(grl), replace = TRUE)
看不懂上面的代码没关系,这个grl画出来是这样子的:
ggplot() + geom_alignment(grl, alpha=.5)
ggtree
定义的facet_plot
可以很容易地将各种geom
画在右侧分面上,并自动调整位置与进化树对齐。但它要求输入的数据是data.frame
,而这里grl
是一个GRangesList
对象,于是我更新了facet_plot
,让它支持geom_alignment
,在这个过程中,我发现了ggbio
的一个bug,并向Michael提交了一个补丁,这个补丁打在了ggbio 1.23.2
上,所以你的ggbio版本必须>=1.23.2
才行。
做为演示,让我们先展示一颗随机树:
library(ggtree)
n <- names(grl) %>% unique %>% length
set.seed(2016-10-25)
tr = rtree(n)
set.seed(2016-10-25)
tr$tip.label = sample(unique(names(grl)), n)
p <- ggtree(tr) + geom_tiplab()
这树画出来,是这样子滴:
好啦,有了facet_plot
的更新和ggbio
的补丁,要画树和alignment
,那是相当容易:
facet_plot(p, 'alignment', grl, geom_alignment,
inherit.aes=FALSE, mapping=aes())
只需要上面这一句,下面这图就生成了:
别看这文章中代码多,其实多半是生成这个演示数据,核心代码来画这个树就是:
#假如我们有一颗树tr, 和一个基因位置信息的数据grl
p <- ggtree(tr) + geom_tiplab()
facet_plot(p, 'alignment', grl, geom_alignment,
inherit.aes=FALSE, mapping=aes())
只此两句,第一句画树,第二句画alignment并对齐。
PS:mapping=aes()在这里是必须的,因为默认的mapping=NULL不被ggbio支持。
PS2:目前我只测试了geom_alignment,它支持GRanges和GRangesList对象,其它的ggbio定义的geom不一定支持,如果发现不支持,欢迎反馈给我,我可以努力让它们基友友情更牢固🤔