From dd850ca3bea72dfc29183ff3f1df515d9f2bf58a Mon Sep 17 00:00:00 2001 From: xiangpin Date: Fri, 4 Jun 2021 17:11:04 +0800 Subject: [PATCH 1/5] optimizing read.mcmctree --- NAMESPACE | 5 ++++ R/MCMCTree.R | 38 +++++++++++++++++++++++++++++-- R/method-as-ultrametric.R | 48 +++++++++++++++++++++++++++++++++++++++ man/as.ultrametric.Rd | 19 ++++++++++++++++ man/read.mcmctree.Rd | 15 ++++++++++-- 5 files changed, 121 insertions(+), 4 deletions(-) create mode 100644 R/method-as-ultrametric.R create mode 100644 man/as.ultrametric.Rd diff --git a/NAMESPACE b/NAMESPACE index 0bcc582..97294a1 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -22,6 +22,8 @@ S3method(as.treedata,phylo4d) S3method(as.treedata,pml) S3method(as.treedata,pvclust) S3method(as.treedata,tbl_df) +S3method(as.ultrametric,phylo) +S3method(as.ultrametric,treedata) S3method(as_tibble,pvclust) S3method(child,phylo) S3method(child,treedata) @@ -59,6 +61,7 @@ export(Ntip) export(ancestor) export(as.phylo) export(as.treedata) +export(as.ultrametric) export(as_tibble) export(child) export(drop.tip) @@ -131,9 +134,11 @@ importFrom(ape,as.AAbin.character) importFrom(ape,as.DNAbin) importFrom(ape,as.DNAbin.character) importFrom(ape,as.phylo) +importFrom(ape,branching.times) importFrom(ape,checkLabel) importFrom(ape,extract.clade) importFrom(ape,is.rooted) +importFrom(ape,is.ultrametric) importFrom(ape,read.FASTA) importFrom(ape,read.nexus) importFrom(ape,read.tree) diff --git a/R/MCMCTree.R b/R/MCMCTree.R index 911d700..67b92d4 100644 --- a/R/MCMCTree.R +++ b/R/MCMCTree.R @@ -2,13 +2,21 @@ ##' ##' @title read.mcmctree ##' @param file the output tree file of MCMCTree +##' @param as.ultrametric logical whether convert the tree +##' to be ultrametric, if it is not ultrametric, default is FALSE. +##' When the tree is ultrametric, branch times will be calculated +##' automatically. +##' @param option integer it should be one of 1 or 2, default is 2. +##' see also the 'option' of 'is.ultrametric' of 'ape'. +##' @param ... additional arguments, see also the arguments of +##' 'is.ultrametric' of 'ape'. ##' @return treedata object ##' @export ##' @examples ##' file <- system.file("extdata/MCMCTree", "mcmctree_output.tree", package="treeio") -##' tr <- read.mcmctree(file) +##' tr <- read.mcmctree(file, option=2) ##' tr -read.mcmctree <- function(file){ +read.mcmctree <- function(file, as.ultrametric = FALSE, option=2, ...){ text <- readLines(file) ind <- grep("^.*tree.*=.*", text, ignore.case=TRUE) text[ind] <- gsub("^.*TREE", "TREE", text[ind], ignore.case=TRUE) @@ -18,10 +26,36 @@ read.mcmctree <- function(file){ obj <- read.beast(file=newfile) if(inherits(obj, "treedata")){ obj@file <- filename(file) + obj <- add.branch.time.mcmctree( + obj=obj, + as.ultrametric = as.ultrametric, + option=option, + ...) }else{ for (i in seq_len(length(obj))){ obj[[i]]@file <- filename(file) + obj[[i]] <- add.branch.time.mcmctree( + obj=obj[[i]], + as.ultrametric = as.ultrametric, + option=option, + ...) } } return(obj) } + +#' @importFrom ape branching.times +#' @importFrom ape is.ultrametric +add.branch.time.mcmctree <- function(obj, as.ultrametric=FALSE, ...){ + if (as.ultrametric && !is.ultrametric(obj@phylo, ...)){ + obj <- as.ultrametric(obj) + } + if (is.ultrametric(obj@phylo, ...)){ + xx <- branching.times(obj@phylo) + ntip <- Ntip(obj) + N <- Nnode(obj, internal.only = FALSE) + xx <- data.frame(node=as.character((ntip+1):N), reltime=as.vector(xx)) + obj@data <- full_join(obj@data, xx, by="node") + } + return(obj) +} diff --git a/R/method-as-ultrametric.R b/R/method-as-ultrametric.R new file mode 100644 index 0000000..893c5eb --- /dev/null +++ b/R/method-as-ultrametric.R @@ -0,0 +1,48 @@ +#' @title as.ultrametric +#' @param tree tree object +#' @param ... additional parameters +#' @return treedata or phylo object +#' @export +as.ultrametric <- function(tree, ...){ + UseMethod("as.ultrametric") +} + +#' @method as.ultrametric phylo +#' @export +## reference +## https://github.com/PuttickMacroevolution/MCMCtreeR/blob/master/R/readMCMCTree.R +as.ultrametric.phylo <- function(tree, ...){ + outer <- tree$edge[, 2] + inner <- tree$edge[, 1] + ntip <- Ntip(tree) + totalPath <- c() + tipindx <- which(outer <= ntip) + for (i in tipindx) { + start <- i + end <- inner[start] + edgeTimes <- tree$edge.length[start] + while (end != inner[1]) { + start <- which(outer == end) + end <- inner[start] + edgeTimes <- c(edgeTimes, tree$edge.length[start]) + } + totalPath <- c(totalPath, sum(edgeTimes)) + } + addLength <- max(totalPath) - totalPath + tree$edge.length[tipindx] <- tree$edge.length[tipindx] + addLength + return (tree) +} + +#' @method as.ultrametric treedata +#' @export +as.ultrametric.treedata <- function(tree, ...){ + tree@phylo <- as.ultrametric(tree=tree@phylo,...) + return (tree) +} + +#' @method as.ultrametric tbl_tree +as.ultrametric.tbl_tree <- function(tree, ...){ + tree <- as.treedata(tree) + tree <- as.ultrametric(tree) + return(tree) +} diff --git a/man/as.ultrametric.Rd b/man/as.ultrametric.Rd new file mode 100644 index 0000000..127ece8 --- /dev/null +++ b/man/as.ultrametric.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/method-as-ultrametric.R +\name{as.ultrametric} +\alias{as.ultrametric} +\title{as.ultrametric} +\usage{ +as.ultrametric(tree, ...) +} +\arguments{ +\item{tree}{tree object} + +\item{...}{additional parameters} +} +\value{ +treedata or phylo object +} +\description{ +as.ultrametric +} diff --git a/man/read.mcmctree.Rd b/man/read.mcmctree.Rd index c318f27..9d561e2 100644 --- a/man/read.mcmctree.Rd +++ b/man/read.mcmctree.Rd @@ -4,10 +4,21 @@ \alias{read.mcmctree} \title{read.mcmctree} \usage{ -read.mcmctree(file) +read.mcmctree(file, as.ultrametric = FALSE, option = 2, ...) } \arguments{ \item{file}{the output tree file of MCMCTree} + +\item{as.ultrametric}{logical whether convert the tree +to be ultrametric, if it is not ultrametric, default is FALSE. +When the tree is ultrametric, branch times will be calculated +automatically.} + +\item{option}{integer it should be one of 1 or 2, default is 2. +see also the 'option' of 'is.ultrametric' of 'ape'.} + +\item{...}{additional arguments, see also the arguments of +'is.ultrametric' of 'ape'.} } \value{ treedata object @@ -17,6 +28,6 @@ read MCMCTree output Tree } \examples{ file <- system.file("extdata/MCMCTree", "mcmctree_output.tree", package="treeio") -tr <- read.mcmctree(file) +tr <- read.mcmctree(file, option=2) tr } From 8046a339a3e4f122790bb1d80e6b585b1c9ddad2 Mon Sep 17 00:00:00 2001 From: xiangpin Date: Fri, 4 Jun 2021 18:14:47 +0800 Subject: [PATCH 2/5] check ultrametric with two methods --- R/MCMCTree.R | 21 ++++++++------------- man/read.mcmctree.Rd | 10 ++-------- 2 files changed, 10 insertions(+), 21 deletions(-) diff --git a/R/MCMCTree.R b/R/MCMCTree.R index 67b92d4..b7045fa 100644 --- a/R/MCMCTree.R +++ b/R/MCMCTree.R @@ -6,17 +6,13 @@ ##' to be ultrametric, if it is not ultrametric, default is FALSE. ##' When the tree is ultrametric, branch times will be calculated ##' automatically. -##' @param option integer it should be one of 1 or 2, default is 2. -##' see also the 'option' of 'is.ultrametric' of 'ape'. -##' @param ... additional arguments, see also the arguments of -##' 'is.ultrametric' of 'ape'. ##' @return treedata object ##' @export ##' @examples ##' file <- system.file("extdata/MCMCTree", "mcmctree_output.tree", package="treeio") -##' tr <- read.mcmctree(file, option=2) +##' tr <- read.mcmctree(file) ##' tr -read.mcmctree <- function(file, as.ultrametric = FALSE, option=2, ...){ +read.mcmctree <- function(file, as.ultrametric = FALSE){ text <- readLines(file) ind <- grep("^.*tree.*=.*", text, ignore.case=TRUE) text[ind] <- gsub("^.*TREE", "TREE", text[ind], ignore.case=TRUE) @@ -29,16 +25,14 @@ read.mcmctree <- function(file, as.ultrametric = FALSE, option=2, ...){ obj <- add.branch.time.mcmctree( obj=obj, as.ultrametric = as.ultrametric, - option=option, - ...) + ) }else{ for (i in seq_len(length(obj))){ obj[[i]]@file <- filename(file) obj[[i]] <- add.branch.time.mcmctree( obj=obj[[i]], - as.ultrametric = as.ultrametric, - option=option, - ...) + as.ultrametric = as.ultrametric + ) } } return(obj) @@ -47,10 +41,11 @@ read.mcmctree <- function(file, as.ultrametric = FALSE, option=2, ...){ #' @importFrom ape branching.times #' @importFrom ape is.ultrametric add.branch.time.mcmctree <- function(obj, as.ultrametric=FALSE, ...){ - if (as.ultrametric && !is.ultrametric(obj@phylo, ...)){ + flag_ultrametric <- !is.ultrametric(obj@phylo, option=2) || !is.ultrametric(obj@phylo) + if (as.ultrametric && flag_ultrametric){ obj <- as.ultrametric(obj) } - if (is.ultrametric(obj@phylo, ...)){ + if (flag_ultrametric){ xx <- branching.times(obj@phylo) ntip <- Ntip(obj) N <- Nnode(obj, internal.only = FALSE) diff --git a/man/read.mcmctree.Rd b/man/read.mcmctree.Rd index 9d561e2..cabf9a0 100644 --- a/man/read.mcmctree.Rd +++ b/man/read.mcmctree.Rd @@ -4,7 +4,7 @@ \alias{read.mcmctree} \title{read.mcmctree} \usage{ -read.mcmctree(file, as.ultrametric = FALSE, option = 2, ...) +read.mcmctree(file, as.ultrametric = FALSE) } \arguments{ \item{file}{the output tree file of MCMCTree} @@ -13,12 +13,6 @@ read.mcmctree(file, as.ultrametric = FALSE, option = 2, ...) to be ultrametric, if it is not ultrametric, default is FALSE. When the tree is ultrametric, branch times will be calculated automatically.} - -\item{option}{integer it should be one of 1 or 2, default is 2. -see also the 'option' of 'is.ultrametric' of 'ape'.} - -\item{...}{additional arguments, see also the arguments of -'is.ultrametric' of 'ape'.} } \value{ treedata object @@ -28,6 +22,6 @@ read MCMCTree output Tree } \examples{ file <- system.file("extdata/MCMCTree", "mcmctree_output.tree", package="treeio") -tr <- read.mcmctree(file, option=2) +tr <- read.mcmctree(file) tr } From 69263feab03da476ffee23257b87f6e163e53726 Mon Sep 17 00:00:00 2001 From: xiangpin Date: Tue, 8 Jun 2021 11:38:39 +0800 Subject: [PATCH 3/5] fix the check of ultrametric --- R/MCMCTree.R | 17 +++++++++-------- man/read.mcmctree.Rd | 4 ++-- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/R/MCMCTree.R b/R/MCMCTree.R index b7045fa..baa64bf 100644 --- a/R/MCMCTree.R +++ b/R/MCMCTree.R @@ -2,7 +2,7 @@ ##' ##' @title read.mcmctree ##' @param file the output tree file of MCMCTree -##' @param as.ultrametric logical whether convert the tree +##' @param force.ultrametric logical whether convert the tree ##' to be ultrametric, if it is not ultrametric, default is FALSE. ##' When the tree is ultrametric, branch times will be calculated ##' automatically. @@ -12,7 +12,7 @@ ##' file <- system.file("extdata/MCMCTree", "mcmctree_output.tree", package="treeio") ##' tr <- read.mcmctree(file) ##' tr -read.mcmctree <- function(file, as.ultrametric = FALSE){ +read.mcmctree <- function(file, force.ultrametric = FALSE){ text <- readLines(file) ind <- grep("^.*tree.*=.*", text, ignore.case=TRUE) text[ind] <- gsub("^.*TREE", "TREE", text[ind], ignore.case=TRUE) @@ -24,14 +24,14 @@ read.mcmctree <- function(file, as.ultrametric = FALSE){ obj@file <- filename(file) obj <- add.branch.time.mcmctree( obj=obj, - as.ultrametric = as.ultrametric, + force.ultrametric = force.ultrametric, ) }else{ for (i in seq_len(length(obj))){ obj[[i]]@file <- filename(file) obj[[i]] <- add.branch.time.mcmctree( obj=obj[[i]], - as.ultrametric = as.ultrametric + force.ultrametric = force.ultrametric ) } } @@ -40,10 +40,11 @@ read.mcmctree <- function(file, as.ultrametric = FALSE){ #' @importFrom ape branching.times #' @importFrom ape is.ultrametric -add.branch.time.mcmctree <- function(obj, as.ultrametric=FALSE, ...){ - flag_ultrametric <- !is.ultrametric(obj@phylo, option=2) || !is.ultrametric(obj@phylo) - if (as.ultrametric && flag_ultrametric){ - obj <- as.ultrametric(obj) +add.branch.time.mcmctree <- function(obj, force.ultrametric=FALSE, ...){ + flag_ultrametric <- is.ultrametric(obj@phylo, option=2) || is.ultrametric(obj@phylo) + if (force.ultrametric && flag_ultrametric){ + message("This tree is not ultrametric, and you has set force.ultrametric to TRUE, so the tree will be convert to ultrametric automatically!") + obj <- force.ultrametric(obj) } if (flag_ultrametric){ xx <- branching.times(obj@phylo) diff --git a/man/read.mcmctree.Rd b/man/read.mcmctree.Rd index cabf9a0..e9b62c9 100644 --- a/man/read.mcmctree.Rd +++ b/man/read.mcmctree.Rd @@ -4,12 +4,12 @@ \alias{read.mcmctree} \title{read.mcmctree} \usage{ -read.mcmctree(file, as.ultrametric = FALSE) +read.mcmctree(file, force.ultrametric = FALSE) } \arguments{ \item{file}{the output tree file of MCMCTree} -\item{as.ultrametric}{logical whether convert the tree +\item{force.ultrametric}{logical whether convert the tree to be ultrametric, if it is not ultrametric, default is FALSE. When the tree is ultrametric, branch times will be calculated automatically.} From b67d8ff760d668051aeaecb34715ce0ae317d22e Mon Sep 17 00:00:00 2001 From: Guangchuang Yu Date: Tue, 8 Jun 2021 16:27:52 +0800 Subject: [PATCH 4/5] Update MCMCTree.R --- R/MCMCTree.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/MCMCTree.R b/R/MCMCTree.R index baa64bf..1902d7a 100644 --- a/R/MCMCTree.R +++ b/R/MCMCTree.R @@ -43,7 +43,7 @@ read.mcmctree <- function(file, force.ultrametric = FALSE){ add.branch.time.mcmctree <- function(obj, force.ultrametric=FALSE, ...){ flag_ultrametric <- is.ultrametric(obj@phylo, option=2) || is.ultrametric(obj@phylo) if (force.ultrametric && flag_ultrametric){ - message("This tree is not ultrametric, and you has set force.ultrametric to TRUE, so the tree will be convert to ultrametric automatically!") + message("This tree is not ultrametric, and you has set force.ultrametric to TRUE, so the tree will be converted to ultrametric automatically!") obj <- force.ultrametric(obj) } if (flag_ultrametric){ From 1ba1f9f4997910f206dd064379349d28832a973a Mon Sep 17 00:00:00 2001 From: xiangpin Date: Tue, 8 Jun 2021 20:00:28 +0800 Subject: [PATCH 5/5] typo for as.ultrametric --- R/MCMCTree.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/MCMCTree.R b/R/MCMCTree.R index baa64bf..4b034fa 100644 --- a/R/MCMCTree.R +++ b/R/MCMCTree.R @@ -44,7 +44,7 @@ add.branch.time.mcmctree <- function(obj, force.ultrametric=FALSE, ...){ flag_ultrametric <- is.ultrametric(obj@phylo, option=2) || is.ultrametric(obj@phylo) if (force.ultrametric && flag_ultrametric){ message("This tree is not ultrametric, and you has set force.ultrametric to TRUE, so the tree will be convert to ultrametric automatically!") - obj <- force.ultrametric(obj) + obj <- as.ultrametric(obj) } if (flag_ultrametric){ xx <- branching.times(obj@phylo)