Browse code

support MCC trees

git-svn-id: https://hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/ggtree@117329 bc3139a8-67e5-0310-9ffc-ced21a209358

g.yu authored on 12/05/2016 07:09:03
Showing 11 changed files

... ...
@@ -20,6 +20,7 @@ S3method(fortify,phylo4d)
20 20
 S3method(fortify,phyloseq)
21 21
 S3method(fortify,r8s)
22 22
 S3method(fortify,raxml)
23
+S3method(print,beastList)
23 24
 export("%<%")
24 25
 export("%<+%")
25 26
 export("%>%")
... ...
@@ -1,3 +1,7 @@
1
+CHANGES IN VERSION 1.5.3
2
+------------------------
3
+ o support reading BEAST MCC trees (multiple trees in one file) via the read.beast function <2016-05-12, Thu>
4
+ 
1 5
 CHANGES IN VERSION 1.5.2
2 6
 ------------------------
3 7
  o add multiplot in ggtreeUtilities vignette <2016-05-12, Thu>
... ...
@@ -12,32 +12,53 @@
12 12
 ##' file <- system.file("extdata/BEAST", "beast_mcc.tree", package="ggtree")
13 13
 ##' read.beast(file)
14 14
 read.beast <- function(file) {
15
+    translation <- read.trans_beast(file)
16
+    treetext <- read.treetext_beast(file)
15 17
     stats <- read.stats_beast(file)
18
+    phylo <- read.nexus(file)
19
+    
20
+    if (length(treetext) == 1) {
21
+        obj <- BEAST(file, treetext, translation, stats, phylo)
22
+    } else {
23
+        obj <- lapply(seq_along(treetext), function(i) {
24
+            BEAST(file, treetext[i], translation, stats[[i]], phylo[[i]])
25
+        })
26
+        class(obj) <- "beastList"
27
+    }
28
+    return(obj)
29
+}
30
+
31
+
32
+BEAST <- function(file, treetext, translation, stats, phylo) {
16 33
     stats$node %<>% gsub("\"*'*", "", .)
17 34
     
18 35
     fields <- sub("_lower|_upper", "", names(stats)) %>% unique
19 36
     fields %<>% `[`(.!="node")
37
+        
38
+    phylo <- remove_quote_in_tree_label(phylo)
39
+        
40
+    obj <- new("beast",
41
+               fields      = fields,
42
+               treetext    = treetext,
43
+               phylo       = phylo,
44
+               translation = translation,
45
+               stats       = stats,
46
+               file        = filename(file)
47
+               )
48
+    return(obj)
49
+}
20 50
 
21
-    phylo <- read.nexus(file)
51
+remove_quote_in_tree_label <- function(phylo) {
22 52
     if (!is.null(phylo$node.label)) {
23 53
         phylo$node.label %<>% gsub("\"*'*", "", .)
24 54
     }
25 55
     if ( !is.null(phylo$tip.label)) {
26 56
         phylo$tip.label %<>% gsub("\"*'*", "", .) 
27 57
     }
28
-    
29
-    new("beast",
30
-        fields      = fields,
31
-        treetext    = read.treetext_beast(file),
32
-        phylo       = phylo,
33
-        translation = read.trans_beast(file),
34
-        stats       = stats,
35
-        file        = filename(file)
36
-        )
58
+    return(phylo)
37 59
 }
38 60
 
39 61
 
40
-
41 62
 ##' @rdname get.fields-methods
42 63
 ##' @exportMethod get.fields
43 64
 setMethod("get.fields", signature(object="beast"),
... ...
@@ -49,24 +70,33 @@ setMethod("get.fields", signature(object="beast"),
49 70
 
50 71
 read.treetext_beast <- function(file) {
51 72
     beast <- readLines(file)
52
-    ii <- grep("tree TREE1\\s+=", beast)
73
+    ## ii <- grep("^tree TREE1\\s+=", beast)
74
+    ii <- grep("^tree ", beast)
53 75
     if (length(ii) == 0) {
54
-        ii <- grep("begin trees;", beast)
76
+        ii <- grep("[Bb]egin trees;", beast)
55 77
     }
56 78
     
57 79
     jj <- grep("[Ee]nd;", beast)
58
-    jj <- jj[jj > ii][1]
59
-    tree <- beast[ii:(jj-1)]
60
-    if (length(tree) > 1) {
61
-        tree <- paste0(tree)
62
-    }
63
-    ## tree %<>% sub("tree TREE1\\s+=\\s+\\[&R\\]\\s+", "", .)
64
-    ## tree %<>% sub("[^(]*", "", .)
65
-    tree %<>% sub("[^=]+=", "", .) %>%
66
-        sub("\\s+\\[&R\\]\\s+", "", .) %>%
67
-            sub("[^(]*", "", .)
80
+    jj <- jj[jj > max(ii)][1]
81
+
82
+    ## tree <- beast[ii:(jj-1)]
83
+    ## if (length(tree) > 1) {
84
+    ##     tree <- paste0(tree)
85
+    ## }
86
+    ## tree %<>% sub("[^=]+=", "", .) %>%
87
+    ##     sub("\\s+\\[&R\\]\\s+", "", .) %>%
88
+    ##         sub("[^(]*", "", .)
89
+    
90
+    jj <- c(ii[-1], jj)
91
+    trees <- sapply(seq_along(ii), function(i) {
92
+        tree <- beast[ii[i]:(jj[i]-1)]
93
+        if (length(tree) > 1) {
94
+            tree <- paste0(tree)
95
+        }
96
+        sub("[^(]*", "", tree)
97
+    })
68 98
 
69
-    return(tree)
99
+    return(trees)
70 100
 }
71 101
 
72 102
 read.trans_beast <- function(file) {
... ...
@@ -89,20 +119,22 @@ read.trans_beast <- function(file) {
89 119
     return(trans)
90 120
 }
91 121
 
122
+
92 123
 read.stats_beast <- function(file) {
93 124
     beast <- readLines(file)
94
-    tree <- read.treetext_beast(file)
125
+    trees <- read.treetext_beast(file)
126
+    if (length(trees) == 1) {
127
+        return(read.stats_beast_internal(beast, trees))
128
+    }
129
+    lapply(trees, read.stats_beast_internal, beast=beast)
130
+}
95 131
 
132
+read.stats_beast_internal <- function(beast, tree) {
96 133
     tree2 <- gsub("\\[[^\\[]*\\]", "", tree)
97 134
     phylo <- read.tree(text = tree2)
98
-    if(is.null(phylo$node.label)) {
99
-        nnode <- phylo$Nnode
100
-        nlab <- paste("X", 1:nnode, sep="")
101
-        for (i in 1:nnode) {
102
-            tree2 <- sub("\\)([:;])", paste0("\\)", nlab[i], "\\1"), tree2)
103
-        }
104
-    }
105 135
 
136
+    tree2 <- add_pseudo_nodelabel(phylo, tree2)
137
+    
106 138
     ## node name corresponding to stats
107 139
     nn <- strsplit(tree2, split=",") %>% unlist %>%
108 140
         strsplit(., split="\\)") %>% unlist %>%
... ...
@@ -244,6 +276,18 @@ read.stats_beast <- function(file) {
244 276
     return(stats3)
245 277
 }
246 278
 
279
+add_pseudo_nodelabel <- function(phylo, treetext) {
280
+    if(is.null(phylo$node.label)) {
281
+        nnode <- phylo$Nnode
282
+        nlab <- paste("X", 1:nnode, sep="")
283
+        for (i in 1:nnode) {
284
+            treetext <- sub("\\)([:;])", paste0("\\)", nlab[i], "\\1"), treetext)
285
+        }
286
+    }
287
+
288
+   return(treetext)
289
+}
290
+
247 291
 
248 292
 ##' @rdname reroot-methods
249 293
 ##' @exportMethod reroot
... ...
@@ -65,7 +65,7 @@ geom_hilight <- function(node, fill="steelblue", alpha=.5, extend=0, extendto=NU
65 65
 stat_hilight <- function(mapping=NULL, data=NULL, geom="rect",
66 66
                          position="identity",  node, 
67 67
                          show.legend=NA, inherit.aes=FALSE,
68
-                        fill, alpha, extend=0, xmax=NULL,
68
+                         fill, alpha, extend=0, extendto=NULL,
69 69
                          ...) {
70 70
     default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent, branch.length=~branch.length)
71 71
     if (is.null(mapping)) {
... ...
@@ -83,8 +83,10 @@ stat_hilight <- function(mapping=NULL, data=NULL, geom="rect",
83 83
         show.legend=show.legend,
84 84
         inherit.aes = inherit.aes,
85 85
         params = list(node=node,
86
-                      fill=fill, alpha=alpha,
87
-                      extend=extend, extendto=extendto,
86
+                      fill=fill,
87
+                      alpha=alpha,
88
+                      extend=extend,
89
+                      extendto=extendto,
88 90
                       ...)
89 91
         )
90 92
 }
... ...
@@ -24,7 +24,9 @@ read.hyphy <- function(nwk, ancseq, tip.fasfile=NULL) {
24 24
     seq       <- anc[(seq.start+1):(seq.end-1)]
25 25
     seq       <- seq[seq != ";"]
26 26
     seq       <- seq[seq != ""]
27
-
27
+    seq       <- gsub(" ", "", seq)
28
+    seq       <- gsub(";", "", seq)
29
+    
28 30
     ## some files may only contains sequences (should have TAXALABELS block that contains seq names).
29 31
     ## some may contains sequence name like phylip format in MATRIX block (no need to have TAXALABELS block).
30 32
     ##
31 33
new file mode 100644
... ...
@@ -0,0 +1,14 @@
1
+##' print information of a list of beast trees
2
+##'
3
+##' 
4
+##' @title print
5
+##' @param x a list of beast object
6
+##' @param ... no used
7
+##' @return message
8
+##' @method print beastList
9
+##' @export
10
+##' @author Guangchuang Yu
11
+print.beastList <- function(x, ...) {
12
+    msg <- paste(length(x), "beast trees")
13
+    cat(msg, "\n")
14
+}
... ...
@@ -1,5 +1,3 @@
1
-
2
-
3 1
 ##' @rdname show-methods
4 2
 ##' @importFrom ape print.phylo
5 3
 ##' @exportMethod show
... ...
@@ -209,6 +209,10 @@ getSubsLabel <- function(seqs, A, B, translate, removeGap) {
209 209
     seqA <- seqs[A]
210 210
     seqB <- seqs[B]
211 211
 
212
+    if (nchar(seqA) != nchar(seqB)) {
213
+        stop("seqA should have equal length to seqB")
214
+    }
215
+    
212 216
     if (translate == TRUE) {
213 217
         AA <- seqA %>% seq2codon %>% codon2AA
214 218
         BB <- seqB %>% seq2codon %>% codon2AA
215 219
new file mode 100644
... ...
@@ -0,0 +1,23 @@
1
+% Generated by roxygen2: do not edit by hand
2
+% Please edit documentation in R/method-print.R
3
+\name{print.beastList}
4
+\alias{print.beastList}
5
+\title{print}
6
+\usage{
7
+\method{print}{beastList}(x, ...)
8
+}
9
+\arguments{
10
+\item{x}{a list of beast object}
11
+
12
+\item{...}{no used}
13
+}
14
+\value{
15
+message
16
+}
17
+\description{
18
+print information of a list of beast trees
19
+}
20
+\author{
21
+Guangchuang Yu
22
+}
23
+
... ...
@@ -6,7 +6,7 @@
6 6
 \usage{
7 7
 stat_hilight(mapping = NULL, data = NULL, geom = "rect",
8 8
   position = "identity", node, show.legend = NA, inherit.aes = FALSE,
9
-  fill, alpha, extend = 0, xmax = NULL, ...)
9
+  fill, alpha, extend = 0, extendto = NULL, ...)
10 10
 }
11 11
 \arguments{
12 12
 \item{mapping}{aes mapping}
... ...
@@ -29,9 +29,9 @@ stat_hilight(mapping = NULL, data = NULL, geom = "rect",
29 29
 
30 30
 \item{extend}{extend xmax of the rectangle}
31 31
 
32
-\item{...}{additional parameter}
33
-
34 32
 \item{extendto}{extend xmax to extendto}
33
+
34
+\item{...}{additional parameter}
35 35
 }
36 36
 \value{
37 37
 layer
... ...
@@ -33,6 +33,7 @@ collapse <- ggtree::collapse
33 33
 expand <- ggtree::expand
34 34
 rotate <- ggtree::rotate
35 35
 flip <- ggtree::flip
36
+MRCA <- ggtree::MRCA
36 37
 ```
37 38
 
38 39
 # Internal node number