#### /spacodiR/R/phylogeneticUtilities.R

http://github.com/eastman/spacodiR
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
```