Nothing
#' @title Check a morphological matrix consistency levels.
#'
#' @description Performs a fast check of the phylogenetic signal in a morphological matrix using parsimony.
#'
#' @param matrix A discrete morphological matrix.
#' @param parsimony Either the parsimony algorithm to be passed to \code{\link[phangorn]{optim.parsimony}} or a parsimony function that can take a \code{\link[phangorn]{phyDat}} object as an input (\code{default = "fitch"}).
#' @param first.tree A list of functions to generate the first most parsimonious tree (default = \code{c(\link[phangorn]{dist.hamming}, \link[phangorn]{NJ})}; see details).
#' @param orig.tree Optional, the input tree to measure the distance between the parsimony and the original tree.
#' @param distance Optional, if orig.tree is provided, the function to use for measuring distance between the trees (default = \code{link[phangorn]{RF.dist}}).
#' @param ... Any additional arguments to be passed to the parsimony algorithm.
#' @param contrast.matrix An optional contrast matrix. By default, the function recognises any character state token as different apart from \code{?} that is treated as all characters.
#' @param verbose Whether to be verbose or not (\code{default = FALSE}).
#'
#' @return
#' Returns the parsimony score (using \code{\link[phangorn]{parsimony}}), the consistency and retention indices (using \code{\link[phangorn]{CI}} and \code{\link[phangorn]{RI}}) from the most parsimonious tree obtained from the matrix.
#' Can also return the topological distance from the original tree if provided.
#'
#' @details
#' \itemize{
#' \item The \code{first.tree} argument must be a list of functions to be used in a cascade to transform the matrix (as a \code{\link[phangorn]{phyDat}} object) into a tree using the functions iteratively.
#' For example the default \code{c(\link[phangorn]{dist.hamming}, \link[phangorn]{NJ})} will apply the following to the matrix: \code{\link[phangorn]{NJ}(\link[phangorn]{dist.hamming}(matrix))}
#' }
#'
#' @examples
#' ## Generating a random tree
#' random_tree <- rcoal(10)
#'
# \dontrun{
#' ## Generating a random matrix
#' random_matrix <- sim.morpho(random_tree, characters = 50, model = "ER",
#' rates = c(rgamma, 1, 1))
#'
#' ## Checking the matrix scores
#' check.morpho(random_matrix, orig.tree = random_tree)
# }
#'
#' @seealso \code{\link{sim.morpho}}, \code{\link{get.contrast.matrix}}, \code{\link[phangorn]{optim.parsimony}}
#'
#' @author Thomas Guillerme
#DEBUG
# stop("DEBUG check.morpho")
# source("sanitizing.R")
# random.tree <- rcoal(10)
# matrix <- sim.morpho(random.tree, characters = 50, model = "ER", rates = c(rgamma, 1, 1))
# orig.tree = random.tree
# parsimony = "fitch"
# first.tree = c(phangorn::dist.hamming, phangorn::NJ)
# distance = phangorn::RF.dist
check.morpho <- function(matrix, orig.tree, parsimony = "fitch", first.tree = c(phangorn::dist.hamming, phangorn::NJ), distance = phangorn::RF.dist, ..., contrast.matrix, verbose = FALSE) {
#SANITIZNG
#matrix
check.class(matrix, "matrix")
#parsimony
if(!is(parsimony, "function")) {
#model is not a sure function
implemented_parsimony <- c("fitch", "sankoff")
if(all(is.na(match(parsimony, implemented_parsimony)))) {
stop.call("", paste0("The parsimony argument must be either a user's function or one of the following: ", paste(implemented_parsimony, collapse = ", ")))
}
#setting the parsimony algorithm
use.optim.parsimony <- TRUE
parsimony.algorithm <- phangorn::optim.parsimony
method <- parsimony
} else {
stop.call("", "User functions not implemented yet for parsimony argument.")
# use.optim.parsimony <- FALSE
# parsimony.algorithm <- phangorn::parsimony
}
#first.tree
if(!is(first.tree, "function")) {
if(any(unlist(lapply(first.tree, class)) != "function")) {
stop.call("", "first.tree argument must be a list of functions to calculate the first tree.")
}
}
#orig.tree
if(!missing(orig.tree)) {
check.class(orig.tree, "phylo")
#must be same size as the matrix
if(any(sort(row.names(matrix)) != sort(orig.tree$tip.label))) {
stop.call("", "Provided orig.tree has not the same taxa as the matrix.")
}
#distance
check.class(distance, "function")
}
#contrast.matrix
if(!missing(contrast.matrix)) {
check.class(contrast.matrix, "matrix")
}
#verbose
check.class(verbose, "logical")
#distance
check.class(distance, "function")
#CHECKING THE MATRIX
#Creating the contrast matrix
if(missing(contrast.matrix)) {
contrast.matrix <- get.contrast.matrix(matrix)
}
#Creating the phyDat object
matrix_phyDat <- phangorn::phyDat(data = matrix, type = "USER", contrast = contrast.matrix)
#Calcualte the first tree
if(is(first.tree, "function")) {
#Calculate the first tree from the phyDat
first_tree <- first.tree(matrix_phyDat)
} else {
#Calculate the operation on the phyDat
first_tree <- first.tree[[1]](matrix_phyDat)
for(operations in 2:length(first.tree)) {
#Transform the first_tree using the other operations (making sure it's binary)
first_tree <- first.tree[[operations]](first_tree)
}
}
#Get the quick and dirty most parsimonious tree
# if(use.optim.parsimony == TRUE) {
verbose.pars <- utils::capture.output(MP_tree <- multi2di(parsimony.algorithm(tree = multi2di(first_tree), data = matrix_phyDat, method = method, ...)))
#verbose.pars <- utils::capture.output(MP_tree <- multi2di(parsimony.algorithm(tree = first_tree, data = matrix_phyDat, method = method))) ; warning("DEBUG")
# } else {
# verbose.pars <- utils::capture.output(MP_tree <- parsimony.algorithm(tree = first_tree, data = matrix_phyDat))
# }
if(verbose != FALSE) {
cat("Most parsimonious tree search:\n")
cat(verbose.pars)
}
#Get the parsimony score
parsimony_score <- phangorn::parsimony(MP_tree, matrix_phyDat)
#Get the CI
consistency_index <- phangorn::CI(MP_tree, matrix_phyDat)
#Get the retention Index
retention_index <- phangorn::RI(MP_tree, matrix_phyDat)
if(!missing(orig.tree)) {
#Get the distance between the trees
tree_distance <- distance(MP_tree, unroot(orig.tree))
#Get the data out vectors (with distance)
data_out <- c(parsimony_score, consistency_index, retention_index, tree_distance)
out_names <- c("Maximum parsimony", "Consistency index", "Retention index", "Robinson-Foulds distance")
} else {
#Get the data out vectors (without distance)
data_out <- c(parsimony_score, consistency_index, retention_index)
out_names <- c("Maximum parsimony", "Consistency index", "Retention index")
}
return(matrix(data = data_out, ncol = 1, dimnames = list(out_names, c(""))))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.