# function reads SIMMAPv1.0 tree into memory as modified {ape} "phylo" object # written by Liam J. Revell 2011 read.simmap.alpha<-function(file="",text){ # check to see if reading from file if(file!=""){ text<-scan(file,sep="\n",what="character") # read from file } tree<-unlist(strsplit(text, NULL)) # split string 'text' Ntip<-1; Nnode<-0; i<-1 while(tree[i]!=";"){ while(tree[i]=="{") while(tree[i]!="}") i<-i+1 # skip maps if(tree[i]==",") Ntip<-Ntip+1 # count tips if(tree[i]=="(") Nnode<-Nnode+1 # count nodes i<-i+1 } tip.label<-vector(mode="character") # create tip label vector edge<-matrix(data=0,Nnode+Ntip-1,2) # create edge matrix edge.length<-rep(0,Nnode+Ntip-1) # create edge.length vector maps<-list() currnode<-Ntip+1; nodecount<-currnode # start with root node i<-1; j<-1; k<-1 while(tree[i]!=";"){ if(tree[i]=="("){ edge[j,1]<-currnode; i<-i+1 # is the next element a label? if(is.na(match(tree[i],c("(",")",",",":",";")))){ l<-1; temp<-vector() while(is.na(match(tree[i],c(",",":",")")))){ temp[l]<-tree[i]; l<-l+1; i<-i+1 } tip.label[k]<-paste(temp,collapse=""); edge[j,2]<-k; k<-k+1 # create tip label # is there a branch length? if(tree[i]==":"){ i<-i+1; l<-1 if(tree[i]=="{"){ maps[[j]]<-vector() i<-i+1 while(tree[i]!="}"){ temp<-vector(); m<-1 while(tree[i]!=","){ temp[m]<-tree[i] i<-i+1; m<-m+1 } state<-paste(temp,collapse="") i<-i+1; temp<-vector(); m<-1 while(tree[i]!=":"&&tree[i]!="}"){ temp[m]<-tree[i] i<-i+1; m<-m+1 } length<-as.numeric(paste(temp,collapse="")) maps[[j]][l]<-length; names(maps[[j]])[l]<-as.character(state) l<-l+1 if(tree[i]==":") i<-i+1 } } edge.length[j]<-sum(maps[[j]]) # create branch length i<-i+1 } } else if(tree[i]=="("){ nodecount<-nodecount+1 # creating a new internal node currnode<-nodecount; edge[j,2]<-currnode; # move to new internal node } j<-j+1; } else if(tree[i]==")"){ i<-i+1 # is there a branch length? if(tree[i]==":"){ i<-i+1; l<-1 if(tree[i]=="{"){ maps[[match(currnode,edge[,2])]]<-vector() i<-i+1 while(tree[i]!="}"){ temp<-vector(); m<-1 while(tree[i]!=","){ temp[m]<-tree[i] i<-i+1; m<-m+1 } state<-paste(temp,collapse="") i<-i+1; temp<-vector(); m<-1 while(tree[i]!=":"&&tree[i]!="}"){ temp[m]<-tree[i] i<-i+1; m<-m+1 } length<-as.numeric(paste(temp,collapse="")) maps[[match(currnode,edge[,2])]][l]<-length; names(maps[[match(currnode,edge[,2])]])[l]<-as.character(state) l<-l+1 if(tree[i]==":") i<-i+1 } } edge.length[match(currnode,edge[,2])]<-sum(maps[[match(currnode,edge[,2])]]) # create branch length i<-i+1 } currnode<-edge[match(currnode,edge[,2]),1] # move down the tree } else if(tree[i]==","){ edge[j,1]<-currnode; i<-i+1 # is the next element a label? if(is.na(match(tree[i],c("(",")",",",":",";")))){ l<-1; temp<-vector() while(is.na(match(tree[i],c(",",":",")")))){ temp[l]<-tree[i]; l<-l+1; i<-i+1 } tip.label[k]<-paste(temp,collapse=""); edge[j,2]<-k; k<-k+1 # create tip label # is there a branch length? if(tree[i]==":"){ i<-i+1; l<-1 if(tree[i]=="{"){ maps[[j]]<-vector() i<-i+1 while(tree[i]!="}"){ temp<-vector(); m<-1 while(tree[i]!=","){ temp[m]<-tree[i] i<-i+1; m<-m+1 } state<-paste(temp,collapse="") i<-i+1; temp<-vector(); m<-1 while(tree[i]!=":"&&tree[i]!="}"){ temp[m]<-tree[i] i<-i+1; m<-m+1 } length<-as.numeric(paste(temp,collapse="")) maps[[j]][l]<-length; names(maps[[j]])[l]<-as.character(state) l<-l+1 if(tree[i]==":") i<-i+1 } } edge.length[j]<-sum(maps[[j]]) # create branch length i<-i+1 } } else if(tree[i]=="("){ nodecount<-nodecount+1 # creating a new internal node currnode<-nodecount; edge[j,2]<-currnode; # move to internal node } j<-j+1 } } # construct the matrix mapped.edge (for backward compatibility allstates<-vector() for(i in 1:nrow(edge)) allstates<-c(allstates,names(maps[[i]])) allstates<-unique(allstates) mapped.edge<-matrix(data=0,length(edge.length),length(allstates),dimnames=list(edge=apply(edge,1,function(x) paste(x,collapse=",")),state=allstates)) for(i in 1:length(maps)) for(j in 1:length(maps[[i]])) mapped.edge[i,names(maps[[i]])[j]]<-mapped.edge[i,names(maps[[i]])[j]]+maps[[i]][j] # assemble into modified "phylo" object phy<-list(edge=edge,Nnode=as.integer(Nnode),tip.label=tip.label,edge.length=edge.length,maps=maps,mapped.edge=mapped.edge) class(phy)<-"phylo" # assign class return(phy) # return object }