2016-09-30 8 views
0

私は、hereと定義された方法を使ってブートストラップを用いて系統樹を構築しようとしています。私が巨大なデータセットを持っているので、各個人の名前をプロットすること(この例のために1,2,3,4、...という名前が付けられています)は有益ではありません。したがって、私は個人の種を反映するためにヒントを色付けしたい。個々の種の情報は種と呼ばれる、分離行列に格納されています個体の種によって系統樹の先端を着色する

それでも
species = data.frame(Ind=c(1 ,2 ,3 ,4 ,5), 
        Spe=c("s1","se1","se2","se2","se3")) 

、私はフィット$ツリー$ラベルにその種を持つ個人の名前を交換しようとしたときに、結果だけがベクトルを先取りしています。おそらくフィットはクラス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 
+0

あなたは、最小限の[再現性の例](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example)を作成する必要があります。サンプルデータを含める(ローカルマシン上のファイルを読み取ることはできません)。コードをRにコピー/ペーストすることができれば、テストして助けが簡単です。質問がプロットに関するものである場合は、質問に関連しないコードをすべて削除します。 – MrFlick

+0

多くの場合、パッケージ関数には、組み込みデータセットを使用して関数の仕組みを示すサンプルコードが含まれています。また、 'dput()'を使って、私が含むリンクごとに変数を再作成することができます。質問自体に使用する最小限のデータを作成することができます。人々があなたを助けやすいようにしたい場合にのみ、私はこれらのことを提案します。 – MrFlick

+0

@MrFlickコードを再現するのに使うことができる最小の "myalignment.fasta"を作成しました。 – j91

答えて

1

をだからあなたの問題は、色ではなくラベルとしてポイントをプロットリストであるならば、あなたは

tiplabels(pch=19, cex=2, 
    col = species$Spe[match(fit$tree$tip.label, species$Ind)]) 

を行うことができます私はちょうど作るためにgynmasticsのビットを行いますラベルが種に一致することを確認してください(この例では同じ順序です)。

enter image description here

関連する問題