查看原文
其他

align genomic features with phylogenetic tree

2017-01-19 Y叔 biobabble

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不一定支持,如果发现不支持,欢迎反馈给我,我可以努力让它们基友友情更牢固🤔

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存