Browse code

Added revised geom_hilight_encircle() function to highlight clades of unrooted trees.

Fixed layoutDaylight ave_change calculation.

Fixed bug in getTreeArcAngles() for cases where the branch root node and origin node are the same.

Added getNodeEuclDistances() function.

Fixed isTip() function to check if nodes has children, instead of checking of the data variable had the isTip field set.

JustGitting authored on 29/06/2017 03:43:04
Showing 112 changed files

... ...
@@ -7,6 +7,8 @@ vignettes/.DS_Store
7 7
 __init__.py
8 8
 __pycache__
9 9
 __init__.pyc
10
-.Rproj.user
10
+*.Rproj.user
11
+*.RData
12
+*.ignore
13
+.web_cache
11 14
 ggtree.Rproj
12
-.RData
... ...
@@ -3,3 +3,7 @@ Justin Silverman
3 3
 + `geom_balance`
4 4
 	- <https://github.com/GuangchuangYu/ggtree/pull/64>
5 5
 
6
+JustGitting
7
+----------------
8
++ `Add daylight layout for unrooted trees`
9
+  - <https://github.com/GuangchuangYu/ggtree/pull/124>
6 10
\ No newline at end of file
... ...
@@ -43,4 +43,4 @@ BugReports: https://github.com/GuangchuangYu/ggtree/issues
43 43
 Packaged: 2014-12-03 08:16:14 UTC; root
44 44
 biocViews: Alignment, Annotation, Clustering, DataImport,
45 45
     MultipleSequenceAlignment, ReproducibleResearch, Software, Visualization
46
-RoxygenNote: 5.0.1
46
+RoxygenNote: 6.0.1
... ...
@@ -27,6 +27,7 @@ export("%<+%")
27 27
 export("%>%")
28 28
 export(.)
29 29
 export(Date2decimal)
30
+export(GeomHilight)
30 31
 export(MRCA)
31 32
 export(StatBalance)
32 33
 export(StatHilight)
... ...
@@ -44,7 +45,9 @@ export(flip)
44 45
 export(geom_aline)
45 46
 export(geom_balance)
46 47
 export(geom_cladelabel)
48
+export(geom_cladelabel2)
47 49
 export(geom_hilight)
50
+export(geom_hilight_encircle)
48 51
 export(geom_label2)
49 52
 export(geom_nodepoint)
50 53
 export(geom_point2)
... ...
@@ -87,6 +90,7 @@ export(scaleClade)
87 90
 export(scale_color)
88 91
 export(scale_x_ggtree)
89 92
 export(stat_balance)
93
+export(stat_chull)
90 94
 export(stat_hilight)
91 95
 export(subview)
92 96
 export(theme_inset)
... ...
@@ -109,6 +113,7 @@ importFrom(ape,ladderize)
109 113
 importFrom(ape,read.tree)
110 114
 importFrom(ape,reorder.phylo)
111 115
 importFrom(ape,write.tree)
116
+importFrom(ggplot2,Geom)
112 117
 importFrom(ggplot2,GeomCurve)
113 118
 importFrom(ggplot2,GeomLabel)
114 119
 importFrom(ggplot2,GeomPoint)
115 120
new file mode 100644
... ...
@@ -0,0 +1 @@
1
+stat_chull.R
... ...
@@ -38,13 +38,20 @@ getMRCA.df <- function(data, tip) {
38 38
 }
39 39
 
40 40
 
41
-getMRCA.df_internal <- function(data, tip1, tip2) {
42
-    node1 <- which(tip1 == data$label | tip1 == data[, "node"])
43
-    node2 <- which(tip2 == data$label | tip2 == data[, "node"])
41
+getMRCA.df_internal <- function(data, tip, anc) {
42
+    node1 <- which(tip == data$label | tip == data[, "node"])
43
+    node2 <- which(anc == data$label | anc == data[, "node"])
44 44
     
45 45
     anc1 <- get.ancestor.df(data, node1)
46 46
     anc2 <- get.ancestor.df(data, node2)
47 47
     
48
+    if(is.null(anc1)){
49
+      print("Warning getMRCA.df_internal(): tip is root")
50
+    }else if(is.null(anc1)){
51
+      print("Warning getMRCA.df_internal(): anc is root")
52
+    }
53
+    
54
+    # Return common ancestors.
48 55
     intersect(c(node1, anc1), c(node2, anc2))[1]
49 56
 }
50 57
 
... ...
@@ -53,7 +60,9 @@ get.ancestor.df <- function(df, node) {
53 60
     pp <- getParent.df(df, node)
54 61
     pp <- pp[pp != 0]
55 62
     if (length(pp) == 0) {
56
-        stop("input node is root...")
63
+      #stop("input node is root...")
64
+      cat("WARNING: input node ",node," is root.")
65
+      return(NULL) # root has no ancestor
57 66
     }
58 67
     i <- 1
59 68
     while(i <= length(pp)) {
... ...
@@ -94,18 +94,18 @@ geom_cladelabel <- function(node, label, offset=0, offset.text=0,
94 94
                                         inherit.aes = inherit.aes, na.rm=na.rm,
95 95
                                         parse = parse, ...)
96 96
         }
97
-
98
-        layer_bar <- stat_cladeBar(node=node, offset=offset, align=align,
99
-                                   size=barsize, color = barcolor,
100
-                                   mapping=mapping, data=data,
101
-                                   position=position, show.legend = show.legend,
102
-                                   inherit.aes = inherit.aes, na.rm=na.rm, ...)
103
-
97
+      
98
+      layer_bar <- stat_cladeBar(node=node, offset=offset, align=align,
99
+                                 size=barsize, color = barcolor,
100
+                                 mapping=mapping, data=data,
101
+                                 position=position, show.legend = show.legend,
102
+                                 inherit.aes = inherit.aes, na.rm=na.rm, ...)
103
+      
104 104
     }
105
-
105
+    
106 106
     list(
107
-       layer_bar,
108
-       layer_text
107
+      layer_bar,
108
+      layer_text
109 109
     )
110 110
 }
111 111
 
... ...
@@ -115,99 +115,100 @@ stat_cladeText <- function(mapping=NULL, data=NULL,
115 115
                            node, label, offset, align, ...,
116 116
                            show.legend=NA, inherit.aes=FALSE,
117 117
                            na.rm=FALSE, parse=FALSE) {
118
-    default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent)
119
-    if (is.null(mapping)) {
120
-        mapping <- default_aes
121
-    } else {
122
-        mapping <- modifyList(mapping, default_aes)
123
-    }
124
-
125
-    layer(stat=StatCladeText,
126
-          data=data,
127
-          mapping=mapping,
128
-          geom=geom,
129
-          position=position,
130
-          show.legend = show.legend,
131
-          inherit.aes = inherit.aes,
132
-          params=list(node=node,
133
-                      label  = label,
134
-                      offset = offset,
135
-                      align  = align,
136
-                      na.rm  = na.rm,
137
-                      parse  = parse,
138
-                      ...),
139
-          check.aes = FALSE
140
-          )
141
-
118
+  default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent)
119
+  if (is.null(mapping)) {
120
+    mapping <- default_aes
121
+  } else {
122
+    mapping <- modifyList(mapping, default_aes)
123
+  }
124
+  
125
+  layer(stat=StatCladeText,
126
+        data=data,
127
+        mapping=mapping,
128
+        geom=geom,
129
+        position=position,
130
+        show.legend = show.legend,
131
+        inherit.aes = inherit.aes,
132
+        params=list(node=node,
133
+                    label  = label,
134
+                    offset = offset,
135
+                    align  = align,
136
+                    na.rm  = na.rm,
137
+                    parse  = parse,
138
+                    ...),
139
+        check.aes = FALSE
140
+  )
141
+  
142 142
 }
143 143
 
144 144
 stat_cladeBar <- function(mapping=NULL, data=NULL,
145 145
                           geom="segment", position="identity",
146 146
                           node, offset, align,  ...,
147 147
                           show.legend=NA, inherit.aes=FALSE, na.rm=FALSE) {
148
-    default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent, xend=~x, yend=~y)
149
-    if (is.null(mapping)) {
150
-        mapping <- default_aes
151
-    } else {
152
-        mapping <- modifyList(mapping, default_aes)
153
-    }
154
-
155
-    layer(stat=StatCladeBar,
156
-          data=data,
157
-          mapping=mapping,
158
-          geom=geom,
159
-          position=position,
160
-          show.legend = show.legend,
161
-          inherit.aes = inherit.aes,
162
-          params=list(node=node,
163
-                      offset=offset,
164
-                      align=align,
165
-                      na.rm=na.rm,
166
-                      ...),
167
-          check.aes = FALSE
168
-          )
148
+  default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent, xend=~x, yend=~y)
149
+  if (is.null(mapping)) {
150
+    mapping <- default_aes
151
+  } else {
152
+    mapping <- modifyList(mapping, default_aes)
153
+  }
154
+  
155
+  layer(stat=StatCladeBar,
156
+        data=data,
157
+        mapping=mapping,
158
+        geom=geom,
159
+        position=position,
160
+        show.legend = show.legend,
161
+        inherit.aes = inherit.aes,
162
+        params=list(node=node,
163
+                    offset=offset,
164
+                    align=align,
165
+                    na.rm=na.rm,
166
+                    ...),
167
+        check.aes = FALSE
168
+  )
169 169
 }
170 170
 
171 171
 StatCladeText <- ggproto("StatCladeText", Stat,
172 172
                          compute_group = function(self, data, scales, params, node, label, offset, align) {
173
-                             df <- get_cladelabel_position(data, node, offset, align, adjustRatio = 1.03)
174
-                             df$y <- mean(c(df$y, df$yend))
175
-                             df$label <- label
176
-                             return(df)
173
+                           df <- get_cladelabel_position(data, node, offset, align, adjustRatio = 1.03)
174
+                           df$y <- mean(c(df$y, df$yend))
175
+                           df$label <- label
176
+                           return(df)
177 177
                          },
178 178
                          required_aes = c("x", "y", "label")
179
-                         )
179
+)
180 180
 
181 181
 
182 182
 
183 183
 StatCladeBar <- ggproto("StatCladBar", Stat,
184 184
                         compute_group = function(self, data, scales, params, node, offset, align) {
185
-                            get_cladelabel_position(data, node, offset, align, adjustRatio=1.02)
185
+                          get_cladelabel_position(data, node, offset, align, adjustRatio=1.02)
186 186
                         },
187 187
                         required_aes = c("x", "y", "xend", "yend")
188
-                        )
188
+)
189 189
 
190 190
 
191 191
 get_cladelabel_position <- function(data, node, offset, align, adjustRatio) {
192
-    df <- get_cladelabel_position_(data, node)
193
-    if (align) {
194
-        mx <- max(data$x, na.rm=TRUE)
195
-    } else {
196
-        mx <- df$x
197
-    }
198
-    mx <- mx * adjustRatio + offset
199
-    data.frame(x=mx, xend=mx, y=df$y, yend=df$yend)
192
+  df <- get_cladelabel_position_(data, node)
193
+  if (align) {
194
+    # Find max x value for all tree nodes so all clade labels align to same position.
195
+    mx <- max(data$x, na.rm=TRUE)
196
+  } else {
197
+    mx <- df$x
198
+  }
199
+  mx <- mx * adjustRatio + offset
200
+  data.frame(x=mx, xend=mx, y=df$y, yend=df$yend)
200 201
 }
201 202
 
202
-
203
+# get x, y and yend of clade region.
203 204
 get_cladelabel_position_ <- function(data, node) {
204
-    sp <- get.offspring.df(data, node)
205
-    sp2 <- c(sp, node)
206
-    sp.df <- data[match(sp2, data$node),]
207
-
208
-    y <- sp.df$y
209
-    y <- y[!is.na(y)]
210
-    mx <- max(sp.df$x, na.rm=TRUE)
211
-    data.frame(x=mx, y=min(y), yend=max(y))
205
+  sp <- get.offspring.df(data, node)
206
+  sp2 <- c(sp, node)
207
+  sp.df <- data[match(sp2, data$node),]
208
+  
209
+  y <- sp.df$y
210
+  y <- y[!is.na(y)]
211
+  mx <- max(sp.df$x, na.rm=TRUE)
212
+  data.frame(x=mx, y=min(y), yend=max(y))
212 213
 }
213 214
 
214 215
new file mode 100644
... ...
@@ -0,0 +1,321 @@
1
+##' annotate a clade with bar and text label
2
+##'
3
+##'
4
+##' @title geom_cladelabel2
5
+##' @param node selected node
6
+##' @param label clade label
7
+##' @param offset offset of bar and text from the clade
8
+##' @param offset.text offset of text from bar
9
+##' @param offset.bar offset of bar from text
10
+##' @param align logical 
11
+##' @param barsize size of bar
12
+##' @param fontsize font size of text
13
+##' @param angle angle of text
14
+##' @param geom one of 'text' or 'label'
15
+##' @param hjust justify text horizontally
16
+##' @param color color for clade & label, of length 1 or 2
17
+##' @param fill fill label background, only work with geom='label'
18
+##' @param family sans by default, can be any supported font
19
+##' @param parse logical, whether parse label
20
+##' @param ... additional parameter
21
+##' @return ggplot layers
22
+##' @export
23
+##' @author JustGitting
24
+geom_cladelabel2 <- function(node, label, offset=0, offset.text=0, offset.bar=0,
25
+                            align=FALSE, barsize=0.5, fontsize=3.88, hjust = 0,
26
+                            geom="text",
27
+                            color = NULL,
28
+                            family="sans", parse=FALSE, ...) {
29
+    mapping <- NULL
30
+    data <- NULL
31
+    position <- "identity"
32
+    show.legend <- NA
33
+    na.rm <- TRUE
34
+    inherit.aes <- FALSE
35
+
36
+    
37
+    # create custom arguments from ellipsis (aka '...') for stat_cladeText2 depending on geom type
38
+    
39
+    # http://ggplot2.tidyverse.org/reference/geom_text.html 
40
+    # geom_label(mapping = NULL, data = NULL, stat = "identity",
41
+    #            position = "identity", ..., parse = FALSE, nudge_x = 0, nudge_y = 0,
42
+    #            label.padding = unit(0.25, "lines"), label.r = unit(0.15, "lines"),
43
+    #            label.size = 0.25, na.rm = FALSE, show.legend = NA,
44
+    #            inherit.aes = TRUE)
45
+    # 
46
+    # geom_text(mapping = NULL, data = NULL, stat = "identity",
47
+    #           position = "identity", ..., parse = FALSE, nudge_x = 0, nudge_y = 0,
48
+    #           check_overlap = FALSE, na.rm = FALSE, show.legend = NA,
49
+    #           inherit.aes = TRUE)
50
+    #
51
+    # Aesthetics: x, y, label, alpha, angle, colour, family, fontface, group, hjust, lineheight, size, vjust
52
+    
53
+    # http://ggplot2.tidyverse.org/reference/geom_segment.html
54
+    # geom_curve(mapping = NULL, data = NULL, stat = "identity",
55
+    #            position = "identity", ..., curvature = 0.5, angle = 90, ncp = 5,
56
+    #            arrow = NULL, lineend = "butt", na.rm = FALSE, show.legend = NA,
57
+    #            inherit.aes = TRUE)
58
+    #
59
+    # Aesthetics:  x, y, xend, yend, alpha, colour, group, linetype, size
60
+    
61
+    # name_mapping = list('oldA'='newA', 'oldB'='newB')
62
+    # data_list = list(oldB=1, oldA=2)
63
+    # names(data_list) = name_mapping[match(names(data_list), names(name_mapping))]
64
+    
65
+    arg_list_geom_label <- c( "nudge_x", "nudge_y", "label.padding", "label.r", "label.size", 
66
+                              "alpha", "angle", "fontface", "group", "lineheight", "size", "vjust", "fill")
67
+    
68
+    arg_list_geom_text <- c( "nudge_x", "nudge_y", "check_overlap", 
69
+                             "alpha", "angle", "fontface", "group", "lineheight", "size", "vjust")
70
+    
71
+    # ignore angle
72
+    arg_list_geom_curve <- c( "curvature", "ncp", "arrow", "lineend", 
73
+                              "alpha", "group", "linetype")
74
+    
75
+    
76
+    # Parse ellipsis to collect parameters for geom_text or geom_label
77
+    ellipsis <- list(...)
78
+    if (geom == "text") {
79
+      args_stat_cladeText2 <- ellipsis[names(ellipsis) %in% arg_list_geom_text]  
80
+    } else {
81
+      args_stat_cladeText2 <- ellipsis[names(ellipsis) %in% arg_list_geom_label]
82
+    }
83
+    
84
+    if (parse == 'emoji') {
85
+      emoji <- get_fun_from_pkg("emojifont", "emoji")
86
+      label <- emoji(label)
87
+      parse <- FALSE
88
+      family <- "EmojiOne"
89
+    }
90
+    
91
+    
92
+    # add parameters to stat_cladeText2 options.
93
+    args_stat_cladeText2$node        <- node
94
+    args_stat_cladeText2$label       <- label
95
+    args_stat_cladeText2$offset      <- offset+offset.text
96
+    args_stat_cladeText2$align       <- align
97
+    args_stat_cladeText2$hjust       <- hjust
98
+    args_stat_cladeText2$size        <- fontsize
99
+    args_stat_cladeText2$family      <- family
100
+    args_stat_cladeText2$mapping     <- mapping
101
+    args_stat_cladeText2$data        <- data
102
+    args_stat_cladeText2$geom        <- geom 
103
+    args_stat_cladeText2$position    <- position
104
+    args_stat_cladeText2$show.legend <- show.legend
105
+    args_stat_cladeText2$inherit.aes <- inherit.aes 
106
+    args_stat_cladeText2$na.rm       <- na.rm
107
+    args_stat_cladeText2$parse       <- parse 
108
+    
109
+    # create arg list of stat_cladeBar2.
110
+    args_stat_cladeBar2 <- ellipsis[names(ellipsis) %in% arg_list_geom_curve]
111
+    
112
+    args_stat_cladeBar2$size        <- barsize
113
+    args_stat_cladeBar2$node        <- node
114
+    args_stat_cladeBar2$offset      <- offset+offset.bar
115
+    args_stat_cladeBar2$align       <- align
116
+    args_stat_cladeBar2$size        <- barsize
117
+    args_stat_cladeBar2$mapping     <- mapping
118
+    args_stat_cladeBar2$data        <- data
119
+    args_stat_cladeBar2$position    <- position
120
+    args_stat_cladeBar2$show.legend <- show.legend
121
+    args_stat_cladeBar2$inherit.aes <- inherit.aes
122
+    args_stat_cladeBar2$na.rm       <- na.rm
123
+    
124
+    
125
+    if (!is.null(color)) {
126
+        if (length(color) > 2) {
127
+          stop("color should be of length 1 or 2")
128
+        }
129
+        if (length(color) == 0) {
130
+          color = NULL
131
+        } else if (length(color) == 1) {
132
+          args_stat_cladeText2$colour <- color
133
+          args_stat_cladeBar2$colour <- color
134
+        } else {
135
+          args_stat_cladeText2$colour <- color[1]
136
+          args_stat_cladeBar2$colour <- color[2]
137
+        }
138
+    }
139
+    
140
+    # print('text opts') # Debug
141
+    # print(args_stat_cladeText2) # Debug
142
+    # print('bar opts') # Debug
143
+    # print(args_stat_cladeBar2) # Debug
144
+
145
+    # create text and bar layers.
146
+    layer_text <- do.call(stat_cladeText2, args_stat_cladeText2)
147
+    layer_bar <- do.call(stat_cladeBar2, args_stat_cladeBar2)
148
+
149
+    list(
150
+      layer_bar,
151
+      layer_text
152
+    )
153
+}
154
+
155
+# Display label at middle angle of clade subtree arc.
156
+stat_cladeText2 <- function(mapping=NULL, data=NULL,
157
+                            geom="text", position="identity",
158
+                            node, label, offset, align, ...,
159
+                            show.legend=NA, inherit.aes=FALSE,
160
+                            na.rm=FALSE, parse=FALSE) {
161
+  # columns from ggplot data data.frame.
162
+  default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent)
163
+  if (is.null(mapping)) {
164
+    mapping <- default_aes
165
+  } else {
166
+    mapping <- modifyList(mapping, default_aes)
167
+  }
168
+  
169
+  layer(stat=StatCladeText2,
170
+        data=data,
171
+        mapping=mapping,
172
+        geom=geom,
173
+        position=position,
174
+        show.legend = show.legend,
175
+        inherit.aes = inherit.aes,
176
+        check.aes = FALSE,
177
+        params=list(node=node,
178
+                    label  = label,
179
+                    offset = offset,
180
+                    align  = align,
181
+                    na.rm  = na.rm,
182
+                    parse  = parse,
183
+                    ...)
184
+        
185
+  )
186
+  
187
+}
188
+      
189
+stat_cladeBar2 <- function(mapping=NULL, data=NULL,
190
+                           geom="curve", position="identity",
191
+                           node, offset, align, ...,
192
+                           show.legend=NA, inherit.aes=FALSE, na.rm=FALSE) {
193
+  default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent, xend=~x, yend=~y)
194
+  if (is.null(mapping)) {
195
+    mapping <- default_aes 
196
+  } else {
197
+    mapping <- modifyList(mapping, default_aes)
198
+  }
199
+  
200
+  layer(stat=StatCladeBar2,
201
+        data=data,
202
+        mapping=mapping, 
203
+        geom=geom,
204
+        position=position,
205
+        show.legend = show.legend,
206
+        inherit.aes = inherit.aes,
207
+        check.aes = FALSE,
208
+        params=list(node=node,
209
+                    offset=offset,
210
+                    align=align,
211
+                    na.rm=na.rm,
212
+                    ...)
213
+        
214
+  )
215
+}
216
+
217
+StatCladeText2 <- ggproto("StatCladeText2", Stat,
218
+                          
219
+                          required_aes = c("x", "y", "label"),
220
+
221
+                          compute_group = function(self, data, scales, params = NULL, node, label, offset, align) {
222
+                            df <- get_cladelabel2_position_label(data, node, offset, align, adjustRatio = 1.2)
223
+
224
+                            # computer_group does not need to return df$label as label is declared in the geom_cladelabel2() function.
225
+                            # The data.frame returned by computer_group() does not override the variables explicitly specified in the geom_cladelabel2()
226
+                            # df$label <- label 
227
+                            
228
+                            if(is.null(params$angle)){
229
+                              df$angle <- df$theta_label * 180
230
+                              if( df$angle > 90 & df$angle < 270){
231
+                                # add 180 to angle so label is easy to ready
232
+                                df$angle <- df$angle + 180
233
+                              }
234
+                            }
235
+                          
236
+                            return(df)
237
+                          }
238
+)
239
+
240
+StatCladeBar2 <- ggproto("StatCladeBar2", Stat,
241
+                        
242
+                        required_aes = c("x", "y", "xend", "yend"),
243
+                        
244
+                        compute_group = function(self, data, scales, params, node, offset, align) {
245
+                          df <- get_cladelabel2_position_bar(data, node, offset, align, adjustRatio=1.1)
246
+                          return(df)
247
+                        }
248
+)
249
+
250
+get_cladelabel2_position_label <- function(data, node, offset, align, adjustRatio) {
251
+  df <- get_cladelabel2_position_(data, node)
252
+  
253
+  if (align) {
254
+    # Find max radius from tree root.
255
+    r <- max(getNodeEuclDistances(data, getRoot.df(data)))
256
+  } else {
257
+    r <- df$r
258
+  }
259
+  
260
+  r <- r * adjustRatio + offset
261
+  
262
+  # Calculate the angle between theta_left and theta_right
263
+  delta <- df$theta_left - df$theta_right
264
+  
265
+  if(delta > 0){
266
+    theta_label <- delta/2 + df$theta_right
267
+  }else if(delta < 0){
268
+    delta_adj <- delta + 2
269
+    theta_label <- delta_adj/2 + df$theta_right
270
+  }else{
271
+    theta_label <- df$theta_left
272
+  }
273
+  
274
+  # correct if theta_label > 360
275
+  if(theta_label > 2){
276
+    theta_label <-  theta_label - 2 
277
+  }
278
+  
279
+  # Calculate the position of the label
280
+  x1 <- r*cospi(theta_label) + data[data$node==node, 'x']
281
+  y1 <- r*sinpi(theta_label) + data[data$node==node, 'y']
282
+  
283
+  data.frame(x=x1, y=y1, theta_label=theta_label)
284
+  
285
+}
286
+
287
+
288
+get_cladelabel2_position_bar <- function(data, node, offset, align, adjustRatio) {
289
+  df <- get_cladelabel2_position_(data, node)
290
+
291
+  if (align) {
292
+    # Find max radius from tree root.
293
+    r <- max(getNodeEuclDistances(data, getRoot.df(data)))
294
+  } else {
295
+    r <- df$r
296
+  }
297
+
298
+  r <- r * adjustRatio + offset
299
+
300
+  # Calculate the left(end) and right(start) points for the arc
301
+  x1 <- r*cospi(df$theta_right) + data[data$node==node, 'x']
302
+  y1 <- r*sinpi(df$theta_right) + data[data$node==node, 'y']
303
+  xend <- r*cospi(df$theta_left) + data[data$node==node, 'x']
304
+  yend <- r*sinpi(df$theta_left) + data[data$node==node, 'y']
305
+  
306
+  data.frame(x=x1, y=y1, xend=xend, yend=yend)  
307
+
308
+}
309
+
310
+# Get clade subtree arc angles and maximum radius from clade node to all other clade nodes.
311
+get_cladelabel2_position_ <- function(data, node) {
312
+  # get left and right angles of the clade subtree.
313
+  subtree <- list( subtree = getSubtree.df(data, node), node = node )
314
+  
315
+  arc <- getTreeArcAngles(data, node, subtree)
316
+  # get max distance from node to clade tips.
317
+  r <- max(getNodeEuclDistances(data[data$node %in% subtree$subtree,], node))
318
+  
319
+  data.frame(r=r, theta_left=as.numeric(arc['left']), theta_right=as.numeric(arc['right']))
320
+}
321
+
... ...
@@ -1,52 +1,51 @@
1 1
 ##' layer of hilight clade with rectangle
2 2
 ##'
3
-##'
4 3
 ##' @title geom_hilight
5
-##' @param node selected node to hilight
6
-##' @param fill color fill
7
-##' @param alpha alpha (transparency)
8
-##' @param extend extend xmax of the rectangle
9
-##' @param extendto extend xmax to extendto
4
+##' @param node selected node to hilight (required)
5
+##' @param fill color fill (default = steelblue)
6
+##' @param alpha alpha transparency, (default = 0.5)
7
+##' @param extend extend xmax of the rectangle (default = 0)
8
+##' @param extendto extend xmax to extendto (default = NULL)
10 9
 ##' @return ggplot2
11 10
 ##' @export
12 11
 ##' @importFrom ggplot2 aes_
13 12
 ##' @importFrom ggplot2 GeomRect
14 13
 ##' @author Guangchuang Yu
15 14
 geom_hilight <- function(node, fill="steelblue", alpha=.5, extend=0, extendto=NULL) {
16
-
17
-
18
-    data = NULL
19
-    stat = "hilight"
20
-    position = "identity"
21
-    show.legend = NA
22
-    na.rm = TRUE
23
-    inherit.aes = FALSE
24
-
25
-    default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent, branch.length=~branch.length)
26
-    mapping <- default_aes
27
-
28
-
29
-    layer(
30
-        stat=StatHilight,
31
-        data = data,
32
-        mapping = mapping,
33
-        geom = GeomRect,
34
-        position = position,
35
-        show.legend=show.legend,
36
-        inherit.aes = inherit.aes,
37
-        params = list(node=node,
38
-                      fill=fill,
39
-                      alpha=alpha,
40
-                      extend=extend,
41
-                      extendto=extendto,
42
-                      na.rm = na.rm),
43
-        check.aes = FALSE
44
-    )
15
+  
16
+  data = NULL
17
+  stat = "hilight"
18
+  position = "identity"
19
+  show.legend = NA
20
+  na.rm = TRUE
21
+  inherit.aes = FALSE
22
+  check.aes = FALSE
23
+  
24
+  default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent, branch.length=~branch.length)
25
+  mapping <- default_aes
26
+  
27
+  l <- layer(
28
+    stat=StatHilight,
29
+    data = data,
30
+    mapping = mapping,
31
+    geom = GeomRect,
32
+    position = position,
33
+    show.legend=show.legend,
34
+    inherit.aes = inherit.aes,
35
+    check.aes = check.aes,
36
+    params = list(node=node,
37
+                  fill=fill,
38
+                  alpha=alpha,
39
+                  extend=extend,
40
+                  extendto=extendto,
41
+                  na.rm = na.rm)
42
+  )
43
+  
44
+  return(l)
45 45
 }
46 46
 
47 47
 ##' stat_hilight
48 48
 ##'
49
-##'
50 49
 ##' @title stat_hilight
51 50
 ##' @param mapping aes mapping
52 51
 ##' @param data data
... ...
@@ -68,28 +67,30 @@ stat_hilight <- function(mapping=NULL, data=NULL, geom="rect",
68 67
                          show.legend=NA, inherit.aes=FALSE,
69 68
                          fill, alpha, extend=0, extendto=NULL,
70 69
                          ...) {
71
-    default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent, branch.length=~branch.length)
72
-    if (is.null(mapping)) {
73
-        mapping <- default_aes
74
-    } else {
75
-        mapping <- modifyList(mapping, default_aes)
76
-    }
77
-
78
-    layer(
79
-        stat=StatHilight,
80
-        data = data,
81
-        mapping = mapping,
82
-        geom = geom,
83
-        position = position,
84
-        show.legend=show.legend,
85
-        inherit.aes = inherit.aes,
86
-        params = list(node=node,
87
-                      fill=fill,
88
-                      alpha=alpha,
89
-                      extend=extend,
90
-                      extendto=extendto,
91
-                      ...)
92
-        )
70
+  
71
+  default_aes <- aes_(x=~x, y=~y, node=~node, parent=~parent, branch.length=~branch.length)
72
+  
73
+  if (is.null(mapping)) {
74
+    mapping <- default_aes
75
+  } else {
76
+    mapping <- modifyList(mapping, default_aes)
77
+  }
78
+  
79
+  layer(
80
+    stat=StatHilight,
81
+    data = data,
82
+    mapping = mapping,
83
+    geom = geom,
84
+    position = position,
85
+    show.legend=show.legend,
86
+    inherit.aes = inherit.aes,
87
+    params = list(node=node,
88
+                  fill=fill,
89
+                  alpha=alpha,
90
+                  extend=extend,
91
+                  extendto=extendto,
92
+                  ...)
93
+  )
93 94
 }
94 95
 
95 96
 ##' StatHilight
... ...
@@ -100,19 +101,20 @@ stat_hilight <- function(mapping=NULL, data=NULL, geom="rect",
100 101
 ##' @export
101 102
 StatHilight <- ggproto("StatHilight", Stat,
102 103
                        compute_group = function(self, data, scales, params, node, extend, extendto) {
103
-                           df <- get_clade_position_(data, node)
104
-                           df$xmax <- df$xmax + extend
105
-                           if (!is.null(extendto) && !is.na(extendto)) {
106
-                               if (extendto < df$xmax) {
107
-                                   warning("extendto is too small, keep the original xmax value...")
108
-                               } else {
109
-                                   df$xmax <- extendto
110
-                               }
104
+                         
105
+                         df <- get_clade_position_(data, node)
106
+                         df$xmax <- df$xmax + extend
107
+                         if (!is.null(extendto) && !is.na(extendto)) {
108
+                           if (extendto < df$xmax) {
109
+                             warning("extendto is too small, keep the original xmax value...")
110
+                           } else {
111
+                             df$xmax <- extendto
111 112
                            }
112
-                           return(df)
113
+                         }
114
+                         return(df)
113 115
                        },
114 116
                        required_aes = c("x", "y", "branch.length")
115
-                       )
117
+)
116 118
 
117 119
 
118 120
 ##' get position of clade (xmin, xmax, ymin, ymax)
... ...
@@ -125,31 +127,31 @@ StatHilight <- ggproto("StatHilight", Stat,
125 127
 ##' @export
126 128
 ##' @author Guangchuang Yu
127 129
 get_clade_position <- function(treeview, node) {
128
-    get_clade_position_(treeview$data, node)
130
+  get_clade_position_(treeview$data, node)
129 131
 }
130 132
 
131 133
 get_clade_position_ <- function(data, node) {
132
-    sp <- tryCatch(get.offspring.df(data, node), error=function(e) NULL)
133
-
134
-    i <- match(node, data$node)
135
-    if (is.null(sp)) {
136
-        ## tip
137
-        sp.df <- data[i,]
138
-    } else {
139
-        sp <- c(sp, node)
140
-        sp.df <- data[match(sp, data$node),]
141
-    }
142
-
143
-    x <- sp.df$x
144
-    y <- sp.df$y
145
-
146
-    if ("branch.length" %in% colnames(data)) {
147
-        xmin <- min(x)-data[i, "branch.length"]/2
148
-    } else {
149
-        xmin <- min(sp.df$branch)
150
-    }
151
-    data.frame(xmin=xmin,
152
-               xmax=max(x),
153
-               ymin=min(y)-0.5,
154
-               ymax=max(y)+0.5)
134
+  sp <- tryCatch(get.offspring.df(data, node), error=function(e) NULL)
135
+  
136
+  i <- match(node, data$node)
137
+  if (is.null(sp)) {
138
+    ## tip
139
+    sp.df <- data[i,]
140
+  } else {
141
+    sp <- c(sp, node)
142
+    sp.df <- data[match(sp, data$node),]
143
+  }
144
+  
145
+  x <- sp.df$x
146
+  y <- sp.df$y
147
+  
148
+  if ("branch.length" %in% colnames(data)) {
149
+    xmin <- min(x)-data[i, "branch.length"]/2
150
+  } else {
151
+    xmin <- min(sp.df$branch)
152
+  }
153
+  data.frame(xmin=xmin,
154
+             xmax=max(x),
155
+             ymin=min(y)-0.5,
156
+             ymax=max(y)+0.5)
155 157
 }
156 158
new file mode 100644
... ...
@@ -0,0 +1,269 @@
1
+# Encircle code originally from:
2
+# https://github.com/hrbrmstr/ggalt/blob/master/R/geom_encircle.r
3
+
4
+# draw_key_hack
5
+draw_key_hack <- function(data, params, size) {
6
+  print('draw_key_hack') ##DEBUG
7
+  data$fill <- alpha(data$fill, data$alpha)
8
+  data$alpha <- 1
9
+  
10
+  grobTree(
11
+    if (!is.na(data$fill)) grid::rectGrob(gp = gpar(col = NA, fill = data$fill)),
12
+    draw_key_path(data, params)
13
+  )
14
+}
15
+
16
+#' GeomHilight
17
+#' @rdname ggtree-ggproto
18
+#' @format NULL
19
+#' @usage NULL
20
+#' @importFrom ggplot2 Geom
21
+#' @export
22
+GeomHilight <- ggproto("GeomHilight", Geom,
23
+                       # Required fields in the 'data' data.frame for draw_group().
24
+                       # Additional fields from the layer param field can be used and will be added to the data column.
25
+                       # eg. required_aes = c("x", "y", "branch.length", "clade_root_node") will add "clade_root_node"
26
+                       # to the 'data' data.frame passed to draw_{panel, group}()
27
+                       required_aes = c("x", "y", "branch.length", "clade_root_node"),
28
+                       # set default aesthetic parameters appended to 'data' data.frame
29
+                       default_aes = aes(colour   = "black",
30
+                                         fill     = "steelblue",
31
+                                         alpha    = 0.5,
32
+                                         expand   = 0.05,
33
+                                         spread   = 0.1,
34
+                                         linetype = 1,
35
+                                         size     = 1,
36
+                                         s_shape  = 0.5,  ## corresponds to default shape in xspline of -0.5
37
+                                         s_open   = FALSE),
38
+                       
39
+                       draw_key = draw_key_hack, ## ???
40
+                       
41
+                       # # Find set of nodes that define the clade.
42
+                       # setup_params = function(data, params) {
43
+                       #   print('setup_params()') ## DEBUG
44
+                       #   if (is.null(params$node)){
45
+                       #     # Assume clade subset is given by user via data = data[subset]
46
+                       #     return(params)
47
+                       #   }
48
+                       #   
49
+                       #   # Find set of child nodes from clade_node.
50
+                       #   clade_node <- 15
51
+                       #   
52
+                       #   #params$clade_root_node <- clade_node
53
+                       #   params
54
+                       # },
55
+                       
56
+                       draw_group = function(data, panel_scales, coord) {
57
+                         # Determine if tree is circular or radial as uses Polar coordinates.
58
+                         #"CoordCartesian" %in% class(coord)
59
+                         #"CoordPolar" %in% class(coord)
60
+                         
61
+                         # Get clade root node and clade node ids.
62
+                         clade_root_node <- data[1,]$clade_root_node
63
+                         
64
+                         # Check if clade parent node exists in data.
65
+                         if( !(clade_root_node %in% data$node) ){
66
+                           cat('ERROR: clade node id (',clade_root_node,') not found in tree data.\n')
67
+                           return(NULL)
68
+                         }
69
+                         
70
+                         clade_ids = ggtree:::getSubtree.df(data, clade_root_node)
71
+                         
72
+                         # Remove non-clade rows.
73
+                         data <- data[data$node %in% clade_ids,]
74
+                         
75
+                         # # Get layout
76
+                         #
77
+                         # layout <- data[1,]$layout
78
+                         # 
79
+                         # If layout is {"rectangular”, “slanted”, “fan”, “circular”, “radial”} then find set of points that define
80
+                         # the retangular region around the clade.
81
+                         # if( layout %in% c('rectangular', 'slanted', 'fan', 'circular', 'radial') ){
82
+                         #   
83
+                         #   # get number of clade nodes.
84
+                         #   n <- nrow(data)
85
+                         #   
86
+                         #   # Find min and max (x,y) coordinates to find rectangle covering the clade.
87
+                         #   X <- data$x
88
+                         #   #Y <- data$y
89
+                         #   
90
+                         #   min_x <- min(X)
91
+                         #   max_x <- max(X)
92
+                         #   #min_y <- min(Y)
93
+                         #   #max_y <- max(Y)
94
+                         #   
95
+                         #   
96
+                         #   # Start with single row
97
+                         #   #data <- data[1,]
98
+                         #   #data <- data[rep(seq_len(nrow(data)), 4), ]
99
+                         #   #data$x <- c(max_x, min_x, min_x, max_x)
100
+                         #   #data$y <- c(min_y, max_y, min_y, max_y)
101
+                         #   
102
+                         #   points_right <- data
103
+                         #   # Update points with bounded box (min and max of X )
104
+                         #   data$x <- min_x
105
+                         #   points_right$x <- max_x
106
+                         #   print('points_right')
107
+                         #   print(points_right)
108
+                         #   
109
+                         #   # Combine left and right extreme points
110
+                         #   data <- rbind(data, points_right)
111
+                         #   
112
+                         #   print('Box data') #DEBUG
113
+                         #   print(data) #DEBUG
114
+                         #                             
115
+                         # }
116
+                         
117
+                         # Create glob
118
+                         glob <- get_glob_encircle(data, panel_scales, coord)
119
+                         
120
+                         return(glob)
121
+                         
122
+                       }
123
+)
124
+
125
+
126
+get_glob_encircle <- function(data, panel_scales, coord){
127
+  coords <- coord$transform(data, panel_scales)
128
+  first_row <- coords[1, , drop = FALSE]
129
+  rownames(first_row) <- NULL ## prevent warning later
130
+  
131
+  m <- lapply(coords[,c("x","y")],mean,na.rm=TRUE)
132
+  ch <- grDevices::chull(coords[c("x","y")])
133
+  
134
+  mkcoords <- function(x,y) {
135
+    data.frame(x,y,first_row[!names(first_row) %in% c("x","y")])
136
+  }
137
+  
138
+  coords <- coords[ch,]
139
+  ## FIXME: using grid:: a lot. importFrom instead?
140
+  
141
+  ## convert from lengths to physical units, for computing *directions*
142
+  cc <- function(x,dir="x")
143
+    grid::convertUnit(grid::unit(x,"native"),"mm",typeFrom="dimension",
144
+                      axisFrom=dir,valueOnly=TRUE)
145
+  
146
+  ## convert back to native (e.g. native + snpc offset)
147
+  cc_inv <- function(x,dir="x")
148
+    grid::convertUnit(x,"native",typeFrom="location",
149
+                      axisFrom=dir,valueOnly=TRUE)
150
+  
151
+  cc_comb <- function(x1,x2,dir="x")
152
+    cc_inv(unit(x1,"native")+unit(x2,"snpc"),dir=dir)
153
+  
154
+  ## find normalized vector: d1 and d2 have $x, $y elements
155
+  normFun <- function(d1,d2) {
156
+    dx <- cc(d1$x-d2$x)
157
+    dy <- cc(d1$y-d2$y)
158
+    r <- sqrt(dx*dx+dy*dy)
159
+    list(x=dx/r,y=dy/r)
160
+  }
161
+  
162
+  if (nrow(coords)==1) {
163
+    ## only one point: make a diamond by spreading points vertically
164
+    ## and horizontally
165
+    coords <- with(coords,
166
+                   mkcoords(
167
+                     c(x,x+spread,x,x-spread),
168
+                     c(y+spread,y,y-spread,y)))
169
+  } else if (nrow(coords)==2) {
170
+    ## only two points: make a diamond by spreading points perpendicularly
171
+    rot <- matrix(c(0,1,-1,0),2)
172
+    dd <- c(rot %*% unlist(normFun(coords[1,],coords[2,])))*
173
+      coords$spread
174
+    coords <- with(coords, {
175
+      ## figure out rotated values, then convert *back* to native units
176
+      ## already in scaled units, so ignore?
177
+      x <- c(x[1],
178
+             m$x+dd[1], ## cc_comb(m$x,dd[1]),
179
+             x[2],
180
+             m$x-dd[1]) ## cc_comb(m$x,-dd[1]))
181
+      y <- c(y[1],
182
+             m$y+dd[2], ## cc_comb(m$y,dd[2],"y"),
183
+             y[2],
184
+             m$y-dd[2]) ## cc_comb(m$y,-dd[2],"y"))
185
+      mkcoords(x,y)
186
+    })
187
+  }
188
+  
189
+  disp <- normFun(coords,m)
190
+  
191
+  ## browser()
192
+  
193
+  gp <- grid::get.gpar()
194
+  pars1 <- c("colour","linetype","alpha","fill","size")
195
+  pars2 <- c("col","lty","alpha","fill","lwd")
196
+  gp[pars2] <- first_row[pars1]
197
+  grid::xsplineGrob(
198
+    with(coords,unit(x,"npc")+disp$x*unit(expand,"snpc")),
199
+    with(coords,unit(y,"npc")+disp$y*unit(expand,"snpc")),
200
+    ## coords$x,
201
+    ## coords$y,
202
+    shape = coords$s_shape-1,  ## kluge!
203
+    open = first_row$s_open,
204
+    gp = gp)
205
+  
206
+}
207
+
208
+
209
+
210
+#' layer of hilight clade with xspline
211
+#'
212
+#' @title geom_hilight_encircle
213
+#' @param node selected node to hilight (required)
214
+#' @param fill colour fill (default = steelblue)
215
+#' @param alpha alpha (transparency) (default = 0.5)
216
+#' @param expand expands the xspline clade region only (default = 0)
217
+#' @param spread spread of shape? (default = 0.1)
218
+#' @param linetype Line type of xspline (default = 1)
219
+#' @param size Size of xspline line (default = 1)
220
+#' @param s_shape Corresponds to shape of xspline (default = 0.5)  
221
+#' @param s_open Boolean switch determines if xspline shape is open or closed. (default = FALSE)
222
+#' @return ggplot2
223
+#' @export
224
+#' @importFrom ggplot2 aes_
225
+geom_hilight_encircle <- function(data = NULL,
226
+                                  node,
227
+                                  mapping     = NULL,
228
+                                  fill        = 'steelblue',
229
+                                  alpha       = 0.5,
230
+                                  expand      = 0, # expand whole hilight region.
231
+                                  ...) {
232
+  
233
+  position    = "identity"
234
+  na.rm       = TRUE
235
+  show.legend = NA
236
+  inherit.aes = FALSE
237
+  check.aes   = FALSE
238
+  
239
+  
240
+  # Select fields(columns) from the ggtree "data" data.frame to be passed to the GeomHilight ggproto object.
241
+  default_aes <- aes_( x=~x, y=~y, node=~node, parent=~parent, branch.length=~branch.length )
242
+  
243
+  if (is.null(mapping)) {
244
+    mapping <- default_aes
245
+  } else {
246
+    mapping <- modifyList(mapping, default_aes)
247
+  }
248
+  
249
+  # create xspline geom for non-uniform trees, e.g. unrooted layout
250
+  l <- layer( 
251
+    geom = GeomHilight,
252
+    stat = "identity",
253
+    mapping = mapping,  
254
+    data = data, 
255
+    position = position, 
256
+    show.legend = show.legend, 
257
+    inherit.aes = inherit.aes,
258
+    check.aes = check.aes,
259
+    params = list(clade_root_node = node, 
260
+                  fill = fill,
261
+                  alpha = alpha,
262
+                  expand = expand,
263
+                  na.rm = na.rm,
264
+                  ...) # Parameters  to geom
265
+  )
266
+  
267
+  return(l)
268
+  
269
+}
0 270
\ No newline at end of file
... ...
@@ -4,7 +4,7 @@
4 4
 ##' @title geom_label2
5 5
 ##' @param mapping the aesthetic mapping
6 6
 ##' @param data A layer specific dataset -
7
-##'             only needed if you want to override he plot defaults.
7
+##'             only needed if you want to override the plot defaults.
8 8
 ##' @param ... other arguments passed on to 'layer'
9 9
 ##' @param position The position adjustment to use for overlapping points on this layer
10 10
 ##' @param family sans by default, can be any supported font
... ...
@@ -77,8 +77,7 @@ geom_label2 <- function(mapping = NULL, data = NULL,
77 77
             label.r = label.r,
78 78
             label.size = label.size,
79 79
             na.rm = na.rm,
80
-            ...
81
-        ),
80
+            ...),
82 81
         check.aes = FALSE
83 82
     )
84 83
 }
... ...
@@ -160,7 +160,7 @@ reroot_node_mapping <- function(tree, tree2) {
160 160
 
161 161
 ##' @importFrom ape reorder.phylo
162 162
 layout.unrooted <- function(tree, branch.length="branch.length", layout.method="equal_angle", ...) {
163
-
163
+  
164 164
     df <- switch(layout.method,
165 165
                  equal_angle = layoutEqualAngle(tree, branch.length),
166 166
                  daylight = layoutDaylight(tree, branch.length)
... ...
@@ -279,8 +279,8 @@ layoutEqualAngle <- function(tree, branch.length ){
279 279
 layoutDaylight <- function( tree, branch.length ){
280 280
 
281 281
     ## How to set optimal
282
-    MAX_COUNT <- 100
283
-    MINIMUM_AVERAGE_ANGLE_CHANGE <- 0.01
282
+    MAX_COUNT <- 5
283
+    MINIMUM_AVERAGE_ANGLE_CHANGE <- 0.05
284 284
 
285 285
 
286 286
     ## Initialize tree.
... ...
@@ -289,7 +289,7 @@ layoutDaylight <- function( tree, branch.length ){
289 289
     ## nodes = get list of nodes in tree_df
290 290
     ## Get list of node id's.
291 291
     ## nodes <- getNodes_by_postorder(tree)
292
-    ## nodes <- GetSubtree.df(tree_df, root)
292
+    ## nodes <- getSubtree.df(tree_df, root)
293 293
 
294 294
     ## Get list of internal nodes
295 295
     ## nodes <- tree_df[tree_df$IsTip != TRUE]$nodes
... ...
@@ -315,12 +315,12 @@ layoutDaylight <- function( tree, branch.length ){
315 315
             result <- applyLayoutDaylight(tree_df, currentNode_id)
316 316
             tree_df <- result$tree
317 317
             total_max <- total_max + result$max_change
318
-
318
+            
319 319
         }
320
+        # Calculate the running average of angle changes.
321
+        ave_change <- total_max / length(nodes) * length(i)
320 322
 
321
-        ave_change <- total_max / length(nodes)
322
-
323
-        cat('Average angle change', ave_change,'\n')
323
+        cat('Average angle change [',i,']', ave_change,'\n')
324 324
 
325 325
         i <- i + 1
326 326
     }
... ...
@@ -372,7 +372,7 @@ layoutDaylight <- function( tree, branch.length ){
372 372
 ##   }
373 373
 ## }
374 374
 ## ```
375
-applyLayoutDaylight <- function(df, node_id ){
375
+applyLayoutDaylight <- function(df, node_id){
376 376
 
377 377
   max_change <- 0.0
378 378
 
... ...
@@ -436,12 +436,12 @@ applyLayoutDaylight <- function(df, node_id ){
436 436
 
437 437
 
438 438
 ##' Find the right (clockwise rotation, angle from +ve x-axis to furthest subtree nodes) and
439
-##' left (anti-clockwise angle from +ve x-axis to subtree)
439
+##' left (anti-clockwise angle from +ve x-axis to subtree) Returning arc angle in [0, 2] (0 to 360) domain.
440 440
 ##'
441 441
 ##' @title getTreeArcAngles
442 442
 ##' @param df tree data.frame
443 443
 ##' @param origin_id node id from which to calculate left and right hand angles of subtree.
444
-##' @param subtree named list of root id of subtree and list of node ids for given subtree.
444
+##' @param subtree named list of root id of subtree (node) and list of node ids for given subtree (subtree).
445 445
 ##' @return named list with right and left angles in range [0,2] i.e 1 = 180 degrees, 1.5 = 270 degrees.
446 446
 getTreeArcAngles <- function(df, origin_id, subtree) {
447 447
   # Initialise variables
... ...
@@ -452,15 +452,66 @@ getTreeArcAngles <- function(df, origin_id, subtree) {
452 452
   # Initialise angle from origin node to parent node.
453 453
   # If subtree_root_id is child of origin_id
454 454
   if( any(subtree_root_id == getChild.df(df, origin_id)) ){
455
-    theta_parent <- getNodeAngle.df(df, origin_id, subtree_root_id)
455
+    # get angle from original node to parent of subtree.
456
+    theta_left <- getNodeAngle.df(df, origin_id, subtree_root_id)
457
+    theta_right <- theta_left
458
+  }else if( subtree_root_id == origin_id){
459
+    # Special case.
460
+    # get angle from parent of subtree to children
461
+    children_ids <- getChild.df(df, subtree_root_id)
462
+    
463
+    if(length(children_ids) == 2){
464
+      # get angles from parent to it's two children.
465
+      theta1 <- getNodeAngle.df(df, origin_id, children_ids[1])
466
+      theta2 <- getNodeAngle.df(df, origin_id, children_ids[2])
467
+      
468
+      delta <- theta1 - theta2
469
+      
470
+  
471
+      # correct delta for points crossing 180/-180 quadrant.
472
+      if(delta > 1){
473
+        delta_adj = delta - 2
474
+      }else if(delta < -1){
475
+        delta_adj = delta + 2
476
+      }else{
477
+        delta_adj <- delta
478
+      }
479
+  
480
+      if(delta_adj >= 0){
481
+        theta_left = theta1
482
+        theta_right = theta2
483
+      }else if(delta_adj < 0){
484
+        theta_left = theta2
485
+        theta_right = theta1
486
+      }
487
+    }else{
488
+      # subtree only has one child node.
489
+      theta_left <- getNodeAngle.df(df, origin_id, children_ids[1])
490
+      theta_right <- theta_left
491
+    }
492
+    
456 493
   }else{
457 494
     # get the real root of df tree to initialise left and right angles.
458
-    theta_parent <- getNodeAngle.df(df, origin_id, getRoot.df(df))
495
+    tree_root <- getRoot.df(df)
496
+    if( !is.na(tree_root) & is.numeric(tree_root) ){
497
+      theta_left <- getNodeAngle.df(df, origin_id, tree_root)
498
+      theta_right <- theta_left
499
+    }else{
500
+      print('ERROR: no root found!')
501
+      theta_left <- NA
502
+    }
503
+
459 504
   }
505
+  
506
+  # no parent angle found.
507
+  if (is.na(theta_left) ){
508
+    return(0)
509
+  }
510
+  
460 511
 
461 512
   # create vector with named columns
462 513
   # left-hand and right-hand angles between origin node and the extremities of the tree nodes.
463
-  arc <- c('left' = theta_parent, 'right' = theta_parent)
514
+  arc <- c('left' = theta_left, 'right' = theta_right)
464 515
 
465 516
   # Subtree has to have 1 or more nodes to compare.
466 517
   if (length(subtree_node_ids) == 0 ){
... ...
@@ -478,8 +529,10 @@ getTreeArcAngles <- function(df, origin_id, subtree) {
478 529
   for( i in seq_along(subtree_node_ids) ){
479 530
     parent_id <- subtree_node_ids[i]
480 531
     # Get angle from origin node to parent node.
481
-    # Skip if parent_id is a tip.
482
-    if(isTip.df(df, parent_id) ){ next }
532
+    # Skip if parent_id is a tip or parent and child node are the same.
533
+    if(origin_id == parent_id | isTip.df(df, parent_id) ){
534
+      next 
535
+    }
483 536
 
484 537
     theta_parent <- getNodeAngle.df(df, origin_id, parent_id)
485 538
 
... ...
@@ -609,7 +662,7 @@ rotateTreePoints.df <- function(df, pivot_node, nodes, angle){
609 662
 ##' @param node_id end node id number
610 663
 ##' @return angle in range [-1, 1], i.e. degrees/180, radians/pi
611 664
 getNodeAngle.df <- function(df, origin_node_id, node_id){
612
-  if(origin_node_id != node_id){
665
+  if( (origin_node_id != node_id) & any(origin_node_id %in% df$node) & any(node_id %in% df$node) ){
613 666
     delta_x <- df[node_id, 'x'] - df[origin_node_id, 'x']
614 667
     delta_y <- df[node_id, 'y'] - df[origin_node_id, 'y']
615 668
     angle <- atan2(delta_y, delta_x) / pi
... ...
@@ -619,6 +672,15 @@ getNodeAngle.df <- function(df, origin_node_id, node_id){
619 672
   }
620 673
 }
621 674
 
675
+euc.dist <- function(x1, x2) sqrt(sum((x1 - x2) ^ 2))
676
+
677
+##' Get the distances from the node to all other nodes in data.frame (including itself if in df)
678
+getNodeEuclDistances <- function(df, node){
679
+  # https://stackoverflow.com/questions/24746892/how-to-calculate-euclidian-distance-between-two-points-defined-by-matrix-contain#24747155
680
+  dist <- NULL
681
+  for(i in 1:nrow(df)) dist[i] <- euc.dist(df[df$node==node, c('x', 'y')], df[i, c('x', 'y')])
682
+  return(dist)
683
+}
622 684
 
623 685
 
624 686
 ##' Get all children of node from tree, including start_node.
... ...
@@ -642,12 +704,13 @@ getSubtree <- function(tree, node){
642 704
 
643 705
 ##' Get all children of node from df tree using breath-first.
644 706
 ##'
645
-##' @title GetSubtree.df
707
+##' @title getSubtree.df
646 708
 ##' @param df tree data.frame
647 709
 ##' @param node id of starting node.
648 710
 ##' @return list of all child node id's from starting node.
649
-GetSubtree.df <- function(df, node){
711
+getSubtree.df <- function(df, node){
650 712
   subtree <- c(node)
713
+  subtree <- subtree[subtree != 0]
651 714
   i <- 1
652 715
   while( i <= length(subtree)){
653 716
     subtree <- c(subtree, getChild.df(df, subtree[i]))
... ...
@@ -661,11 +724,11 @@ GetSubtree.df <- function(df, node){
661 724
 ##' Get all subtrees of specified node. This includes all ancestors and relatives of node and
662 725
 ##' return named list of subtrees.
663 726
 ##'
664
-##' @title GetSubtreeUnrooted
727
+##' @title getSubtreeUnrooted
665 728
 ##' @param tree ape phylo tree object
666 729
 ##' @param node is the tree node id from which the subtrees are derived.
667 730
 ##' @return named list of subtrees with the root id of subtree and list of node id's making up subtree.
668
-GetSubtreeUnrooted <- function(tree, node){
731
+getSubtreeUnrooted <- function(tree, node){
669 732
   # if node leaf, return nothing.
670 733
   if( isTip(tree, node) ){
671 734
     # return NA
... ...
@@ -704,7 +767,7 @@ GetSubtreeUnrooted <- function(tree, node){
704 767
 
705 768
 ##' Get all subtrees of node, as well as remaining branches of parent (ie, rest of tree structure as subtree)
706 769
 ##' return named list of subtrees with list name as starting node id.
707
-##' @title GetSubtreeUnrooted
770
+##' @title getSubtreeUnrooted
708 771
 ##' @param df tree data.frame
709 772
 ##' @param node is the tree node id from which the subtrees are derived.
710 773
 ##' @return named list of subtrees with the root id of subtree and list of node id's making up subtree.
... ...
@@ -726,7 +789,7 @@ getSubtreeUnrooted.df <- function(df, node){
726 789
   remaining_nodes <- setdiff(remaining_nodes, node)
727 790
 
728 791
   for( child in children_ids ){
729
-    subtree <- GetSubtree.df(df, child)
792
+    subtree <- getSubtree.df(df, child)
730 793
     # Append subtree nodes to list if more than 1 node in subtree (i.e. not a tip)
731 794
     #if(length(subtree) >= 2){
732 795
       subtrees[[length(subtrees)+1]] <- list( node = child, subtree = subtree)
... ...
@@ -751,7 +814,13 @@ getSubtreeUnrooted.df <- function(df, node){
751 814
 
752 815
 
753 816
 getRoot.df <- function(df, node){
817
+  
754 818
   root <- which(is.na(df$parent))
819
+  # Check if root was found.
820
+  if(length(root) == 0){
821
+    # Alternatively, root can self reference, eg node = 10, parent = 10
822
+    root <- unlist(apply(df, 1, function(x){ if(x['node'] == x['parent']){ x['node'] } }))
823
+  }
755 824
   return(root)
756 825