Hello Phylofolks,
I've been using ape::phymltest() for a while now on a project and I just
realised there is a potential glitch in the function. The reason is that
the order of the models in the 'mdl' object inside the function is
different from the order in .phymltest.model, what causes the cmd object to
mess up.
The following MWE should do the job of showing my point:
### begin MWE
library(ape)
.phymltest.model <-
c("JC69", "JC69+I", "JC69+G", "JC69+I+G",
"K80", "K80+I", "K80+G", "K80+I+G",
"F81", "F81+I", "F81+G", "F81+I+G",
"F84", "F84+I", "F84+G", "F84+I+G",
"HKY85", "HKY85+I", "HKY85+G", "HKY85+I+G",
"TN93", "TN93+I", "TN93+G", "TN93+I+G",
"GTR", "GTR+I", "GTR+G", "GTR+I+G")
phymltest.test <- function(seqfile, format = "interleaved", itree = NULL,
exclude = NULL, execname = NULL, append = TRUE)
{
N <- length(.phymltest.model)
format <- match.arg(format, c("interleaved", "sequential"))
fmt <- rep("", N)
if (format != "interleaved") fmt[] <- "-q"
boot <- rep("-b 0", N) # to avoid any testing
mdl <- paste("-m", rep(c("JC69", "K80", "F81","HKY85","F84", "TN93",
"GTR"), each = 4))
tstv <- rep("-t e", N) # ignored by PhyML with JC69 or F81
inv <- rep(c("", "-v e"), length.out = N)
## no need to use the -c option of PhyML (4 categories by default if '-a
e' is set):
alpha <- rep(rep(c("-c 1", "-a e"), each = 2), length.out = N)
tree <- rep("", N)
cmd <- paste(execname, "-i", seqfile, fmt, boot, mdl, tstv, inv, alpha,
tree, "--append ")
imod <- 1:N
if (!is.null(exclude)) imod <- imod[!.phymltest.model %in% exclude]
for (i in imod) print(cmd[i])
}
phymltest.test(seqfile = "seqname",
exclude = c("K80", "K80+I", "K80+G", "K80+I+G",
"F81", "F81+I", "F81+G", "F81+I+G",
"F84", "F84+I", "F84+G", "F84+I+G",
"TN93", "TN93+I", "TN93+G", "TN93+I+G"),
execname = "execname")
1] "execname -i seqname -b 0 -m JC69 -t e -c 1 --append "
[1] "execname -i seqname -b 0 -m JC69 -t e -v e -c 1 --append "
[1] "execname -i seqname -b 0 -m JC69 -t e -a e --append "
[1] "execname -i seqname -b 0 -m JC69 -t e -v e -a e --append "
[1] "execname -i seqname -b 0 -m F84 -t e -c 1 --append "
[1] "execname -i seqname -b 0 -m F84 -t e -v e -c 1 --append "
[1] "execname -i seqname -b 0 -m F84 -t e -a e --append "
[1] "execname -i seqname -b 0 -m F84 -t e -v e -a e --append "
[1] "execname -i seqname -b 0 -m GTR -t e -c 1 --append "
[1] "execname -i seqname -b 0 -m GTR -t e -v e -c 1 --append "
[1] "execname -i seqname -b 0 -m GTR -t e -a e --append "
[1] "execname -i seqname -b 0 -m GTR -t e -v e -a e --append "
## end MWE
I think it's easily solved by replacing mdl with paste("-m", rep(c("JC69",
"K80", "F81", "F84""HKY85", "TN93", "GTR"), each = 4))
Anyways, hope to clarify this minor issue.
Cheers,
Luiz
