文章/答案/技术大牛

发布

社区首页 >专栏 >R语言之系统进化树的美化

R语言之系统进化树的美化

一粒沙

发布于 2019-07-31 06:40:35

5.5K00

代码可运行

文章被收录于专栏：R语言交流中心R语言交流中心

运行总次数：0

代码可运行

百度百科对进化树的定义是：在生物学中，用来表示物种之间的进化关系。生物分类学家和进化论者根据各类生物间的亲缘关系的远近，把各类生物安置在有分枝的树状的图表上，简明地表示生物的进化历程和亲缘关系。在进化树上每个叶子结点代表一个物种，如果每一条边都被赋予一个适当的权值，那么两个叶子结点之间的最短距离就可以表示相应的两个物种之间的差异程度。同时有很多算法应运而生主要包括：贝叶斯法（Bayesian），最大似然法（Maximum likelihood，ML），最大简约法（Maximum parsimony，MP），邻接法（Neighbor-Joining，NJ），最小进化法（Minimum Evolution，ME），类平均法（UPGMA）。与此同时相对应的软件也出现，下图总结来源于网络：

我们今天就不一一介绍树状图的形成，如果实在没操作过那可以参考下面的创建进化树数据的实例：

mafft --auto ggtree.fasta > ggtree_aligned.fasta##序列的比对
./FastTree ggtree_ aligned.fasta >ggtree_resource.tree ###最大似然法进化树构建

接下来就是我们今天的主题如何对获得进化树数据进行美化。需要用到R语言的包ggtree。

包的安装我们就是不赘述了，请参考bioconductor安装方式。

首先我们看数据的导入。所用到的函数是read.tree。

主要的参数：

file和text参数是相互补充的数据。如果设置了file就不用管text，如果传入自己的文本数据那就会忽略file。

Skip参数主要是在读入数据时需要忽略的前多少行。同时去掉后面多少字符的则是coment.char,它会忽略此处设置的单个字符串以后的数据。

keep.multi 主要是面对多个进化树时，可以设置为true，并且利用tree.names进行选择对应的树。

实例：

单个树的文件构建：

s <-"owls(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3, Tyto_alba:13.5);"
cat(s, file = "ex.tre", sep ="\n")
tree.owls <-read.tree("ex.tre")
str(tree.owls)

s <-"owls(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3, Tyto_alba:13.5);\nowls1(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3,Tyto_alba:13.5);"
cat(s, file = "ex.tre", sep ="\n")
tree.owls <-read.tree("ex.tre", keep.multi = TRUE)
str(tree.owls)

当然有时候我们直接从给我们的其他软件直接生成的tree文件，就可以直接利用file参数读进去，那个就不赘述了。

接下来我们看下其绘制进化树图的主函数ggtree。

主要参数：

Tr主要是进化树的数据

Layout主要是展现的样式具体的我们引用下官网的例子：

Ladderize主要是指的是否对进化树进行阶梯状重构，同时结合right进行设置阶梯化的方向。

Open.angle主要是对于“fan”布局角度设置。

当然其他的一些修饰是完全和ggplot2一致的，比如color，linetype等。

接下来就是对构建好的基础系统进化树进行美化，添加图层，主要会用到下面的函数，其语法类似于ggplot2：

我们介绍其中几个比较常用的：

theme_tree2 添加X轴。

实例：

set.seed(2017-02-16)
tree <- rtree(50)
ggtree(tree) + theme_tree2()

geom_point2 注释想要注释的样本信息。

实例：

ggtree(tree) +theme_tree2()+geom_point2(aes(subset=(node == 21)), size=5, shape=23,fill="steelblue")

geom_point 注释每一个节点。

实例：

ggtree(tree)+theme_tree2()+geom_point(aes(shape=isTip, color=isTip), size=3)

geom_treescale 添加距离标尺。

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip), size=3)+geom_treescale(x=0,y=45, width=0.2, color='red')

geom_nodepoint 注释节点信息，可以在节点生成圆形的阴影。

实例：

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip),size=3)+geom_treescale(x=0, y=45, width=0.2,color='red')+geom_nodepoint(color="#b5e521", alpha=1/4, size=10)

geom_tippoint 对每个样本进行标注。

实例：

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip),size=3)+geom_treescale(x=0, y=45, width=0.2,color='red')+geom_nodepoint(color="#b5e521", alpha=1/4,size=10)+geom_tippoint(color="#FDAC4F", shape=8, size=3)

geom_tiplab 标注每一个样本的名称。如果是圆形的图像需要设置aes(angle=angle)。

实例：

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip),size=3)+geom_treescale(x=0, y=45, width=0.2,color='red')+geom_nodepoint(color="#b5e521", alpha=1/4,size=10)+geom_tippoint(color="#FDAC4F", shape=8,size=3)+geom_tiplab(size=3, color="purple")

geom_text2 注释每个节点的编号。

实例：

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip),size=3)+geom_treescale(x=0, y=45, width=0.2, color='red')+geom_nodepoint(color="#b5e521",alpha=1/4, size=10)+geom_tippoint(color="#FDAC4F", shape=8,size=3)+geom_tiplab(size=3,color="purple")+geom_text2(aes(subset=!isTip, label=node), hjust=-.3)

geom_cladelabel 注释某一个节点区域。

实例：

ggtree(tree) + theme_tree2()+geom_point(aes(shape=isTip,color=isTip), size=3)+geom_treescale(x=0, y=45, width=0.2,color='red')+geom_nodepoint(color="#b5e521", alpha=1/4,size=10)+geom_tippoint(color="#FDAC4F", shape=8,size=3)+geom_tiplab(size=3, color="purple")+geom_text2(aes(subset=!isTip,label=node), hjust=-.3)+geom_cladelabel(node=84, label="test label")+geom_cladelabel(node=59, label="another clade")

geom_strip 可以更加美化的进行色块的直接注释。

实例：

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip),size=3)+geom_treescale(x=0, y=45, width=0.2,color='red')+geom_nodepoint(color="#b5e521", alpha=1/4,size=10)+geom_tippoint(color="#FDAC4F", shape=8, size=3)+geom_tiplab(size=3,color="purple")+geom_text2(aes(subset=!isTip, label=node),hjust=-.3)+geom_cladelabel(node=84, label="test label")+geom_cladelabel(node=59, label="another clade")+geom_strip(57, 59,barsize=2, color='red') + geom_strip(87, 88, barsize=2, color='blue')

geom_hilight 注释某一个节点的区域。

实例：

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip),size=3)+geom_treescale(x=0, y=45, width=0.2,color='red')+geom_nodepoint(color="#b5e521", alpha=1/4,size=10)+geom_tippoint(color="#FDAC4F", shape=8, size=3)+geom_tiplab(size=3,color="purple")+geom_text2(aes(subset=!isTip, label=node),hjust=-.3)+geom_hilight(node=59, fill="darkgreen", alpha=.6)

geom_balance 注释某个节点后相邻子代的两代区域。

实例：

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip),size=3)+geom_treescale(x=0, y=45, width=0.2,color='red')+geom_nodepoint(color="#b5e521", alpha=1/4,size=10)+geom_tippoint(color="#FDAC4F", shape=8, size=3)+geom_tiplab(size=3,color="purple")+geom_text2(aes(subset=!isTip, label=node),hjust=-.3)+geom_balance(node=95, fill='darkgreen', color='white', alpha=0.6,extend=1)

Theme 进行legend注释设置。

实例：

ggtree(tree) +theme_tree2()+geom_point(aes(shape=isTip, color=isTip), size=3)+geom_treescale(x=0,y=45, width=0.2, color='red')+geom_nodepoint(color="#b5e521",alpha=1/4, size=10)+geom_tippoint(color="#FDAC4F", shape=8,size=3)+geom_tiplab(size=3,color="purple")+geom_text2(aes(subset=!isTip, label=node), hjust=-.3)+geom_balance(node=95,fill='darkgreen', color='white', alpha=0.6,extend=1)+theme(legend.position="right")

然后介绍一个可以对进化树局部放大的函数gzoom。其主要的参数就是focus，他可以支持各种对label筛选的结果。实例：

library("ape")
data(chiroptera)
library("ggtree")
gzoom(chiroptera,grep("Plecotus", chiroptera$tip.label))

以上是基础的绘图，那么我们还可以进行美化，我们就直接借用官网的实例：

groupInfo <- split(chiroptera$tip.label,gsub("_\\w+", "", chiroptera$tip.label))
chiroptera <- groupOTU(chiroptera,groupInfo)
p <- ggtree(chiroptera,aes(color=group)) + geom_tiplab() + xlim(NA, 23)
gzoom(p, grep("Plecotus",chiroptera$tip.label), xmax_adjust=2)

我们也可以进行多个系统进化树的同时展示，需要用到函数facet_wrap，示例如下：

trees <- lapply(c(10, 20, 40), rtree)
class(trees) <- "multiPhylo"
ggtree(trees) + facet_wrap(~.id,scale="free") + geom_tiplab()

那么接下来我们看下更加复杂的多图像可视化，首先是如何将每个样本对应的其他信息以热图形式组合展示。具体我们不挨个剖析了，看下官方给的实例：

在运行实例前，我们需要先加载ape和ggplot2两个包：

beast_file <-system.file("examples/MCC_FluA_H3.tree", package="ggtree")
beast_tree <- read.beast(beast_file)
 
genotype_file <-system.file("examples/Genotype.txt", package="ggtree")
genotype <- read.table(genotype_file,sep="\t", stringsAsFactor=F)
colnames(genotype) <-sub("\\.$", "", colnames(genotype))
p <- ggtree(beast_tree,mrsd="2013-01-01") + geom_treescale(x=2008, y=1, offset=2)
p <- p + geom_tiplab(size=2)
gheatmap(p, genotype, offset=5, width=0.5,font.size=3, colnames_angle=-45, hjust=0) +
scale_fill_manual(breaks=c("HuH3N2","pdm", "trig"), values=c("steelblue","firebrick", "darkgreen"))

现在有个问题就是，感觉X轴并不是那么协调，我们也可以进行进一步美化：

p <- ggtree(beast_tree,mrsd="2013-01-01") + geom_tiplab(size=2,, linesize=.5) +theme_tree2()
pp <- (p +scale_y_continuous(expand=c(0, 0.3))) %>%
   gheatmap(genotype, offset=8, width=0.6, colnames=FALSE) %>%
       scale_x_ggtree()
pp +theme(legend.position="right")

我们这是系统进化树可视化，那么制定离不了序列。当然此包也引入了序列的可视化展示：

fasta <-system.file("examples/FluA_H3_AA.fas", package="ggtree")
msaplot(ggtree(beast_tree), fasta) +theme(legend.position="right")

最后，我们来看个更加复杂的多图的展示，还是要回到我们最初的组合函数faset_plot:

tr <- rtree(30)
 
d1 <- data.frame(id=tr$tip.label,val=rnorm(30, sd=3))
p <- ggtree(tr)
 
p2 <- facet_plot(p,panel="dot", data=d1, geom=geom_point, aes(x=val), color='firebrick')
d2 <- data.frame(id=tr$tip.label,value=abs(rnorm(30, mean=100, sd=50)))
 
facet_plot(p2, panel='bar', data=d2,geom=geom_segment, aes(x=0, xend=value, y=y, yend=y), size=3,color='steelblue') + theme_tree2()