-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathaot.R
129 lines (122 loc) · 4.38 KB
/
aot.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
#' aot
#'
#' AOT conducts univariate and bivariate tests of phylogenetic signal and
#' trait correlations, respectively, and node-level analyses of trait
#' means and diversification.
#'
#' @export
#' @param traits (data.frame/character) trait data.frame or path to
#' traits file. required. See Details.
#' @template phylo
#' @param randomizations (numeric) number of randomizations. Default: 999
#' @param trait_contrasts (numeric) Specify which trait should be used as 'x'
#' variable for contrasts. Default: 1
#' @param ebl_unstconst (logical) Use equal branch lengths and unstandardized
#' contrasts. Default: `FALSE`
#' @return a list of data.frames
#' @details See [phylocomr-inputs] for expected input formats
#' @section Taxon name case:
#' In the `traits` table, if you're passing in a file, the names in the
#' first column must be all lowercase; if not, we'll lowercase them for you.
#' If you pass in a data.frame, we'll
#' lowercase them for your. All phylo tip/node labels are also lowercased
#' to avoid any casing problems
#' @examples \dontrun{
#' traits_file <- system.file("examples/traits_aot", package = "phylocomr")
#' phylo_file <- system.file("examples/phylo_aot", package = "phylocomr")
#'
#' # from data.frame
#' traitsdf_file <- system.file("examples/traits_aot_df",
#' package = "phylocomr")
#' traits <- read.table(text = readLines(traitsdf_file), header = TRUE,
#' stringsAsFactors = FALSE)
#' phylo_str <- readLines(phylo_file)
#' (res <- ph_aot(traits, phylo = phylo_str))
#'
#' # from files
#' traits_str <- paste0(readLines(traits_file), collapse = "\n")
#' traits_file2 <- tempfile()
#' cat(traits_str, file = traits_file2, sep = '\n')
#' phylo_file2 <- tempfile()
#' cat(phylo_str, file = phylo_file2, sep = '\n')
#' (res <- ph_aot(traits_file2, phylo_file2))
#' }
ph_aot <- function(traits, phylo, randomizations = 999, trait_contrasts = 1,
ebl_unstconst = FALSE) {
assert(traits, c("data.frame", "character"))
assert(phylo, c("phylo", "character"))
assert(randomizations, c('integer', 'numeric'))
assert(trait_contrasts, c('integer', 'numeric'))
assert(ebl_unstconst, 'logical')
stopifnot(class(traits) %in% c('data.frame', 'character'))
tfile <- tempfile("trait_")
if (inherits(traits, "data.frame")) {
if (colnames(traits)[1] != "name") {
stop("first column name in `traits` must be `name`", call. = FALSE)
}
if (inherits(traits[,1], c("character", "factor")))
traits[,1] <- tolower(traits[,1])
utils::write.table(prep_traits(traits), file = tfile,
quote = FALSE, row.names = FALSE)
} else {
stopifnot(file.exists(traits))
zz <- tolower(readLines(traits))
cat(zz, file = tfile, sep = "\n")
}
phylo <- phylo_check(phylo)
cdir <- getwd()
bdir <- dirname(tfile)
stopifnot(bdir == dirname(phylo))
setwd(bdir)
on.exit(setwd(cdir))
out <- suppressWarnings(
phylocom(c(
"aotf",
"-t", basename(tfile),
"-f", basename(phylo),
"-r", randomizations,
"-x", trait_contrasts,
if (ebl_unstconst) "-e "
), intern = TRUE)
)
phylocom_error(out)
out <- strsplit(out, split = "\n")[[1]]
list(
trait_conservatism = {
astbl(utils::read.delim(text = out[
seq(from = grep("Trait conservatism by node", out) + 1,
to = grep("Output of independent contrast", out) - 2)
], header = TRUE, stringsAsFactors = FALSE, sep = "\t"))
},
independent_contrasts = {
astbl(utils::read.delim(text = out[
seq(from = grep("Output of independent contrast", out) + 1,
to = grep("Phylogenetic signal", out) - 2)
], header = TRUE, stringsAsFactors = FALSE, sep = "\t"))
},
phylogenetic_signal = {
astbl(utils::read.delim(text = out[
seq(from = grep("Phylogenetic signal", out) + 1,
to = grep("Independent contrast correlations", out) - 2)
], header = TRUE, stringsAsFactors = FALSE, sep = "\t"))
},
ind_contrast_corr = {
astbl(utils::read.delim(text = out[
seq(from = grep("Independent contrast correlations", out) + 1,
to = length(out))
], header = TRUE, stringsAsFactors = FALSE, sep = "\t"))
}
)
}
prep_traits <- function(x) {
hd <- names(x)
x <- rbind(hd, x)
types <- unname(apply(x[,-1], 2, function(z) {
if (length(unique(z)) == 2) {
0
} else {
3
}
}))
stats::setNames(x, c("type", types))
}