Muller fish plots to OncoSimul

OncoSimulR/DESCRIPTION

@@ -29,7 +29,7 @@ License: GPL (>= 3)
 URL: https://github.com/rdiaz02/OncoSimul, https://popmodels.cancercontrol.cancer.gov/gsr/packages/oncosimulr/
 BugReports: https://github.com/rdiaz02/OncoSimul/issues
 Depends: R (>= 3.3.0)
-Imports: Rcpp (>= 0.12.4), parallel, data.table, graph, Rgraphviz, gtools, igraph, methods, RColorBrewer, grDevices, car, dplyr, smatr, ggplot2, ggrepel, nem
+Imports: Rcpp (>= 0.12.4), parallel, data.table, graph, Rgraphviz, gtools, igraph, methods, RColorBrewer, grDevices, car, dplyr, smatr, ggplot2, ggrepel, nem, ggmuller
 Suggests: BiocStyle, knitr, Oncotree, testthat (>= 1.0.0), rmarkdown, bookdown, pander
 LinkingTo: Rcpp
 VignetteBuilder: knitr

OncoSimulR/NAMESPACE

@@ -64,6 +64,7 @@ importFrom("dplyr", "full_join", "left_join", "right_join", "%>%", "mutate",
 importFrom("smatr", "ma") ## for major axis regression in some tests
 importFrom("car", "linearHypothesis")
 importFrom("nem", "transitive.reduction")
+importFrom("ggmuller", "get_Muller_df", "Muller_plot")
 ## importFrom("slam", "simple_triplet_zero_matrix", ## "colapply_simple_triplet_matrix",
 ##            "col_sums")
 

OncoSimulR/R/OncoSimulR.R

@@ -1006,6 +1006,54 @@ plot.oncosimulpop <- function(x, ask = TRUE,
 
 ## }
 
+plot_muller <- function(onco_simul_obj, 
+                        wild_type_legend = "wt", 
+                        cutoff = 0.0,
+                        xlab = "Time", 
+                        ylab = "Frequency",
+                        ...) {
+
+    # Testing input:
+    if(!("oncosimul2" %in% class(onco_simul_obj))) {
+      stop("Input must be an oncosimul2 object")
+    }
+
+    if(class(wild_type_legend) != "character"){
+      stop("wild_type_legend must be a character")
+    }
+
+    # Fill adjacency matrix
+    adjacency_matrix <- data.frame(onco_simul_obj$other$LOD_DF)
+    levels(adjacency_matrix$parent)[1] <- wild_type_legend
+
+    population_evolution <- as.data.frame(onco_simul_obj$pops.by.time)
+
+    onco_simul_obj$GenotypesLabels[1] <- wild_type_legend
+    genotypes_labels <- c("Time", onco_simul_obj$GenotypesLabels)
+    colnames(population_evolution) <- genotypes_labels
+
+    # Fill generation vector
+    Generation <- rep(population_evolution[ ,1], length(colnames(population_evolution))-1)
+
+    # Fill population vector
+    Population <- c()
+    for(i in 2:length(colnames(population_evolution))) {
+      Population <- c(Population, population_evolution[ ,i])
+    }
+    Population <- as.numeric(Population)
+
+    # Fill identity vector
+    Identity <- c()
+    for(i in 2:length(colnames(population_evolution))) {
+      Identity <- c(Identity, rep(colnames(population_evolution)[i], 
+                    length(rownames(population_evolution))))
+    }
+
+    muller_df <- get_Muller_df(adjacency_matrix, 
+                               data.frame(Identity, Population, Generation),
+                               cutoff)
+    Muller_plot(muller_df, xlab = xlab, ylab = ylab, ...)
+}
 
 plot.oncosimul <- function(x,
                            show = "drivers", 
@@ -1017,8 +1065,8 @@ plot.oncosimul <- function(x,
                            lwdClone = 0.9,
                            ltyDrivers = 1,
                            lwdDrivers = 3,
-                           xlab = "Time units",
-                           ylab = "Number of cells",
+                           xlab = ifelse(type == "fish", "Generation", "Time units"),
+                           ylab = ifelse(type == "fish", "Genotype Frequency", "Number of cells"),
                            plotClones = TRUE,
                            plotDrivers = TRUE,
                            addtot = FALSE,
@@ -1039,133 +1087,149 @@ plot.oncosimul <- function(x,
                            vrange = c(0.8, 1),
                            breakSortColors = "oe",
                            legend.ncols = "auto",
+                           fish_wild_type_legend = "wt",
+                           fish_cutoff = 0.0,
                            ...
                            ) {
 
 
-    if(!(type %in% c("stacked", "stream", "line")))
-        stop("Type of plot unknown: it must be one of",
-             "stacked, stream or line")
-
-    if(!(show %in% c("genotypes", "drivers")))
-        stop("show must be one of ",
-             "genotypes or drivers")
-
-    if(!(breakSortColors %in% c("oe", "distave", "random")))
-        stop("breakSortColors must be one of ",
-             "oe, distave, or random")
-
+    if(!(type %in% c("stacked", "stream", "line", "fish")))
+        stop("Type of plot unknown: it must be one of ",
+             "stacked, stream, line or fish")
+
+    if(type == "fish") {
+        if(thinData)
+            x <- thin.pop.data(x, keep = thinData.keep, min.keep = thinData.min)
+
+        if(!is.null(xlim))
+            x <- xlim.pop.data(x, xlim)
+
+        plot_muller(x,
+                    wild_type_legend = fish_wild_type_legend, 
+                    cutoff = fish_cutoff,
+                    xlab = xlab, 
+                    ylab = ylab,
+                    ...)
+    } else {
+        if(!(show %in% c("genotypes", "drivers")))
+            stop("show must be one of ",
+                 "genotypes or drivers")
 
-
-    colauto <- FALSE
-    if(col == "auto" && (type == "line") && (show == "drivers"))
-        col <- c(8, "orange", 6:1)
-    if(col == "auto" && (show == "genotypes")) {
-        ## For categorical data, I find Dark2, Paired, or Set1 to work best.
-        col <- colorRampPalette(brewer.pal(8, "Dark2"))(ncol(x$pops.by.time) - 1)
-        colauto <- TRUE
-    }
+        if(!(breakSortColors %in% c("oe", "distave", "random")))
+            stop("breakSortColors must be one of ",
+                 "oe, distave, or random")
 
-    if(show == "genotypes") {
-        plotDrivers <- FALSE
-        plotClones <- TRUE
-    }
 
-    if(thinData)
-        x <- thin.pop.data(x, keep = thinData.keep, min.keep = thinData.min)
-
-    if(!is.null(xlim))
-        x <- xlim.pop.data(x, xlim)
+        colauto <- FALSE
+        if(col == "auto" && (type == "line") && (show == "drivers"))
+            col <- c(8, "orange", 6:1)
+        if(col == "auto" && (show == "genotypes")) {
+            ## For categorical data, I find Dark2, Paired, or Set1 to work best.
+            col <- colorRampPalette(brewer.pal(8, "Dark2"))(ncol(x$pops.by.time) - 1)
+            colauto <- TRUE
+        }
+
+        if(show == "genotypes") {
+            plotDrivers <- FALSE
+            plotClones <- TRUE
+        }
+
+        if(thinData)
+            x <- thin.pop.data(x, keep = thinData.keep, min.keep = thinData.min)
 
-    ## For genotypes, ndr is now the genotypes.  Actually, ndr is now just
-    ## a sequence 1:(ncol(y) - 1)
-
-    ## The user will want to change the colors, like a colorRamp, etc. Or
-    ## rainbow.
-
-    ## genotypes and line, always call plotDrivers0
-    if(show == "drivers") {
-        if(!inherits(x, "oncosimul2"))
-            ndr <- colSums(x$Genotypes[1:x$NumDrivers, , drop = FALSE])
-        else {
-            ndr <- colSums(x$Genotypes[x$Drivers, , drop = FALSE])
+        if(!is.null(xlim))
+            x <- xlim.pop.data(x, xlim)
+
+        ## For genotypes, ndr is now the genotypes.  Actually, ndr is now just
+        ## a sequence 1:(ncol(y) - 1)
+
+        ## The user will want to change the colors, like a colorRamp, etc. Or
+        ## rainbow.
+
+        ## genotypes and line, always call plotDrivers0
+        if(show == "drivers") {
+            if(!inherits(x, "oncosimul2"))
+                ndr <- colSums(x$Genotypes[1:x$NumDrivers, , drop = FALSE])
+            else {
+                ndr <- colSums(x$Genotypes[x$Drivers, , drop = FALSE])
+            }
+        } else { ## show we are showing genotypes
+            ndr <- 1:(ncol(x$pops.by.time) - 1)
+        }
+
+        if((type == "line") && is.null(ylim)) {
+            if(log %in% c("y", "xy", "yx") )
+                ylim <- c(1, max(apply(x$pops.by.time[, -1, drop = FALSE], 1, sum)))
+            else
+                ylim <- c(0, max(apply(x$pops.by.time[, -1, drop = FALSE], 1, sum)))
+        }
+        if(plotDiversity) {
+            oppd <- par(fig = c(0, 1, 0.8, 1))
+            m1 <- par()$mar
+            m <- m1
+            m[c(1, 3)] <- c(0, 0.7)
+            op <- par(mar = m )
+            plotShannon(x)
+            par(op)
+            m1[c(3)] <- 0.2
+            op <- par(mar = m1)
+            par(fig = c(0, 1, 0, 0.8), new = TRUE)  
         }
-    } else { ## show we are showing genotypes
-        ndr <- 1:(ncol(x$pops.by.time) - 1)
-    }
 
-    if((type == "line") && is.null(ylim)) {
-        if(log %in% c("y", "xy", "yx") )
-            ylim <- c(1, max(apply(x$pops.by.time[, -1, drop = FALSE], 1, sum)))
-        else
-            ylim <- c(0, max(apply(x$pops.by.time[, -1, drop = FALSE], 1, sum)))
-    }
-    if(plotDiversity) {
-        oppd <- par(fig = c(0, 1, 0.8, 1))
-        m1 <- par()$mar
-        m <- m1
-        m[c(1, 3)] <- c(0, 0.7)
-        op <- par(mar = m )
-        plotShannon(x)
-        par(op)
-        m1[c(3)] <- 0.2
-        op <- par(mar = m1)
-        par(fig = c(0, 1, 0, 0.8), new = TRUE)  
-    }
-
-    ## Shows its history: plotClones makes plotDrivers0 unneeded with
-    ## stacked and stream plots. But now so with line plot.
-    ## When showing genotypes, plotDrivers0 with line only used for
-    ## showing the legend.
-    if(plotClones) {
-        plotClonesSt(x,
-                     ndr = ndr,
-                     show = show,
-                     na.subs = TRUE,
-                     log = log,
-                     lwd = lwdClone,
-                     lty = ifelse(show == "drivers", ltyClone, ltyDrivers),
-                     col = col, 
-                     order.method = order.method,
-                     stream.center = stream.center,
-                     stream.frac.rand = stream.frac.rand,
-                     stream.spar = stream.spar,
-                     border = border,
-                     srange = srange,
-                     vrange = vrange,
-                     type = type,
-                     breakSortColors = breakSortColors,
-                     colauto = colauto,
-                     legend.ncols = legend.ncols,
-                     lwdStackedStream = lwdStackedStream,
-                     xlab = xlab,
-                     ylab = ylab,
-                     ylim = ylim,
-                     xlim = xlim,
-                     ...)
-    }
-
-    if(plotClones && plotDrivers && (type == "line"))
-        par(new = TRUE)
+        ## Shows its history: plotClones makes plotDrivers0 unneeded with
+        ## stacked and stream plots. But now so with line plot.
+        ## When showing genotypes, plotDrivers0 with line only used for
+        ## showing the legend.
+        if(plotClones) {
+            plotClonesSt(x,
+                         ndr = ndr,
+                         show = show,
+                         na.subs = TRUE,
+                         log = log,
+                         lwd = lwdClone,
+                         lty = ifelse(show == "drivers", ltyClone, ltyDrivers),
+                         col = col, 
+                         order.method = order.method,
+                         stream.center = stream.center,
+                         stream.frac.rand = stream.frac.rand,
+                         stream.spar = stream.spar,
+                         border = border,
+                         srange = srange,
+                         vrange = vrange,
+                         type = type,
+                         breakSortColors = breakSortColors,
+                         colauto = colauto,
+                         legend.ncols = legend.ncols,
+                         lwdStackedStream = lwdStackedStream,
+                         xlab = xlab,
+                         ylab = ylab,
+                         ylim = ylim,
+                         xlim = xlim,
+                         ...)
+        }
 
-    if( plotDrivers && (type == "line") ) {
-        plotDrivers0(x,
-                     ndr,
-                     timescale = 1,
-                     trim.no.drivers = FALSE,
-                     xlab = "", ylab = "",
-                     lwd = lwdDrivers,
-                     lty = ltyDrivers,
-                     col = col, 
-                     addtot = addtot,
-                     addtotlwd = addtotlwd,
-                     log = log, ylim = ylim,
-                     xlim = xlim,
-                     legend.ncols = legend.ncols,
-                     ...)
-    }
-    if(plotDiversity) {
-        par(oppd)
+        if(plotClones && plotDrivers && (type == "line"))
+            par(new = TRUE)
+
+        if( plotDrivers && (type == "line") ) {
+            plotDrivers0(x,
+                         ndr,
+                         timescale = 1,
+                         trim.no.drivers = FALSE,
+                         xlab = "", ylab = "",
+                         lwd = lwdDrivers,
+                         lty = ltyDrivers,
+                         col = col, 
+                         addtot = addtot,
+                         addtotlwd = addtotlwd,
+                         log = log, ylim = ylim,
+                         xlim = xlim,
+                         legend.ncols = legend.ncols,
+                         ...)
+        }
+        if(plotDiversity) {
+            par(oppd)
+        }
     }
 
 }