-
Notifications
You must be signed in to change notification settings - Fork 1
/
tree_time.R
158 lines (138 loc) · 5.12 KB
/
tree_time.R
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
156
157
158
#' Plot Phylogenies and Time
#'
#' These two functions plot a tree built from the sequences available
#' in the \linkS4class{obkData} object together with the time of
#' sampling of each sequence (also taken from the same object).
#'
#' @param x an \linkS4class{obkData} object.
#' @return nothing, the results are plotted on the current graphical device.
#' @author Emmanuel Paradis
#'
#' @details The first function \code{tree_time_1} plots the NJ tree
#' growing upwards, without the branch length information, and the
#' time axis is represented on top with from left to right. The second
#' function \code{tree_time_2} plots the NJ tree in the usual way,
#' from left to right, with the branch length information, and the
#' time axis is represented on bottom also from left to right.
#'
#' Coloured lines show the link between the tips of the tree (i.e.,
#' unique sequences) and the time of sampling of the individuals. The
#' colour scale goes from blue (oldest dates) to red (youngest dates).
#'
#' The functions scan the data for the different genes and plots the
#' output successively (the user is asked to type enter).
#'
#' @note The neighbor-joigning (NJ) method reconstructs unrooted trees
#' whereas the present functions represent them as rooted.
#'
#' @examples
#' \dontrun{
#' require(OutbreakTools)
#' data(ToyOutbreak)
#' tree_time_1(ToyOutbreak)
#' tree_time_2(ToyOutbreak)
#' }
tree_time_1 <- function(x)
{
require(pegas)
DNA <- x@dna@dna
LOCI <- names(DNA)
META <- x@dna@meta
cat("okbData object has", length(LOCI), "genetic data set(s)\n")
for (i in seq_along(LOCI)) {
cat("Processing data from", LOCI[i], "\n")
dna <- DNA[[i]]
h <- haplotype(dna)
nh <- nrow(h)
attr(h, "index") <- lapply(attr(h, "index"), function(x) rownames(dna)[x])
names(attr(h, "index")) <- rownames(h) <- paste0("uniqseqID", seq_len(nh))
phy <- nj(dist.dna(h, "N", pairwise.deletion = TRUE))
n <- Ntip(phy)
sel <- META$locus == LOCI[i]
submeta <- META[sel, ]
dates <- as.integer(submeta$date)
range.dates <- range(dates)
footrans <- function(x) (n - 1) * (x - range.dates[1])/diff(range.dates) + 1
transfdates <- footrans(dates)
x1 <- transfdates
names(x1) <- rownames(submeta)
x2 <- x1
cat("Type enter to plot\n")
readLines(n = 1)
plot(phy, "p", FALSE, show.tip.label = FALSE, direction = "u",
font = 1, label.offset = 1, node.depth = 2, y.lim = 40)
psr <- par("usr")
nms <- rownames(submeta)
index <- attr(h, "index")
lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
for (j in 1:n) {
s <- index[[phy$tip.label[j]]]
x2[s] <- lastPP$xx[j]
}
y2 <- psr[4] - 0.1
y1 <- max(lastPP$yy) + 0.1
segments(x2, y1, x1, y2, lwd = 1, col = rgb(x1/n, 0, 1 - x1/n, alpha = .5))
at <- pretty(dates)
axis(3, at = footrans(at), labels = as.Date(at, origin = "1970-01-01"))
}
}
#old code:
#
#nms <- rownames(submeta)
#index <- attr(h, "index")
#lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
#for (i in 1:n) {
# s <- index[[phy$tip.label[i]]]
# x2[s] <- lastPP$xx[i]
#}
#
#y1 <- psr[3]
#y2 <- 5
#segments(x1, y1, x2, y2, lwd = .5)
#at <- pretty(dates)
#axis(1, at = footrans(at), labels = as.Date(at, origin = "1970-01-01"))
# o <- order(dates)
# phy2 <- rotateConstr(phy, rownames(submeta)[o])
tree_time_2 <- function(x)
{
require(pegas)
DNA <- x@dna@dna
LOCI <- names(DNA)
META <- x@dna@meta
cat("okbData object has", length(LOCI), "genetic data set(s)\n")
for (i in seq_along(LOCI)) {
cat("Processing data from", LOCI[i], "...\n")
dna <- DNA[[i]]
h <- haplotype(dna)
nh <- nrow(h)
attr(h, "index") <- lapply(attr(h, "index"), function(x) rownames(dna)[x])
names(attr(h, "index")) <- rownames(h) <- paste0("uniqseqID", seq_len(nh))
phy <- nj(dist.dna(h, "N", pairwise.deletion = TRUE))
n <- Ntip(phy)
sel <- META$locus == LOCI[i]
submeta <- META[sel, ]
dates <- as.integer(submeta$date)
range.dates <- range(dates)
cat("Type enter to plot\n")
readLines(n = 1)
plot(phy, "p", show.tip.label = FALSE, font = 1, label.offset = 1, node.depth = 2, y.lim = c(-100, n), edge.width = 2)
psr <- par("usr")
footrans2 <- function(x) psr[2] * (x - range.dates[1])/diff(range.dates)
x1 <- footrans2(dates)
y1 <- psr[3]
lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
index <- attr(h, "index")
x2 <- y2 <- numeric(length(dates))
names(x2) <- names(y2) <- rownames(submeta)
for (j in 1:n) {
s <- index[[phy$tip.label[j]]]
x2[s] <- lastPP$xx[j]
y2[s] <- lastPP$yy[j]
}
segments(x1, y1, x2, y2, lwd = 1, col = rgb(x1/max(x1), 0, 1 - x1/max(x1), alpha = .5))
at <- pretty(dates)
axis(1, at = footrans2(at), labels = as.Date(at, origin = "1970-01-01"))
add.scale.bar(psr[2] - 3, psr[4] - 5, 1)
mtext(LOCI[i])
}
}