Come comprimere i rami in un albero filogenetico dall'etichetta nei loro nodi o foglie?

Question

Ho costruito un albero filogenetico per una famiglia proteica che può essere suddivisa in diversi gruppi, classificandoli ognuno per tipo di recettore o tipo di risposta. I nodi nell'albero sono etichettati come il tipo di recettore.

Nell'albero filogenetico posso vedere che le proteine ​​che appartengono allo stesso gruppo o tipo di recettore si sono raggruppate insieme negli stessi rami. Quindi vorrei raggruppare questi rami che hanno etichette in comune, raggruppandoli in base a un determinato elenco di parole chiave.

Il comando sarebbe qualcosa di simile a questo:

./collapse_tree_by_label -f phylogenetic_tree.newick -l list_of_labels_to_collapse.txt -o collapsed_tree.eps (o pdf)

mio list_of_labels_to_collapse.txt sarebbe come questo : A B C D

mio albero Newick sarebbe come questo: (A_1: 0.05, A_2: 0.03, A_3: 0.2, A_4: 0.1): 0,9, (((B_1: 0.05, q_2 : 0.02, B_3: 0.04): 0,6, (C_1: 0.6, C_2: 0.08): 0. 7): 0.5, (D_1: 0.3, D_2: 0.4, D_3: 0.5, D_4: 0.7, D_5: 0.4): 1.2)

L'immagine in uscita senza collassare è così:

L'uscita immagine collassamento dovrebbe essere come questa (collapsed_tree.eps): http://i.stack.imgur.com/TLXd0.png

la larghezza dei triangoli dovrebbe rappresentare la lunghezza dei rami, e l'alto dei triangoli deve rappresentare il numero di nodi nel ramo.

Ho giocato con il pacchetto "scimmia" in R. Sono stato in grado di tracciare un albero filogenetico, ma non riesco ancora a capire come crollare i rami con parole chiave nelle etichette:

require("ape") 

Questo caricherà l'albero:

cat("((A_1:0.05,A_2:0.03,A_3:0.2,A_4:0.1):0.9,(((B_1:0.05,B_2:0.02,B_3:0.04):0.6,(C_1:0.6,C_2:0.08):0.7):0.5,(D_1:0.3,D_2:0.4,D_3:0.5,D_4:0.7,D_5:0.4):1.2):0.5);", file = "ex.tre", sep = "\n") 
tree.test <- read.tree("ex.tre") 

qui dovrebbe essere il codice a crollare

Ciò tracciare l'albero:

01.235.164,106 mila

plot(tree.test) 

Answer 1

L'albero così come è memorizzato in R ha già i suggerimenti memorizzati come politestie. È solo questione di tracciare l'albero con triangoli che rappresentano le politomie.

Non v'è alcuna funzione in ape per fare questo, che io sono a conoscenza, ma se si fa confusione con la funzione di tracciato un po 'si può tirare fuori

# Step 1: make edges for descendent nodes invisible in plot: 
groups <- c("A", "B", "C", "D") 
group_edges <- numeric(0) 
for(group in groups){ 
    group_edges <- c(group_edges,getMRCA(tree.test,tree.test$tip.label[grepl(group, tree.test$tip.label)])) 
} 
edge.width <- rep(1, nrow(tree.test$edge)) 
edge.width[tree.test$edge[,1] %in% group_edges ] <- 0 


# Step 2: plot the tree with the hidden edges 
plot(tree.test, show.tip.label = F, edge.width = edge.width) 

# Step 3: add triangles 
add_polytomy_triangle <- function(phy, group){ 
    root <- length(phy$tip.label)+1 
    group_node_labels <- phy$tip.label[grepl(group, phy$tip.label)] 
    group_nodes <- which(phy$tip.label %in% group_node_labels) 
    group_mrca <- getMRCA(phy,group_nodes) 

    tip_coord1 <- c(dist.nodes(phy)[root, group_nodes[1]], group_nodes[1]) 
    tip_coord2 <- c(dist.nodes(phy)[root, group_nodes[1]], group_nodes[length(group_nodes)]) 
    node_coord <- c(dist.nodes(phy)[root, group_mrca], mean(c(tip_coord1[2], tip_coord2[2]))) 

    xcoords <- c(tip_coord1[1], tip_coord2[1], node_coord[1]) 
    ycoords <- c(tip_coord1[2], tip_coord2[2], node_coord[2]) 
    polygon(xcoords, ycoords) 
} 

Poi basta avere a scorrere i gruppi per aggiungere i triangoli

for(group in groups){ 
    add_polytomy_triangle(tree.test, group) 
} 

Answer 2

Penso che la sceneggiatura stia finalmente facendo ciò che volevo. Dalla risposta fornita da @CactusWoman, ho modificato un po 'il codice in modo che lo script proverà a trovare l'MRCA che rappresenta il ramo più grande che corrisponde al mio modello di ricerca.Questo risolve il problema di non unire rami non politiomici o di comprimere l'intero albero perché un nodo corrispondente era erroneamente fuori dal ramo corretto.

Inoltre, ho incluso un parametro che rappresenta il limite per il rapporto di abbondanza del motivo in un determinato ramo, quindi possiamo selezionare e comprimere/raggruppare i rami che hanno almeno il 90% dei suggerimenti corrispondenti al modello di ricerca, per esempio.

library(geiger) 
library(phylobase) 
library(ape) 

#functions 
find_best_mrca <- function(phy, group, threshold){ 

    group_matches <- phy$tip.label[grepl(group, phy$tip.label, ignore.case=TRUE)] 
    group_mrca <- getMRCA(phy,phy$tip.label[grepl(group, phy$tip.label, ignore.case=TRUE)]) 
    group_leaves <- tips(phy, group_mrca) 
    match_ratio <- length(group_matches)/length(group_leaves) 

     if(match_ratio < threshold){ 

      #start searching for children nodes that have more than 95% of descendants matching to the search pattern 
      mrca_children <- descendants(as(phy,"phylo4"), group_mrca, type="all") 
      i <- 1 
      new_ratios <- NULL 
      nleaves <- NULL 
      names(mrca_children) <- NULL 

      for(new_mrca in mrca_children){ 
       child_leaves <- tips(tree.test, new_mrca) 
       child_matches <- grep(group, child_leaves, ignore.case=TRUE) 
       new_ratios[i] <- length(child_matches)/length(child_leaves) 
       nleaves[i] <- length(tips(phy, new_mrca)) 
       i <- i+1 
      } 



      match_result <- data.frame(mrca_children, new_ratios, nleaves) 


      match_result_sorted <- match_result[order(-match_result$nleaves,match_result$new_ratios),] 
      found <- numeric(0); 

      print(match_result_sorted) 

      for(line in 1:nrow(match_result_sorted)){ 
       if(match_result_sorted$ new_ratios[line]>=threshold){ 
        return(match_result_sorted$mrca_children[line]) 
        found <- 1 
       } 

      } 

      if(found==0){return(found)} 
     }else{return(group_mrca)} 




} 

add_triangle <- function(phy, group,phylo_plot){ 

    group_node_labels <- phy$tip.label[grepl(group, phy$tip.label)] 
    group_mrca <- getMRCA(phy,group_node_labels) 
    group_nodes <- descendants(as(tree.test,"phylo4"), group_mrca, type="tips") 
    names(group_nodes) <- NULL 

    x<-phylo_plot$xx 
    y<-phylo_plot$yy 


    x1 <- max(x[group_nodes]) 
    x2 <-max(x[group_nodes]) 
    x3 <- x[group_mrca] 

    y1 <- min(y[group_nodes]) 
    y2 <- max(y[group_nodes]) 
    y3 <- y[group_mrca] 

    xcoords <- c(x1,x2,x3) 
    ycoords <- c(y1,y2,y3) 

    polygon(xcoords, ycoords) 

    return(c(x2,y3)) 

} 



#main 

    cat("((A_1:0.05,E_2:0.03,A_3:0.2,A_4:0.1,A_5:0.1,A_6:0.1,A_7:0.35,A_8:0.4,A_9:01,A_10:0.2):0.9,((((B_1:0.05,B_2:0.05):0.5,B_3:0.02,B_4:0.04):0.6,(C_1:0.6,C_2:0.08):0.7):0.5,(D_1:0.3,D_2:0.4,D_3:0.5,D_4:0.7,D_5:0.4):1.2):0.5);", file = "ex.tre", sep = "\n") 
tree.test <- read.tree("ex.tre") 


# Step 1: Find the best MRCA that matches to the keywords or search patten 

groups <- c("A", "B|C", "D") 
group_labels <- groups 

group_edges <- numeric(0) 
edge.width <- rep(1, nrow(tree.test$edge)) 
count <- 1 


for(group in groups){ 

    best_mrca <- find_best_mrca(tree.test, group, 0.90) 

    group_leaves <- tips(tree.test, best_mrca) 

    groups[count] <- paste(group_leaves, collapse="|") 
    group_edges <- c(group_edges,best_mrca) 

    #Step2: Remove the edges of the branches that will be collapsed, so they become invisible 
    edge.width[tree.test$edge[,1] %in% c(group_edges[count],descendants(as(tree.test,"phylo4"), group_edges[count], type="all")) ] <- 0 
    count = count +1 

} 


#Step 3: plot the tree hiding the branches that will be collapsed/grouped 

last_plot.phylo <- plot(tree.test, show.tip.label = F, edge.width = edge.width) 

#And save a copy of the plot so we can extract the xy coordinates of the nodes 
#To get the x & y coordinates of a plotted tree created using plot.phylo 
#or plotTree, we can steal from inside tiplabels: 
last_phylo_plot<-get("last_plot.phylo",envir=.PlotPhyloEnv) 

#Step 4: Add triangles and labels to the collapsed nodes 
for(i in 1:length(groups)){ 

    text_coords <- add_triangle(tree.test, groups[i],last_phylo_plot) 

    text(text_coords[1],text_coords[2],labels=group_labels[i], pos=4) 

} 

Answer 3

Ho anche stato alla ricerca di questo tipo di strumento per secoli, non tanto per crollare gruppi categoriche, ma per collassare i nodi interni sulla base di un valore di supporto numerico.

La funzione di2multi nel pacchetto ape può comprimere i nodi in polytomies, ma al momento può farlo solo in base alla soglia della lunghezza del ramo. Ecco un adattamento approssimativo che consente invece il collasso di una soglia del valore di supporto nodo (soglia predefinita = 0,5).

Utilizzare a proprio rischio, ma funziona per me sul mio albero bayesiano radicato.

di2multi4node <- function (phy, tol = 0.5) 
    # Adapted di2multi function from the ape package to plot polytomies 
    # based on numeric node support values 
    # (di2multi does this based on edge lengths) 
    # Needs adjustment for unrooted trees as currently skips the first edge 
{ 
    if (is.null(phy$edge.length)) 
    stop("the tree has no branch length") 
    if (is.na(as.numeric(phy$node.label[2]))) 
    stop("node labels can't be converted to numeric values") 
    if (is.null(phy$node.label)) 
    stop("the tree has no node labels") 
    ind <- which(phy$edge[, 2] > length(phy$tip.label))[as.numeric(phy$node.label[2:length(phy$node.label)]) < tol] 
    n <- length(ind) 
    if (!n) 
    return(phy) 
    foo <- function(ancestor, des2del) { 
    wh <- which(phy$edge[, 1] == des2del) 
    for (k in wh) { 
     if (phy$edge[k, 2] %in% node2del) 
     foo(ancestor, phy$edge[k, 2]) 
     else phy$edge[k, 1] <<- ancestor 
    } 
    } 
    node2del <- phy$edge[ind, 2] 
    anc <- phy$edge[ind, 1] 
    for (i in 1:n) { 
    if (anc[i] %in% node2del) 
     next 
    foo(anc[i], node2del[i]) 
    } 
    phy$edge <- phy$edge[-ind, ] 
    phy$edge.length <- phy$edge.length[-ind] 
    phy$Nnode <- phy$Nnode - n 
    sel <- phy$edge > min(node2del) 
    for (i in which(sel)) phy$edge[i] <- phy$edge[i] - sum(node2del < 
                  phy$edge[i]) 
    if (!is.null(phy$node.label)) 
    phy$node.label <- phy$node.label[-(node2del - length(phy$tip.label))] 
    phy 
} 

Answer 4

Questa è la mia risposta basata su phytools::phylo.toBackbone funzione, vedere http://blog.phytools.org/2013/09/even-more-on-plotting-subtrees-as.html e http://blog.phytools.org/2013/10/finding-edge-lengths-of-all-terminal.html. Innanzitutto, carica la funzione alla fine del codice.

library(ape) 
library(phytools) #phylo.toBackbone 
library(phangorn) 

cat("((A_1:0.05,E_2:0.03,A_3:0.2,A_4:0.1,A_5:0.1,A_6:0.1,A_7:0.35,A_8:0.4,A_9:01,A_10:0.2):0.9,((((B_1:0.05,B_2:0.05):0.5,B_3:0.02,B_4:0.04):0.6,(C_1:0.6,C_2:0.08):0.7):0.5,(D_1:0.3,D_2:0.4,D_3:0.5,D_4:0.7,D_5:0.4):1.2):0.5);", file = "ex.tre", sep = "\n") 

phy <- read.tree("ex.tre") 
groups <- c("A", "B|C", "D") 

backboneoftree<-makebackbone(groups,phy) 
# tip.label clade.label N  depth 
# 1  A_1   A 10 0.2481818 
# 2  B_1   B|C 6 0.9400000 
# 3  D_1   D 5 0.4600000 

par(mfrow=c(2,2), mar=c(0,2,2,0)) 
plot(phy, main="Complete tree") 
plot(backboneoftree) 

makebackbone<-function(groupings,phy){ 
    listofspecies<-phy$tip.label 
    listtopreserve<-list() 
    lengthofclades<-list() 
    meandistnode<-list() 
    newedgelengths<-list() 
    for (i in 1:length(groupings)){ 
    groupings<-groups 
    bestmrca<-getMRCA(phy,grep(groupings[i], phy$tip.label)) 
    mrcatips<-phy$tip.label[unlist(phangorn::Descendants(phy,bestmrca, type="tips"))] 
    listtopreserve[i]<- mrcatips[1] 
    meandistnode[i]<- mean(dist.nodes(phy)[unlist(lapply(mrcatips, 
    function(x) grep(x, phy$tip.label))),bestmrca]) 
    lengthofclades[i]<-length(mrcatips) 
    provtree<-drop.tip(phy,mrcatips, trim.internal=F, subtree = T) 
    n3<-length(provtree$tip.label) 
    newedgelengths[i]<-setNames(provtree$edge.length[sapply(1:n3,function(x,y) 
     which(y==x),y=provtree$edge[,2])], 
     provtree$tip.label)[provtree$tip.label[grep("tips",provtree$tip.label)] ] 
    } 
    newtree<-drop.tip(phy,setdiff(listofspecies,unlist(listtopreserve)), 
        trim.internal = T) 
    n<-length(newtree$tip.label) 
    newtree$edge.length[sapply(1:n,function(x,y) 
    which(y==x),y=newtree$edge[,2])]<-unlist(newedgelengths)+unlist(meandistnode) 
    trans<-data.frame(tip.label=newtree$tip.label,clade.label=groupings, 
        N=unlist(lengthofclades), depth=unlist(meandistnode)) 
    rownames(trans)<-NULL 
    print(trans) 
    backboneoftree<-phytools::phylo.toBackbone(newtree,trans) 
    return(backboneoftree) 
} 

EDIT: Non ho provato questo, ma potrebbe essere un'altra risposta: "Script e la funzione di trasformare i rami di un albero di punta, vale a dire lo spessore o triangoli, con il larghezza sia correlazione con determinati parametri (ad esempio, il numero di specie del clade) (tip.branches.R)" http://www.sysbot.biologie.uni-muenchen.de/en/people/cusimano/use_r.html http://www.sysbot.biologie.uni-muenchen.de/en/people/cusimano/tip.branches.R