2016-09-09 10 views
1

おそらくタイトルをより良く言えるかもしれませんが、私は系統樹(クレードに1人のメンバーがいても)のクレードを崩壊したいと思っています。 foo "を表示し、その特定のクレードからドロップされたヒントの数を数え、チップラベルが35個の枝を表示するブランチを作成します。系統樹の位置を維持しながらチップラベルでクレードを縮める

カウント部分は簡単です。しかし、私が使用するとき

drop.tip(rooted.tree,tip=which(rooted.tree$tip.label=='foo'),subtree=TRUE) 

ドロップされたヒントはツリー内でその位置を維持しません。むしろ、彼らはすべて最後にグループ化されています(しかし、適切に数えられます)。とにかくチップラベルでクレードを崩壊させ、木の上にその位置を維持することはありますか?

答えて

0

あなたが持っている個体数はどれくらいですか?この関数はトリックを行う必要があります。

collapse_identical_tips <- function(phy,tip_label){ 
    matching_tips <- which(phy$tip.label==tip_label) 
    nt <- length(phy$tip.label) # number of tips in tree 
    nm <- length(matching_tips) # Number of tips matching the label 
    keep <- numeric(nm) 

    cur_tip <- 1 
    while(cur_tip<=nm){ 
    if(cur_tip == nm){ 
     keep[cur_tip] <- 1 
     break 
    } 
    next_tip <- cur_tip + 1 
    mrca_ <- getMRCA(phy,c(matching_tips[cur_tip],matching_tips[next_tip])) 
    descendants <- getDescendants(phy, mrca_) 
    descendant_tips <- descendants[descendants<=nt] 
    if(all(descendant_tips %in% matching_tips)){ 
     keep[cur_tip] <- 1 
     cur_tip <- cur_tip + length(descendant_tips) 
    }else{ 
     keep[cur_tip] <- 1 
     cur_tip <- cur_tip + 1 
    } 
    } 
    to_drop <- matching_tips[!keep] 
    new_phy <- drop.tip(phy,to_drop) 
    return(new_phy) 
} 

テストそれを:

tc <- "(A:3.135206161,(B:2.757615245,(((C:0.5796267872,((foo:0.1917981792,(foo:0.08246947568,foo:0.08246947568):0.1093287035):0.2328473818,(D:0.3107268924,E:0.3107268924):0.1139186686):0.1549812262):0.3387382152,F:0.9183650024):1.172666972,(((G:0.02437174382,H:0.02437174382):0.4727952475,foo:0.4971669913):0.8701228492,(foo:1.124632261,(foo:0.6503599778,foo:0.6503599778):0.4742722831):0.2426575797):0.7237421344):0.6665832698):0.3775909166);" 
phy <- read.tree(textConnection(tc)) 

plot(phy) 

enter image description here

plot(collapse_identical_tips(phy,"foo")) 

enter image description here

+0

あなたの方法を動作します。しかし、それを使用して、私はどのように崩壊した後の 'foos'を数えるか分からない。たとえば、ヒントを「2 foo、1 foo、1 foo、2 foo's」としたいと思います。どのようにしてコードを編集してこれを行うか、それをどうやってやっていくかについてのアイディアを教えてください。 –

関連する問題