/spacodiR/R/phylogeneticUtilities.R
R | 148 lines | 111 code | 27 blank | 10 comment | 29 complexity | 68f20a7ea80c162859cf6ef68bafe361 MD5 | raw file
1#general statistical function for comparing two vectors on non-independent samples 2#author: JM EASTMAN 2010; updated 03.06.2011 to remove NA 3 4resamp.test <- 5function(obs=obs, exp=exp, mu=0, iter=10000, two.tailed=FALSE, na.rm=TRUE){ 6 if(na.rm) { 7 oe=lapply(list(obs,exp), function(x) return(x[!is.na(x)])) 8 obs=oe[[1]] 9 exp=oe[[2]] 10 } 11 O=as.numeric(obs)[sample(1:length(obs), iter, replace=TRUE)] 12 E=as.numeric(exp)[sample(1:length(exp), iter, replace=TRUE)] 13 14 result=c(O-(E+mu)) 15 p=round(sum(result>=0)/iter, digits=nchar(iter)) 16 q=1-p 17 18 if(two.tailed) { 19 res=list(diffs=result, 2*(c(p,q)[which(c(p,q)==min(c(p,q)))])) 20 } else { 21 res=list(diffs=result, greater=p, lesser=q) 22 } 23 return(res) 24} 25 26#general phylogenetic utility for returning the first (usually, unless a polytomy exists) two descendants the supplied node 27#author: JM EASTMAN 2010 28 29get.desc.of.node <- 30function(node, phy) { 31 return(phy$edge[which(phy$edge[,1]==node),2]) 32} 33 34 35#general phylogenetic utility for returning all descendants (listed as given in phy$edge[,2]) of a node (excluding the supplied node) 36#author: JM EASTMAN 2010, based on GEIGER:::node.leaves() 37 38get.descendants.of.node <- 39function (node, phy, tips=FALSE) 40{ 41 node <- as.numeric(node) 42 n <- Ntip(phy) 43 if(node <= n) return(NULL) 44 l <- c() 45 d <- get.desc.of.node(node, phy) 46 for (j in d) { 47 l <- c(l, j) 48 if (j > n) { 49 l <- c(l, get.descendants.of.node(j, phy)) 50 } 51 } 52 if(tips) return(l[l<=n]) else return(l) 53} 54 55 56#determines if a value falls within a range [min,max] 57#author: JM Eastman 2010 58 59withinrange <- 60function(x, min, max) { 61 a=sign(x-min) 62 b=sign(x-max) 63 if(abs(a+b)==2) return(FALSE) else return(TRUE) 64} 65 66 67 68 69###### 70phy.deresolve=function(phy, time.range=c(0,0), relative=TRUE) 71{ 72 if (is.null(phy$edge.length)) 73 stop("The tree does not appear to have branch lengths") 74 if (class(phy)!="phylo") 75 stop("The tree does not appear to be a valid phylo object") 76 if(length(time.range)>2) 77 stop("Cannot interpret the range of time with more than two elements") 78 if(length(time.range)==1) 79 time.range=c(0,time.range) 80 time.range=time.range[order(time.range)] 81 bb <- branching.times(phy) 82 if(relative) bb=bb/max(bb) 83 inr=sapply(bb, function(x) withinrange(x, time.range[1], time.range[2])) 84 ind <- as.numeric(names(bb[inr])) 85 if(any(ind==Ntip(phy)+1)) ind=ind[-which(ind==Ntip(phy)+1)] 86 n <- length(ind) 87 if (!n) { 88 return(phy) 89 } else { 90 ind.tmp = match(ind, phy$edge[,2]) 91 ind = ind.tmp[!is.na(ind.tmp)] 92 } 93 orig.edge=phy$edge 94 orig.phy=phy 95 ntips=Ntip(phy) 96 reedge <- function(ancestor, des.to.drop) { 97 wh <- which(phy$edge[, 1] == des.to.drop) 98 dd <- which(orig.edge[, 2] == des.to.drop) 99 dropped.branch <- phy$edge.length[dd] 100 d.d <- c(get.desc.of.node(des.to.drop, orig.phy)) 101 if(length(d.d)) phy$edge.length[match(d.d, orig.edge[,2])]<<-phy$edge.length[match(d.d, orig.edge[,2])]+dropped.branch 102 103 for (k in wh) { 104 if (phy$edge[k, 2] %in% node.to.drop) { 105 reedge(ancestor, phy$edge[k, 2]) 106 } else { 107 phy$edge[k, 1] <<- ancestor 108 109 } 110 } 111 } 112 113 node.to.drop <- phy$edge[ind, 2] 114 anc <- phy$edge[ind, 1] 115 for (i in 1:n) { 116 if (anc[i] %in% node.to.drop) 117 next 118 reedge(anc[i], node.to.drop[i]) 119 } 120 phy$edge <- phy$edge[-ind, ] 121 phy$edge.length <- phy$edge.length[-ind] 122 phy$Nnode <- phy$Nnode - n 123 sel <- phy$edge > min(node.to.drop) 124 for (i in which(sel)) phy$edge[i] <- phy$edge[i] - sum(node.to.drop < phy$edge[i]) 125 if (!is.null(phy$node.label)) 126 phy$node.label <- phy$node.label[-(node.to.drop - length(phy$tip.label))] 127 phy 128} 129 130###### 131phy.nodetimes <- 132function(phy, time.range=c(0,0), proportion=TRUE) { 133 N=Ntip(phy) 134 phy$node.label=NULL 135 xx = c(rep(0,N), branching.times(phy)) 136 names(xx)[1:N]=1:N 137 if(proportion) xx=xx/max(branching.times(phy)) 138 tt=sapply(xx,function(x) withinrange(x,min(time.range),max(time.range))) 139 return(xx[tt]) 140} 141 142 143 144 145 146 147 148