【发布时间】:2016-09-30 14:32:10
【问题描述】:
我正在尝试使用定义的方法here 构建一个系统发育树。由于我有一个庞大的数据集,因此绘制每个人的姓名(在此示例中命名为 1、2、3、4、...)并不提供信息。因此,我希望对提示进行着色以反映个人的物种。个体物种信息存储在一个称为物种的分离矩阵中:
species = data.frame(Ind=c(1 ,2 ,3 ,4 ,5 ),
Spe=c("s1","se1","se2","se2","se3"))
尽管如此,每当我尝试用 fit$tree$label 中的物种替换个体的名称时,唯一的结果就是清空向量。可能这是因为 fit 是 pml 类的对象,它与我的策略产生冲突,但我不确定。有没有办法根据个别物种为树的尖端着色?我的代码现在看起来像这样:
library(phangorn)
#Create a tree from data in fasta format
dat = read.phyDat(file = "myalignment.fasta", format ="fasta")
tree <- pratchet(dat) # parsimony tree
mt <- modelTest(dat, tree=tree, multicore=TRUE)
mt[order(mt$AICc),]
bestmodel <- mt$Model[which.min(mt$AICc)]
env = attr(mt, "env")
fitStart = eval(get("GTR+G+I", env), env)
fit = optim.pml(fitStart, rearrangement = "stochastic", optGamma=TRUE, optInv=TRUE, model="GTR")
bs = bootstrap.pml(fit, bs=100, optNni=TRUE, multicore=TRUE)
#Replace the names with the species...
fit$tree$tip.label <- species[which(species[,1] == fit$tree$tip.label),2]
#If I print fit$tree$tip.label here, the output is factor(0)
#...and create the tree with colored tips
plotBS(midpoint(fit$tree), bs, p = 50, type="p", show.tip.label = FALSE)
tiplabels(pch=19, col = as.factor(fit$tree$tip.label), adj = 2.5, cex = 2)
可重复的示例
使用名称“myalignment.fasta”保存以下内容并运行上面的代码。它应该创建一个玩具示例来玩:
>1
AACCAGGAGAAAATTAA
>2
AAAAA---GAAAATTAA
>3
ACACAGGAGAAAATTAA
>4
AACCTTGAGAAAATTAT
>5
CCTGAGGAGAAAATTAA
【问题讨论】:
-
您应该创建一个最小的reproducible example。包括示例数据(我们无法读取您本地计算机上的文件)。如果我们可以将代码复制/粘贴到 R 中,则更容易测试和提供帮助。如果问题是关于绘图的,请删除任何与问题无关的代码。
-
通常包函数包含示例代码,这些代码使用内置数据集来演示函数的工作方式。此外,
dput()可用于根据我包含的链接重新创建变量。您可以创建最少的数据以用于问题本身。如果你想让人们更容易地帮助你,我只建议这些东西。 -
@MrFlick 我创建了一个最小的“myalignment.fasta”,您可以使用它来重现代码。
标签: r tree bioinformatics phylogeny ape