1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
|
#'
#' Phylogenetic tree exploration
#'
#' Compares phylogenetic trees using a choice of metrics / measures, and maps their pairwise distances into a small number of dimensions for easy visualisation and identification of clusters.
#'
#' @param x an object of the class multiPhylo
#' @param method the method for summarising the tree as a vector.
#' Choose from:
#' \code{treeVec} (default) the Kendall Colijn metric vector
#' \code{BHV} the Billera, Holmes Vogtmann metric using \code{dist.multiPhylo} from package \code{distory}
#' \code{KF} the Kuhner Felsenstein metric (branch score distance) using \code{KF.dist} from package \code{phangorn} (Note: this considers the trees as unrooted)
#' \code{RF} the Robinson Foulds metric using \code{RF.dist} from package \code{phangorn} (Note: this considers the trees as unrooted and issues a corresponding warning)
#' \code{wRF} the weighted Robinson Foulds metric using \code{wRF.dist} from package \code{phangorn} (Note: this considers the trees as unrooted and issues a corresponding warning)
#' \code{nNodes} the Steel & Penny tip-tip path difference metric, (topological, ignoring branch lengths), using \code{path.dist} from package \code{phangorn} (Note: this considers the trees as unrooted)
#' \code{patristic} the Steel & Penny tip-tip path difference metric, using branch lengths, calling \code{path.dist} from package \code{phangorn} (Note: this considers the trees as unrooted)
#' others inherited from \code{distTips} in \code{adephylo}:
#' \itemize{
#' \item \code{Abouheif}: performs Abouheif's test. See Pavoine et al. (2008) and \code{adephylo}.
#' \item \code{sumDD}: sum of direct descendants of all nodes on the path, related to Abouheif's test. See \code{adephylo}.
#' }
#' @param nf the number of principal components to retain
#' @param lambda a number in [0,1] which specifies the extent to which topology (default, with lambda=0) or branch lengths (lambda=1) are emphasised in the Kendall Colijn metric.
#' @param return.tree.vectors option to also return the tree vectors. Note that this can use a lot of memory so defaults to \code{FALSE}.
#' @param processors value (default 1) to be passed to mcmapply specifying the number of cores to use. Must be 1 on Windows (see \code{mcmapply} for more details).
#' @param ... further arguments to be passed to \code{method}.
#'
#' @author Thibaut Jombart \email{thibautjombart@@gmail.com}
#' @author Michelle Kendall \email{michelle.louise.kendall@@gmail.com}
#'
#' @import ape
#' @importFrom ade4 dudi.pco cailliez
#' @importFrom adephylo distTips
#' @importFrom distory dist.multiPhylo
#' @importFrom fields rdist
#' @importFrom phangorn KF.dist
#' @importFrom phangorn path.dist
#' @importFrom phangorn RF.dist
#' @importFrom phangorn wRF.dist
#' @importFrom parallel mcmapply
#'
#' @examples
#'
#' ## generate list of trees
#' x <- rmtree(10, 20)
#' names(x) <- paste("tree", 1:10, sep = "")
#'
#' ## use treespace
#' res <- treespace(x, nf=3)
#' table.paint(as.matrix(res$D))
#' scatter(res$pco)
#'
#' data(woodmiceTrees)
#' woodmiceDists <- treespace(woodmiceTrees,nf=3)
#' plot(woodmiceDists$pco$li[,1],woodmiceDists$pco$li[,2])
#' woodmicedf <- woodmiceDists$pco$li
#' if(require(ggplot2)){
#' woodmiceplot <- ggplot(woodmicedf, aes(x=A1, y=A2)) # create plot
#' woodmiceplot + geom_density2d(colour="gray80") + # contour lines
#' geom_point(size=6, shape=1, colour="gray50") + # grey edges
#' geom_point(size=6, alpha=0.2, colour="navy") + # transparent blue points
#' xlab("") + ylab("") + theme_bw(base_family="") # remove axis labels and grey background
#' }
#'
#' \dontrun{
#' if(require(rgl)){
#' plot3d(woodmicedf[,1], woodmicedf[,2], woodmicedf[,3], type="s", size=1.5,
#' col="navy", alpha=0.5, xlab="", ylab="", zlab="")
#' }
#' }
#'
#'
#' @export
treespace <- function(x, method="treeVec", nf=NULL, lambda=0, return.tree.vectors=FALSE, processors=1, ...){
## CHECKS ##
if(!inherits(x, "multiPhylo")) stop("x should be a multiphylo object")
num_trees <- length(x) # number of trees
## fix potential bug with input of two trees
if(num_trees<3) {
stop("treespace expects at least three trees. The function treeDist is suitable for comparing two trees.")
}
# check for user supplying invalid options (these gave unhelpful error messages before)
dots <- list(...)
if(!is.null(dots$return.lambda.function)) stop("return.lambda.function is not compatible with treespace. Consider using multiDist instead.")
if(!is.null(dots$save.memory)) stop("save.memory is not compatible with treespace. Consider using multiDist instead.")
# make name labels well defined
if(is.null(names(x))) names(x) <- 1:num_trees
else if(length(unique(names(x)))!=num_trees){
warning("duplicates detected in tree labels - using generic names")
names(x) <- 1:num_trees
}
lab <- names(x)
# check all trees have same tip labels
for (i in 1:num_trees) {
if (!setequal(x[[i]]$tip.label,x[[1]]$tip.label)) {
stop(paste0("Tree ",lab[[i]]," has different tip labels from the first tree."))
}
}
## GET DISTANCES BETWEEN TREES, according to method ##
## get summary vectors then compute pairwise distances ##
if (method=="treeVec") {
df <- t(mcmapply(treeVec, x, lambda=lambda, MoreArgs=dots, mc.cores=processors))
## get pairwise Euclidean distances ##
D <- as.dist(rdist(df))
}
else if(method %in% c("Abouheif","sumDD")){
df <- t(mcmapply(adephylo::distTips, x, method=method, MoreArgs=dots, mc.cores=processors))
## get pairwise Euclidean distances ##
D <- as.dist(rdist(df))
}
else if(method=="patristic"){
D <- path.dist(x, use.weight=TRUE)
}
else if(method=="nNodes"){
D <- path.dist(x, use.weight=FALSE)
}
else if(method=="RF"){
D <- RF.dist(x)
## make the distance Euclidean
D <- ade4::cailliez(D, print=FALSE)
}
else if(method=="wRF"){
D <- wRF.dist(x)
## make the distance Euclidean
D <- ade4::cailliez(D, print=FALSE)
}
else if(method=="KF"){
D <- KF.dist(x)
}
else if(method=="BHV"){
D <- dist.multiPhylo(x)
## make the distance Euclidean
D <- ade4::cailliez(D, print=FALSE)
}
## restore labels
attr(D,"Labels") <- lab
## perform PCoA/MDS ##
pco <- dudi.pco(D, scannf=is.null(nf), nf=nf)
## BUILD RESULT AND RETURN ##
if (return.tree.vectors==TRUE) {
out <- list(D=D, pco=pco, vectors=df)
}
else {
out <- list(D=D, pco=pco)
}
return(out)
} # end treespace
|
