## ##' read HYPHY output
## ##'
## ##'
## ##' @title read.hyphy
## ##' @param nwk tree file in nwk format, one of hyphy output
## ##' @param ancseq ancestral sequence file in nexus format,
## ##' one of hyphy output
## ##' @param tip.fasfile tip sequence file
## ##' @return A hyphy object
## ## @importFrom Biostrings readBStringSet
## ## @importFrom Biostrings toString
## ##' @export
## ##' @author Guangchuang Yu \url{http://ygc.name}
## ##' @examples
## ##' nwk <- system.file("extdata/HYPHY", "labelledtree.tree", package="ggtree")
## ##' ancseq <- system.file("extdata/HYPHY", "ancseq.nex", package="ggtree")
## ##' read.hyphy(nwk, ancseq)
## read.hyphy <- function(nwk, ancseq, tip.fasfile=NULL) {
## anc <- scan(ancseq, what="", sep="\n", quiet=TRUE)
## end <- grep("END;", anc, ignore.case=TRUE)
## seq.start <- grep("MATRIX", anc, ignore.case=TRUE)
## seq.end <- end[end > seq.start][1]
## seq <- anc[(seq.start+1):(seq.end-1)]
## seq <- seq[seq != ";"]
## seq <- seq[seq != ""]
## seq <- gsub(" ", "", seq)
## seq <- gsub(";", "", seq)
## ## some files may only contains sequences (should have TAXALABELS block that contains seq names).
## ## some may contains sequence name like phylip format in MATRIX block (no need to have TAXALABELS block).
## ##
## ## extract sequence name if available
## if (all(grepl("\\s+", seq))) {
## ## if contains blank space, may contains seq name
## sn <- gsub("(\\w*)\\s.*", "\\1", seq)
## }
## seq <- gsub("\\w*\\s+", "", seq)
## label.start <- grep("TAXLABELS", anc, ignore.case = TRUE)
## if (length(label.start) == 0) {
## if (all(sn == "")) {
## stop("taxa labels is not available...")
## }
## label <- sn
## } else {
## label.end <- end[end > label.start][1]
## label <- anc[(label.start+1):(label.end-1)]
## label <- sub("^\t+", "", label)
## label <- sub("\\s*;$", "", label)
## label <- unlist(strsplit(label, split="\\s+"))
## label <- gsub("'|\"", "", label)
## }
## names(seq) <- label
## tr <- read.tree(nwk)
## nl <- tr$node.label
## ## root node may missing, which was supposed to be 'Node1'
## ##
## ## from a user's file, which is 'Node0', but it seems the file is not from the output of HYPHY.
## ##
## ## I am not sure. But it's safe to use "label[!label %in% nl]" instead of just assign it to "Node1".
## ##
## ## nl[nl == ""] <- "Node1"
## nl[nl == ""] <- label[!label %in% nl]
## tr$node.label <- nl
## type <- get_seqtype(seq)
## fields <- "subs"
## if (type == "NT") {
## fields <- c(fields, "AA_subs")
## }
## res <- new("hyphy",
## fields = fields,
## treetext = scan(nwk, what='', quiet=TRUE),
## phylo = tr,
## seq_type = type,
## ancseq = seq,
## tree.file = filename(nwk),
## ancseq.file = ancseq
## )
## if ( !is.null(tip.fasfile) ) {
## readBStringSet <- get_fun_from_pkg("Biostrings", "readBStringSet")
## toString <- get_fun_from_pkg("Biostrings", "toString")
## tip_seq <- readBStringSet(tip.fasfile)
## nn <- names(tip_seq)
## tip_seq <- sapply(seq_along(tip_seq), function(i) {
## toString(tip_seq[i])
## })
## names(tip_seq) <- nn
## res@tip_seq <- tip_seq
## res@tip.fasfile <- tip.fasfile
## }
## set.hyphy_(res)
## }
## ##' @rdname groupOTU-methods
## ##' @exportMethod groupOTU
## setMethod("groupOTU", signature(object="hyphy"),
## function(object, focus, group_name="group") {
## groupOTU_(object, focus, group_name)
## }
## )
## ##' @rdname groupClade-methods
## ##' @exportMethod groupClade
## setMethod("groupClade", signature(object="hyphy"),
## function(object, node, group_name="group") {
## groupClade_(object, node, group_name)
## }
## )
## ##' @rdname scale_color-methods
## ##' @exportMethod scale_color
## setMethod("scale_color", signature(object="hyphy"),
## function(object, by, ...) {
## scale_color_(object, by, ...)
## })
## ##' @rdname gzoom-methods
## ##' @exportMethod gzoom
## setMethod("gzoom", signature(object="hyphy"),
## function(object, focus, subtree=FALSE, widths=c(.3, .7)) {
## gzoom.phylo(get.tree(object), focus, subtree, widths)
## })
## ##' @rdname show-methods
## ##' @exportMethod show
## setMethod("show", signature(object = "hyphy"),
## function(object) {
## cat("'hyphy' S4 object that stored information of \n\t",
## paste0("'", object@tree.file, "'"))
## if (length(object@tip_seq) == 0) {
## cat(paste0("and '", object@ancseq.file, "'"), ".\n")
## } else {
## cat(paste0(", \n\t'", object@ancseq.file, "'"),
## paste0("and \n\t'", object@tip.fasfile, "'."),
## "\n\n")
## }
## cat("...@ tree:")
## print.phylo(get.tree(object))
## cat("\nwith the following features available:\n")
## cat("\t", paste0("'",
## paste(get.fields(object), collapse="',\t'"),
## "'."),
## "\n")
## })
## ##' @rdname get.tree-methods
## ##' @exportMethod get.tree
## ##' @examples
## ##' nwk <- system.file("extdata/HYPHY", "labelledtree.tree", package="ggtree")
## ##' ancseq <- system.file("extdata/HYPHY", "ancseq.nex", package="ggtree")
## ##' hy <- read.hyphy(nwk, ancseq)
## ##' get.tree(hy)
## setMethod("get.tree", signature(object = "hyphy"),
## function(object) {
## object@phylo
## }
## )
## ##' @rdname get.fields-methods
## ##' @exportMethod get.fields
## setMethod("get.fields", signature(object = "hyphy"),
## function(object, ...) {
## if(length(object@tip_seq) == 0) {
## warning("tip sequence not available...\n")
## } else {
## get.fields.tree(object)
## }
## })
## ##' @rdname get.subs-methods
## ##' @exportMethod get.subs
## ##' @examples
## ##' nwk <- system.file("extdata/HYPHY", "labelledtree.tree", package="ggtree")
## ##' ancseq <- system.file("extdata/HYPHY", "ancseq.nex", package="ggtree")
## ##' tipfas <- system.file("extdata", "pa.fas", package="ggtree")
## ##' hy <- read.hyphy(nwk, ancseq, tipfas)
## ##' get.subs(hy, type="AA_subs")
## setMethod("get.subs", signature(object="hyphy"),
## function(object, type, ...) {
## if (length(object@tip_seq) == 0) {
## stop("tip sequence not available...\n")
## }
## if (type == "subs") {
## return(object@subs)
## } else {
## return(object@AA_subs)
## }
## })
## set.hyphy_ <- function(object, ...) {
## if (!is(object, "hyphy")) {
## stop("object should be an instance of 'hyphy'")
## }
## if (length(object@tip_seq) == 0) {
## return(object)
## }
## types <- get.fields(object)
## seqs <- c(object@tip_seq, object@ancseq)
## for (type in types) {
## if (type == "subs") {
## translate <- FALSE
## } else {
## translate <- TRUE
## }
## subs <- get.subs_(object@phylo, seqs, translate, ...)
## if (type == "subs") {
## object@subs <- subs
## } else {
## object@AA_subs <- subs
## }
## }
## return(object)
## }
