@@ -1,7 +1,7 @@

 Package: ggtree

 Type: Package

 Title: a phylogenetic tree viewer for different types of tree annotations

-Version: 1.1.1

+Version: 1.1.2

 Author: Guangchuang Yu and Tommy Tsan-Yuk Lam

 Maintainer: Guangchuang Yu <guangchuangyu@gmail.com>

 Description: ggtree extends the ggplot2 plotting system which implemented the

@@ -44,6 +44,7 @@ export(groupClade)

 export(groupOTU)

 export(gzoom)

 export(hilight)

+export(merge_tree)

 export(phylopic)

 export(plot)

 export(read.baseml)

@@ -112,7 +113,9 @@ importFrom(ggplot2,geom_segment)

 importFrom(ggplot2,geom_text)

 importFrom(ggplot2,geom_tile)

 importFrom(ggplot2,ggplot)

+importFrom(ggplot2,labs)

 importFrom(ggplot2,scale_color_manual)

+importFrom(ggplot2,scale_fill_discrete)

 importFrom(ggplot2,scale_fill_gradient)

 importFrom(ggplot2,scale_x_continuous)

 importFrom(ggplot2,scale_x_reverse)

@@ -1,3 +1,19 @@

+CHANGES IN VERSION 1.1.2

+------------------------

+ o update vignette to add example of merge_tree <2015-04-29, Wed>

+ o in addition to parsing beast time scale tree in XXX_year[\\.\\d]*, now supports XXX/year[\\.\\d]* <2015-04-29, Wed>

+ o add examples folder in inst that contains sample data <2015-04-29, Wed>

+ o update gplot, now rowname of heatmap will not be displayed <2015-04-28, Tue>

+ o add line break if substitution longer than 50 character <2015-04-28, Tue>

+ o upport calculating branch for time scale tree <2015-04-28, Tue>

+ o remove parsing tip sequence from mlb and mlc file <2015-04-28, Tue>

+ o remove tip.fasfile in read.paml_rst for rstfile already contains tip sequence <2015-04-28, Tue>

+ o scale_color accepts user specific interval and output contains 'scale' attribute that can be used for adding legend <2015-04-28, Tue>

+ o extend fortify methods to support additional fields <2015-04-28, Tue>

+ o extend get.fields methods to support additional fields <2015-04-28, Tue>

+ o extend tree class to support additional info by merging two tree <2015-04-28, Tue>

+ o implement merge_tree function to merge two tree objects into one <2015-04-28, Tue>

+ 

 CHANGES IN VERSION 1.1.1

 ------------------------

  o minor bug fixed in extracting node ID of rst file <2015-04-27, Mon>

@@ -18,6 +18,7 @@ setOldClass("phylo")

 ##' @slot tip.fasfile fasta file of tip sequences

 ##' @slot tree.file tree file

 ##' @slot ancseq.file ancestral sequence file, nexus format

+##' @slot extraInfo extra information

 ##' @exportClass hyphy

 ##' @author Guangchuang Yu \url{http://ygc.name}

 ##' @seealso \linkS4class{paml_rst}

@@ -34,7 +35,8 @@ setClass("hyphy",

              tip_seq     = "character",

              tip.fasfile = "character",

              tree.file   = "character",

-             ancseq.file = "character"

+             ancseq.file = "character",

+             extraInfo   = "data.frame"

              )

          )

@@ -60,8 +62,8 @@ setClass("hyphy",

 ##' @slot joint_subs sequence substitutions based on joint_ancseq

 ##' @slot marginal_AA_subs Amino acid sequence substitutions based on marginal_ancseq

 ##' @slot joint_AA_subs Amino acid sequence substitutions based on joint_ancseq

-##' @slot tip.fasfile fasta file of tip sequences

 ##' @slot rstfile rst file

+##' @slot extraInfo extra information

 ##' @exportClass paml_rst

 ##' @author Guangchuang Yu \url{http://ygc.name}

 ##' @seealso \linkS4class{codeml} \linkS4class{codeml_mlc}

@@ -79,8 +81,8 @@ setClass("paml_rst",

              joint_subs       = "data.frame",

              marginal_AA_subs = "data.frame",

              joint_AA_subs    = "data.frame",

-             tip.fasfile      = "character",

-             rstfile          = "character"

+             rstfile          = "character",

+             extraInfo        = "data.frame"

          )

          )

@@ -94,9 +96,8 @@ setClass("paml_rst",

 ##' @slot treetext tree text

 ##' @slot phylo phylo object

 ##' @slot dNdS dN dS information

-##' @slot seq_type one of "NT" and "AA"

-##' @slot tip_seq sequences of tips

 ##' @slot mlcfile mlc file

+##' @slot extraInfo extra information

 ##' @exportClass codeml_mlc

 ##' @author Guangchuang Yu

 ##' @seealso \linkS4class{paml_rst} \linkS4class{codeml}

@@ -107,9 +108,10 @@ setClass("codeml_mlc",

              treetext   = "character",

              phylo      = "phylo",

              dNdS       = "matrix",

-             seq_type   = "character",

-             tip_seq    = "character",

-             mlcfile    = "character"

+             ## seq_type   = "character",

+             ## tip_seq    = "character",

+             mlcfile    = "character",

+             extraInfo  = "data.frame"

              )

          )

@@ -121,13 +123,15 @@ setClass("codeml_mlc",

 ##' @docType class

 ##' @slot mlc A \code{code_mlc} object

 ##' @slot rst A \code{paml_rst} object

+##' @slot extraInfo extra information

 ##' @exportClass codeml

 ##' @seealso \linkS4class{codeml_mlc} \linkS4class{paml_rst}

 ##' @keywords codeml

 setClass("codeml",

          representation = representation(

-             mlc = "codeml_mlc",

-             rst = "paml_rst"

+             mlc       = "codeml_mlc",

+             rst       = "paml_rst",

+             extraInfo = "data.frame"

              )

          )

@@ -151,6 +155,7 @@ setClass("codeml",

 ##' @slot version version

 ##' @slot metadata metadata

 ##' @slot file jplace file

+##' @slot extraInfo extra information

 ##' @exportClass jplace

 ##' @author Guangchuang Yu \url{http://ygc.name}

 ##' @seealso \code{\link{show}} \code{\link{get.tree}}

@@ -164,7 +169,8 @@ setClass("jplace",

              placements = "data.frame",

              version    = "numeric",

              metadata   = "list",

-             file       = "character"

+             file       = "character",

+             extraInfo  = "data.frame"

              )

          )

@@ -183,6 +189,7 @@ setClass("jplace",

 ##' @slot translation tip number to name translation in beast file

 ##' @slot stats beast statistics

 ##' @slot file beast file, nexus format

+##' @slot extraInfo extra information

 ##' @exportClass beast

 ##' @author Guangchuang Yu \url{http://ygc.name}

 ##' @seealso \code{\link{show}} \code{\link{get.fields}}

@@ -195,6 +202,7 @@ setClass("beast",

              phylo       = "phylo",

              translation = "matrix",

              stats       = "data.frame",

-             file        = "character"

+             file        = "character",

+             extraInfo   = "data.frame"

              )

          )

@@ -96,7 +96,7 @@ setMethod("show", signature(object = "beast"),

               cat("...@ tree: ")

               print.phylo(get.tree(object))                  

               cat("\nwith the following features available:\n")

-              print_fields(object)              

+              print_fields(object)

           })

 ##' @rdname groupOTU-methods

@@ -150,7 +150,7 @@ setMethod("get.tree", signature(object="beast"),

 ##' @exportMethod get.fields

 setMethod("get.fields", signature(object="beast"),

           function(object, ...) {

-              object@fields

+              get.fields.tree(object)

           }

           )

@@ -12,9 +12,9 @@

 ##' mlcfile <- system.file("extdata/PAML_Codeml", "mlc", package="ggtree")

 ##' read.codeml(rstfile, mlcfile) 

 read.codeml <- function(rstfile, mlcfile) {

-    rst = read.paml_rst(rstfile)

-    mlc = read.codeml_mlc(mlcfile)

-    rst@tip_seq <- mlc@tip_seq

+    rst <- read.paml_rst(rstfile)

+    mlc <- read.codeml_mlc(mlcfile)

+    ## rst@tip_seq <- mlc@tip_seq

     new("codeml",

         rst = set.paml_rst_(rst),

         mlc = mlc

@@ -99,10 +99,8 @@ setMethod("get.subs", signature(object = "codeml"),

 ##' @rdname get.fields-methods

 ##' @exportMethod get.fields

 setMethod("get.fields", signature(object="codeml"),

-          function(object) {

-              fields <- c(get.fields(object@rst),

-                          get.fields(object@mlc))

-              return(unique(fields))

+          function(object, ...) {

+              get.fields.tree(object)

           }

           )

@@ -10,7 +10,7 @@

 ##' mlcfile <- system.file("extdata/PAML_Codeml", "mlc", package="ggtree")

 ##' read.codeml_mlc(mlcfile)

 read.codeml_mlc <- function(mlcfile) {

-    tip_seq <- read.tip_seq_mlc(mlcfile)

+    ## tip_seq <- read.tip_seq_mlc(mlcfile)

     dNdS <- read.dnds_mlc(mlcfile)

     new("codeml_mlc",

@@ -18,8 +18,8 @@ read.codeml_mlc <- function(mlcfile) {

         treetext = read.treetext_paml_mlc(mlcfile),

         phylo    = read.phylo_paml_mlc(mlcfile),

         dNdS     = dNdS,

-        seq_type = get_seqtype(tip_seq),

-        tip_seq  = tip_seq,

+        ## seq_type = get_seqtype(tip_seq),

+        ## tip_seq  = tip_seq,

         mlcfile  = mlcfile)

 }

@@ -79,7 +79,7 @@ setMethod("show", signature(object = "codeml_mlc"),

 ##' @exportMethod get.fields

 setMethod("get.fields", signature(object = "codeml_mlc"),

           function(object) {

-              object@fields

+              get.fields.tree(object)

           })

 ##' @rdname plot-methods

deleted file mode 100644

@@ -1,79 +0,0 @@

-

-merge_baseml_hyphy <- function(baseml, hyphy) {

-    ml_phylo <- baseml@mlc@phylo

-    hy_phylo <- hyphy@phylo

-    phylo_merge_check(ml_phylo, hy_phylo)

-}

-

-

-phylo_merge_check <- function(phy1, phy2) {

-    if (Ntip(phy1) != Ntip(phy2)) {

-        stop("Number of tip is not eqaul...")

-    }

-    if ( ! all(phy1$tip.label %in% phy2$tip.label)) {

-        stop("tip label is not consistent...")

-    }

-

-    if (phy1$Nnode != phy2$Nnode) {

-        stop("Number of node is not eqaul...")

-    }

-    if (!is.null(phy1$node.label) && !is.null(phy2$node.label)) {

-        n <- getNodeNum(phy1)

-        node <- (Ntip(phy1)+1):n

-        if (all(phy1$node.label == node) || all(phy2$node.label == node)) {

-            ## nothing happend.

-        } else if (any(phy1$node.label != phy2$node.label)){

-            stop("node label is not consistent...")

-        }

-    }

-    

-    degree1 <- table(phy1$edge[,1]) %>% sort

-    degree2 <- table(phy2$edge[,1]) %>% sort

-    if (length(degree1) != length(degree2)) {

-        stop("node degree is not consistent...")

-    } else if (!all(degree1 == degree2)) {

-        stop("node degree is not consistent...")

-    }

-}

-

-

-node_mapper <- function(phy1, phy2) {

-   df1 <- fortify(phy1)

-   df2 <- fortify(phy2)

-

-   df1$node2 <- NA

-   df1$visited <- FALSE

-   ntip <- Ntip(phy1)

-   df1$node2[1:ntip] <- match(phy1$tip.label, phy2$tip.label)

- 

-   root <- getRoot(phy1)

-   currentNode <- 1:ntip

-   while(any(df1$visited == FALSE)) {

-       pNode <- c()

-       for (i in currentNode) {

-           df1[i, "visited"] <- TRUE

-           j <- df1[i, "node2"]

-           ip <- df1[i, "parent"]

-           if (df1[ip, "visited"] == FALSE) {

-               jp <- df2[j, "parent"]

-               df1[ip, "node2"] <- jp

-               df1[ip, "visited"] <- TRUE

-           }

-           pNode <- c(pNode, ip)

-       }

-       currentNode <- unique(pNode)

-   }

-

-   if (df1[root, "node"] != df1[root, "node2"]) {

-       stop("phylogenies not compatible...")

-   }

-   

-   

-   node <- df1$node

-   names(node) <- df1$node2

-   phy2$edge <- cbind(node[as.character(phy2$edge[,1])],

-                      node[as.character(phy2$edge[,2])])

-   phy2$tip.label <- phy1$tip.label

-

-   reorder(phy2)

-}

@@ -22,9 +22,10 @@

 ##' colnames(d) <- paste0("G", 1:4)

 ##' gplot(p, d, low="green", high="red")

 gplot <- function(p, data, low="green", high="red", widths=c(0.5, 0.5), font.size=14) {

-    ## p <- p + scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0))

-    p1 <- p + theme(panel.margin=unit(0, "null"))

-    p1 <- p1 + theme(plot.margin = unit(c(1, -1, 1.5, 1), "lines"))

+    ## p <- p + scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0.6))

+    p1 <- p + scale_y_continuous(expand = c(0, 0.6))

+    ## p1 <- p + theme(panel.margin=unit(0, "null"))

+    ## p1 <- p1 + theme(plot.margin = unit(c(1, -1, 1.5, 1), "lines"))

     p2 <- gplot.heatmap(p, data, low, high, font.size)

     grid.arrange(p1, p2, ncol=2, widths=widths)

     invisible(list(p1=p1, p2=p2))

@@ -33,8 +34,10 @@ gplot <- function(p, data, low="green", high="red", widths=c(0.5, 0.5), font.siz

 ##' @importFrom grid unit

 ##' @importFrom ggplot2 scale_fill_gradient

+##' @importFrom ggplot2 scale_fill_discrete

 ##' @importFrom ggplot2 element_text

 ##' @importFrom ggplot2 geom_tile

+##' @importFrom ggplot2 labs

 ##' @importFrom reshape2 melt

 gplot.heatmap <- function(p, data, low, high, font.size) {

     isTip <- x <- Var1 <- Var2 <- value <- NULL

@@ -46,25 +49,34 @@ gplot.heatmap <- function(p, data, low, high, font.size) {

     df=df[df$isTip,]

     dd$Var1 <- factor(dd$Var1, levels = df$label[order(df$y)])

+    if (any(dd$value == "")) {

+        dd$value[dd$value == ""] <- NA

+    }

     p2 <- ggplot(dd, aes(Var2, Var1, fill=value))+geom_tile(color="black")

-    p2 <- p2 + scale_fill_gradient(low=low, high=high)

+    if (is(dd$value,"numeric")) {

+        p2 <- p2 + scale_fill_gradient(low=low, high=high, na.value="white")

+    } else {

+        p2 <- p2 + scale_fill_discrete(na.value="white")

+    }

+    

     p2 <- p2+xlab("")+ylab("")

     p2 <- p2+theme_tree2() + theme(axis.ticks.x = element_blank(),

                                    axis.line.x=element_blank())

     ## p1 <- p1 + theme(axis.text.x = element_text(size = font.size))

     p2 <- p2 + theme(axis.ticks.margin = unit(0, "lines")) 

     p2 <- p2 + theme(axis.text.x = element_text(size = font.size))

-    p2 <- p2 + theme(axis.text.y = element_text(size=font.size))

+    ## p2 <- p2 + theme(axis.text.y = element_text(size=font.size))

     ## plot.margin   margin around entire plot (unit with the sizes of the top, right, bottom, and left margins) 

     ## units can be given in "lines" or  something more specific like "cm"...

     p2 <- p2 + theme(panel.margin=unit(0, "null"))

-    p2 <- p2 + theme(plot.margin = unit(c(1, 1, 1, -0.5), "lines"))

+    p2 <- p2 + theme(plot.margin = unit(c(1, 1, .5, -0.5), "lines"))

     p2 <- p2 + theme(legend.position = "right")

-

+    ## p2 <- p2 + labs(fill="")

+    

     return(p2)

 }

@@ -159,11 +159,11 @@ setMethod("get.tree", signature(object = "hyphy"),

 ##' @rdname get.fields-methods

 ##' @exportMethod get.fields

 setMethod("get.fields", signature(object = "hyphy"),

-          function(object) {

+          function(object, ...) {

               if(length(object@tip_seq) == 0) {

                   warning("tip sequence not available...\n")

               } else {

-                  object@fields

+                  get.fields.tree(object)

               }

           })

@@ -159,7 +159,7 @@ setMethod("get.treetext", signature(object = "jplace"),

 ##' get.fields(jp)

 setMethod("get.fields", signature(object = "jplace"),

           function(object, ...) {

-              object@fields

+              get.fields.tree(object)

           }

           )

new file mode 100644

@@ -0,0 +1,119 @@

+##' merge two tree object

+##'

+##' 

+##' @title merge_tree

+##' @param obj1 tree object 1

+##' @param obj2 tree object 2

+##' @return tree object

+##' @importFrom magrittr %<>%

+##' @importFrom ape Ntip

+##' @export

+##' @author Guangchuang Yu

+merge_tree <- function(obj1, obj2) {

+    

+    if (has.slot(obj1, "extraInfo") == FALSE) {

+        stop("input tree object is not supported...")

+    }

+    

+    if ((is.tree(obj1) & is.tree(obj2)) == FALSE) {

+        stop("input should be tree objects...")

+    }

+

+    tr1 <- get.tree(obj1)

+    tr2 <- get.tree(obj2)

+

+    if (getNodeNum(tr1) != getNodeNum(tr2)) {

+        stop("number of nodes not equals...")

+    }

+

+    if (Ntip(tr1) != Ntip(tr2)) {

+        stop("number of tips not equals...")

+    }

+    

+    if (all(tr1$tip.label %in% tr2$tip.label) == FALSE) {

+        stop("tip names not match...")

+    }

+

+

+    ## order tip.label in tr2 as in tr1

+    ## mapping corresponding ID

+    idx <- match(tr2$tip.label, tr1$tip.label)

+    tr2$edge[match(1:Ntip(tr2), tr2$edge[,2]), 2] <- idx

+    tr2$tip.label <- tr1$tip.label

+

+    node_map <- list()

+    node_map$from %<>% c(1:Ntip(tr2))

+    node_map$to %<>% c(idx)

+

+    root <- getRoot(tr1)

+    root.2 <- getRoot(tr2)

+    tr2$edge[tr2$edge[,1] == root.2, 1] <-  root

+

+    node_map$from %<>% c(root.2)

+    node_map$to %<>% c(root)

+

+    

+    currentNode <- 1:Ntip(tr1)

+    while(length(currentNode)) {

+        p1 <- sapply(currentNode, getParent, tr=tr1)

+        p2 <- sapply(currentNode, getParent, tr=tr2)

+

+        if (length(p1) != length(p2)) {

+            stop("trees are not identical...")

+        }

+

+        jj <- match(p2, tr2$edge[,1])

+        if (length(jj)) {

+            notNA <- which(!is.na(jj))

+            jj <- jj[notNA]

+        }

+        if (length(jj)) {

+            tr2$edge[jj,1] <- p1[notNA]

+        }

+

+        

+        ii <- match(p2, tr2$edge[,2])

+        if (length(ii)) {

+            notNA <- which(!is.na(ii))

+            ii <- ii[notNA]

+        }

+        if (length(ii)) {

+            tr2$edge[ii,2] <- p1[notNA]

+        }

+

+        node_map$from %<>% c(p2)

+        node_map$to %<>% c(p1)

+        

+        ## parent of root will return 0, which is in-valid node ID

+        currentNode <- unique(p1[p1 != 0])

+    }

+

+    if ( any(tr2$edge != tr2$edge) ) {

+        stop("trees are not identical...")

+    }

+    

+    node_map.df <- do.call("cbind", node_map)

+    node_map.df <- unique(node_map.df)

+    node_map.df <- node_map.df[node_map.df[,1] != 0,]

+    i <- order(node_map.df[,1], decreasing = FALSE)

+    node_map.df <- node_map.df[i,]

+

+    info2 <- fortify(obj2)

+    info2$node <- node_map.df[info2$node, 2]

+    info2$parent <- node_map.df[info2$parent, 2]

+

+    cn <- colnames(info2)

+    i <- match(c("x", "y", "isTip", "label", "branch", "branch.length"), cn)

+    i <- i[!is.na(i)]

+    info2 <- info2[, -i]

+

+    extraInfo <- obj1@extraInfo

+    if (nrow(extraInfo) == 0) {

+        obj1@extraInfo <- info2

+    } else {

+        info <- merge(extraInfo, info2, by.x =c("node", "parent"), by.y = c("node", "parent"))

+        obj1@extraInfo <- info

+    }

+    

+    return(obj1)

+}

@@ -13,7 +13,7 @@

 ##' read.baseml(rstfile, mlbfile)

 read.baseml <- function(rstfile, mlbfile) {

     res <- read.paml_rst(rstfile)

-    res@tip_seq <- read.tip_seq_mlb(mlbfile)

+    ## res@tip_seq <- read.tip_seq_mlb(mlbfile)

     set.paml_rst_(res)

 }

@@ -24,14 +24,13 @@ read.baseml <- function(rstfile, mlbfile) {

 ##' @importFrom Biostrings toString

 ##' @title read.paml_rst

 ##' @param rstfile rst file

-##' @param tip.fasfile fasta file of tips

 ##' @return A \code{paml_rst} object

 ##' @export

 ##' @author Guangchuang Yu \url{http://ygc.name}

 ##' @examples

 ##' rstfile <- system.file("extdata/PAML_Baseml", "rst", package="ggtree")

 ##' read.paml_rst(rstfile)

-read.paml_rst <- function(rstfile, tip.fasfile = NULL) {

+read.paml_rst <- function(rstfile) {

     ms <- read.ancseq_paml_rst(rstfile, by="Marginal")

     phylo <- read.phylo_paml_rst(rstfile)

     ## class(phylo) <- "list"

@@ -49,14 +48,16 @@ read.paml_rst <- function(rstfile, tip.fasfile = NULL) {

                joint_ancseq    = read.ancseq_paml_rst(rstfile, by = "Joint"),

                rstfile = rstfile

                )

-    if (!is.null(tip.fasfile)) {

-        seqs <- readBStringSet(tip.fasfile)

-        tip_seq <- sapply(seq_along(seqs), function(i) {

-            toString(seqs[i])

-        })

-        res@tip_seq <- tip_seq

-        res@tip.fasfile <- tip.fasfile

-    }

+    ## if (!is.null(tip.fasfile)) {

+    ##     seqs <- readBStringSet(tip.fasfile)

+    ##     tip_seq <- sapply(seq_along(seqs), function(i) {

+    ##         toString(seqs[i])

+    ##     })

+    ##     res@tip_seq <- tip_seq

+    ##     res@tip.fasfile <- tip.fasfile

+    ## }

+    res@tip_seq <- ms[names(ms) %in% phylo$tip.label]

+    

     set.paml_rst_(res)

 }

@@ -107,12 +108,12 @@ setMethod("get.tipseq", signature(object="paml_rst"),

 setMethod("show", signature(object = "paml_rst"),

           function(object) {

               cat("'paml_rst' S4 object that stored information of\n\t",

-                  paste0("'", object@rstfile, "'"))

-              if (length(object@tip.fasfile) != 0) {

-                  cat(paste0(" and \n\t'", object@tip.fasfile, "'.\n\n"))

-              } else {

-                  cat(".\n\n")

-              }

+                  paste0("'", object@rstfile, "'.\n\n"))

+              ## if (length(object@tip.fasfile) != 0) {

+              ##     cat(paste0(" and \n\t'", object@tip.fasfile, "'.\n\n"))

+              ## } else {

+              ##     cat(".\n\n")

+              ## }

               cat("...@ tree:")

               print.phylo(get.tree(object))                  

               cat("\nwith the following features available:\n")

@@ -129,7 +130,7 @@ setMethod("get.fields", signature(object = "paml_rst"),

               if (length(object@tip_seq) == 0) {

                   warning("tip sequence not available...\n")

               } else {

-                  object@fields

+                  get.fields.tree(object)

               }

           }

           )

@@ -1,15 +1,30 @@

 ##' @importFrom colorspace rainbow_hcl

-scale_color_ <- function(phylo, by, low=NULL, high=NULL, na.color=NULL, default.color="grey") {

+scale_color_ <- function(phylo, by, low=NULL, high=NULL, na.color=NULL, default.color="darkgrey", interval=NULL) {

+    df <- fortify(phylo)    

+    vals <- df[, by]

+

+    MIN=min(vals, na.rm=TRUE)

+    MAX=max(vals, na.rm=TRUE)

+

+    if (is.null(interval)) {

+        interval <- seq(MIN, MAX, length.out=100)

+    }

+    n <- length(interval)

+    

     if (!is.null(low) & ! is.null(high)) {

-        cols <- color_scale(c(low, high))

+        cols <- color_scale(low, high, n)

     } else {

-        cols <- rainbow_hcl(100)

+        cols <- rainbow_hcl(n)

     }

-    df <- fortify(phylo)

+

+    ## if (by == "dN_vs_dS") {

+    ##     interval <- seq(0, 1.5, length.out = 100)

+    ## } 

+

+    idx <- getIdx(vals, MIN=MIN, MAX=MAX, interval=interval)

+    interval <- attr(idx, "interval")

-    vals <- df[, by]

-    idx <- sapply(vals, getIdx, min(vals, na.rm=TRUE), max(vals, na.rm=TRUE))

     df$color <- cols[idx]

     tree <- get.tree(phylo)

@@ -38,7 +53,10 @@ scale_color_ <- function(phylo, by, low=NULL, high=NULL, na.color=NULL, default.

     }

     ## cols[is.na(cols)] <- "grey"

-    return(df$color)

+    color <- df$color

+

+    attr(color, "scale") <- list(interval=interval, color=cols)

+    return(color)

 }

 groupClade_ <- function(object, node) {

@@ -708,6 +726,27 @@ getYcoord_scale_numeric <- function(tr, df, yscale, ...) {

     return(df)

 }

+.assign_child_status <- function(tr, df, variable) {

+    yy <- df[, variable]

+    na.idx <- which(is.na(yy))

+    if (length(na.idx) > 0) {

+        tree <- get.tree(tr)

+        nodes <- rev(getNodes_by_postorder(tree))

+        for (curNode in nodes) {

+            parent <- getParent(tree, curNode)

+            if (parent == 0) { ## already reach root

+                next

+            }

+            idx <- which(is.na(yy[parent]))

+            if (length(idx) > 0) {

+                yy[parent[idx]] <- yy[curNode]

+            }

+        }

+    }

+    df[, variable] <- yy

+    return(df)

+}

+

 getYcoord_scale_category <- function(tr, df, yscale, yscale_mapping=NULL, ...) {

     if (is.null(yscale_mapping)) {

@@ -719,7 +758,10 @@ getYcoord_scale_category <- function(tr, df, yscale, yscale_mapping=NULL, ...) {

         stop("yscale_mapping should be a named numeric vector...")

     }

+    ## assign to parent status is more prefer...

     df <- .assign_parent_status(tr, df, yscale)

+    df <- .assign_child_status(tr, df, yscale)

+    

     yy <- df[, yscale]

     ## y <- as.numeric(factor(yy))

@@ -193,14 +193,15 @@ fortify.beast <- function(model, data,

     if (time_scale) {

         df <- scaleX_by_time(df)

     }

-    

-    return(df)

+

+    append_extraInfo(df, model)

 }

 scaleX_by_time <- function(df) {

-    time <- with(df, gsub(".*_(\\d+\\.*\\d+)$", "\\1", label[isTip])) %>% as.numeric

+    time <- with(df, gsub(".*[_/]{1}(\\d+\\.*\\d+)$", "\\1", label[isTip])) %>% as.numeric

     latest <- which.max(time)

     df$x <- df$x + time[latest] - max(df$x)

+    df$branch <- (df[df$parent, "x"] + df[, "x"])/2

     return(df)

 }

@@ -243,7 +244,8 @@ fortify.codeml <- function(model, data,

     res <- merge_phylo_anno.codeml_mlc(df, dNdS, ndigits)

     df <- merge_phylo_anno.paml_rst(res, model@rst)

-    scaleY(phylo, df, yscale, layout, ...)

+    df <- scaleY(phylo, df, yscale, layout, ...)

+    append_extraInfo(df, model)

 }

@@ -266,7 +268,8 @@ fortify.codeml_mlc <- function(model, data,

     dNdS <- model@dNdS

     df <- merge_phylo_anno.codeml_mlc(df, dNdS, ndigits)

-    scaleY(phylo, df, yscale, layout, ...)

+    df <- scaleY(phylo, df, yscale, layout, ...)

+    append_extraInfo(df, model)

 }

 merge_phylo_anno.codeml_mlc <- function(df, dNdS, ndigits = NULL) {

@@ -320,7 +323,8 @@ fortify.paml_rst <- function(model, data, layout = "phylogram", yscale="none",

                              ladderize=TRUE, right=FALSE, ...) {

     df <- fortify.phylo(model@phylo, data, layout, ladderize, right, ...)

     df <- merge_phylo_anno.paml_rst(df, model)

-    scaleY(model@phylo, df, yscale, layout, ...)

+    df <- scaleY(model@phylo, df, yscale, layout, ...)

+    append_extraInfo(df, model)

 }

 merge_phylo_anno.paml_rst <- function(df, model) {

@@ -349,7 +353,9 @@ fortify.jplace <- function(model, data,

     df <- df %add2% place

-    scaleY(model@phylo, df, yscale, layout, ...)

+    df <- scaleY(model@phylo, df, yscale, layout, ...)

+

+    append_extraInfo(df, model)    

 }

 scaleY <- function(phylo, df, yscale, layout, ...) {

@@ -1,5 +1,66 @@

+has.slot <- function(object, slotName) {

+    if (!isS4(object)) {

+        return(FALSE)

+    }

+    

+    slot <- tryCatch(slot(object, slotName), error=function(e) NULL)

+    ! is.null(slot)

+}

+

+has.extraInfo <- function(object) {

+    if (!is.tree(object)) {

+        return(FALSE)

+    }

+

+    if (! has.slot(object, "extraInfo")) {

+        return(FALSE)

+    }

+

+    extraInfo <- object@extraInfo

+

+    if (nrow(extraInfo) > 0) {

+        return(TRUE)

+    }

+

+    return(FALSE)        

+}

+

+append_extraInfo <- function(df, object) {

+    if (has.extraInfo(object)) {

+        info <- object@extraInfo

+        res <- merge(df, info, by.x=c("node", "parent"), by.y=c("node", "parent"))

+    } else {

+        res <- df

+    }

+

+    i <- order(res$node, decreasing = FALSE)

+    res <- res[i,]

+    return(res)

+}

+

+get.fields.tree <- function(object) {

+    if (is(object, "codeml")) {

+        fields <- c(get.fields(object@rst),

+                    get.fields(object@mlc))

+        fields <- unique(fields)

+    } else {

+        fields <- object@fields

+    }

+    

+    if (has.slot(object, "extraInfo")) {

+        extraInfo <- object@extraInfo

+        if (nrow(extraInfo) > 0) {

+            cn <- colnames(extraInfo)

+            i <- match(c("x", "y", "isTip", "node", "parent", "label", "branch", "branch.length"), cn)

+            i <- i[!is.na(i)]

+            fields %<>% c(cn[-i])

+        }

+    }

+    return(fields)

+}

+

 print_fields <- function(object, len=5) {

-    fields <- get.fields(object)

+    fields <- get.fields(object)    

     n <- length(fields)

     i <- floor(n/len)

     for (j in 0:i) {

@@ -60,6 +121,12 @@ get.subs_ <- function(tree, fasta, translate=TRUE, removeGap=TRUE) {

         if (is.null(res)) {

             return('')

         }

+        if (nchar(res) > 50) {

+            idx <- gregexpr("/", res)[[1]]

+            i <- idx[floor(length(idx)/2)]

+            res <- paste0(substring(res, 1, i-1), "\n", substring(res, i+1))

+        }

+        

         return(res)

     })

@@ -209,6 +276,7 @@ is.tree <- function(x) {

                         "baseml",

                         "paml_rst",

                         "baseml_mlc",

+                        "codeml",

                         "hyphy",

                         "beast")

         ) {

@@ -219,21 +287,27 @@ is.tree <- function(x) {

-color_scale <- function(c1="grey", c2="red") {

+color_scale <- function(c1="grey", c2="red", n=100) {

     pal <- colorRampPalette(c(c1, c2))

-    colors <- pal(100)

+    colors <- pal(n)

     return(colors)

 }

-getIdx <- function(v, MIN, MAX) {

+getIdx <- function(v, MIN, MAX, interval=NULL) {

+    res <- sapply(v, getIdx_internal, MIN=MIN, MAX=MAX, interval=interval)

+    attr(res, "interval") <- interval

+    return(res)

+}

+

+getIdx_internal <- function(v, MIN, MAX, interval=NULL) {

     if (is.na(v)) {