天际线绘图的错误条？

Question

我创建了一个共识树( newick格式)，并使用ape包中的函数对其进行了重新分支。

有没有人知道是否可以将误差条添加到天际线图中？

干杯!输入：

>write.tree(ccc)
[1] "(13.1,43.1,28.1,21.1,20.1,4.1,37.1,23.1,29.1,15.1,36.1,40.1,26.1,16.1,33.1,10.1,18.1,6.1,11.1,34.1,2.1,38.1,35.1,8.1,42.1,22.1,5.1,27.1,39.1,17.1,30.1,31.1,41.1,7.1,25.1,3.1,12.1,14.1,1.1,24.1,9.1,32.1,19.1);"

分支树和最终输入：

> write.tree(new_tree)
[1] "   (13.1:1.244090638,43.1:0.1755118382,28.1:4.442071925,21.1:4.866247957,20.1:3.863557777,4.1:0.08206463186,37.1:1.180065627,23.1:4.778251834,29.1:4.65377018,15.1:3.726316827,36.1:3.799621998,40.1:1.707243007,26.1:3.655532246,16.1:2.436848865,33.1:1.315581813,10.1:0.5973169219,18.1:2.561718575,6.1:1.409437305,11.1:4.000977813,34.1:2.769218051,2.1:1.143412833,38.1:2.636477078,35.1:3.435094416,8.1:0.3714309109,42.1:4.429772993,22.1:0.2805716533,5.1:2.991251546,27.1:3.269689628,39.1:1.192241808,17.1:1.866855259,30.1:1.363407132,31.1:3.150236446,41.1:4.079005712,7.1:2.695786237,25.1:2.867184281,3.1:3.351797838,12.1:1.958669939,14.1:3.119812493,1.1:4.864444738,24.1:1.734493783,9.1:3.283970281,32.1:3.055829309,19.1:0.8193949074);"

require(ape)
new_tree<-compute.brlen(ccc,runif,min=0,max=5)
sk1 <- skyline(new_tree)   
sk2 <- skyline(new_tree, 0.0119) 
plot(sk1, show.years=TRUE, subst.rate=0.0023, present.year = 2011, col="red")
lines(sk2,  show.years=TRUE, subst.rate=0.0023, present.year = 2011)
legend(.15,500, c("classic", "generalized"), col="red",lty=1)

Answer 1

您可以从plotrix包中使用plotCI。追加plot.skyline树来实现错误栏会很有趣……

不管怎样，这是我用来生成下面图片的代码。我使用的是ape包附带的另一个数据集。你得处理好细节。

library(ape)
library(plotrix)
data(bird.orders)

new_tree <- compute.brlen(bird.orders, runif, min = 0, max = 5)
sk1 <- skyline(new_tree)   

subst.rate <- 0.0023
year <- 2011
m <- sk1$population.size
lm <- length(m)

plot(sk1, show.years=TRUE, subst.rate = subst.rate, present.year = year, col="red")

plotCI(x = (-c(0, sk1$time))/subst.rate + year, y = c(m, m[lm]),
        uiw = runif(length(sk1$time)+1, min = 1, max = 30), # some random deviation
        #liw = runif(length(sk1$time)+1, min = 1, max = 30), # this one can be missing, see ?plotCI
        add = TRUE,
        pt.bg = NA
)

legend(.15,500, c("classic", "generalized"), col="red",lty=1)

​

​