-
Notifications
You must be signed in to change notification settings - Fork 56
/
consensus.edges.R
31 lines (30 loc) · 1.15 KB
/
consensus.edges.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
## compute consensus edge lengths from a set of trees given (or not) a consensus topology
## written by Liam J. Revell 2016
consensus.edges<-function(trees,method=c("mean.edge","least.squares"),...){
if(hasArg(consensus.tree)) tree<-list(...)$consensus.tree
else tree<-consensus(trees,p=0.5)
if(hasArg(if.absent)) if.absent<-list(...)$if.absent
else if.absent<-"zero"
N<-length(trees)
if(method[1]=="mean.edge"){
M<-lapply(trees,function(x,y) rbind(matchLabels(y,x),matchNodes(y,x)),y=tree)
nodes<-M[[1]][,1]
edge.length<-vector(mode="numeric",length=length(nodes))
for(i in 2:length(nodes)){
ii<-which(tree$edge[,2]==nodes[i])
n.absent<-0
for(j in 1:N){
edge.length[ii]<-edge.length[ii] +
if(!is.na(M[[j]][i,2])) trees[[j]]$edge.length[which(trees[[j]]$edge[,2]==M[[j]][i,2])]/N
else 0
if(is.na(M[[j]][i,2])) n.absent<-n.absent+1
}
if(if.absent=="ignore") edge.length[ii]<-edge.length[ii]*N/(N-n.absent)
}
tree$edge.length<-edge.length
} else if(method[1]=="least.squares"){
D<-Reduce('+',lapply(trees,function(x,t) cophenetic(x)[t,t],t=tree$tip.label))/N
tree<-nnls.tree(D,tree=tree,rooted=is.rooted(tree))
}
tree
}