[Picante-commits] r212 - in pkg: R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Feb 9 22:42:43 CET 2010
Author: skembel
Date: 2010-02-09 22:42:42 +0100 (Tue, 09 Feb 2010)
New Revision: 212
Modified:
pkg/R/phylosor.R
pkg/man/phylosor.Rd
pkg/man/phylosor.rnd.Rd
Log:
Fix phylosor and phylosor.rnd to allow non-ultrametric trees
Modified: pkg/R/phylosor.R
===================================================================
--- pkg/R/phylosor.R 2010-02-08 22:57:56 UTC (rev 211)
+++ pkg/R/phylosor.R 2010-02-09 21:42:42 UTC (rev 212)
@@ -1,32 +1,44 @@
-phylosor <- function(samp,tree)
+phylosor <- function (samp, tree)
{
-
- #If phylo has no given branch lengths
- if(is.null(tree$edge.length)) {
- stop("Tree has no branch lengths, cannot compute pd")
+ if (is.null(tree$edge.length)) {
+ stop("Tree has no branch lengths, cannot compute pd")
}
- # Make sure that the species line up
- tree<-prune.sample(samp,tree)
- samp<-samp[,tree$tip.label]
-
- s=nrow(samp)
- phylodist=matrix(NA,s,s)
- rownames(phylodist)=rownames(samp)
- colnames(phylodist)=rownames(samp)
-
- for (l in 1:(s-1))
- {
- pdl=.pdshort(samp[l,],tree)
- for (k in (l+1):s)
- {
- pdk=.pdshort(samp[k,],tree)
- pdtot=.pdshort((samp[l,]+samp[k,]),tree)
- pdsharedlk=pdl+pdk-pdtot
- phylodist[k,l]=2*pdsharedlk/(pdl+pdk)
+ if (!is.rooted(tree)) {
+ stop("Rooted phylogeny required for phylosor calculation")
+ }
+
+ s <- nrow(samp)
+ phylodist <- matrix(NA, s, s)
+ rownames(phylodist) <- rownames(samp)
+ colnames(phylodist) <- rownames(samp)
+
+ samp_comb<-matrix(NA,s*(s-1)/2,ncol(samp))
+ colnames(samp_comb)<-colnames(samp)
+
+ i<-1
+ for (l in 1:(s - 1))
+ {
+ for (k in (l + 1):s)
+ {samp_comb[i,]<-samp[l, ] + samp[k, ]
+ i<-i+1}
+ }
+
+ pdsamp<-pd(samp, tree)
+ pdsamp_comb<-pd(samp_comb, tree)
+
+ i<-1
+ for (l in 1:(s - 1)) {
+ pdl <- pdsamp[l,"PD"]
+ for (k in (l + 1):s) {
+ pdk <- pdsamp[k,"PD"]
+ pdcomb <- pdsamp_comb[i,"PD"]
+ pdsharedlk <- pdl + pdk - pdcomb
+ phylodist[k, l] = 2 * pdsharedlk/(pdl + pdk)
+ i<-i+1
}
}
- return(as.dist(phylodist))
+ return(as.dist(phylodist))
}
phylosor.rnd <- function(samp, tree, cstSor=TRUE, null.model=c("taxa.labels", "frequency", "richness", "independentswap", "trialswap"), runs=999, iterations=1000)
@@ -97,54 +109,12 @@
{
sampr=samp
colnames(sampr)=sample(colnames(samp))
- pdl=.pdshort(sampr[l,],tree)
- pdk=.pdshort(sampr[k,],tree)
- pdtot=.pdshort((sampr[l,]+sampr[k,]),tree)
+ pdl=pd(sampr[l,,drop=FALSE],tree)$PD
+ pdk=pd(sampr[k,,drop=FALSE],tree)$PD
+ pdtot=pd((sampr[l,,drop=FALSE]+sampr[k,,drop=FALSE]),tree)$PD
pdsharedlk=pdl+pdk-pdtot
phylodist[k,l]=2*pdsharedlk/(pdl+pdk)
- }
- }
- return(as.dist(phylodist))
+ }
+ }
+ return(as.dist(phylodist))
}
-
-#############################################################################################
-
-.pdshort=function(comm,tree,keep.root=TRUE)
-{
-
- if (!is.rooted(tree) || !is.ultrametric(tree)) {
- stop("Rooted ultrametric tree required for phylosor calculation")
- }
-
- nbspecies <- length(comm)
- species <- names(comm)
-
- present <- species[comm>0] #species in sample
- treeabsent <- tree$tip.label[which(!(tree$tip.label %in% present))]
-
- if(length(present)==0)
- {
- #no species present
- PD<-0
- }
- else if(length(present)==1)
- {
- #one species present - PD = length from root to that tip
- PD <- node.age(tree)$ages[ which(tree$edge[,2] ==
- which(tree$tip.label==present))]
- }
- else if(length(treeabsent)==0)
- {
- #all species in tree present in community
- PD <- sum(tree$edge.length)
- }
- else
- {
- #subset of tree species present in community
- sub.tree<-drop.tip(tree,treeabsent)
- sub.tree.depth <- max(node.age(sub.tree)$ages)
- orig.tree.depth <- max(node.age(tree)$ages)
- PD<-sum(sub.tree$edge.length) + (orig.tree.depth-sub.tree.depth)
- }
- return(PD)
-}
Modified: pkg/man/phylosor.Rd
===================================================================
--- pkg/man/phylosor.Rd 2010-02-08 22:57:56 UTC (rev 211)
+++ pkg/man/phylosor.Rd 2010-02-09 21:42:42 UTC (rev 212)
@@ -9,14 +9,14 @@
\arguments{
\item{samp}{ Community data matrix }
- \item{tree}{ Object of class phylo}
+ \item{tree}{ Object of class phylo - a rooted phylogeny}
}
\value{A distance object of the PhyloSor index of similarity between communities, the fraction of PD (branch-length) shared between two samples}
\references{Bryant, J.B., Lamanna, C., Morlon, H., Kerkhoff, A.J., Enquist, B.J., Green, J.L. 2008. Microbes on mountainsides: Contrasting elevational patterns of bacterial and plant diversity. Proceedings of the National Academy of Sciences 105 Supplement 1: 11505-11511}
\author{ Helene Morlon <morlon.helene at gmail.com> and Steven Kembel <skembel at uoregon.edu> }
\seealso{\code{\link{phylosor.rnd}}, \code{\link{pd}}}
\note{
-The root of the supplied tree is included in calculations of PhyloSor. The supplied tree must be rooted and ultrametric. Single-species samples will be assigned a PD value equal to the distance from the root to the present.
+The root of the supplied tree is included in calculations of PhyloSor. The supplied tree must be rooted. Single-species samples will be assigned a PD value equal to the distance from the root to the present.
}
\section{Warning }{
The phylosor of all samples will include the branch length connecting taxa in those samples and the root of the supplied tree. The root of the supplied tree may not be spanned by any taxa in the sample. If you want the root of your tree to correspond to the most recent ancestor of the taxa actually present in your sample, you should prune the tree before running \code{phylosor}:
Modified: pkg/man/phylosor.rnd.Rd
===================================================================
--- pkg/man/phylosor.rnd.Rd 2010-02-08 22:57:56 UTC (rev 211)
+++ pkg/man/phylosor.rnd.Rd 2010-02-09 21:42:42 UTC (rev 212)
@@ -13,7 +13,7 @@
\arguments{
\item{ samp }{ Community data matrix }
- \item{ tree }{ Object of class phylo}
+ \item{ tree }{ Object of class phylo - a rooted phylogeny}
\item{ cstSor }{ TRUE if the Sorensen similarity should be kept constant across communities. FALSE otherwise }
\item{ null.model }{ Null model to use (see Details section) }
\item{ runs }{ Number of randomizations }
More information about the Picante-commits
mailing list