#### ggtree可以直接读取的数据格式

Newick (via ape)
Nexus (via ape)
New Hampshire eXtended format (NHX) （http://home.cc.umanitoba.ca/~psgendb/doc/atv/NHX.pdf）
Jplace （http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0031009）
Phylip

以及以下软件输出的数据:

BEAST2
EPA3
HYPHY5
PAML1
PHYLDOG6
pplacer4
r8s7
RAxML8
RevBayes9

####################################################
用以下函数读取
read.beast()       ## for parsing output of BEASE
read.codeml()      ## for parsing output of CODEML (rst and mlc files)
read.codeml_mlc()  ## for parsing mlc file (output of CODEML)
read.hyphy()       ## for parsing output of HYPHY
read.jplace()      ## for parsing jplace file including output from EPA and pplacer
read.nhx()         ## for parsing NHX file including output from PHYLODOG and RevBayes
read.paml_rst()    ## for parsing rst file (output of BASEML and CODEML)
read.r8s()         ## for parsing output of r8s
read.raxml()       ## for parsing output of RAxML

#####################################################
###### 定义的 S4 Classes #####################################
apeBootstrap ##for bootstrap analysis of ape::boot.phylo()10, output of apeBoot() defined in ggtree
beast        ##for storing output of read.beast()
codeml       ## for storing output of read.codeml()
codeml_mlc   ##for storing output of read.codeml_mlc()
hyphy        ## for storing output of read.hyphy()
jplace       ## for storing output of read.jplace()
nhx          ## for storing output of read.nhx()
paml_rst     ## for rst file obtained by PAML, including BASEML and CODEML.
phangorn     ## for storing ancestral sequences inferred by R package phangorn11, output of phyPML defined in ggtree
r8s          ## for storing output of read.r8s()
raxml        ## for storing output of read.raxml()

### 其他支持的S4 类
###

phylo, 
multiPhylo (defined by ape10), 
phylo4 (defined by phylobase) 
obkData (defined in OutbreakTools) and 
phyloseq (defined in phyloseq).

### 获取任何能够在进化树上显示的参数
get.fields

### 显示进化树
ggtree(tree_object)

### 将任何格式的进化树转换为 phylo格式
get.tree

### 界定分类单元
groupOTU()
groupClade()

### 将进化树转换为 data.frame
fortify

############################################
#########3 举例： BEAST

library(ggtree)
file <- system.file("extdata/BEAST", "beast_mcc.tree", package="ggtree")
beast <- read.beast(file)
beast

get.fields(beast)
ggtree(beast, ndigits = 2, branch.length = 'none') + geom_text(aes(x = branch, label =length_0.95_HPD), vjust = -.5, color = 'firebrick')

######### 导入的数据通过 fortify函数能够转换为 data.frame
beast_data <- fortify(beast)
head(beast_data)


################ PAML
brstfile <- system.file("extdata/PAML_Baseml", "rst", package="ggtree")
brst <- read.paml_rst(brstfile)
brst

p <- ggtree(brst) +
   geom_text(aes(x = branch, label = marginal_AA_subs), vjust=-.5, color='steelblue') +
   theme_tree2()
print(p)

################# CODEML
crstfile <- system.file("extdata/PAML_Codeml", "rst", package="ggtree")
crst <- read.paml_rst(crstfile)
crst

##### %<% 类似于格式刷， 将某一进化树绘制时的参数， 转移到要绘制的进化树中
p %<% crst

################ 
#### CODEML mlc文件
mlcfile <- system.file("extdata/PAML_Codeml", "mlc", package="ggtree")
mlc <- read.codeml_mlc(mlcfile)
mlc

ggtree(mlc) + geom_text(aes(x=branch, label=dN_vs_dS), color='blue', vjust=-.2)

#### 获取有用的信息
get.fields(mlc)

#### 指定需要调整的枝长

ggtree(mlc, branch.length = "dN_vs_dS", aes(color=dN_vs_dS)) +
   scale_color_continuous(name='dN/dS', limits=c(0, 1.5),
                          oob=scales::squish, low="darkgreen", high="red")+
   theme_tree2(legend.position=c(.9, .5))

#### CODEML 的 rst以及 mlc文件

ml <- read.codeml(crstfile, mlcfile)
ml
#### 
head(fortify(ml))

#### 读取 HYPHY文件并绘制进化树

nwk <- system.file("extdata/HYPHY", "labelledtree.tree", package="ggtree")
ancseq <- system.file("extdata/HYPHY", "ancseq.nex", package="ggtree")
tipfas <- system.file("extdata", "pa.fas", package="ggtree")
hy <- read.hyphy(nwk, ancseq, tipfas)
hy

ggtree(hy) + geom_text(aes(x=branch, label=AA_subs), vjust=-.5)


##### 读取r8s生成的文件并绘制进化树
r8s <- read.r8s(system.file("extdata/r8s", "H3_r8s_output.log", package="ggtree"))
ggtree(r8s, branch.length="TREE", mrsd="2014-01-01") + theme_tree2()

ggtree(get.tree(r8s), aes(color=.id)) + facet_wrap(~.id, scales="free_x")

##### 读取RAxML生成的bootstrap文件并绘制进化树

raxml_file <- system.file("extdata/RAxML", "RAxML_bipartitionsBranchLabels.H3",package="ggtree")
raxml <- read.raxml(raxml_file)
ggtree(raxml) + geom_label(aes(label=bootstrap, fill=bootstrap)) + geom_tiplab() +
   scale_fill_continuous(low='darkgreen', high='red') + theme_tree2(legend.position='right')

##### 读取NHX New Hampshire eXtended文件并绘制进化树
nhxfile <- system.file("extdata", "ADH.nhx", package="ggtree")
nhx <- read.nhx(nhxfile)
ggtree(nhx) + geom_tiplab() + geom_point(aes(color=S), size=5, alpha=.5) + 
   theme(legend.position="right") +
   geom_text(aes(label=branch.length, x=branch), vjust=-.5) + 
   xlim(NA, 0.3)

##### 读取PHYLIP文件并绘制进化树以及序列比对图
phyfile <- system.file("extdata", "sample.phy", package="ggtree")
phylip <- read.phylip(phyfile)
phylip

ggtree(phylip) + geom_tiplab()

msaplot(phylip, offset=1)

##### jplace格式数据同时保存自定义数据

jpf <- system.file("extdata/sample.jplace",  package="ggtree")
jp <- read.jplace(jpf)
print(jp)
## get only best hit
get.placements(jp, by="best")

## get all placement
get.placements(jp, by="all")

########################################################################
############ 进化树合并 ### 分类单元匹配 ##################################
############ 作者讲的不太清楚 
beast_file <- system.file("examples/MCC_FluA_H3.tree", package="ggtree")
beast_tree <- read.beast(beast_file)

rst_file <- system.file("examples/rst", package="ggtree")
mlc_file <- system.file("examples/mlc", package="ggtree")
codeml_tree <- read.codeml(rst_file, mlc_file)

merged_tree <- merge_tree(beast_tree, codeml_tree)

merged_tree

#################################################################
#### 读取 jplace 文件， 并将保存在csv文件种的性状信息

tree <- system.file("extdata", "pa.nwk", package="ggtree")
data <- read.csv(system.file("extdata", "pa_subs.csv", package="ggtree"),stringsAsFactor=FALSE)
print(tree)

outfile <- tempfile()
write.jplace(tree, data, outfile)
jp <- read.jplace(outfile)

ggtree(jp) + 
   geom_text(aes(x=branch, label=subs), color="purple", vjust=-1, size=3) + 
   geom_text(aes(label=gc), color="steelblue", hjust=-.6, size=3) +
   geom_tiplab()

#############################################################################
########## 绘制进化树
library("ggtree")
nwk <- system.file("extdata", "sample.nwk", package="ggtree")
tree <- read.tree(nwk)

ggplot(tree, aes(x, y)) + geom_tree() + theme_tree()
ggtree(tree, color="firebrick", size=1, linetype="dotted")

### ggtree会自动将进化树转换为 ladderized， 如果要取消 ladderized， 则应该改为 ladderize
ggtree(tree, ladderize=FALSE)

### 不要枝长， 只查
ggtree(tree, branch.length="none")

########## 进化树显示的方式 ############

矩形 rectangular (by default)
偏斜 slanted
扇形/圆形 fan or circular

ggtree(tree) + ggtitle("(Phylogram) rectangular layout")
ggtree(tree, layout="slanted") + ggtitle("(Phylogram) slanted layout")
ggtree(tree, layout="circular") + ggtitle("(Phylogram) circular layout")

###### 只显示拓扑结构
### 矩形
ggtree(tree, branch.length='none') + ggtitle("(Cladogram) rectangular layout")

### 偏斜
ggtree(tree, layout="slanted", branch.length='none') + ggtitle("(Cladogram) slanted layout")

### 圆形
ggtree(tree, layout="circular", branch.length="none") + ggtitle("(Cladogram) circular layout")

##### 无根树
ggtree(tree, layout="unrooted") + ggtitle("unrooted layout")

#######################################################################################
########### Ultrametric Tree 绘图会自动添加标尺 ###################
ggtree(tree) + geom_treescale()

#### 更改标尺的长短颜色， 粗细 
ggtree(tree)+geom_treescale(x=0, y=12, width=6, color='red')

#################### 
ggtree(tree)+geom_treescale(fontsize=8, linesize=2, offset=-1)

###### 添加时间轴 依靠 theme_tree2()
ggtree(tree) + theme_tree2()

#########################################################################################
###### 节点以及分类单元的标识， 根据是否为 tiplabel以及node分别显示
ggtree(tree)+geom_point(aes(shape=isTip, color=isTip), size=3)

#### Node Point参数 以及 tippoint分别显示
p <- ggtree(tree) + geom_nodepoint(color="#b5e521", alpha=1/4, size=10)
p + geom_tippoint(color="#FDAC4F", shape=8, size=3)

#### Tiplabel的显示
p + geom_tiplab(size=3, color="purple")

#### 圆形进化树的标签
ggtree(tree, layout="circular") + geom_tiplab(aes(angle=angle), color='blue')

#### 进化树tip label显示到branch上
p + geom_tiplab(aes(x=branch), size=3, color="purple", vjust=-0.3)

#### %<% 若之前已经有ggtree object完成了进化树的绘制， 那么可以用 %>% 将前一个进化树的属性转移到当前进化树中。 %<% 的作用相当于word中的格式刷。 
p %<% rtree(50)

#################################################################
######### ggtree中有两个主题 theme， theme_tree2() 和 theme_tree()
### theme_tree2 与 后者的区别在于 增加了时间坐标轴。 
### 无坐标轴
ggtree(rtree(30), color="red") + theme_tree()

### 有坐标轴
ggtree(rtree(30), color="red") + theme_tree2()

### 改变背景颜色
ggtree(rtree(30), color="red") + theme_tree("steelblue")

### 
ggtree(rtree(20), color="white") + theme_tree("#b5e521")

#################################################################
############# 显示多个进化树 ######################################
trees <- lapply(c(10, 20, 40), rtree)
class(trees) <- "multiPhylo"
ggtree(trees) + facet_wrap(~.id, scale="free") + geom_tiplab()

#################################################################
btrees <- read.tree(system.file("extdata/RAxML", "RAxML_bootstrap.H3", package="ggtree"))
ggtree(btrees) + facet_wrap(~.id, ncol=10)


#################################################################
########### 显示bootstrap进化树作为背景 #################################
p <- ggtree(btrees, layout="rectangular",   color="lightblue", alpha=.3)

best_tree <- read.tree(system.file("extdata/RAxML", "RAxML_bipartitionsBranchLabels.H3",package="ggtree"))
df <- fortify(best_tree, branch.length = 'none')
p + geom_tree(data=df, color='firebrick')

################################################################################
############# 进化树的操作 ##########################################

#### 对进化树操作之前， 必须要定位进化树内部的节点， 每一个节点都有一个顺序号
#### ggtree提供了显示内部节点号的函数

nwk <- system.file("extdata", "sample.nwk", package="ggtree")
tree <- read.tree(nwk)
ggtree(tree) + geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + geom_tiplab()

#### 用MRCA显示节点编号
MRCA(tree, tip=c('A', 'E'))

########### 显示某一个节点下的整个类群 clade
p <- ggtree(tree)
viewClade(p + geom_tiplab(), node=21)

########### 基于内部节点号， 划定clade， 按照clade划分颜色
tree <- groupClade(tree, node=21)
ggtree(tree, aes(color=group, linetype=group))

##########  不同group分不同颜色显示， 显示tiplabel时， 不同的组别用不同的颜色表示
tree <- groupClade(tree, node=c(21, 17))
ggtree(tree, aes(color=group, linetype=group)) + geom_tiplab(aes(subset=(group==2)))

########### 手工指定类群的分组，并且表示为不同的颜色 
cls <- list(c1=c("A", "B", "C", "D", "E"),
           c2=c("F", "G", "H"),
           c3=c("L", "K", "I", "J"),
           c4="M")

#### s输入对象未 ggtree时
tree <- groupOTU(tree, cls)
library("colorspace")
ggtree(tree, aes(color=group, linetype=group)) + geom_tiplab() +
    scale_color_manual(values=c("black", rainbow_hcl(4))) + theme(legend.position="right")

#### 输入对象为
p <- ggtree(tree)
groupOTU(p, LETTERS[1:5]) + aes(color=group) + geom_tiplab() +scale_color_manual(values=c("black", "firebrick"))


##### groupOTU划分组
library("ape")
data(chiroptera)
groupInfo <- split(chiroptera$tip.label, gsub("_\w+", "", chiroptera$tip.label))
chiroptera <- groupOTU(chiroptera, groupInfo)
ggtree(chiroptera, aes(color=group), layout='circular') + geom_tiplab(size=1,aes(angle=angle))


##### 类群折叠与展开
cp <- ggtree(tree) %>% collapse(node=21)

### 折叠节点号为21的类群
cp + geom_point2(aes(subset=(node == 21)), size=5, shape=23, fill="steelblue")

### 展开
cp %>% expand(node=21)

### 折叠以及展开多个节点
p1 <- ggtree(tree)
p2 <- collapse(p1, 21) + geom_point2(aes(subset=(node==21)), size=5, shape=23, fill="blue")
p3 <- collapse(p2, 17) + geom_point2(aes(subset=(node==17)), size=5, shape=23, fill="red")
p4 <- expand(p3, 17)
p5 <- expand(p4, 21)

multiplot(p1, p2, p3, p4, p5, ncol=5)


#### 类群的缩放
multiplot(ggtree(tree) + geom_hilight(21, "steelblue"),
         ggtree(tree) %>% scaleClade(21, scale=0.3) + geom_hilight(21, "steelblue"),
         ncol=2)

#### 指定节点编号， 高亮某一个类群
ggtree(tree) + geom_hilight(21, "steelblue")

#### 类群clade旋转180度 (此例子不能运行)
tree <- groupClade(tree, c(21, 17))
### p <- ggtree(tree, aes(color=group)) + scale_color_manual(values=c("black", "firebrick", "steelblue"))
#### p2 <- rotate(p, 21) %>% rotate(p, 17)
multiplot(p, p2, ncol=2)

########################################################################
##### 两个分类群位置互换
multiplot(p, flip(p, 17, 21), ncol=2)

############################################################################
############################################################################
##### 进化树
beast_file <- system.file("examples/MCC_FluA_H3.tree", package="ggtree")
beast_tree <- read.beast(beast_file)
beast_tree


##### 绘制进化树
p1 <- ggtree(beast_tree, mrsd='2013-01-01') + theme_tree2() +
   ggtitle("Divergence time")
p2 <- ggtree(beast_tree, branch.length = 'rate') + theme_tree2() +
   ggtitle("Substitution rate")
multiplot(p1, p2, ncol=2)

###### 局部放大某一个分类群 clade
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)

###### 各分支按照条件设色
ggtree(beast_tree, aes(color = rate)) +
   scale_color_continuous(low='darkgreen', high='red') +
   theme(legend.position="right")

##### 标注某一个clade， 在右侧画竖线
set.seed(2015-12-21)
tree = rtree(30)
p <- ggtree(tree) + xlim(NA, 6)

p + geom_cladelabel(node=45, label="test label") + 
   geom_cladelabel(node=34, label="another clade") 

##### 右侧标注用的竖线可以对齐
p + geom_cladelabel(node=45, label="test label", align=TRUE, offset=.5) +
   geom_cladelabel(node=34, label="another clade", align=TRUE, offset=.5)

##### 改变标注竖线的颜色
p+geom_cladelabel(node=45, label="test label", align=T, color='red') +
   geom_cladelabel(node=34, label="another clade", align=T, color='blue')

##### 旋转标注类群的文字
p+geom_cladelabel(node=45, label="test label", align=T, angle=270, hjust='center',offset.text=.5) + geom_cladelabel(node=34, label="another clade", align=T, angle=45)     

##### 改变标注竖线的粗细， 以及文字的大小
p+geom_cladelabel(node=45, label="test label", align=T, angle=270, hjust='center',offset.text=.5, barsize=1.5) +
   geom_cladelabel(node=34, label="another clade", align=T, angle=45, fontsize=8) 

##### 为文字添加边框， 并设置文本框的背景颜色
p+ geom_cladelabel(node=34, label="another clade", align=T, geom='label', fill='lightblue') 

########################################################################
################# 对类群进行高亮显示 ######################################
nwk <- system.file("extdata", "sample.nwk", package="ggtree")
tree <- read.tree(nwk)
ggtree(tree) + geom_hilight(node=21, fill="steelblue", alpha=.6) +
   geom_hilight(node=17, fill="darkgreen", alpha=.6)

#### 圆形树的类群高亮
ggtree(tree, layout="circular") + geom_hilight(node=21, fill="steelblue", alpha=.6) +
   geom_hilight(node=23, fill="darkgreen", alpha=.6)

#########################################################################
#### 标注 bootstrap数值 
#### 用ape生成一个NJ树， 并且进行bootstrap， 将bootstrap的结果在NJ树的节点上显示。 

library(ape)
data(woodmouse)
d <- dist.dna(woodmouse)
tr <- nj(d)
bp <- boot.phylo(tr, woodmouse, function(xx) nj(dist.dna(xx)))

##### 
tree <- apeBoot(tr, bp)  ### merge phylo and output of boot.phylo to 'apeBootstrap' object
ggtree(tree) + geom_label(aes(label=bootstrap)) + geom_tiplab()


##### 对phangorn生成进化树的标注
library(phangorn)
treefile <- system.file("extdata", "pa.nwk", package="ggtree")
tre <- read.tree(treefile)
tipseqfile <- system.file("extdata", "pa.fas", package="ggtree")
tipseq <- read.phyDat(tipseqfile,format="fasta")
fit <- pml(tre, tipseq, k=4)
#### fit是设定条件生成的最优进化树
fit <- optim.pml(fit, optNni=FALSE, optBf=T, optQ=T,
                optInv=T, optGamma=T, optEdge=TRUE,
                optRooted=FALSE, model = "GTR")

#### 将pml格式的数据转换为 ggtree能够识别的类型
phangorn <- phyPML(fit, type="ml")
ggtree(phangorn) + geom_text(aes(x=branch, label=AA_subs, vjust=-.5))

#### ggtree可以读取如下软件生成的进化树， 并且同时保存这些软件输出结果关于进化速率等相应的信息
BEAST
EPA
HYPHY
PAML
PHYLDOG
pplacer
r8s
RAxML
RevBayes

########################################################################
######## 将性状信息与进化树匹配绘图 ###############

nwk <- system.file("extdata", "sample.nwk", package="ggtree")
tree <- read.tree(nwk)
p <- ggtree(tree)

##### 生成随机数据
dd <- data.frame(taxa  = LETTERS[1:13], 
                place = c(rep("GZ", 5), rep("HK", 3), rep("CZ", 4), NA),
                value = round(abs(rnorm(13, mean=70, sd=10)), digits=1))
## you don't need to order the data
## data was reshuffled just for demonstration
dd <- dd[sample(1:13, 13), ]
row.names(dd) <- NULL

##### 数据框匹配字符
##### %<+% 运算符的要求
##### 左侧为进化树
##### 右侧为 data.frame, 要求是该进化树的第一列种的物种名必须能与进化树的tip.label匹配上。 

p <- p %<+% dd + geom_tiplab(aes(color=place)) + 
      geom_tippoint(aes(size = value, shape = place, color = place), alpha = 0.25)
      
p + theme(legend.position="right")

##### 将采集地点以及性状值写在枝长上
p + geom_text(aes(color=place, label=place), hjust=1, vjust=-0.4, size=3) +
 geom_text(aes(color=place, label=value), hjust=1, vjust=1.4, size=3)

########################################################################
########################################################################
########################################################################
####################### 在进化树右侧绘制热图 ###############################

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)
head(genotype)

p <- ggtree(beast_tree, mrsd="2013-01-01") + geom_treescale(x=2008, y=1)
p <- p + geom_tiplab(size=3)
gheatmap(p, genotype, offset = 2, width=0.5)

########################################################################
#### scale_x_ggtree() 保证右侧矩阵的文字说明能正常显示。 
p <- ggtree(beast_tree, mrsd="2013-01-01") + geom_tiplab(size=3, align=TRUE) + theme_tree2()
pp <- (p + scale_y_continuous(expand=c(0, 0.3))) %>%
   gheatmap(genotype, offset=4, width=0.5, colnames=FALSE) %>%
       scale_x_ggtree()
pp + theme(legend.position="right")

######################## msaplot ###### 
##### 进化树 以及 fasta比对矩阵 
fasta <- system.file("examples/FluA_H3_AA.fas", package="ggtree")
msaplot(ggtree(beast_tree), fasta) 

######################## 扇形的 msaplot #####
msaplot(ggtree(beast_tree), fasta, window=c(150, 200)) + coord_polar(theta='y')

######################## 在进化树的总图上高亮显示某一个clade： geom_hilight
######################## 并且在右侧的大图中放大该clade viewClade subview
set.seed(2016-01-04)
tr <- rtree(30)
tr <- groupClade(tr, node=45)
p <- ggtree(tr, aes(color=group)) + geom_tippoint()
p1 <- p + geom_hilight(node=45)
p2 <- viewClade(p, node=45) + geom_tiplab()
subview(p2, p1+theme_transparent(), x=2.3, y=28.5)


#############################################################################
########## 用inset为内部节点添加直方图， 以显示每个节点在性状上所占比例  ##############
#############################################################################

set.seed(2015-12-31)
tr <- rtree(15)
p <- ggtree(tr)

a <- runif(14, 0, 0.33)
b <- runif(14, 0, 0.33)
c <- runif(14, 0, 0.33)
d <- 1 - a - b - c
dat <- data.frame(a=a, b=b, c=c, d=d)
## input data should have a column of `node` that store the node number 
### 数据必须有一列名为node， 并且值为节点数
dat$node <- 15+1:14

## cols parameter indicate which columns store stats (a, b, c and d in this example)
## cols表示dat的1:4列储存了用于绘图的数据
bars <- nodebar(dat, cols=1:4)

inset(p, bars)

## 柱状图的大小和颜色
inset(p, bars, width=.03, height=.06)

### 柱状图的排列方式， 改为并排排列
bars2 <- nodebar(dat, cols=1:4, position='dodge',
                color=c(a='blue', b='red', c='green', d='cyan'))
p2 <- inset(p, bars2, x='branch', width=.03, vjust=-.3)
print(p2)


##############################################################
### 节点上的饼图 

pies <- nodepie(dat, cols=1:4, alpha=.6)
inset(p, pies) 

inset(p, pies, hjust=-.06) ### 饼图向右平移

##############################################################
#### 同时用 pie和bar标注内部节点

pies_and_bars <- bars2
pies_and_bars[9:14] <- pies[9:14]
inset(p, pies_and_bars)

##############################################################
##### 用 inset函数用tip添加对应的箱线图 boxplot， 以展示数据分布

d <- lapply(1:15, rnorm, n=100)
ylim <- range(unlist(d))
bx <- lapply(d, function(y) {
   dd <- data.frame(y=y)
   ggplot(dd, aes(x=1, y=y)) + geom_boxplot() + ylim(ylim) + theme_inset()
})
names(bx) <- 1:15
inset(p, bx, width=.03, height=.1, hjust=-.05)


#########################################################################
### 混合标注， 未能成功
### p2 <- inset(p, bars2, x='branch', width=.03, vjust=-.4)
### p2 <- inset(p2, pies, x='branch', vjust=.4)
### bx2 <- lapply(bx, function(g) g+coord_flip())
### inset(p2, bx2, width=.2, height=.03, vjust=.04, hjust=p2datadatax[1:15]-4) + xlim(NA, 4.5)


########################################################################
###### ggplot统计与进化树绘制 ################ 
tr <- rtree(30)
df <- fortify(tr)
df$tipstats <- NA
d1 <- df
d2 <- df
d2$tipstats[d2$isTip] <- abs(rnorm(30))
d1$panel <- 'Tree'
d2$panel <- 'Stats'
d1$panel <- factor(d1$panel, levels=c("Tree", "Stats"))
d2$panel <- factor(d2$panel, levels=c("Tree", "Stats"))

p <- ggplot(mapping=aes(x=x, y=y)) + facet_grid(.~panel, scale="free_x") + theme_tree2()
p+geom_tree(data=d1) + geom_point(data=d2, aes(x=tipstats)) 

########################################################################
##### 从 phylopic数据库上下载剪影图像， 并标在进化树上

pp <- ggtree(tree) %>% phylopic("79ad5f09-cf21-4c89-8e7d-0c82a00ce728", color="steelblue",alpha = .3)
print(pp)

#### 首先生成一套模拟数据simdat, 通过Volcano数据随机提取1000个点, 用raster的extract函数提取

library(raster)
rrr <- raster(t(volcano), xmn=0, xmx= 870 , ymn=0, ymx=610)
x11()
plot(rrr, col = terrain.colors(100))

set.seed(12345)
loc <- data.frame(x = round(runif(1000, min = 0, max = 870), 2), y = round(runif(1000,min = 0, max = 610), 2))
value <- extract(rrr, loc)
simdat <- data.frame(loc, alt = value) ### The simulated data


#### Krig准备工作开始, 导入程序包
library(sp)
library(gstat)

#### simdat需要转换为sp格式, 需要指定坐标
coordinates(simdat) <- ~x+y
### bubble(simdat, "alt")
spplot(simdat, "alt") ### 查看每个点的值 

#### 创建克里格所用的网格, 20m方格的中心位置
xgrid <- seq(ceiling(min(simdat$x)), floor(max(simdat$x)), by = 20) 
ygrid <- seq(ceiling(min(simdat$y)), floor(max(simdat$y)), by = 20) 
basexy <- expand.grid(xgrid, ygrid)
#### plot(y ~ x, basexy)
colnames(basexy) <- c("x", "y")
coordinates(basexy) <- ~x+y
gridded(basexy) <- TRUE

### 因Kriging需要基于 Variogram以及相应的模型.Variogram横轴是距离, 纵轴是方差. 做模型拟合时, 应先查看 Variogram. Variogram根据点与点之间的距离, 分成若干段, 并计算每一段的方差, 再通过拟合Variogram模型, 获得Krig所需要的参数. 拟合Variogram模型时, 需要提供模型的初始值,再进行迭代优化, 否则很难获得模型的准确参数. 

######### 若所进行的是普通克里格插值, 则需要插值的变量只与自身取样点的距离相关,用公式表示为 alt~1, 这里的~1是variogram的规定 (参见 ?variogram). 
vgm1 <- variogram(alt~1, simdat)
plot(vgm1, plot.numbers = TRUE) 

### variogram的模型常包括
### Circular
### Spherical
### Exponential
### Gaussian
### Linear
### gstat程序包提供的模型可以通过 vgm() 查看

m <- fit.variogram(vgm1,vgm(.01,"Sph",100,300))
plot(vgm1, model=m)

### 设定vgm模型的初始值
m <- vgm(.01,"Sph",300, 0.2)

### 进行克里格
krige_res <- krige(alt~1, simdat, basexy, model = m)

### 查看克里格插值的结果
spplot(krige_res, zcol = "var1.pred", main = "Predictions of altitude based on the randomly sampled data", col.regions = terrain.colors(100))