Browse code

move source code to tidytree

guangchuang yu authored on 07/12/2017 08:55:59
Showing 4 changed files

... ...
@@ -15,7 +15,7 @@ Description: 'ggtree' extends the 'ggplot2' plotting system which implemented
15 15
 Depends:
16 16
     R (>= 3.4.0),
17 17
     ggplot2 (>= 2.2.0),
18
-    treeio (>= 1.3.2)
18
+    treeio (>= 1.3.3)
19 19
 Imports:
20 20
     ape,
21 21
     dplyr,
22 22
deleted file mode 100644
... ...
@@ -1,64 +0,0 @@
1
-## not used.
2
-## we don't guess mrsd from tip labels
3
-## user should provide mrsd parameter in ggtree if they want to use time-scaled tree
4
-scaleX_by_time <- function(df, as.Date=FALSE) {
5
-    time <- with(df, gsub(".*[_/]{1}(\\d+\\.*\\d+)$", "\\1", label[isTip])) %>% as.numeric
6
-    latest <- which.max(time)
7
-
8
-    scaleX_by_time_from_mrsd(df, decimal2Date(time[latest]), as.Date)
9
-}
10
-
11
-
12
-scaleX_by_time_from_mrsd <- function(df, mrsd, as.Date) {
13
-    mrsd %<>% as.Date
14
-    date <- Date2decimal(mrsd)
15
-
16
-    df$x <- df$x + date - max(df$x)
17
-    df$branch <- with(df, (x[match(parent, node)] + x)/2)
18
-    ## df$branch <- (df[df$parent, "x"] + df[, "x"])/2
19
-
20
-    if (as.Date) {
21
-        df$x <- decimal2Date(df$x)
22
-        df$branch <- decimal2Date(df$branch)
23
-    }
24
-
25
-    return(df)
26
-}
27
-
28
-
29
-
30
-##' convert Date to decimal format, eg "2014-05-05" to "2014.34"
31
-##'
32
-##'
33
-##' @title Date2decimal
34
-##' @param x Date
35
-##' @return numeric
36
-##' @export
37
-##' @author Guangchuang Yu
38
-Date2decimal <- function(x) {
39
-    if (is(x, "numeric")) {
40
-        return(x)
41
-    }
42
-
43
-    if (is(x, "character")) {
44
-        x <- as.Date(x)
45
-    }
46
-    year <- format(x, "%Y")
47
-    y <- x - as.Date(paste0(year, "-01-01"))
48
-    as.numeric(year) + as.numeric(y)/365
49
-}
50
-
51
-##' convert decimal format to Date, eg "2014.34" to "2014-05-05"
52
-##'
53
-##'
54
-##' @title decimal2Date
55
-##' @param x numerical number, eg 2014.34
56
-##' @return Date
57
-##' @export
58
-##' @author Guangchuang Yu
59
-decimal2Date <- function(x) {
60
-    date <- as.Date(paste0(floor(x), "-01-01"))
61
-    date + as.numeric(sub("^\\d+", "0", x)) * 365
62
-}
63
-
64
-
... ...
@@ -347,29 +347,6 @@ rm.singleton.newick <- function(nwk, outfile = NULL) {
347 347
 ## fortify.hyphy <- fortify.paml_rst
348 348
 
349 349
 
350
-scaleY <- function(phylo, df, yscale, layout, ...) {
351
-    if (yscale == "none") {
352
-        return(df)
353
-    }
354
-    if (! yscale %in% colnames(df)) {
355
-        warning("yscale is not available...\n")
356
-        return(df)
357
-    }
358
-    if (is.numeric(df[[yscale]])) {
359
-        y <- getYcoord_scale_numeric(phylo, df, yscale, ...)
360
-        ## if (order.y) {
361
-        ##     y <- getYcoord_scale2(phylo, df, yscale)
362
-        ## } else {
363
-        ##     y <- getYcoord_scale(phylo, df, yscale)
364
-        ## }
365
-    } else {
366
-        y <- getYcoord_scale_category(phylo, df, yscale, ...)
367
-    }
368
-
369
-    df[, "y"] <- y
370
-
371
-    return(df)
372
-}
373 350
 
374 351
 
375 352
 ##' @method fortify phylo4
... ...
@@ -1,221 +1,9 @@
1
-##' @importFrom ggplot2 fortify
2
-##' @method fortify treedata
3
-##' @export
4
-fortify.treedata <- function(model, data,
5
-                             layout        = "rectangular",
6
-                             yscale        = "none",
7
-                             ladderize     = TRUE,
8
-                             right         = FALSE,
9
-                             branch.length = "branch.length",
10
-                             mrsd          = NULL,
11
-                             as.Date       = FALSE, ...) {
12
-
13
-    model <- set_branch_length(model, branch.length)
14
-
15
-    fortify.phylo(model, data,
16
-                  layout        = layout,
17
-                  yscale        = yscale,
18
-                  ladderize     = ladderize,
19
-                  right         = right,
20
-                  branch.length = branch.length,
21
-                  mrsd          = mrsd,
22
-                  as.Date       = as.Date, ...)
23
-}
24
-
25
-##' @importFrom ape ladderize
26
-##' @method fortify phylo
27
-##' @export
28
-fortify.phylo <- function(model, data,
29
-                          layout        = "rectangular",
30
-                          ladderize     = TRUE,
31
-                          right         = FALSE,
32
-                          branch.length = "branch.length",
33
-                          mrsd          = NULL,
34
-                          as.Date       = FALSE,
35
-                          yscale        = "none",
36
-                          ...) {
37
-
38
-    x <- as.phylo(model) ## reorder.phylo(get.tree(model), "postorder")
39
-    if (ladderize == TRUE) {
40
-        x <- ladderize(x, right=right)
41
-    }
42
-
43
-    if (! is.null(x$edge.length)) {
44
-        if (anyNA(x$edge.length)) {
45
-            warning("'edge.length' contains NA values...\n## setting 'edge.length' to NULL automatically when plotting the tree...")
46
-            x$edge.length <- NULL
47
-        }
48
-    }
49
-
50
-    if (is.null(x$edge.length) || branch.length == "none") {
51
-        xpos <- getXcoord_no_length(x)
52
-    } else {
53
-        xpos <- getXcoord(x)
54
-    }
55
-
56
-    ypos <- getYcoord(x)
57
-    N <- Nnode(x, internal.only=FALSE)
58
-    xypos <- data_frame(node=1:N, x=xpos, y=ypos)
59
-
60
-    df <- as_data_frame(model)
61
-
62
-    res <- full_join(df, xypos, by = "node")
63
-
64
-    ## add branch mid position
65
-    res <- calculate_branch_mid(res)
66
-
67
-    if (!is.null(mrsd)) {
68
-        res <- scaleX_by_time_from_mrsd(res, mrsd, as.Date)
69
-    }
70
-
71
-    if (layout == "slanted") {
72
-        res <- add_angle_slanted(res)
73
-    } else {
74
-        ## angle for all layout, if 'rectangular', user use coord_polar, can still use angle
75
-        res <- calculate_angle(res)
76
-    }
77
-    scaleY(as.phylo(model), res, yscale, layout, ...)
78
-}
79
-
80
-
81
-##' @method as_data_frame treedata
82
-##' @importFrom tibble as_data_frame
83
-##' @export
84
-##' @importFrom treeio Nnode
85
-##' @importFrom treeio Ntip
86
-as_data_frame.treedata <- function(x, ...) {
87
-    res <- as_data_frame(x@phylo)
88
-    tree_anno <- as_data_frame(get_tree_data(x))
89
-    if (nrow(tree_anno) > 0) {
90
-        by <- "node"
91
-        tree_anno$node <- as.integer(tree_anno$node)
92
-        if ("parent" %in% colnames(tree_anno)) {
93
-            by <- c(by, "parent")
94
-            tree_anno$parent <- as.integer(tree_anno$parent)
95
-        }
96
-
97
-        res <- full_join(res, tree_anno, by=by)
98
-    }
99
-    return(res)
100
-}
101
-
102
-##' @method as_data_frame phylo
103
-##' @export
104
-##' @importFrom tibble data_frame
105
-##' @importFrom dplyr full_join
106
-as_data_frame.phylo <- function(x, ...) {
107
-    phylo <- x
108
-    ntip <- Ntip(phylo)
109
-    N <- Nnode(phylo, internal.only=FALSE)
110
-
111
-    tip.label <- phylo[["tip.label"]]
112
-    res <- as_data_frame(phylo[["edge"]])
113
-    colnames(res) <- c("parent", "node")
114
-    if (!is.null(phylo$edge.length))
115
-        res$branch.length <- phylo$edge.length
116
-
117
-    label <- rep(NA, N)
118
-    label[1:ntip] <- tip.label
119
-    if ( !is.null(phylo$node.label) ) {
120
-        label[(ntip+1):N] <- phylo$node.label
121
-    }
122
-    isTip <- rep(FALSE, N)
123
-    isTip[1:ntip] <- TRUE
124
-
125
-    label.df <- data_frame(node=1:N, label=label, isTip = isTip)
126
-    res <- full_join(res, label.df, by='node')
127
-
128
-    idx <- is.na(res$parent)
129
-    res$parent[idx] <- res$node[idx]
130
-
131
-    res <- res[order(res$node),]
132
-    aa <- names(attributes(phylo))
133
-    group <- aa[ ! aa %in% c("names", "class", "order", "reroot", "node_map")]
134
-    if (length(group) > 0) {
135
-        for (group_ in group) {
136
-            ## groupOTU & groupClade
137
-            group_info <- attr(phylo, group_)
138
-            if (length(group_info) == nrow(res)) {
139
-                res[[group_]] <- group_info
140
-            }
141
-        }
142
-    }
143
-    return(res)
144
-}
145
-
146
-get_tree_data <- function(tree_object) {
147
-    tree_anno <- tree_object@data
148
-    extraInfo <- tree_object@extraInfo
149
-
150
-    if (nrow(tree_anno) == 0) {
151
-        return(extraInfo)
152
-    }
153
-    if (nrow(extraInfo) == 0) {
154
-        return(tree_anno)
155
-    }
156
-
157
-    full_join(tree_anno, extraInfo, by = "node")
158
-}
159
-
160
-
161
-##' convert tip or node label(s) to internal node number
162
-##'
163
-##'
164
-##' @title nodeid
165
-##' @param x tree object or graphic object return by ggtree
166
-##' @param label tip or node label(s)
167
-##' @return internal node number
168
-##' @export
169
-##' @author Guangchuang Yu
170
-nodeid <- function(x, label) {
171
-    if (is(x, "gg"))
172
-        return(nodeid.gg(x, label))
173
-
174
-    nodeid.tree(x, label)
175
-}
176
-
177
-nodeid.tree <- function(tree, label) {
178
-    tr <- get.tree(tree)
179
-    lab <- c(tr$tip.label, tr$node.label)
180
-    match(label, lab)
181
-}
182
-
183
-nodeid.gg <- function(p, label) {
184
-    p$data$node[match(label, p$data$label)]
185
-}
186 1
 
187 2
 
188
-reroot_node_mapping <- function(tree, tree2) {
189
-    root <- getRoot(tree)
190
-
191
-    node_map <- data.frame(from=1:getNodeNum(tree), to=NA, visited=FALSE)
192
-    node_map[1:Ntip(tree), 2] <- match(tree$tip.label, tree2$tip.label)
193
-    node_map[1:Ntip(tree), 3] <- TRUE
194 3
 
195
-    node_map[root, 2] <- root
196
-    node_map[root, 3] <- TRUE
197 4
 
198
-    node <- rev(tree$edge[,2])
199
-    for (k in node) {
200
-        ip <- getParent(tree, k)
201
-        if (node_map[ip, "visited"])
202
-            next
203 5
 
204
-        cc <- getChild(tree, ip)
205
-        node2 <- node_map[cc,2]
206
-        if (anyNA(node2)) {
207
-            node <- c(node, k)
208
-            next
209
-        }
210 6
 
211
-        to <- unique(sapply(node2, getParent, tr=tree2))
212
-        to <- to[! to %in% node_map[,2]]
213
-        node_map[ip, 2] <- to
214
-        node_map[ip, 3] <- TRUE
215
-    }
216
-    node_map <- node_map[, -3]
217
-    return(node_map)
218
-}
219 7
 
220 8
 
221 9
 
... ...
@@ -899,166 +687,12 @@ getRoot.df <- function(df, node){
899 687
   return(root)
900 688
 }
901 689
 
902
-##' Get parent node id of child node.
903
-##'
904
-##' @title getParent.df
905
-##' @param df tree data.frame
906
-##' @param node is the node id of child in tree.
907
-##' @return integer node id of parent
908
-getParent.df <- function(df, node) {
909
-    i <- which(df$node == node)
910
-    parent_id <- df$parent[i]
911
-    if (parent_id == node | is.na(parent_id)) {
912
-        ## root node
913
-        return(0)
914
-    }
915
-    return(parent_id)
916
-}
917
-
918
-getAncestor.df <- function(df, node) {
919
-    anc <- getParent.df(df, node)
920
-    anc <- anc[anc != 0]
921
-    if (length(anc) == 0) {
922
-        # stop("selected node is root...")
923
-      return(0)
924
-    }
925
-    i <- 1
926
-    while(i<= length(anc)) {
927
-        anc <- c(anc, getParent.df(df, anc[i]))
928
-        anc <- anc[anc != 0]
929
-        i <- i+1
930
-    }
931
-    return(anc)
932
-}
933
-
934
-
935
-##' Get list of child node id numbers of parent node
936
-##'
937
-##' @title getChild.df
938
-##' @param df tree data.frame
939
-##' @param node is the node id of child in tree.
940
-##' @return list of child node ids of parent
941
-getChild.df <- function(df, node) {
942
-    i <- which(df$parent == node)
943
-    if (length(i) == 0) {
944
-        return(0) # it has no children, hence tip node.
945
-    }
946
-    res <- df$node[i]
947
-    res <- res[res != node] ## node may root
948
-    return(res)
949
-}
950
-
951
-get.offspring.df <- function(df, node) {
952
-    sp <- getChild.df(df, node)
953
-    sp <- sp[sp != 0] # Remove root node.
954
-    if (length(sp) == 0) {
955
-        #stop("input node is a tip...")
956
-      return(0)
957
-    }
958
-
959
-    i <- 1
960
-    while(i <= length(sp)) {
961
-        sp <- c(sp, getChild.df(df, sp[i]))
962
-        sp <- sp[sp != 0]
963
-        i <- i + 1
964
-    }
965
-    return(sp)
966
-}
967
-
968
-
969
-##' extract offspring tips
970
-##'
971
-##'
972
-##' @title get.offspring.tip
973
-##' @param tr tree
974
-##' @param node node
975
-##' @return tip label
976
-##' @author ygc
977
-##' @importFrom ape extract.clade
978
-##' @export
979
-get.offspring.tip <- function(tr, node) {
980
-    if ( ! node %in% tr$edge[,1]) {
981
-        ## return itself
982
-        return(tr$tip.label[node])
983
-    }
984
-    clade <- extract.clade(tr, node)
985
-    clade$tip.label
986
-}
987
-
988 690
 
989
-## ##' calculate total number of nodes
990
-## ##'
991
-## ##'
992
-## ##' @title getNodeNum
993
-## ##' @param tr phylo object
994
-## ##' @return number
995
-## ##' @author Guangchuang Yu
996
-## ##' @export
997
-## getNodeNum <- function(tr) {
998
-##     Ntip <- length(tr[["tip.label"]])
999
-##     Nnode <- tr[["Nnode"]]
1000
-##     ## total nodes
1001
-##     N <- Ntip + Nnode
1002
-##     return(N)
1003
-## }
1004
-getParent <- function(tr, node) {
1005
-    if ( node == getRoot(tr) )
1006
-        return(0)
1007
-    edge <- tr[["edge"]]
1008
-    parent <- edge[,1]
1009
-    child <- edge[,2]
1010
-    res <- parent[child == node]
1011
-    if (length(res) == 0) {
1012
-        stop("cannot found parent node...")
1013
-    }
1014
-    if (length(res) > 1) {
1015
-        stop("multiple parent found...")
1016
-    }
1017
-    return(res)
1018
-}
1019 691
 
1020
-getChild <- function(tr, node) {
1021
-    # Get edge matrix from phylo object.
1022
-    edge <- tr[["edge"]]
1023
-    # Select all rows that match "node".
1024
-    res <- edge[edge[,1] == node, 2]
1025
-    ## if (length(res) == 0) {
1026
-    ##     ## is a tip
1027
-    ##     return(NA)
1028
-    ## }
1029
-    return(res)
1030
-}
1031
-
1032
-getSibling <- function(tr, node) {
1033
-    root <- getRoot(tr)
1034
-    if (node == root) {
1035
-        return(NA)
1036
-    }
1037 692
 
1038
-    parent <- getParent(tr, node)
1039
-    child <- getChild(tr, parent)
1040
-    sib <- child[child != node]
1041
-    return(sib)
1042
-}
1043 693
 
1044 694
 
1045
-getAncestor <- function(tr, node) {
1046
-    root <- getRoot(tr)
1047
-    if (node == root) {
1048
-        return(NA)
1049
-    }
1050
-    parent <- getParent(tr, node)
1051
-    res <- parent
1052
-    while(parent != root) {
1053
-        parent <- getParent(tr, parent)
1054
-        res <- c(res, parent)
1055
-    }
1056
-    return(res)
1057
-}
1058 695
 
1059
-isRoot <- function(tr, node) {
1060
-    getRoot(tr) == node
1061
-}
1062 696
 
1063 697
 isTip <- function(tr, node) {
1064 698
   children_ids <- getChild(tr, node)
... ...
@@ -1082,109 +716,6 @@ isTip.df <- function(df, node) {
1082 716
 }
1083 717
 
1084 718
 
1085
-getNodeName <- function(tr) {
1086
-    if (is.null(tr$node.label)) {
1087
-        n <- length(tr$tip.label)
1088
-        nl <- (n + 1):(2 * n - 2)
1089
-        nl <- as.character(nl)
1090
-    }
1091
-    else {
1092
-        nl <- tr$node.label
1093
-    }
1094
-    nodeName <- c(tr$tip.label, nl)
1095
-    return(nodeName)
1096
-}
1097
-
1098
-## ##' get the root number
1099
-## ##'
1100
-## ##'
1101
-## ##' @title getRoot
1102
-## ##' @param tr phylo object
1103
-## ##' @return root number
1104
-## ##' @export
1105
-## ##' @author Guangchuang Yu
1106
-## getRoot <- function(tr) {
1107
-##     edge <- tr[["edge"]]
1108
-##     ## 1st col is parent,
1109
-##     ## 2nd col is child,
1110
-##     if (!is.null(attr(tr, "order")) && attr(tr, "order") == "postorder")
1111
-##         return(edge[nrow(edge), 1])
1112
-
1113
-##     parent <- unique(edge[,1])
1114
-##     child <- unique(edge[,2])
1115
-##     ## the node that has no parent should be the root
1116
-##     root <- parent[ ! parent %in% child ]
1117
-##     if (length(root) > 1) {
1118
-##         stop("multiple roots founded...")
1119
-##     }
1120
-##     return(root)
1121
-## }
1122
-get.trunk <- function(tr) {
1123
-    root <- getRoot(tr)
1124
-    path_length <- sapply(1:(root-1), function(x) get.path_length(tr, root, x))
1125
-    i <- which.max(path_length)
1126
-    return(get.path(tr, root, i))
1127
-}
1128
-
1129
-##' path from start node to end node
1130
-##'
1131
-##'
1132
-##' @title get.path
1133
-##' @param phylo phylo object
1134
-##' @param from start node
1135
-##' @param to end node
1136
-##' @return node vectot
1137
-##' @export
1138
-##' @author Guangchuang Yu
1139
-get.path <- function(phylo, from, to) {
1140
-    anc_from <- getAncestor(phylo, from)
1141
-    anc_from <- c(from, anc_from)
1142
-    anc_to <- getAncestor(phylo, to)
1143
-    anc_to <- c(to, anc_to)
1144
-    mrca <- intersect(anc_from, anc_to)[1]
1145
-
1146
-    i <- which(anc_from == mrca)
1147
-    j <- which(anc_to == mrca)
1148
-
1149
-    path <- c(anc_from[1:i], rev(anc_to[1:(j-1)]))
1150
-    return(path)
1151
-}
1152
-
1153
-get.path_length <- function(phylo, from, to, weight=NULL) {
1154
-    path <- get.path(phylo, from, to)
1155
-    if (is.null(weight)) {
1156
-        return(length(path)-1)
1157
-    }
1158
-
1159
-    df <- fortify(phylo)
1160
-    if ( ! (weight %in% colnames(df))) {
1161
-        stop("weight should be one of numerical attributes of the tree...")
1162
-    }
1163
-
1164
-    res <- 0
1165
-
1166
-    get_edge_index <- function(df, from, to) {
1167
-        which((df[,1] == from | df[,2] == from) &
1168
-                  (df[,1] == to | df[,2] == to))
1169
-    }
1170
-
1171
-    for(i in 1:(length(path)-1)) {
1172
-        ee <- get_edge_index(df, path[i], path[i+1])
1173
-        res <- res + df[ee, weight]
1174
-    }
1175
-
1176
-    return(res)
1177
-}
1178
-
1179
-getNodes_by_postorder <- function(tree) {
1180
-  tree <- reorder.phylo(tree, "postorder")
1181
-    unique(rev(as.vector(t(tree$edge[,c(2,1)]))))
1182
-}
1183
-
1184
-#getNodes_by_postorder.df <- function(df) {
1185
-    #tree <- reorder.phylo(tree, "postorder")
1186
-    #unique(rev(as.vector(t(tree$edge[,c(2,1)]))))
1187
-#}
1188 719
 
1189 720
 ##' Get the nodes of tree from root in breadth-first order.
1190 721
 ##'
... ...
@@ -1225,516 +756,11 @@ getNodesBreadthFirst.df <- function(df){
1225 756
 }
1226 757
 
1227 758
 
1228
-getXcoord2 <- function(x, root, parent, child, len, start=0, rev=FALSE) {
1229
-    x[root] <- start
1230
-    x[-root] <- NA  ## only root is set to start, by default 0
1231
-
1232
-    currentNode <- root
1233
-    direction <- 1
1234
-    if (rev == TRUE) {
1235
-        direction <- -1
1236
-    }
1237
-
1238
-    while(anyNA(x)) {
1239
-        idx <- which(parent %in% currentNode)
1240
-        newNode <- child[idx]
1241
-        x[newNode] <- x[parent[idx]]+len[idx] * direction
1242
-        currentNode <- newNode
1243
-    }
1244
-
1245
-    return(x)
1246
-}
1247
-
1248
-getXcoord_no_length <- function(tr) {
1249
-    edge <- tr$edge
1250
-    parent <- edge[,1]
1251
-    child <- edge[,2]
1252
-    root <- getRoot(tr)
1253
-
1254
-    len <- tr$edge.length
1255
-
1256
-    N <- getNodeNum(tr)
1257
-    x <- numeric(N)
1258
-    ntip <- Ntip(tr)
1259
-    currentNode <- 1:ntip
1260
-    x[-currentNode] <- NA
1261
-
1262
-    cl <- split(child, parent)
1263
-    child_list <- list()
1264
-    child_list[as.numeric(names(cl))] <- cl
1265
-
1266
-    while(anyNA(x)) {
1267
-        idx <- match(currentNode, child)
1268
-        pNode <- parent[idx]
1269
-        ## child number table
1270
-        p1 <- table(parent[parent %in% pNode])
1271
-        p2 <- table(pNode)
1272
-        np <- names(p2)
1273
-        i <- p1[np] == p2
1274
-        newNode <- as.numeric(np[i])
1275
-
1276
-        exclude <- rep(NA, max(child))
1277
-        for (j in newNode) {
1278
-            x[j] <- min(x[child_list[[j]]]) - 1
1279
-            exclude[child_list[[j]]] <- child_list[[j]]
1280
-        }
1281
-        exclude <- exclude[!is.na(exclude)]
1282 759
 
1283
-        ## currentNode %<>% `[`(!(. %in% exclude))
1284
-        ## currentNode %<>% c(., newNode) %>% unique
1285
-        currentNode <- currentNode[!currentNode %in% exclude]
1286
-        currentNode <- unique(c(currentNode, newNode))
1287 760
 
1288
-    }
1289
-    x <- x - min(x)
1290
-    return(x)
1291
-}
1292 761
 
1293 762
 
1294
-getXcoord <- function(tr) {
1295
-    edge <- tr$edge
1296
-    parent <- edge[,1]
1297
-    child <- edge[,2]
1298
-    root <- getRoot(tr)
1299 763
 
1300
-    len <- tr$edge.length
1301 764
 
1302
-    N <- getNodeNum(tr)
1303
-    x <- numeric(N)
1304
-    x <- getXcoord2(x, root, parent, child, len)
1305
-    return(x)
1306
-}
1307 765
 
1308
-getXYcoord_slanted <- function(tr) {
1309 766
 
1310
-    edge <- tr$edge
1311
-    parent <- edge[,1]
1312
-    child <- edge[,2]
1313
-    root <- getRoot(tr)
1314
-
1315
-    N <- getNodeNum(tr)
1316
-    len <- tr$edge.length
1317
-    y <- getYcoord(tr, step=min(len)/2)
1318
-
1319
-    len <- sqrt(len^2 - (y[parent]-y[child])^2)
1320
-    x <- numeric(N)
1321
-    x <- getXcoord2(x, root, parent, child, len)
1322
-    res <- data.frame(x=x, y=y)
1323
-    return(res)
1324
-}
1325
-
1326
-
1327
-## @importFrom magrittr %>%
1328
-##' @importFrom magrittr equals
1329
-getYcoord <- function(tr, step=1) {
1330
-    Ntip <- length(tr[["tip.label"]])
1331
-    N <- getNodeNum(tr)
1332
-
1333
-    edge <- tr[["edge"]]
1334
-    parent <- edge[,1]
1335
-    child <- edge[,2]
1336
-
1337
-    cl <- split(child, parent)
1338
-    child_list <- list()
1339
-    child_list[as.numeric(names(cl))] <- cl
1340
-
1341
-    y <- numeric(N)
1342
-    tip.idx <- child[child <= Ntip]
1343
-    y[tip.idx] <- 1:Ntip * step
1344
-    y[-tip.idx] <- NA
1345
-
1346
-    ## use lookup table
1347
-    pvec <- integer(max(tr$edge))
1348
-    pvec[child] = parent
1349
-
1350
-    currentNode <- 1:Ntip
1351
-    while(anyNA(y)) {
1352
-        ## pNode <- unique(parent[child %in% currentNode])
1353
-        pNode <- unique(pvec[currentNode])
1354
-
1355
-        ## piping of magrittr is slower than nested function call.
1356
-        ## pipeR is fastest, may consider to use pipeR
1357
-        ##
1358
-        ## child %in% currentNode %>% which %>% parent[.] %>% unique
1359
-        ## idx <- sapply(pNode, function(i) all(child[parent == i] %in% currentNode))
1360
-        idx <- sapply(pNode, function(i) all(child_list[[i]] %in% currentNode))
1361
-        newNode <- pNode[idx]
1362
-
1363
-        y[newNode] <- sapply(newNode, function(i) {
1364
-            mean(y[child_list[[i]]], na.rm=TRUE)
1365
-            ##child[parent == i] %>% y[.] %>% mean(na.rm=TRUE)
1366
-        })
1367
-
1368
-        currentNode <- c(currentNode[!currentNode %in% unlist(child_list[newNode])], newNode)
1369
-        ## currentNode <- c(currentNode[!currentNode %in% child[parent %in% newNode]], newNode)
1370
-        ## parent %in% newNode %>% child[.] %>%
1371
-        ##     `%in%`(currentNode, .) %>% `!` %>%
1372
-        ##         currentNode[.] %>% c(., newNode)
1373
-    }
1374
-
1375
-    return(y)
1376
-}
1377
-
1378
-
1379
-getYcoord_scale <- function(tr, df, yscale) {
1380
-
1381
-    N <- getNodeNum(tr)
1382
-    y <- numeric(N)
1383
-
1384
-    root <- getRoot(tr)
1385
-    y[root] <- 0
1386
-    y[-root] <- NA
1387
-
1388
-    edge <- tr$edge
1389
-    parent <- edge[,1]
1390
-    child <- edge[,2]
1391
-
1392
-    currentNodes <- root
1393
-    while(anyNA(y)) {
1394
-        newNodes <- c()
1395
-        for (currentNode in currentNodes) {
1396
-            idx <- which(parent %in% currentNode)
1397
-            newNode <- child[idx]
1398
-            direction <- -1
1399
-            for (i in seq_along(newNode)) {
1400
-                y[newNode[i]] <- y[currentNode] + df[newNode[i], yscale] * direction
1401
-                direction <- -1 * direction
1402
-            }
1403
-            newNodes <- c(newNodes, newNode)
1404
-        }
1405
-        currentNodes <- unique(newNodes)
1406
-    }
1407
-    if (min(y) < 0) {
1408
-        y <- y + abs(min(y))
1409
-    }
1410
-    return(y)
1411
-}
1412
-
1413
-getYcoord_scale2 <- function(tr, df, yscale) {
1414
-    root <- getRoot(tr)
1415
-
1416
-    pathLength <- sapply(1:length(tr$tip.label), function(i) {
1417
-        get.path_length(tr, i, root, yscale)
1418
-    })
1419
-
1420
-    ordered_tip <- order(pathLength, decreasing = TRUE)
1421
-    ii <- 1
1422
-    ntip <- length(ordered_tip)
1423
-    while(ii < ntip) {
1424
-        sib <- getSibling(tr, ordered_tip[ii])
1425
-        if (length(sib) == 0) {
1426
-            ii <- ii + 1
1427
-            next
1428
-        }
1429
-        jj <- which(ordered_tip %in% sib)
1430
-        if (length(jj) == 0) {
1431
-            ii <- ii + 1
1432
-            next
1433
-        }
1434
-        sib <- ordered_tip[jj]
1435
-        ordered_tip <- ordered_tip[-jj]
1436
-        nn <- length(sib)
1437
-        if (ii < length(ordered_tip)) {
1438
-            ordered_tip <- c(ordered_tip[1:ii],sib, ordered_tip[(ii+1):length(ordered_tip)])
1439
-        } else {
1440
-            ordered_tip <- c(ordered_tip[1:ii],sib)
1441
-        }
1442
-
1443
-        ii <- ii + nn + 1
1444
-    }
1445
-
1446
-
1447
-    long_branch <- getAncestor(tr, ordered_tip[1]) %>% rev
1448
-    long_branch <- c(long_branch, ordered_tip[1])
1449
-
1450
-    N <- getNodeNum(tr)
1451
-    y <- numeric(N)
1452
-
1453
-    y[root] <- 0
1454
-    y[-root] <- NA
1455
-
1456
-    ## yy <- df[, yscale]
1457
-    ## yy[is.na(yy)] <- 0
1458
-
1459
-    for (i in 2:length(long_branch)) {
1460
-        y[long_branch[i]] <- y[long_branch[i-1]] + df[long_branch[i], yscale]
1461
-    }
1462
-
1463
-    parent <- df[, "parent"]
1464
-    child <- df[, "node"]
1465
-
1466
-    currentNodes <- root
1467
-    while(anyNA(y)) {
1468
-        newNodes <- c()
1469
-        for (currentNode in currentNodes) {
1470
-            idx <- which(parent %in% currentNode)
1471
-            newNode <- child[idx]
1472
-            newNode <- c(newNode[! newNode %in% ordered_tip],
1473
-                         rev(ordered_tip[ordered_tip %in% newNode]))
1474
-            direction <- -1
1475
-            for (i in seq_along(newNode)) {
1476
-                if (is.na(y[newNode[i]])) {
1477
-                    y[newNode[i]] <- y[currentNode] + df[newNode[i], yscale] * direction
1478
-                    direction <- -1 * direction
1479
-                }
1480
-            }
1481
-            newNodes <- c(newNodes, newNode)
1482
-        }
1483
-        currentNodes <- unique(newNodes)
1484
-    }
1485
-    if (min(y) < 0) {
1486
-        y <- y + abs(min(y))
1487
-    }
1488
-    return(y)
1489
-}
1490
-
1491
-getYcoord_scale_numeric <- function(tr, df, yscale, ...) {
1492
-    df <- .assign_parent_status(tr, df, yscale)
1493
-    df <- .assign_child_status(tr, df, yscale)
1494
-
1495
-    y <- df[, yscale]
1496
-
1497
-    if (anyNA(y)) {
1498
-        warning("NA found in y scale mapping, all were setting to 0")
1499
-        y[is.na(y)] <- 0
1500
-    }
1501
-
1502
-    return(y)
1503
-}
1504
-
1505
-.assign_parent_status <- function(tr, df, variable) {
1506
-    yy <- df[[variable]]
1507
-    na.idx <- which(is.na(yy))
1508
-    if (length(na.idx) > 0) {
1509
-        tree <- get.tree(tr)
1510
-        nodes <- getNodes_by_postorder(tree)
1511
-        for (curNode in nodes) {
1512
-            children <- getChild(tree, curNode)
1513
-            if (length(children) == 0) {
1514
-                next
1515
-            }
1516
-            idx <- which(is.na(yy[children]))
1517
-            if (length(idx) > 0) {
1518
-                yy[children[idx]] <- yy[curNode]
1519
-            }
1520
-        }
1521
-    }
1522
-    df[, variable] <- yy
1523
-    return(df)
1524
-}
1525
-
1526
-.assign_child_status <- function(tr, df, variable, yscale_mapping=NULL) {
1527
-    yy <- df[[variable]]
1528
-    if (!is.null(yscale_mapping)) {
1529
-        yy <- yscale_mapping[yy]
1530
-    }
1531
-
1532
-    na.idx <- which(is.na(yy))
1533
-    if (length(na.idx) > 0) {
1534
-        tree <- get.tree(tr)
1535
-        nodes <- rev(getNodes_by_postorder(tree))
1536
-        for (curNode in nodes) {
1537
-            parent <- getParent(tree, curNode)
1538
-            if (parent == 0) { ## already reach root
1539
-                next
1540
-            }
1541
-            idx <- which(is.na(yy[parent]))
1542
-            if (length(idx) > 0) {
1543
-                child <- getChild(tree, parent)
1544
-                yy[parent[idx]] <- mean(yy[child], na.rm=TRUE)
1545
-            }
1546
-        }
1547
-    }
1548
-    df[, variable] <- yy
1549
-    return(df)
1550
-}
1551
-
1552
-
1553
-getYcoord_scale_category <- function(tr, df, yscale, yscale_mapping=NULL, ...) {
1554
-    if (is.null(yscale_mapping)) {
1555
-        stop("yscale is category variable, user should provide yscale_mapping,
1556
-             which is a named vector, to convert yscale to numberical values...")
1557
-    }
1558
-    if (! is(yscale_mapping, "numeric") ||
1559
-        is.null(names(yscale_mapping))) {
1560
-        stop("yscale_mapping should be a named numeric vector...")
1561
-    }
1562
-
1563
-    if (yscale == "label") {
1564
-        yy <- df[[yscale]]
1565
-        ii <- which(is.na(yy))
1566
-        if (length(ii)) {
1567
-            df[ii, yscale] <- df[ii, "node"]
1568
-        }
1569
-    }
1570
-
1571
-    ## assign to parent status is more prefer...
1572
-    df <- .assign_parent_status(tr, df, yscale)
1573
-    df <- .assign_child_status(tr, df, yscale, yscale_mapping)
1574
-
1575
-    y <- df[[yscale]]
1576
-
1577
-    if (anyNA(y)) {
1578
-        warning("NA found in y scale mapping, all were setting to 0")
1579
-        y[is.na(y)] <- 0
1580
-    }
1581
-    return(y)
1582
-}
1583
-
1584
-
1585
-add_angle_slanted <- function(res) {
1586
-    x <- res[["x"]]
1587
-    y <- res[["y"]]
1588
-    dy <- (y - y[match(res$parent, res$node)]) / diff(range(y))
1589
-    dx <- (x - x[match(res$parent, res$node)]) / diff(range(x))
1590
-    theta <- atan(dy/dx)
1591
-    theta[is.na(theta)] <- 0 ## root node
1592
-    res$angle <- theta/pi * 180
1593
-
1594
-    branch.y <- (y[match(res$parent, res$node)] + y)/2
1595
-    idx <- is.na(branch.y)
1596
-    branch.y[idx] <- y[idx]
1597
-    res[, "branch.y"] <- branch.y
1598
-    return(res)
1599
-}
1600
-
1601
-calculate_branch_mid <- function(res) {
1602
-    res$branch <- with(res, (x[match(parent, node)] + x)/2)
1603
-    if (!is.null(res$branch.length)) {
1604
-        res$branch.length[is.na(res$branch.length)] <- 0
1605
-    }
1606
-    res$branch[is.na(res$branch)] <- 0
1607
-    return(res)
1608
-}
1609
-
1610
-
1611
-set_branch_length <- function(tree_object, branch.length) {
1612
-    if (branch.length == "branch.length") {
1613
-        return(tree_object)
1614
-    } else if (branch.length == "none") {
1615
-        tree_object@phylo$edge.length <- NULL
1616
-        return(tree_object)
1617
-    }
1618
-
1619
-    if (is(tree_object, "phylo")) {
1620
-        return(tree_object)
1621
-    }
1622
-
1623
-    tree_anno <- get_tree_data(tree_object)
1624
-
1625
-    if (is(tree_anno, "matrix"))
1626
-        tree_anno <- as.data.frame(tree_anno)
1627
-
1628
-    phylo <- get.tree(tree_object)
1629
-
1630
-    cn <- colnames(tree_anno)
1631
-    cn <- cn[!cn %in% c('node', 'parent')]
1632
-
1633
-    length <- match.arg(branch.length, cn)
1634
-
1635
-    if (all(is.na(as.numeric(tree_anno[[length]])))) {
1636
-        stop("branch.length should be numerical attributes...")
1637
-    }
1638
-
1639
-    edge <- as.data.frame(phylo$edge)
1640
-    colnames(edge) <- c("parent", "node")
1641
-
1642
-    dd <- merge(edge, tree_anno,
1643
-                by  = "node",
1644
-                all.x = TRUE)
1645
-    dd <- dd[match(edge$node, dd$node),]
1646
-    len <- unlist(dd[[length]])
1647
-    len <- as.numeric(len)
1648
-    len[is.na(len)] <- 0
1649
-
1650
-    phylo$edge.length <- len
1651
-
1652
-    tree_object@phylo <- phylo
1653
-    return(tree_object)
1654
-}
1655
-
1656
-## set_branch_length <- function(tree_object, branch.length) {
1657
-##     if (is(tree_object, "phylo4d")) {
1658
-##         phylo <- as.phylo.phylo4(tree_object)
1659
-##         d <- tree_object@data
1660
-##         tree_anno <- data.frame(node=rownames(d), d)
1661
-##     } else {
1662
-##         phylo <- get.tree(tree_object)
1663
-##     }
1664
-
1665
-##     if (branch.length %in%  c("branch.length", "none")) {
1666
-##         return(phylo)
1667
-##     }
1668
-
1669
-##     ## if (is(tree_object, "codeml")) {
1670
-##     ##     tree_anno <- tree_object@mlc@dNdS
1671
-##     ## } else
1672
-
1673
-##     if (is(tree_object, "codeml_mlc")) {
1674
-##         tree_anno <- tree_object@dNdS
1675
-##     } else if (is(tree_object, "beast")) {
1676
-##         tree_anno <- tree_object@stats
1677
-##     }
1678
-
1679
-##     if (has.extraInfo(tree_object)) {
1680
-##         tree_anno <- merge(tree_anno, tree_object@extraInfo, by.x="node", by.y="node")
1681
-##     }
1682
-##     cn <- colnames(tree_anno)
1683
-##     cn <- cn[!cn %in% c('node', 'parent')]
1684
-
1685
-##     length <- match.arg(branch.length, cn)
1686
-
1687
-##     if (all(is.na(as.numeric(tree_anno[, length])))) {
1688
-##         stop("branch.length should be numerical attributes...")
1689
-##     }
1690
-
1691
-##     edge <- as.data.frame(phylo$edge)
1692
-##     colnames(edge) <- c("parent", "node")
1693
-
1694
-##     dd <- merge(edge, tree_anno,
1695
-##                 by.x  = "node",
1696
-##                 by.y  = "node",
1697
-##                 all.x = TRUE)
1698
-##     dd <- dd[match(edge$node, dd$node),]
1699
-##     len <- unlist(dd[, length])
1700
-##     len <- as.numeric(len)
1701
-##     len[is.na(len)] <- 0
1702
-
1703
-##     phylo$edge.length <- len
1704
-
1705
-##     return(phylo)
1706
-## }
1707
-
1708
-
1709
-re_assign_ycoord_df <- function(df, currentNode) {
1710
-    while(anyNA(df$y)) {
1711
-        pNode <- with(df, parent[match(currentNode, node)]) %>% unique
1712
-        idx <- sapply(pNode, function(i) with(df, all(node[parent == i & parent != node] %in% currentNode)))
1713
-        newNode <- pNode[idx]
1714
-        ## newNode <- newNode[is.na(df[match(newNode, df$node), "y"])]
1715
-
1716
-        df[match(newNode, df$node), "y"] <- sapply(newNode, function(i) {
1717
-            with(df, mean(y[parent == i], na.rm = TRUE))
1718
-        })
1719
-        traced_node <- as.vector(sapply(newNode, function(i) with(df, node[parent == i])))
1720
-        currentNode <- c(currentNode[! currentNode %in% traced_node], newNode)
1721
-    }
1722
-    return(df)
1723
-}
1724
-
1725
-
1726
-## ##' test whether input object is produced by ggtree function
1727
-## ##'
1728
-## ##'
1729
-## ##' @title is.ggtree
1730
-## ##' @param x object
1731
-## ##' @return TRUE or FALSE
1732
-## ##' @export
1733
-## ##' @author guangchuang yu
1734
-## is.ggtree <- function(x) inherits(x, 'ggtree')
1735
-
1736
-
1737
-calculate_angle <- function(data) {
1738
-    data$angle <- 360/(diff(range(data$y)) + 1) * data$y
1739
-    return(data)
1740
-}