/spacodiR/R/phylogeneticUtilities.R

  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