@@ -18,10 +18,8 @@ Imports:

     ape,

     grDevices,

     grid,

-    jsonlite,

     magrittr,

     methods,

-    stats4,

     tidyr,

     utils

 Suggests:

@@ -8,9 +8,7 @@ S3method(fortify,beast)

 S3method(fortify,codeml)

 S3method(fortify,codeml_mlc)

 S3method(fortify,hyphy)

-S3method(fortify,jplace)

 S3method(fortify,multiPhylo)

-S3method(fortify,nhx)

 S3method(fortify,obkData)

 S3method(fortify,paml_rst)

 S3method(fortify,phangorn)

@@ -73,8 +71,6 @@ export(ggtree)

 export(gheatmap)

 export(gzoom)

 export(inset)

-export(mask)

-export(merge_tree)

 export(msaplot)

 export(multiplot)

 export(nodebar)

@@ -82,9 +78,6 @@ export(nodeid)

 export(nodepie)

 export(open_tree)

 export(phylopic)

-export(plot)

-export(raxml2nwk)

-export(read.phyloT)

 export(reroot)

 export(rescale_tree)

 export(revts)

@@ -105,7 +98,6 @@ export(xlim_tree)

 exportMethods(groupClade)

 exportMethods(groupOTU)

 exportMethods(gzoom)

-exportMethods(plot)

 exportMethods(reroot)

 importFrom(ape,di2multi)

 importFrom(ape,extract.clade)

@@ -183,7 +175,8 @@ importFrom(magrittr,"%<>%")

 importFrom(magrittr,"%>%")

 importFrom(magrittr,add)

 importFrom(magrittr,equals)

-importFrom(stats4,plot)

+importFrom(methods,is)

+importFrom(methods,missingArg)

 importFrom(tidyr,gather)

 importFrom(treeio,as.treedata)

 importFrom(treeio,groupClade)

@@ -1,6 +1,6 @@

 ## setOldClass("phylo")

 ## setOldClass("multiPhylo")

-## setOldClass("ggtree")

+setOldClass("ggtree")

 ## setClassUnion("phyloOrmultiPhylo", c("phylo", "multiPhylo"))

@@ -7,20 +7,20 @@ as.binary <- function(tree, ...) {

     UseMethod("as.binary")

 }

-##' plot method generics

-##'

-##'

-##' @docType methods

-##' @name plot

-##' @rdname plot-methods

-##' @title plot method

-##' @param x object

-##' @param ... Additional argument list

-##' @return plot

-##' @importFrom stats4 plot

-##' @export

-if ( !isGeneric("plot") )

-	setGeneric("plot", function(x, ...) standardGeneric("plot"))

+## ##' plot method generics

+## ##'

+## ##'

+## ##' @docType methods

+## ##' @name plot

+## ##' @rdname plot-methods

+## ##' @title plot method

+## ##' @param x object

+## ##' @param ... Additional argument list

+## ##' @return plot

+## ##' @importFrom stats4 plot

+## ##' @export

+## if ( !isGeneric("plot") )

+## 	setGeneric("plot", function(x, ...) standardGeneric("plot"))

 ##' @docType methods

 ##' @name reroot

@@ -1,6 +1,6 @@

 ##' add insets in a tree

 ##'

-##' 

+##'

 ##' @title inset

 ##' @param tree_view tree view

 ##' @param insets a list of ggplot objects, named by node number

@@ -22,7 +22,7 @@ inset <- function(tree_view, insets, width=0.1, height=0.1, hjust=0, vjust=0, x=

         xx <- df$branch

     }

     yy <- df$y

-    

+

     xx <- xx - hjust

     yy <- yy - vjust

@@ -38,7 +38,7 @@ inset <- function(tree_view, insets, width=0.1, height=0.1, hjust=0, vjust=0, x=

 ##' generate a list of bar charts for results of ancestral state reconstruction

 ##'

-##' 

+##'

 ##' @title nodebar

 ##' @param position position of bar, one of 'stack' and 'dodge'

 ##' @inheritParams nodepie

@@ -52,7 +52,7 @@ nodebar <- function(data, cols, color, alpha=1, position="stack") {

         stop("data should have a column 'node'...")

     }

     type <- value <- NULL

-    

+

     ldf <- gather(data, type, value, cols) %>% split(., .$node)

     bars <- lapply(ldf, function(df) ggplot(df, aes_(x=1, y=~value, fill=~type)) +

                                      geom_bar(stat='identity', alpha=alpha, position=position) +

@@ -69,7 +69,7 @@ nodebar <- function(data, cols, color, alpha=1, position="stack") {

 ##' generate a list of pie charts for results of ancestral stat reconstruction

 ##'

-##' 

+##'

 ##' @title nodepie

 ##' @param data a data.frame of stats with an additional column of node number

 ##' @param cols column of stats

@@ -91,11 +91,12 @@ nodepie <- function(data, cols, color, alpha=1) {

 }

+##' @importFrom methods missingArg

 ggpie <- function(data, y, fill, color, alpha=1) {

     p <- ggplot(data, aes_(x=1, y=y, fill=fill)) +

         geom_bar(stat='identity', alpha=alpha) +

         coord_polar(theta='y') + theme_inset()

-    

+

     if (missingArg(color) || is.null(color) || is.na(color)) {

         ## do nothing

     } else {

deleted file mode 100644

@@ -1,122 +0,0 @@

-##' merge two tree object

-##'

-##'

-##' @title merge_tree

-##' @param obj1 tree object 1

-##' @param obj2 tree object 2

-##' @return tree object

-##' @importFrom magrittr %<>%

-##' @export

-##' @author Guangchuang Yu

-merge_tree <- function(obj1, obj2) {

-    ##

-    ## INFO:

-    ## ape::all.equal.phylo can be used to test equal phylo topology.

-    ##

-

-    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", "angle"), 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)

-}

@@ -52,6 +52,7 @@ as.binary.phylo <- function(tree, ...) {

 ##' @importFrom magrittr %<>%

 ##' @importFrom magrittr add

 ##' @importFrom ape write.tree

+##' @importFrom ape read.tree

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

 rm.singleton.newick <- function(nwk, outfile = NULL) {

     tree <- readLines(nwk)

@@ -341,21 +342,21 @@ fortify.phangorn <- fortify.paml_rst

 fortify.hyphy <- fortify.paml_rst

-##' @method fortify jplace

-##' @importFrom ape read.tree

-##' @export

-fortify.jplace <- function(model, data,

-                           layout="rectangular", yscale="none",

-                           ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {

-    df <- get.treeinfo(model, layout, ladderize, right, mrsd=mrsd, ...)

-    place <- get.placements(model, by="best")

+## ##' @method fortify jplace

+## ##' @importFrom ape read.tree

+## ##' @export

+## fortify.jplace <- function(model, data,

+##                            layout="rectangular", yscale="none",

+##                            ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {

+##     df <- get.treeinfo(model, layout, ladderize, right, mrsd=mrsd, ...)

+##     place <- get.placements(model, by="best")

-    df <- df %add2% place

+##     df <- df %add2% place

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

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

-    append_extraInfo(df, model)

-}

+##     append_extraInfo(df, model)

+## }

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

     if (yscale == "none") {

@@ -568,14 +569,14 @@ as.data.frame.phylo_ <- function(x, layout="rectangular",

     return(res)

 }

-##' @method fortify nhx

-##' @export

-fortify.nhx <- function(model, data, layout= "rectangular",

-                        ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {

-    df <- fortify(get.tree(model), layout=layout, ladderize=ladderize, right=right, mrsd=mrsd, ...)

-    df <- merge(df, model@nhx_tags, by.x="node", by.y="node", all.x=TRUE)

-    append_extraInfo(df, model)

-}

+## ##' @method fortify nhx

+## ##' @export

+## fortify.nhx <- function(model, data, layout= "rectangular",

+##                         ladderize=TRUE, right=FALSE, mrsd=NULL, ...) {

+##     df <- fortify(get.tree(model), layout=layout, ladderize=ladderize, right=right, mrsd=mrsd, ...)

+##     df <- merge(df, model@nhx_tags, by.x="node", by.y="node", all.x=TRUE)

+##     append_extraInfo(df, model)

+## }

 ##' @method fortify raxml

@@ -13,11 +13,18 @@

 ##           }

 ##           )

+##' groupClade method for ggtree object

+##'

+##'

 ##' @name groupClade

 ##' @title groupClade method

 ##' @rdname groupClade-methods

+##' @param object ggtree object

+##' @param node internal node number

+##' @param group_name name of the group

 ##' @importFrom treeio groupClade

 ##' @exportMethod groupClade

+##' @aliases groupClade,ggtree-method

 setMethod("groupClade", signature(object="ggtree"),

           function(object, node, group_name) {

               groupClade.ggtree(object, node, group_name)

@@ -33,11 +33,19 @@

 ##           )

+##' groupOTU method for ggtree object

+##'

+##'

 ##' @name groupOTU

 ##' @title groupOTU method

 ##' @rdname groupOTU-methods

+##' @param object ggtree object

+##' @param focus OTU to focus

+##' @param group_name name of the group

+##' @param ... additional parameters

 ##' @importFrom treeio groupOTU

 ##' @exportMethod groupOTU

+##' @aliases groupOTU,ggtree-method

 setMethod("groupOTU", signature(object="ggtree"),

           function(object, focus, group_name="group", ...) {

               groupOTU.ggtree(object, focus, group_name, ...)

@@ -172,14 +180,14 @@ setMethod("groupOTU", signature(object="ggtree"),

 ##     return(phy)

 ## }

-groupOTU_ <- function(object, focus, group_name, ...) {

-    if (is(object, "phylo")) {

-        object <- groupOTU.phylo(object, focus, group_name, ...)

-    } else {

-        object@phylo <- groupOTU.phylo(get.tree(object), focus, group_name, ...)

-    }

-    return(object)

-}

+## groupOTU_ <- function(object, focus, group_name, ...) {

+##     if (is(object, "phylo")) {

+##         object <- groupOTU.phylo(object, focus, group_name, ...)

+##     } else {

+##         object@phylo <- groupOTU.phylo(get.tree(object), focus, group_name, ...)

+##     }

+##     return(object)

+## }

 groupOTU.ggtree <- function(object, focus, group_name, ...) {

@@ -53,6 +53,7 @@ gzoom.ggtree <- function(tree_view, focus, widths=c(.3, .7), xmax_adjust=0) {

 ##' @rdname gzoom-methods

 ##' @exportMethod gzoom

 ##' @param xmax_adjust adjust xmax (xlim[2])

+##' @aliases gzoom,ggtree-method

 setMethod("gzoom", signature(object="ggtree"),

           function(object, focus, widths=c(.3, .7), xmax_adjust=0) {

               gzoom.ggtree(object, focus, widths, xmax_adjust)

@@ -96,6 +96,7 @@

 ##' @param p tree view

 ##' @param data data.frame

 ##' @return updated data.frame

+##' @importFrom methods is

 ##' @export

 ##' @author Guangchuang Yu

 `%+>%` <- function(p, data) {

@@ -1,394 +1,394 @@

-##' site mask

-##'

-##'

-##' @title mask

-##' @param tree_object tree object

-##' @param field selected field

-##' @param site site

-##' @param mask_site if TRUE, site will be masked.

-##'                  if FALSE, selected site will not be masked, while other sites will be masked.

-##' @return updated tree object

-##' @export

-##' @author Guangchuang Yu

-mask <- function(tree_object, field, site, mask_site=FALSE) {

-    has_field <- has.field(tree_object, field)

-    if (has_field == FALSE) {

-        stop("'field' is not available in 'tree_object'...")

-    }

-

-    has_slot <- attr(has_field, "has_slot")

-    is_codeml <- attr(has_field, "is_codeml")

-

-    if (has_slot) {

-        if (is_codeml) {

-            field_info <- slot(tree_object@rst, field)

-        } else {

-            field_info <- slot(tree_object, field)

-        }

-        field_data <- field_info[,2]

-    } else {

-        field_data <- tree_object@extraInfo[, field]

-    }

-

-    field_data <- sapply(field_data, gsub, pattern="\n", replacement="/")

-

-    x <- field_data[field_data != ""]

-    x <- x[!is.na(x)]

-    pos <- strsplit(x, " / ") %>% unlist %>%

-        gsub("^[a-zA-Z]+", "", . ) %>%

-        gsub("[a-zA-Z]\\s*$", "", .) %>%

-        as.numeric

-

-    if (mask_site == FALSE) {

-        pos2 <- 1:max(pos)

-        pos2 <- pos2[-site]

-        site <- pos2

-    }

-

-    site <- site[site %in% pos]

-

-    for (i in seq_along(field_data)) {

-        if (is.na(field_data[i]))

-            next

-        for (j in seq_along(site)) {

-            pattern <- paste0("/*\\s*[a-zA-Z]", site[j], "[a-zA-Z]\\s*")

-            field_data[i] <- gsub(pattern, "",  field_data[i])

-        }

-        field_data[i] <- gsub("^/\\s", "", field_data[i]) %>% .add_new_line

-    }

-

-    if (has_slot) {

-        field_info[,2] <- field_data

-        if (is_codeml) {

-            slot(tree_object@rst, field) <- field_info

-        } else {

-            slot(tree_object, field) <- field_info

-        }

-    } else {

-        tree_object@extraInfo[, field] <- field_data

-    }

-    tree_object

-}

-

-

-read.tip_seq_mlc <- function(mlcfile) {

-    info <- getPhyInfo(mlcfile)

-    mlc <- readLines(mlcfile)

-    ## remove blank lines

-    blk <- grep("^\\s*$", mlc)

-    if (length(blk) > 0) {

-        mlc <- mlc[-blk]

-    }

-

-    seqs <- mlc[2:(info$num+1)]

-    seqs <- gsub("\\s+", "", seqs)

-    wd <- info$width

-    res <- sapply(seqs, function(x) substring(x, (base::nchar(x) - wd + 1), base::nchar(x)))

-    nn <- sapply(seqs, function(x) substring(x, 1, (base::nchar(x) - wd)))

-    names(res) <- nn

-    return(res)

-}

-

-read.tip_seq_mlb <- read.tip_seq_mlc

-

-read.dnds_mlc <- function(mlcfile) {

-    mlc <- readLines(mlcfile)

-    i <- grep("dN & dS for each branch", mlc)

-    j <- grep("tree length for dN", mlc)

-

-    mlc <- mlc[i:j]

-    hi <- grep("dN/dS", mlc)

-    cn <- strsplit(mlc[hi], " ") %>% unlist %>% `[`(nzchar(.))

-

-    ii <- grep("\\d+\\.\\.\\d+", mlc)

-    info <- mlc[ii]

-    info %<>% sub("^\\s+", "", .)

-    info %<>% sub("\\s+$", "", .)

-    res <- t(sapply(info, function(x) {

-        y <- unlist(strsplit(x, "\\s+"))

-        edge <- unlist(strsplit(y[1], "\\.\\."))

-        yy <- c(edge, y[-1])

-        as.numeric(yy)

-    }))

-

-    row.names(res) <- NULL

-    colnames(res) <- c("parent", "node", cn[-1])

-    colnames(res) <- gsub("\\*", "_x_", colnames(res))

-    colnames(res) <- gsub("\\/", "_vs_", colnames(res))

-    return(res)

-}

-

-read.treetext_paml_mlc <- function(mlcfile) {

-    read.treetext_paml(mlcfile, "mlc")

-}

-

-read.treetext_paml_rst<- function(rstfile) {

-    read.treetext_paml(rstfile, "rst")

-}

-

-read.treetext_paml <- function(file, by) {

-    ## works fine with baseml and codeml

-    x <- readLines(file)

-    tr.idx <- get_tree_index_paml(x)

-

-    if ( by == "rst" ) {

-        ii <- 1

-    } else if (by == "mlc") {

-        ii <- 3

-    } else {

-        stop("_by_ should be one of 'rst' or 'mlc'")

-    }

-

-    return(x[tr.idx][ii])

-}

-

-read.phylo_paml_mlc <- function(mlcfile) {

-    parent <- node <- label <- NULL

-

-    mlc <- readLines(mlcfile)

-    edge <- get_tree_edge_paml(mlc)

-

-    tr.idx <- get_tree_index_paml(mlc)

-    tr2 <- read.tree(text=mlc[tr.idx[2]])

-    tr3 <- read.tree(text=mlc[tr.idx[3]])

-

-    treeinfo <- as.data.frame(edge)

-    colnames(treeinfo) <- c("parent", "node")

-    treeinfo$length <- NA

-    treeinfo$label <- NA

-    treeinfo$isTip <- FALSE

-    ntip <- Ntip(tr3)

-    ii <- match(tr2$tip.label, treeinfo[, "node"])

-    treeinfo[ii, "label"] <- tr3$tip.label

-    treeinfo[ii, "isTip"] <- TRUE

-    ## jj <- match(1:ntip, tr3$edge[,2])

-    treeinfo[ii, "length"] <- tr2$edge.length[ii] ##tr3$edge.length[jj]

-

-    root <- getRoot(tr3) ## always == (Ntip(tr3) + 1)

-    currentNode <- treeinfo$label[ii]

-    treeinfo.tr3 <- fortify(tr3)

-

-    treeinfo$visited <- FALSE

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

-        pNode <- c()

-        for( kk in currentNode ) {

-            i <- which(treeinfo$label == kk)

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

-            j <- which(treeinfo.tr3$label == kk)

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

-            if (ip != root) {

-                ii <- which(treeinfo[, "node"] == ip)

-                if (treeinfo$visited[ii] == FALSE) {

-                    jp <- treeinfo.tr3[j, "parent"]

-                    jj <- which(treeinfo.tr3[, "node"] == jp)

-                    treeinfo[ii, "label"] <- as.character(ip)

-                    treeinfo.tr3[jj, "label"] <- as.character(ip)

-                    treeinfo[ii, "length"] <- treeinfo.tr3[jj, "branch.length"]

-                    pNode <- c(pNode, ip)

-                }

-                treeinfo[ii, "visited"] <- TRUE

-            }

-

-        }

-        currentNode <- unique(pNode)

-    }

-

-    phylo <- with(treeinfo,

-                  list(edge= cbind(as.numeric(parent),

-                           as.numeric(node)),

-                       edge.length = length,

-                       tip.label = label[order(node)][1:ntip],

-                       Nnode = tr3$Nnode,

-                       node.label = c(root, label[order(node)][-c(1:ntip)])

-                       )

-                  )

-    class(phylo) <- "phylo"

-    phylo <- reorder.phylo(phylo, "cladewise")

-    return(phylo)

-}

-

-##' @importFrom ape reorder.phylo

-read.phylo_paml_rst <- function(rstfile) {

-    parent <- node <- label <- NULL

-

-    ## works fine with baseml and codeml

-    rst <- readLines(rstfile)

-    tr.idx <- get_tree_index_paml(rst)

-

-    tr1 <- read.tree(text=rst[tr.idx][1])

-    tr3 <- read.tree(text=rst[tr.idx][3])

-

-    edge <- get_tree_edge_paml(rst)

-

-    label=c(tr3$tip.label, tr3$node.label)

-    root <- getRoot(tr3)

-    label %<>% `[`(. != root)

-

-    node.length <- data.frame(label=label,

-                              length=tr1$edge.length)

-

-    ## node.length$node <- sub("_\\w+", "", node.length$label

-    node.length$node <- gsub("^(\\d+)_.*", "\\1", node.length$label)

-    node.length$label %<>% sub("\\d+_", "", .)

-

-    edge <- as.data.frame(edge)

-    colnames(edge) <- c("parent", "node")

-

-    treeinfo <- merge(edge, node.length, by.x="node", by.y="node")

-    edge2 <- treeinfo[, c("parent", "node")]

-    edge2 %<>% as.matrix

-

-    ntip <- Ntip(tr3)

-

-    phylo <- with(treeinfo,

-                  list(edge= cbind(as.numeric(parent),

-                      as.numeric(node)),

-                       edge.length = length,

-                       tip.label = label[order(node)][1:ntip],

-                       Nnode = tr3$Nnode,

-                       node.label = c(root, label[order(node)][-c(1:ntip)])

-                       )

-                  )

-

-    class(phylo) <- "phylo"

-    phylo <- reorder.phylo(phylo, "cladewise")

-

-    return(phylo)

-}

-

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

-    ## works fine with baseml and codeml

-    rst <- readLines(rstfile)

-

-    by <- match.arg(by, c("Marginal", "Joint"))

-    query <- paste(by, "reconstruction of ancestral sequences")

-    idx <- grep(query, rst)

-    if(length(idx) == 0) {

-        ## in some paml setting, joint_ancseq are not available.

-        return("")

-    }

-    si <- grep("reconstructed sequences", rst)

-    idx <- si[which.min(abs(si-idx))]

-

-    nl <- strsplit(rst[idx+2], split=" ") %>% unlist %<>% `[`(nzchar(.))

-    N <- as.numeric(nl[1])

-    seq.leng <- as.numeric(nl[2])

-

-    seqs <- rst[(idx+4):(idx+3+N)]

-

-    seq.name <- character(N)

-    res <- character(N)

-    for (i in 1:N) {

-        ss <- gsub(" ", "", seqs[i])

-        nn <- base::nchar(ss)

-        res[i] <- substring(ss, nn-seq.leng+1,nn)

-        seq.name[i] <- substring(ss, 1, nn-seq.leng)

-    }

-    seq.name <- sub("\\w+#", "", seq.name)

-    names(res) <- seq.name

-

-    return(res)

-}

-

-

-get_tree_index_paml <- function(paml) {

-    grep("\\)[ \\.0-9]*;", paml)

-}

-

-get_tree_edge_index_paml <- function(paml) {

-    grep("\\d+\\.\\.\\d+", paml)

-}

-

-get_tree_edge_paml <- function(paml) {

-    tr.idx <- get_tree_index_paml(paml)

-

-    edge.idx <- get_tree_edge_index_paml(paml)

-    edge.idx <- edge.idx[edge.idx < tr.idx[3]]

-

-    nodeNum <- strsplit(paml[edge.idx], split="\\.\\.") %>%

-        unlist %>% strsplit(split="[[:space:]]") %>% unlist

-

-    nodeNum %<>% `[`(nzchar(.))

-

-    edge <- matrix(as.numeric(nodeNum), ncol=2, byrow = TRUE)

-

-    return(edge)

-}

-

-set.paml_rst_ <- function(object) {

-    if (!is(object, "paml_rst")) {

-        stop("object should be an instance of 'paml_rst'")

-    }

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

-        return(object)

-    }

-

-    types <- get.fields(object)

-    for (type in types) {

-        value <- subs_paml_rst(object, type)

-        if (all(is.na(value)))

-            next

-

-        if (type == "marginal_subs") {

-            object@marginal_subs <- value

-        } else if (type == "marginal_AA_subs") {

-            object@marginal_AA_subs <- value

-        } else if (type == "joint_subs") {

-            object@joint_subs <- value

-        } else if (type == "joint_AA_subs") {

-            object@joint_AA_subs <- value

-        }

-    }

-    return(object)

-}

-

-

-get.subs_paml_rst <- function(object, type) {

-    if (!is(object, "paml_rst")) {

-        stop("object should be an instance of 'paml_rst'")

-    }

-    if (type == "marginal_subs") {

-        res <- object@marginal_subs

-    } else if (type == "marginal_AA_subs") {

-        res <- object@marginal_AA_subs

-    } else if (type == "joint_subs") {

-        res <- object@joint_subs

-    } else if (type == "joint_AA_subs") {

-        res <- object@joint_AA_subs

-    } else {

-        stop("type should be one of 'marginal_subs',

-                             'marginal_AA_subs', 'joint_subs' or 'joint_AA_subs'. ")

-    }

-    return(res)

-}

-

-subs_paml_rst <- function(x, type, ...) {

-    if (class(x) != "paml_rst") {

-        stop("x should be an object of paml_rst...")

-    }

-    seqs <- x@tip_seq

-    if (length(seqs) == 0) {

-        stop("tip sequences is not available...")

-    }

-    if (type %in% c("marginal_subs", "marginal_AA_subs")) {

-        ancseq <- x@marginal_ancseq

-        ## seqs <- c(seqs, x@marginal_ancseq)

-    } else if (type %in% c("joint_subs", "joint_AA_subs")){

-        ancseq <- x@joint_ancseq

-        ## seqs <- c(seqs, x@joint_ancseq)

-    } else {

-        stop("type should be one of 'marginal_subs',

-                             'marginal_AA_subs', 'joint_subs' or 'joint_AA_subs'. ")

-    }

-    if( type %in% c("marginal_subs", "joint_subs")) {

-        translate <- FALSE

-    } else {

-        translate <- TRUE

-    }

-

-    if (all(ancseq == "")) {

-        return(NA)

-    }

-    seqs <- c(seqs, ancseq)

-    get.subs_(x@phylo, seqs, translate=translate, ...)

-}

+## ##' site mask

+## ##'

+## ##'

+## ##' @title mask

+## ##' @param tree_object tree object

+## ##' @param field selected field

+## ##' @param site site

+## ##' @param mask_site if TRUE, site will be masked.

+## ##'                  if FALSE, selected site will not be masked, while other sites will be masked.

+## ##' @return updated tree object

+## ##' @export

+## ##' @author Guangchuang Yu

+## mask <- function(tree_object, field, site, mask_site=FALSE) {

+##     has_field <- has.field(tree_object, field)

+##     if (has_field == FALSE) {

+##         stop("'field' is not available in 'tree_object'...")

+##     }

+

+##     has_slot <- attr(has_field, "has_slot")

+##     is_codeml <- attr(has_field, "is_codeml")

+

+##     if (has_slot) {

+##         if (is_codeml) {

+##             field_info <- slot(tree_object@rst, field)

+##         } else {

+##             field_info <- slot(tree_object, field)

+##         }

+##         field_data <- field_info[,2]

+##     } else {

+##         field_data <- tree_object@extraInfo[, field]

+##     }

+

+##     field_data <- sapply(field_data, gsub, pattern="\n", replacement="/")

+

+##     x <- field_data[field_data != ""]

+##     x <- x[!is.na(x)]

+##     pos <- strsplit(x, " / ") %>% unlist %>%

+##         gsub("^[a-zA-Z]+", "", . ) %>%

+##         gsub("[a-zA-Z]\\s*$", "", .) %>%

+##         as.numeric

+

+##     if (mask_site == FALSE) {

+##         pos2 <- 1:max(pos)

+##         pos2 <- pos2[-site]

+##         site <- pos2

+##     }