Merge pull request #87 from uclahs-cds/danknight-angle-calculation

Reimplement angle calculation

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: CancerEvolutionVisualization
 Title: Publication Quality Phylogenetic Tree Plots
 Version: 2.0.0
-Date: 2023-10=20
+Date: 2023-11-09
 Authors@R: c(
 	person("Paul Boutros", role = "cre", email = "[email protected]"),
 	person("Adriana Salcedo", role = "aut"),

NEWS

@@ -1,12 +1,18 @@
-CancerEvolutionVisualization 2.0.0 2023-10-20 (Helena Winata, Dan Knight)
+CancerEvolutionVisualization 2.0.0 2023-11-09
 
 ADDED
+* Support for specifying tree angles in either radians or degrees using
+  an optional "angle" column
 * Generic functions to generate accompanying heatmaps
 * Option to specify tree node colours with "node.col" column
 * Option to specify tree node border colour, width, and line-type with 
   "border.col", "border.width", and "border.type" columns
 * Option ot specify tree node label colour with "node.label.col" column
 
+UPDATE
+* Reimplemented tree angle calculations
+* Fixed lopsided radial tree bug
+
 REMOVED
 * "node.col" parameter to SRCGrob. (Node colour only customizable through
   tree input data.frame.)

R/SRCGrob.R

@@ -25,6 +25,7 @@ SRCGrob <- function(
     colour.scheme = CancerEvolutionVisualization::colours,
     draw.nodes = TRUE,
     add.normal = FALSE,
+    use.radians = FALSE,
     normal.cex = 1,
     sig.shape = 3,
     label.nodes = TRUE,
@@ -47,6 +48,7 @@ SRCGrob <- function(
         tree,
         node.text,
         colour.scheme = colour.scheme,
+        use.radians = use.radians,
         default.node.colour = node.col
         );
 

R/add.segs.R

@@ -149,8 +149,6 @@ get.seg.coords <- function(
             }
         );
 
-    tree.segs.adjusted <- tree.segs.adjusted[which(tree.segs.adjusted$basey != tree.segs.adjusted$tipy), ];
-
     if (length(grep('length',colnames(tree))) == 4) {
         second.tree.segs <- tree.segs;
         second.tree.segs$tipy <- second.tree.segs$basey + second.tree.segs$length2.c * cos(second.tree.segs$angle);

R/angles.R

@@ -0,0 +1,105 @@
+calculate.angles.radial <- function(v, tree, spread, total.angle) {
+    root.node.id <- v$id[[1]];
+    node.ids <- c(root.node.id);
+
+    total.angle <- total.angle * spread;
+    angles <- numeric(nrow(tree));
+
+    while (length(node.ids) > 0) {
+        # "Pops" next element in FIFO queue node.ids
+        current.node.id <- node.ids[1];
+        node.ids <- node.ids[-1];
+
+        parent.id <- tree$parent[tree$tip == current.node.id];
+
+        if (parent.id == -1) {
+            tree$angle[tree$tip == current.node.id] <- 0;
+            }
+
+        child.ids <- tree$tip[tree$parent == current.node.id & !is.na(tree$parent)];
+        num.children <- length(child.ids);
+
+        if (num.children > 0) {
+            parent.angle <- tree$angle[tree$tip == current.node.id];
+            child.weight <- assign.weight(current.node.id, v);
+
+            start.angle <- parent.angle - (total.angle) * (num.children > 1) / 2;
+            num.slices <- max(num.children - 1, 1);
+            angle.increment <- total.angle / num.slices;
+
+            for (i in seq_along(child.ids)) {
+                child.id <- child.ids[i];
+                angle <- start.angle + (i - 1) * (angle.increment);
+                angles[tree$tip == child.id] <- angle;
+                }
+
+            # Appending to end of queue for breadth-first traversal
+            node.ids <- append(node.ids, child.ids);
+            }
+        }
+
+    angles <- override.angles(tree, v, angles);
+    return(angles);
+    }
+
+calculate.angles.fixed <- function(v, tree, fixed.angle) {
+    angles <- numeric(nrow(tree));
+    node.ids <- c(v$id[[1]]);
+
+    while (length(node.ids) > 0) {
+        # "Pops" next element in FIFO queue node.ids
+        current.node.id <- node.ids[1];
+        node.ids <- node.ids[-1];
+
+        child.ids <- tree$tip[tree$parent == current.node.id & !is.na(tree$parent)];
+        num.children <- length(child.ids);
+        if (num.children > 0) {
+            # Safe to hardcode temporarily. This will only ever apply to
+            # cases with 0, 1, or 2 children. 3+ will use radial calculation.
+
+            # In future, I would like to remove this fixed angle calculation entirely.
+            # It would be ideal to handle all calculations in the same way, and
+            # rely more on user defined spread and explicit angle overrides.
+            child.angles <- if (num.children == 1) c(0) else c(-1, 1) * fixed.angle;
+
+            for (i in seq_along(child.ids)) {
+                child.id <- child.ids[i];
+                angle <- child.angles[i];
+                angles[tree$tip == child.id] <- angle;
+                }
+            }
+
+        # Appending to end of queue for breadth-first traversal
+        node.ids <- append(node.ids, child.ids);
+        }
+
+    angles <- override.angles(tree, v, angles);
+    return(angles);
+    }
+
+override.angles <- function(tree, v, angles) {
+    if (is.null(v$angle)) {
+        v$angle <- NA;
+        }
+
+    angle.index <- as.list(v$angle);
+
+    names(angle.index) <- v$id;
+    angle.index <- angle.index[!is.na(angle.index)];
+
+    angles <- apply(
+        data.frame(tip = tree$tip, angle = angles),
+        MARGIN = 1,
+        FUN = function(x) {
+            node.id <- as.character(x['tip']);
+            angle.override <- angle.index[[node.id]];
+            if (!is.null(angle.override)) {
+                angle.override <- as.numeric(angle.override);
+                }
+            angle <- if (is.null(angle.override) || is.na(angle.override)) x['angle'] else angle.override;
+            return(as.numeric(angle));
+            }
+        );
+
+    return(angles);
+    }

R/calculate.clone.polygons.R

@@ -330,27 +330,26 @@ compute.clones <- function(
 	root <- v[!is.na(v$parent) & v$parent == -1, ];
 	v <- v[is.na(v$parent) | v$parent != -1, ];
 	v <- rbind(root, v);
+	v <- count.leaves.per.node(v);
+
+	if (no.ccf) {
+	    tree$angle <- if ((is.null(fixed.angle) && nrow(v) > 6) || any(table(v$parent) > 2)) {
+    		tau <- -(pi / 2.5);
+    		vi <- v[v$parent == -1, ];
+    		calculate.angles.radial(v, tree, spread, abs(tau));
+	    } else {
+	        calculate.angles.fixed(v, tree, fixed.angle);
+	        }
+	    tmp <- position.nodes(v, tree, extra.len);
+
+	    clone.env <-  new.env(parent = emptyenv());
+	    clone.env$v <- tmp$v;
+	    clone.env$tree <- tmp$tree;
+
+	    return(clone.env);
+	    }
 
-	if (no.ccf && (is.null(fixed.angle) && nrow(v) > 6) || any(table(v$parent) > 2)) {
-		v <- count.leaves.per.node(v);
-		tmp <-  position.nodes.node.radiusial(v, tree, extra.len, spread);
-		clone.env <-  new.env(parent = emptyenv());
-		clone.env$v <- tmp$v;
-		clone.env$tree <- tmp$tree;
-		return(clone.env);
-	} else if (no.ccf && !is.null(fixed.angle)) {
-		#position nodes fixed angle
-		clone.env <- position.nodes.fixed(
-		    v,
-		    tree,
-		    fixed.angle = fixed.angle,
-		    len = extra.len
-		    );
-
-		return(clone.env)
-	} else{
-		v <- position.clones(v, tree, wid);
-    	}
+	v <- position.clones(v, tree, wid);
 
 	v$x <- v$y <- v$len <- 0;
 	len <- extra.len;

R/position.clones.R

@@ -80,62 +80,6 @@ position.clones <- function(v, tree, wid) {
     return(v);
     }
 
-position.nodes.fixed <- function(v, tree, fixed.angle, len) {
-    for (i in seq_along(v$id)) {
-        vi <- v[i, ];
-
-        if (!is.na(vi$parent) && vi$parent == -1) {
-            # If root the clone extends the full width of the plot
-            x0 <- 0;
-            y0 <- tree$length[tree$parent == -1];
-            len0 <- len + y0;
-        } else {
-            # Parent not root -- not trunk clone
-            par <- v[v$id == vi$parent, ];
-
-            #get parent clone
-            siblings <- v[which(v$parent == par$id),]
-
-            if (nrow(siblings) == 1) {
-                parent.angle <- 0;
-            } else if (nrow(siblings) == 2) {
-                if (any(siblings$x > par$x)) {
-                    parent.angle <- -(fixed.angle);
-                } else {
-                    parent.angle <- fixed.angle;
-                }
-            } else if (nrow(siblings) == 3) {
-                if (any(siblings$x > par$x)) {
-                    parent.angle <- -(fixed.angle);
-                } else if (any(siblings$x < par$x)) {
-                  parent.angle <- fixed.angle;
-                } else {
-                    parent.angle <- 0;
-                    }
-                }
-
-            r <- tree$length[which(tree$parent == par$id & tree$tip == vi$id)];
-            x.shift <- r * sin(parent.angle);
-            x0 <- par$x + x.shift;
-            y.shift <- r * cos(parent.angle);
-            y0 <- par$y + y.shift;
-            len0 <- par$len + y.shift;
-
-            tree$angle[which(tree$parent == par$id & tree$tip == vi$id)] <- parent.angle;
-            }
-
-      v[i,]$len <- len0;
-      v[i,]$y <- y0;
-      v[i,]$x <- x0;
-      }
-
-    clone.env <-  new.env(parent = emptyenv());
-    clone.env$v <- v;
-    clone.env$tree <- tree;
-
-    return(clone.env);
-    }
-
 position.clones.no.vaf <- function(v, wid, spread = TRUE) {
     v$y.mid <- v$y1 <- v$y2 <- 0;
 
@@ -206,3 +150,31 @@ position.clones.no.vaf <- function(v, wid, spread = TRUE) {
 
     return(v);
     }
+
+reposition.clones <- function(tree, v) {
+    for (i in seq_along(v$id)) {
+        vi <- v[i, ];
+        angle <- tree$angle[tree$tip == vi$id];
+
+        if (!is.na(vi$parent) && vi$parent == -1) {
+            x0 <- 0;
+            y0 <- tree$length[tree$parent == -1];
+            len0 <- 0;
+        } else {
+            par <- v[v$id == vi$parent, ];
+
+            r <- tree$length[tree$tip == vi$id];
+            x.shift <- r * sin(angle);
+            x0 <- par$x + x.shift;
+            y.shift <- r * cos(angle);
+            y0 <- par$y + y.shift;
+            len0 <- par$len + y.shift;
+            }
+
+        v[i,]$len <- len0;
+        v[i,]$y <- y0;
+        v[i,]$x <- x0;
+        }
+
+    return(v);
+    }

R/position.nodes.radial.R

@@ -33,60 +33,31 @@ assign.weight <- function(node,v, extra.len, spread) {
 	return(node.weight);
     }
 
-position.nodes.node.radiusial <- function(v, tree, extra.len, spread = 1) {
-	w <- spread * pi;
+position.nodes <- function(v, tree, extra.len) {
 	xpos <- ypos <- 0;
-	tau <- -(pi / 2.5);
 	vi <- v[v$parent == -1, ];
 
-	preorder.traversal <- function(
-	    node = NULL,
-	    tree = NULL,
-	    w = NULL,
-	    tau = NULL,
-	    eta = NULL,
-	    spread = 1
-	    ) {
-
+	preorder.traversal <- function(node, tree) {
 		vi <- v[v$id == node, ];
 		d <- tree$length[tree$tip == vi$id & tree$parent == vi$parent];
 
 		if (vi$parent != -1) {
-			v$x[v$id == vi$id] <<- v$x[v$id == vi$parent] + d * sin(tau + w / 2);
-			v$y[v$id == vi$id] <<- v$y[v$id == vi$parent] + d * cos(tau + w / 2);
-			tree$angle[tree$tip == vi$id & tree$parent == vi$parent] <<- tau + w / 2;
+		    angle <- tree$angle[tree$tip == vi$id & tree$parent == vi$parent];
+			v$x[v$id == vi$id] <<- v$x[v$id == vi$parent] + d * sin(angle);
+			v$y[v$id == vi$id] <<- v$y[v$id == vi$parent] + d * cos(angle);
 		} else {
-			v$x[v$id == vi$id] <<- 0;
+		    v$x[v$id == vi$id] <<- 0;
 			v$y[v$id == vi$id] <<- d;
-			tree$angle[tree$tip == vi$id & tree$parent == vi$parent] <<- 0;
 		    }
 
-		eta <- tau;
-
 		for (child in v$id[v$parent == vi$id]) {
-			child.weight <- assign.weight(child, v);
-			w <- child.weight * spread * pi;
-			tau <- eta;
-			eta <- eta + w;
-
-			preorder.traversal(
-			    node = child,
-			    tree = tree,
-			    w = w,
-			    tau = tau,
-			    eta = eta,
-			    spread = spread
-			    );
+			preorder.traversal(node = child, tree = tree);
 		    }
 	     }
 
-	preorder.traversal(
-	    node = 1,
-	    tree = tree,
-	    w = w,
-	    tau = tau,
-	    spread = spread
-	    );
+	preorder.traversal(node = 1, tree = tree);
+
+	v <- reposition.clones(tree, v);
 
 	v$len <- sapply(
 	    v$y,