Merge reworked dynplot2

.Rbuildignore

@@ -9,3 +9,5 @@
 ^\.readme_files
 ^README
 ^tests
+^LICENSE\.md$
+^data-raw

DESCRIPTION

@@ -1,53 +1,61 @@
 Package: dynplot2
 Type: Package
-Title: Grammar of graphics for single-cell data
-Version: 1.0.0
+Title: Modular trajectory plotting
+Version: 1.1.0
 Authors@R: c(
   person("Robrecht", "Cannoodt", email = "[email protected]", role = c("aut")),
   person("Wouter", "Saelens", email = "[email protected]", role = c("aut", "cre")))
-Description: Package that create ggplot2 plots of single-cell datasets.
-License: GPL-3
+Description: Finetune visualisations of trajectories using ggplot2's grammar of graphics principles.
+License: MIT + file LICENSE
 Encoding: UTF-8
 LazyData: true
-RoxygenNote: 6.1.1
+RoxygenNote: 7.1.1
 SystemRequirements: libudunits2-dev (deb) or udunits2-devel (rpm)
+URL: https://github.com/dynverse/dynplot2
+BugReports: https://github.com/dynverse/dynplot2/issues
 Depends:
   R (>= 3.0.0)
 Imports:
   assertthat,
   cowplot,
   dplyr,
+  digest,
   dynutils (>= 1.0.2),
   dynfeature (>= 1.0.0),
   dyndimred (>= 1.0.0),
   dynwrap (>= 1.0.0),
   GA,
   ggbeeswarm,
+  ggnewscale,
   ggplot2 (>= 3.0),
   ggraph,
   ggrepel,
+  grDevices,
   igraph,
   MASS,
   magrittr,
+  Matrix,
   methods,
   purrr,
-  rje,
-  shades,
+  reticulate,
+  rlang,
   RColorBrewer,
   reshape2,
+  scvelo,
+  shades,
+  stats,
+  stringr,
   testthat,
   tibble,
   tidyr,
   tidygraph,
+  utils,
   vipor,
   viridis
 Suggests:
   knitr,
   rmarkdown
-Remotes:  
-    dynverse/dynutils@master,
-    dynverse/dyndimred@master,
-    dynverse/dynwrap@master,
-    dynverse/dynfeature@master
+Remotes:
+  dynverse/scvelo@devel
 VignetteBuilder: knitr
 Roxygen: list(markdown = TRUE)

LICENSE

@@ -0,0 +1,2 @@
+YEAR: 2014-2020
+COPYRIGHT HOLDER: Robrecht Cannoodt, Wouter Saelens

LICENSE.md

@@ -0,0 +1,21 @@
+# MIT License
+
+Copyright (c) 2014-2020 Robrecht Cannoodt, Wouter Saelens
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

NAMESPACE

@@ -1,36 +1,55 @@
 # Generated by roxygen2: do not edit by hand
 
 S3method(ggplot_build,dynplot)
+export(color_milestone_percentages)
+export(define_milestone_colors)
 export(dynplot)
+export(dynplot_dendro)
+export(dynplot_dimred)
+export(dynplot_graph)
+export(dynplot_onedim)
 export(empty_plot)
+export(facet_grid_data)
+export(facet_wrap_data)
 export(geom_cell_hex)
 export(geom_cell_point)
-export(geom_expression_raster)
+export(geom_expression_tile)
+export(geom_feature_text)
 export(geom_milestone_label)
 export(geom_milestone_point)
 export(geom_trajectory_connection)
 export(geom_trajectory_divergence)
 export(geom_trajectory_segments)
 export(geom_velocity_arrow)
+export(geom_velocity_stream)
 export(get_features)
 export(get_milestone_palette_names)
-export(layout_dendro)
-export(layout_dimred)
-export(layout_graph)
 export(layout_heatmap)
 export(layout_modules)
-export(layout_onedim)
+export(milestone_palette)
+export(milestone_percentage_breaks)
 export(new_scale)
-export(new_scale_fillcolour)
-export(scale_expression_fillcolour)
-export(scale_milestones_fillcolour)
-export(select_expression)
-export(shadow)
+export(new_scale_color)
+export(new_scale_colour)
+export(new_scale_fill)
+export(scale_expression_color)
+export(scale_expression_colour)
+export(scale_expression_fill)
+export(scale_milestones_color)
+export(scale_milestones_colour)
+export(scale_milestones_fill)
+export(scale_velocity_color)
+export(scale_velocity_fill)
+export(scale_x_heatmap)
+export(select_feature_expression)
+export(select_feature_velocity)
+export(shadow_defaults)
 export(stat_velocity_cells)
 export(stat_velocity_grid)
 export(stat_velocity_random)
+export(stat_velocity_stream)
 export(theme_clean)
-export(theme_graph)
+export(theme_dynplot)
 import(dplyr)
 import(dyndimred)
 import(dynutils)
@@ -41,18 +60,29 @@ import(purrr)
 import(tibble)
 import(tidyr)
 importFrom(GA,ga)
+importFrom(Matrix,t)
 importFrom(RColorBrewer,brewer.pal)
 importFrom(assertthat,assert_that)
+importFrom(grDevices,col2rgb)
 importFrom(grDevices,rainbow)
+importFrom(grDevices,rgb)
 importFrom(magrittr,"%$%")
 importFrom(magrittr,"%<>%")
 importFrom(magrittr,set_colnames)
 importFrom(magrittr,set_rownames)
-importFrom(rje,cubeHelix)
+importFrom(rlang,"%|%")
+importFrom(rlang,eval_tidy)
+importFrom(rlang,quo)
+importFrom(rlang,quo_text)
 importFrom(shades,hue)
 importFrom(stats,as.dendrogram)
 importFrom(stats,as.dist)
+importFrom(stats,cor)
 importFrom(stats,dist)
 importFrom(stats,hclust)
+importFrom(stats,kmeans)
 importFrom(stats,runif)
-importFrom(stringr,str_detect)
+importFrom(stringr,str_glue)
+importFrom(stringr,str_split)
+importFrom(utils,head)
+importFrom(utils,tail)

inst/NEWS.md → NEWS.md

@@ -1,3 +1,7 @@
+# dynplot2 1.0.0 (unreleased)
+
+ * Spin-off from dynplot, dynplot2 provides a more modular interface to visualize trajectories
+
 # dynplot 1.0.0 (28-03-2019)
 
  * Initial release of dynplot

R/dynplot.R

@@ -1,62 +1,97 @@
 #' Create a dynplot with a specified layout
 #'
-#' @param dataset A dynwrap dataset object, typically containing a trajectory
-#' @param layout A `layout_*` function from dynplot, such as [layout_dimred()] or [layout_graph()]
+#' It's recommended you use one of the layout-specific dynplot functions such as
+#' [dynplot_dendro()], [dynplot_dimred()], [dynplot_graph()] or [dynplot_onedim()]
+#' to start creating a new dynplot visualisation. Only use this function if you
+#' know how dynplot layouts work.
 #'
-#' @return A ggplot2 object, with the processed data in `plot$data` and `attr(plot$data, "data")`
+#' @param dataset A dynwrap dataset object, typically containing a trajectory.
+#' @param trajectory The trajectory dynwrap object, if available.
+#' @param layout A layout list object. Can contain data frames `cell_positions`, `milestone_positions`, `edge_positions`, `segment_positions`, `segment_progressions`.
+#'
+#' @return A ggplot2 object, with the processed data in `plot$data` and `attr(plot$data, "data")`.
+#'
+#' @seealso [dynplot_dendro()], [dynplot_dimred()], [dynplot_graph()] or [dynplot_onedim()]
 #'
 #' @examples
-#' toy <- dyntoy::generate_dataset()
-#' dynplot(dataset) +
-#'   geom_cell_point(aes(colour = select_expression("G1"))) +
-#'   scale_expression_fillcolour() +
-#'   new_scale_fillcolour() +
+#' library(ggplot2)
+#' dynplot_dimred(example_bifurcating) +
+#'   geom_cell_point(aes(colour = select_feature_expression("G1", d = .data))) +
+#'   scale_expression_colour() +
+#'   new_scale_colour() +
 #'   geom_trajectory_segments(aes(colour = edge_id))
 #'
-#' @import dynwrap
-#' @import ggplot2
-#'
 #' @export
 dynplot <- function(
   dataset,
-  layout = layout_dimred(dataset)
+  trajectory,
+  layout
 ) {
   data <- list(dataset = dataset)
 
   # cell info ---------------------------------------------------------------
-  milestone_id_levels <- dataset$milestone_ids
-
   cell_info <-
     bind_cols(
       dataset$cell_info %||% tibble(cell_id = dataset$cell_ids)
     ) %>%
     left_join(layout$cell_positions, "cell_id")
 
+  # add trajectory cell info
+  if (!is.null(trajectory$cell_info)) {
+    cell_info <- left_join(
+      cell_info,
+      trajectory$cell_info[c("cell_id", setdiff(colnames(trajectory$cell_info), colnames(cell_info)))],
+      "cell_id"
+    )
+  }
+
+  # add pseudotime info if available
+  if (!is.null(trajectory$pseudotime)) {
+    cell_info <- left_join(
+      cell_info,
+      enframe(trajectory$pseudotime, "cell_id", "pseudotime"),
+      "cell_id"
+    )
+  }
+
+  # grouping ----------------------------------------------------------------
+  if ("grouping" %in% names(dataset)) {
+    cell_info <- left_join(
+      cell_info,
+      dataset$grouping %>% enframe("cell_id", "group_id"),
+      "cell_id"
+    )
+  }
+
   # trajectory --------------------------------------------------------------
-  if (dynwrap::is_wrapper_with_trajectory(dataset)) {
+  if (dynwrap::is_wrapper_with_trajectory(trajectory)) {
+    milestone_id_levels <- trajectory$milestone_ids
+
     # add milestone percentages to cell info
-    cell_info_milestone_percentages <- dataset$milestone_percentages %>%
-      mutate(milestone_id = factor(milestone_id, milestone_id_levels)) %>%
-      nest(-cell_id, .key = "milestone_percentages") %>%
+    cell_info_milestone_percentages <-
+      trajectory$milestone_percentages %>%
+      mutate(milestone_id = factor(.data$milestone_id, milestone_id_levels)) %>%
+      nest(milestone_percentages = c(.data$milestone_id, .data$percentage)) %>%
       deframe()
     cell_info$milestone_percentages <- unname(cell_info_milestone_percentages[cell_info$cell_id])
 
     # milestone info
     milestone_info <- tibble(
-      milestone_id = dataset$milestone_ids
+      milestone_id = trajectory$milestone_ids
     ) %>%
       mutate(
-        label = milestone_id
+        labelling = trajectory$milestone_labelling[.data$milestone_id] %||% NA_character_,
+        label = .data$labelling %|% .data$milestone_id
       ) %>%
       left_join(layout$milestone_positions, "milestone_id") %>%
       mutate(milestone_id = factor(milestone_id, milestone_id_levels))
 
     # milestone network
     edge_info <-
       bind_cols(
-        dataset$milestone_network
+        trajectory$milestone_network
       ) %>%
-      mutate(edge_id = paste0(from, "->", to), label = edge_id) %>%
+      mutate(edge_id = paste0(.data$from, "->", .data$to), label = .data$edge_id) %>%
       left_join(layout$edge_positions, c("from", "to"))
 
     data <- c(data, lst(
@@ -67,23 +102,24 @@ dynplot <- function(
     if ("segment_progressions" %in% names(layout)) {
       # segment info (produced by layout)
       segment_info <- layout$segment_progressions %>%
-        mutate(edge_id = paste0(from, "->", to)) %>%
-        arrange(edge_id, percentage) %>%
-        left_join(layout$segment_positions, "point_id")
+        mutate(edge_id = paste0(.data$from, "->", .data$to)) %>%
+        arrange(.data$edge_id, .data$percentage) %>%
+        left_join(layout$segment_positions, "point_id") %>%
+        left_join(trajectory$milestone_network %>% select(-.data$length), c("from", "to"))
 
       # get milestone percentages of segments from progressions
       segment_milestone_percentages <- convert_progressions_to_milestone_percentages(
-        cell_ids = dataset$cell_ids,
-        milestone_ids = dataset$milestone_ids,
-        milestone_network = dataset$milestone_network,
-        progressions = segment_info %>% mutate(cell_id = point_id)
+        cell_ids = trajectory$cell_ids,
+        milestone_ids = trajectory$milestone_ids,
+        milestone_network = trajectory$milestone_network,
+        progressions = segment_info %>% mutate(cell_id = .data$point_id)
       ) %>%
         mutate(
-          point_id = cell_id,
-          milestone_id = factor(milestone_id, milestone_id_levels)
+          point_id = .data$cell_id,
+          milestone_id = factor(.data$milestone_id, milestone_id_levels)
         ) %>%
-        select(-cell_id) %>%
-        nest(-point_id, .key = "milestone_percentages") %>%
+        select(-.data$cell_id) %>%
+        nest(milestone_percentages = c(.data$milestone_id, .data$percentage)) %>%
         deframe()
       segment_info$milestone_percentages <- segment_milestone_percentages[segment_info$point_id]
 
@@ -103,18 +139,25 @@ dynplot <- function(
 
   # features ----------------------------------------------------------------
   if ("feature_positions" %in% names(layout)) {
-    data$feature_info <- layout$feature_positions
+    data$feature_info <- left_join(
+      layout$feature_positions,
+      dataset$feature_info,
+      "feature_id"
+    )
+  } else {
+    data$feature_info <- dataset$feature_info
   }
 
+
   # finalise ----------------------------------------------------------------
   data$cell_info <- cell_info
   attr(cell_info, "data") <- data
 
 
   # plot --------------------------------------------------------------------
-  envir <- parent.frame()
-  p <- ggplot(data = cell_info, environment = envir) +
-    theme_graph()
+  p <-
+    ggplot(data = cell_info) +
+    theme_dynplot()
   class(p) <- c("dynplot", class(p))
   p
 }
@@ -128,7 +171,6 @@ aesIntersect <- function(aes1, aes2) {
 
 
 
-#' @importFrom stringr str_detect
 #' @export
 ggplot_build.dynplot <- function(plot) {
   # do some checks for aesthetics
@@ -142,7 +184,7 @@ ggplot_build.dynplot <- function(plot) {
     names()
 
   milestone_percentage_aesthetics_covered <- plot$scales$scales %>%
-    keep(~any(str_detect(class(.), "^ScaleMilestone"))) %>%
+    keep(~any(grepl("^ScaleMilestone", class(.)))) %>%
     map(~.$aesthetics) %>%
     unlist()
 

R/layout_dendro.R → R/dynplot_dendro.R

@@ -1,12 +1,39 @@
-#' Dendrogram layout of a trajectory
+#' Plot a trajectory as a dendrogram
 #'
-#' @inheritParams dynwrap::common_param
-#' @param diag_offset The x-offset (percentage of the edge lenghts) between milestones
+#' @inheritParams dynplot
+#' @param diag_offset The x-offset (percentage of the edge lengths) between milestones
+#' @param y_offset The y-offset (percentage of the edge lengths) between milestones
 #'
-#' @keywords plot_trajectory
+#' @return A ggplot2 object, with the processed data in `plot$data` and `attr(plot$data, "data")`.
+#'
+#' @seealso [dynplot_dendro()], [dynplot_dimred()], [dynplot_graph()] or [dynplot_onedim()]
 #'
 #' @export
-layout_dendro <- function(trajectory, diag_offset = 0.05) {
+#' @examples
+#' library(ggplot2)
+#' data(example_bifurcating)
+#' dataset <- example_bifurcating
+#'
+#' dynplot_dendro(dataset) +
+#'   geom_trajectory_segments() +
+#'   geom_trajectory_connection() +
+#'   geom_cell_point(size = 2, colour = "black") +
+#'   geom_cell_point(aes(colour = milestone_percentages), size = 1.8) +
+#'   geom_milestone_label(aes(fill = milestone_id)) +
+#'   scale_milestones_fill() +
+#'   scale_milestones_colour()
+#'
+#' dynplot_dendro(dataset) +
+#'   geom_trajectory_segments(aes(color = milestone_percentages), size = 5, color = "#333333") +
+#'   geom_trajectory_connection() +
+#'   scale_milestones_colour() +
+#'   new_scale_colour() +
+#'   geom_cell_point(aes(colour = select_feature_expression("G1", d = .data))) +
+#'   scale_expression_colour() +
+#'   new_scale_colour() +
+#'   geom_milestone_label(aes(fill = milestone_id)) +
+#'   scale_milestones_fill()
+dynplot_dendro <- function(dataset, trajectory = dataset, diag_offset = 0.05, y_offset = 0.2) {
   # root if necessary
   if ("root_milestone_id" %in% names(trajectory)) {
     root <- trajectory$root_milestone_id
@@ -46,7 +73,8 @@ layout_dendro <- function(trajectory, diag_offset = 0.05) {
         comp_1 = milestone_positions$comp_1[match(milestone_network_to$from, milestone_positions$node_id)] + milestone_network_to$length + diag_offset,
         comp_2 = map_dbl(milestone_network_to$to, ~mean(leaves_comp_2[descendants[[.]]])),
         parent_node_id = milestone_network_to$from,
-        edge_id = milestone_network_to$edge_id
+        edge_id = milestone_network_to$edge_id,
+        milestone_id = milestone_network_to$to
       )
     )
 
@@ -65,13 +93,14 @@ layout_dendro <- function(trajectory, diag_offset = 0.05) {
     mutate(
       child_node_id = node_id,
       comp_1 = milestone_positions_to$comp_1[match(parent_node_id, milestone_positions_to$node_id)] + diag_offset,
+      milestone_id = parent_node_id,
       node_id = paste0(parent_node_id, "-", node_id)
     )
 
   # combine positions
   milestone_positions <- bind_rows(
-    milestone_positions_to %>% mutate(node_type = "milestone"),
-    milestone_positions_from %>% mutate(node_type = "fake_milestone")
+    milestone_positions_to %>% mutate(node_type = "milestone") %>% mutate(hjust = 0),
+    milestone_positions_from %>% mutate(node_type = "fake_milestone") %>% mutate(hjust = 1)
   )
 
   # now generate network between milestones
@@ -86,6 +115,8 @@ layout_dendro <- function(trajectory, diag_offset = 0.05) {
   connection_positions <- tibble(
     from = milestone_positions_from$parent_node_id,
     to = milestone_positions_from$node_id
+  ) %>% filter(
+    from != root
   ) %>%
     left_join(
       milestone_positions %>% select(node_id, comp_1, comp_2) %>% rename_all(~paste0(., "_from")),
@@ -125,10 +156,16 @@ layout_dendro <- function(trajectory, diag_offset = 0.05) {
       comp_2 = comp_2_from
     )
 
+  # add quasirandom
+  cell_positions <- cell_positions %>%
+    mutate(
+      comp_2 = comp_2 + vipor::offsetX(comp_1, edge_id, method = "quasirandom", width = y_offset)
+    )
+
   # clean up milestone positions & edges
   milestone_positions <- milestone_positions %>%
-    filter(node_type == "milestone") %>%
-    select(milestone_id = node_id, comp_1, comp_2)
+    # filter(node_type == "milestone") %>%
+    select(milestone_id, comp_1, comp_2, hjust)
 
   edge_positions <- edge_positions %>%
     select(from, to, comp_1_from, comp_2_from, comp_1_to, comp_2_to)
@@ -139,7 +176,7 @@ layout_dendro <- function(trajectory, diag_offset = 0.05) {
   segment_progressions <- segments$segment_progressions
   segment_positions <- segments$segment_positions
 
-  lst(
+  layout <- lst(
     milestone_positions = milestone_positions %>% rename_dimred_xy(),
     edge_positions = edge_positions %>% rename_dimred_xy(),
     cell_positions = cell_positions %>% rename_dimred_xy(),
@@ -148,4 +185,9 @@ layout_dendro <- function(trajectory, diag_offset = 0.05) {
     connection_positions = connection_positions %>% rename_dimred_xy()
   )
 
+  dynplot(
+    dataset = dataset,
+    trajectory = trajectory,
+    layout = layout
+  )
 }

R/dynplot_dimred.R

@@ -0,0 +1,123 @@
+#' Plot a trajectory based on a dimensionality reduction
+#'
+#' @inheritParams dynplot
+#' @param dimred A dimensionality reduction matrix of the cells. Default is `dataset$dimred`.
+#'
+#' @return A ggplot2 object, with the processed data in `plot$data` and `attr(plot$data, "data")`.
+#'
+#' @seealso [dynplot_dendro()], [dynplot_dimred()], [dynplot_graph()] or [dynplot_onedim()]
+#'
+#' @export
+#' @examples
+#' library(ggplot2)
+#' data(example_bifurcating)
+#' dataset <- example_bifurcating
+#'
+#' dynplot_dimred(dataset) +
+#'   geom_trajectory_segments() +
+#'   geom_cell_point(size = 2, colour = "black") +
+#'   geom_cell_point(aes(colour = milestone_percentages), size = 1.8) +
+#'   geom_milestone_label(aes(fill = milestone_id)) +
+#'   scale_milestones_color() +
+#'   scale_milestones_fill()
+#'
+#' dynplot_dimred(dataset) +
+#'   geom_cell_point(aes(colour = select_feature_expression("G1", d = .data))) +
+#'   scale_expression_color() +
+#'   new_scale_color() +
+#'   geom_trajectory_segments(aes(colour = milestone_percentages), size = 2) +
+#'   geom_milestone_label(aes(fill = milestone_id)) +
+#'   scale_milestones_color() +
+#'   scale_milestones_fill()
+dynplot_dimred <- function(dataset, trajectory = dataset, dimred = dataset$dimred) {
+  layout <- list()
+
+  if (is.null(dimred)) {
+    message("No dimred specified, calculating it")
+    dimred <- dyndimred::dimred_landmark_mds(dynwrap::get_expression(dataset), ndim = 2, distance_method = "spearman")
+  }
+
+  # if the dimred originated from the trajectory AND this trajectory already contains a trajectory dimred, we don't calculate the projection. In all other cases we do
+  trajectory_dimred <- trajectory$dimred # store this in a variable for pryr::address
+
+  recalculate_traj_dimred <- !(
+    dynwrap::is_wrapper_with_dimred(trajectory)
+    &&
+      !is.null(trajectory_dimred)
+    &&
+      identical(trajectory_dimred[1, ], dimred[1, ], )
+    &&
+      all(c("dimred_segment_points", "dimred_segment_progressions") %in% names(trajectory))
+  )
+
+  if (!is.null(dataset$dimred_future)) {
+    # check whether the dimred_future has to be recalculated, based on the dimred_digest attribute
+    if(!is.null(attr(dataset$dimred_future, "dimred_digest")) && attr(dataset$dimred_future, "dimred_digest") == digest::digest(dimred, "md5")) {
+      dimred_future <- dataset$dimred_future
+    } else {
+      message("Embedding velocity in dimensionality reduction")
+      dimred_future <- scvelo::embed_velocity(dataset, dimred)
+    }
+    dimred <- cbind(
+      dimred,
+      dimred_future %>% {set_colnames(., paste0(colnames(.), "_future"))}
+    )
+  }
+  cell_positions <- dimred %>%
+    rename_dimred_xy() %>%
+    as.data.frame() %>%
+    rownames_to_column("cell_id")
+
+  assert_that(cell_positions$cell_id %all_in% dataset$cell_ids)
+  layout$cell_positions <- cell_positions
+
+  # trajectory --------------------------------------------------------------
+  if (dynwrap::is_wrapper_with_trajectory(trajectory)) {
+    # trajectory dimred
+    if (!recalculate_traj_dimred) {
+      traj_dimred <- trajectory
+    } else {
+      message("Projecting trajectory onto dimensionality reduction")
+      traj_dimred <- trajectory %>% dynwrap::project_trajectory(dimred)
+    }
+
+    # milestone positions
+    milestone_positions <- as.data.frame(traj_dimred$dimred_milestones[trajectory$milestone_ids, , drop = FALSE]) %>%
+      rename_dimred_xy() %>%
+      as.data.frame() %>%
+      rownames_to_column("milestone_id")
+
+    edge_positions <- trajectory$edge_positions
+    segment_positions <- trajectory$edge_positions
+    segment_progressions <- trajectory$edge_positions
+
+    # trajectory edge positions
+    edge_positions <- trajectory$milestone_network %>%
+      select(from, to) %>%
+      left_join(milestone_positions %>% rename_all(~paste0(., "_from")), c("from" = "milestone_id_from")) %>%
+      left_join(milestone_positions %>% rename_all(~paste0(., "_to")), c("to" = "milestone_id_to"))
+
+    # trajectory segment positions
+    segment_positions <- traj_dimred$dimred_segment_points %>%
+      rename_dimred_xy() %>%
+      as.data.frame() %>%
+      rownames_to_column("point_id")
+
+    segment_progressions <- traj_dimred$dimred_segment_progressions %>%
+      mutate(point_id = segment_positions$point_id)
+
+    # add to layout
+    layout <- c(layout, lst(
+      milestone_positions,
+      edge_positions,
+      segment_positions,
+      segment_progressions
+    ))
+  }
+
+  dynplot(
+    dataset = dataset,
+    trajectory = trajectory,
+    layout = layout
+  )
+}

R/dynplot_graph.R

@@ -0,0 +1,47 @@
+#' Plot a trajectory as a graph
+#'
+#' @inheritParams dynplot
+#'
+#' @return A ggplot2 object, with the processed data in `plot$data` and `attr(plot$data, "data")`.
+#'
+#' @seealso [dynplot_dendro()], [dynplot_dimred()], [dynplot_graph()] or [dynplot_onedim()]
+#'
+#' @export
+#' @examples
+#' library(ggplot2)
+#' data(example_bifurcating)
+#' dataset <- example_bifurcating
+#'
+#' dynplot_graph(dataset) +
+#'   geom_trajectory_divergence() +
+#'   geom_trajectory_segments(size = 2, color = "#333333", arrow_size = 0.5) +
+#'   geom_cell_point(colour = "black", size = 3) +
+#'   geom_cell_point(aes(colour = milestone_percentages), size = 2.8) +
+#'   geom_milestone_label(aes(fill = milestone_id)) +
+#'   scale_milestones_fill() +
+#'   scale_milestones_color()
+dynplot_graph <- function(dataset, trajectory = dataset) {
+  assert_that(dynwrap::is_wrapper_with_trajectory(dataset))
+  trajectory_dimred <- dynwrap::calculate_trajectory_dimred(dataset)
+
+  segments <- calculate_segments_from_edges(trajectory_dimred$edge_positions)
+
+  segment_progressions <- segments$segment_progressions %>% rename_dimred_xy()
+  segment_positions <- segments$segment_positions %>% rename_dimred_xy()
+
+  layout <- lst(
+    cell_positions = trajectory_dimred$cell_positions %>% rename_dimred_xy(),
+    milestone_positions = trajectory_dimred$milestone_positions %>% rename_dimred_xy(),
+    edge_positions = trajectory_dimred$edge_positions %>% rename_dimred_xy(),
+    segment_progressions,
+    segment_positions,
+    divergence_edge_positions = trajectory_dimred$divergence_edge_positions %>% rename_dimred_xy(),
+    divergence_polygon_positions = trajectory_dimred$divergence_polygon_positions %>% rename_dimred_xy()
+  )
+
+  dynplot(
+    dataset = dataset,
+    trajectory = trajectory,
+    layout = layout
+  )
+}

R/dynplot_onedim.R

@@ -0,0 +1,68 @@
+#' Plot a trajectory as a graph
+#'
+#' @inheritParams dynplot
+#'
+#' @return A ggplot2 object, with the processed data in `plot$data` and `attr(plot$data, "data")`.
+#'
+#' @seealso [dynplot_dendro()], [dynplot_dimred()], [dynplot_graph()] or [dynplot_onedim()]
+#'
+#' @export
+#' @examples
+#' library(ggplot2)
+#' data(example_bifurcating)
+#' dataset <- example_bifurcating
+#'
+#' dynplot_onedim(dataset, equal_cell_width = TRUE, margin = 0.05) +
+#'   geom_trajectory_segments() +
+#'   geom_trajectory_connection() +
+#'   geom_cell_point(aes(colour = select_feature_expression("G2", d = .data))) +
+#'   scale_expression_colour() +
+#'   new_scale_colour() +
+#'   geom_milestone_label(aes(y = -0.1, hjust = as.integer(type == "end")))
+#'
+#' dynplot_onedim(dataset, equal_cell_width = TRUE, margin = 0.05) +
+#'   geom_trajectory_segments(aes(color = milestone_percentages), size = 1, color = "#333333") +
+#'   geom_trajectory_connection() +
+#'   scale_milestones_colour() +
+#'   new_scale_colour() +
+#'   geom_cell_point(aes(colour = select_feature_expression("G2", d = .data))) +
+#'   scale_expression_colour() +
+#'   new_scale_colour() +
+#'   geom_milestone_label(aes(y = -0.1, fill = milestone_id, hjust = as.integer(type == "end"))) +
+#'   scale_milestones_fill()
+dynplot_onedim <- function(dataset, trajectory = dataset, margin = 0.02, equal_cell_width = TRUE) {
+  assert_that(dynwrap::is_wrapper_with_trajectory(dataset))
+
+  # reorder
+  dataset$milestone_network <- optimize_order(dataset$milestone_network)
+
+  # linearise
+  linearised <- linearise_trajectory(
+    dataset,
+    margin = margin,
+    equal_cell_width = equal_cell_width
+  )
+
+  # calculate positions of connections
+  connection_positions <- calculate_connections(linearised) %>% rename_dimred_xy()
+
+  segments <- calculate_segments_from_edges(linearised$edge_positions)
+
+  segment_progressions <- segments$segment_progressions %>% rename_dimred_xy()
+  segment_positions <- segments$segment_positions %>% rename_dimred_xy()
+
+  layout <- lst(
+    cell_positions = linearised$cell_positions %>% rename_dimred_xy(),
+    milestone_positions = linearised$milestone_positions %>% rename_dimred_xy(),
+    edge_positions = linearised$edge_positions %>% rename_dimred_xy(),
+    segment_progressions,
+    segment_positions,
+    connection_positions
+  )
+
+  dynplot(
+    dataset = dataset,
+    trajectory = trajectory,
+    layout = layout
+  )
+}

R/facet_grid_data.R

@@ -0,0 +1,93 @@
+#' Grid facetting based on data
+#'
+#' @param ... Data to facet along.
+#' @importFrom stringr str_glue
+#'
+#' @examples
+#' library(ggplot2)
+#' cutoff <- function(x, cutoff) {
+#'   x > cutoff
+#' }
+#'
+#' d <- data.frame(
+#'   x = 1:100,
+#'   y = sqrt(1:100)
+#' )
+#'
+#' ggplot(d) +
+#'   geom_text(aes(x, y, color = x > cutoff, label = hello)) +
+#'   facet_grid_data(hello = c("a", "b", "c"), cutoff = c(30, -100, 10))
+#'
+#' @export
+facet_grid_data <- function(...) {
+  dots <- list(...)
+
+  # make sure at least row or columns are provided
+  assert_that(length(dots) %in% c(1, 2))
+
+  # generate the layout
+  if (names(dots)[[1]] != "") {
+    # rows was provided
+    row_name <- names(dots)[[1]]
+    row_data <- tibble(
+      !!row_name := dots[[1]],
+      ROW = seq_along(dots[[1]])
+    )
+    rows <-  vars(!!row_name := !!as.name(row_name))
+  } else {
+    row_name <- rows <- NULL
+    row_data <- tibble(ROW = 1)
+  }
+
+  if (length(dots) == 2 && names(dots)[[2]] != "")  {
+    # cols was provided
+    col_name <- names(dots)[[2]]
+    col_data <- tibble(
+      !!col_name := dots[[2]],
+      COL = seq_along(dots[[2]])
+    )
+    cols <-  vars(!!col_name := !!as.name(col_name))
+  } else {
+    col_name <- cols <-NULL
+    col_data <- tibble(COL = 1)
+  }
+
+  layout <- crossing(col_data, row_data) %>%
+    mutate(PANEL = row_number(), SCALE_X = 1, SCALE_Y = 1) %>%
+    as.data.frame() # convert to dataframe to avoid warnings when running unique(layout[names(params$cols)])
+
+  # generate the labeller
+  labeller <- function(x) {
+    if (col_name %in% colnames(x)) {
+      x %>%
+        mutate_at(vars(!!col_name), ~str_glue("{col_name} = {.}"))
+    } else if (row_name %in% colnames(x)) {
+      x %>%
+        mutate_at(vars(!!row_name), ~str_glue("{row_name} = {.}"))
+    }
+  }
+
+  ggproto(
+    "FacetGridData",
+    facet_grid(rows = rows, cols = cols, labeller = labeller),
+    compute_layout = function(self, data, params) {
+      layout
+    },
+    map_data = function(self, data, layout, params) {
+      results <- pmap(layout, function(PANEL, ROW, COL, SCALE_X, SCALE_Y, ...) {
+        variable_values <- list(...)
+
+        for (variable in names(variable_values)) {
+          data[[variable]] <- variable_values[[variable]]
+        }
+
+        data$PANEL <- PANEL
+
+        data
+      })
+      result <- bind_rows(results)
+
+      result
+    }
+  )
+}

R/facet_wrap_data.R

@@ -0,0 +1,80 @@
+#' Wrap facetting based on data
+#'
+#' @importFrom stringr str_glue
+#'
+#' @param ... Data to facet along.
+#'
+#' @examples
+#' library(ggplot2)
+#' cutoff <- function(x, cutoff) {
+#'   x > cutoff
+#' }
+#'
+#' d <- data.frame(
+#'   x = 1:100,
+#'   y = sqrt(1:100)
+#' )
+#'
+#' ggplot(d) +
+#'   geom_text(aes(x, y, color = x > cutoff, label = cutoff)) +
+#'   facet_wrap_data(cutoff = c(1, 5, 15))
+#'
+#' @export
+facet_wrap_data <- function(..., nrow = NULL, ncol = 4, labeller = NULL) {
+  dots <- list(...)
+
+  # make sure that one dots was provided
+  assert_that(length(dots) == 1)
+
+  # prepare the facet data
+  facet_name <- names(dots)[[1]]
+  facet_data <- tibble(
+    !!facet_name := dots[[1]]
+  )
+  facets <-  vars(!!as.name(facet_name))
+
+  # calculate the layout (number of rows and columns)
+  if (is.null(nrow) && is.null(ncol)) {
+    ncol <- 4
+  }
+  if (is.null(nrow)) {
+    nrow <- ceiling(nrow(facet_data)/ncol)
+  }
+  if (is.null(ncol)) {
+    ncol <- ceiling(nrow(facet_data)/ncol)
+  }
+
+  # generate the layout
+  layout <- facet_data %>%
+    mutate(PANEL = row_number(), SCALE_X = 1, SCALE_Y = 1, COL = ((PANEL - 1) %% ncol) + 1, ROW = ceiling(PANEL / ncol))
+
+  # generate the labeller
+  labeller <- function(x) {
+    x %>%
+      mutate_at(vars(!!facet_name), ~str_glue("{facet_name} = {.}"))
+  }
+
+  ggproto(
+    "FacetWrapData",
+    facet_wrap(facets = facets, labeller = labeller),
+    compute_layout = function(self, data, params) {
+      layout
+    },
+    map_data = function(self, data, layout, params) {
+      results <- pmap(layout, function(PANEL, ROW, COL, SCALE_X, SCALE_Y, ...) {
+        variable_values <- list(...)
+
+        for (variable in names(variable_values)) {
+          data[[variable]] <- variable_values[[variable]]
+        }
+
+        data$PANEL <- PANEL
+
+        data
+      })
+      result <- bind_rows(results)
+
+      result
+    }
+  )
+}

R/geom_cell_.R

@@ -4,6 +4,12 @@ GeomCellPoint <- ggproto(
   default_aes = aesIntersect(GeomPoint$default_aes, aes(color = "grey80"))
 )
 
+GeomHexPoint <- ggproto(
+  "GeomHexPoint",
+  GeomHex,
+  default_aes = aesIntersect(GeomHex$default_aes, aes(color = "grey80"))
+)
+
 #' Plotting cells
 #'
 #' @param mapping Set of aesthetic mappings created by aes().
@@ -37,19 +43,102 @@ geom_cell_point <- function(
   )
 }
 
+
+StatCellHex <- ggproto(
+  "StatCellHex",
+  Stat,
+  compute_group = function(
+    data,
+    scales,
+    binwidth = NULL,
+    bins = 30,
+    na.rm = FALSE
+  ) {
+    ggplot2:::try_require("hexbin", "stat_binhex")
+
+    # parts of this code were taken from ggplot2 stat_binhex
+
+    binwidth <- binwidth %||% ggplot2:::hex_binwidth(bins, scales)
+
+    # calculate bounds and bins
+    if (length(binwidth) == 1) {
+      binwidth <- rep(binwidth, 2)
+    }
+    xbnds <- ggplot2:::hex_bounds(data$x, binwidth[1])
+    xbins <- diff(xbnds)/binwidth[1]
+    ybnds <- ggplot2:::hex_bounds(data$y, binwidth[2])
+    ybins <- diff(ybnds)/binwidth[2]
+    hb <- hexbin::hexbin(data$x, xbnds = xbnds, xbins = xbins, data$y,
+                         ybnds = ybnds, shape = ybins/xbins, IDs = TRUE)
+
+    out_coords <- bind_cols(hexbin::hcell2xy(hb)) %>%
+      mutate(hexagon_id = hb@cell)
+
+    # get count and density
+    out_count <- data %>%
+      mutate(hexagon_id = hb@cID) %>%
+      group_by(hexagon_id) %>%
+      summarise(
+        count = n()
+      ) %>%
+      mutate(
+        density = count / max(count)
+      )
+
+    # summarise all numeric variables to first value
+    out_data <- data %>%
+      select(-x, -y, -group, -PANEL) %>%
+      mutate(hexagon_id = hb@cID) %>%
+      group_by(hexagon_id) %>%
+      slice(1) %>%
+      ungroup()
+
+    out <- left_join(out_coords, out_data, "hexagon_id") %>%
+      left_join(out_count, "hexagon_id")
+
+    out
+  }
+)
+
+
 #' @rdname geom_cell
 #' @export
-geom_cell_hex <- function() {
+geom_cell_hex <- function(
+  mapping = NULL,
+  data = construct_get_cell_info(),
+  bins = 100,
+  ...,
+  show.legend = NA
+) {
+  mapping <- aesIntersect(mapping, aes_(x=~x, y=~y))
 
+  layer(
+    data = data,
+    mapping = mapping,
+    stat = StatCellHex,
+    geom = GeomHexPoint,
+    position = "identity",
+    show.legend = show.legend,
+    inherit.aes = FALSE,
+    params = list(
+      na.rm = FALSE,
+      bins = bins,
+      ...
+    )
+  )
 }
 
 construct_get_cell_info <- function() {
   function(data) {
     # first parse the mapping to know what to put inside the cell info
-    walk(mapping, function(mapping_element) {
-      assign("data", data, envir = environment(mapping$colour))
-    })
+    out <- attr(data, "data")$cell_info
+
+    # create expression, counts and velocity_vector
+    data_env <- new.env(parent = emptyenv())
+    assign("data", data, data_env)
 
-    attr(data, "data")$cell_info
+    d <- list(data_env)
+    out$d <- d
+    out
   }
 }

R/geom_cell_contour.R

@@ -0,0 +1,97 @@
+GeomCellContour <- ggproto(
+  "GeomCellContour",
+  GeomPolygon,
+  default_aes = aesIntersect(aes(alpha = 0.2), GeomPolygon$default_aes)
+)
+
+StatCellContour <- ggproto(
+  "StatCellContour",
+  StatDensity2d,
+  setup_params = function(data, params) {
+    xlims <- c(min(data$x), max(data$x))
+    ylims <- c(min(data$y), max(data$y))
+
+    # determine bandwidth
+    assert_that(!is.null(params$relative_bandwidth))
+    params$bandwidth_y <- diff(ylims) * params$relative_bandwidth
+    params$bandwidth_x <- diff(xlims) * params$relative_bandwidth
+
+    # determine padding
+    assert_that(!is.null(params$padding))
+
+    xpad <- diff(xlims) * params$padding
+    ypad <- diff(ylims) * params$padding
+
+    xlims <- xlims + c(-xpad, xpad)
+    ylims <- ylims + c(-ypad, ypad)
+
+    params$xlims <- xlims
+    params$ylims <- ylims
+
+    params
+  },
+  compute_group = function(data, scales, bandwidth_x, bandwidth_y, xlims, ylims, na.rm = F, resolution = 200, relative_density_cutoff = 0.2, relative_bandwidth = NULL, padding = NULL) {
+    density <- MASS::kde2d(data$x, data$y, h = c(bandwidth_x, bandwidth_y), lims = c(xlims, ylims), n = resolution)
+    df <- expand.grid(x = density$x, y = density$y)
+    df$group <- data$group[1]
+    df$z <- as.vector(density$z)
+
+    density_cutoff <- min(df$z) + (max(df$z) - min(df$z)) * relative_density_cutoff
+
+    # output2 <- ggplot2:::contour_lines(df, breaks = density_cutoff, complete = FALSE)
+
+    lines <- isoband::isolines(density$x, density$y, t(density$z), levels = density_cutoff)[[1]]
+
+    # plot <- qplot(data$x, data$y)
+    # plot + geom_point(aes(lines$x, lines$y))
+
+    # ggplot(df) +
+    #   geom_tile(aes(x, y, color = z)) +
+    #   geom_point(aes(x, y), data = data)
+
+    output2 <- tibble(
+      y = lines$y,
+      x = lines$x
+    )
+    # browser()
+
+    output2
+  }
+)
+
+#' Plot contour around cells based
+#'
+#' @param relative_density_cutoff At whtat level of density the contour should be drawn, should be between 0 and 1
+#' @param relative_bandwidth Bandwidth calculated relative to the x and y limits of the points, should be between 0 and 1
+#' @param resolution The higher, the more accurate the polygon will be drawn at the cost of longer computing/drawing time
+#' @param padding How much padding to add to the limits, to avoid the contour to be drawn outside the plot
+geom_cell_contour <- function(
+  mapping = NULL,
+  data = construct_get_cell_info(),
+  relative_density_cutoff = 0.2,
+  relative_bandwidth = 0.2,
+  padding = 1,
+  resolution = 200,
+  ...,
+  show.legend = NA
+) {
+  mapping <- aesIntersect(mapping, aes_(x=~x, y=~y, group=~group_id))
+
+  layer(
+    data = data,
+    mapping = mapping,
+    stat = StatCellContour,
+    geom = GeomCellContour,
+    position = "identity",
+    show.legend = show.legend,
+    inherit.aes = FALSE,
+    params = lst(
+      na.rm = FALSE,
+      resolution = resolution,
+      padding = padding,
+      relative_bandwidth = relative_bandwidth,
+      relative_density_cutoff = relative_density_cutoff,
+      ...
+    )
+  )
+}

R/geom_cell_contour_label.R

@@ -0,0 +1,64 @@
+GeomCellContourLabel <- ggproto(
+  "GeomCellContourLabel",
+  ggrepel::GeomLabelRepel
+)
+
+StatCellContourLabel <- ggproto(
+  "StatCellContourLabel",
+  StatDensity2d,
+  setup_params = function(data, params) {
+    xlims <- c(min(data$x), max(data$x))
+    ylims <- c(min(data$y), max(data$y))
+
+    # determine bandwidth
+    assert_that(!is.null(params$relative_bandwidth))
+    params$bandwidth_y <- diff(ylims) * params$relative_bandwidth
+    params$bandwidth_x <- diff(xlims) * params$relative_bandwidth
+
+    params$xlims <- xlims
+    params$ylims <- ylims
+
+    params
+  },
+  compute_group = function(data, scales, bandwidth_x, bandwidth_y, xlims, ylims, na.rm = F, relative_bandwidth = NULL) {
+    density <- MASS::kde2d(data$x, data$y, h = c(bandwidth_x, bandwidth_y), lims = c(xlims, ylims), n = 100)
+    df <- expand.grid(x = density$x, y = density$y)
+    df$group <- data$group[1]
+    df$z <- as.vector(density$z)
+
+    row <- df[which.max(df$z),]
+    row <- cbind(row, data[1, !colnames(data) %in% c("x", "y")])
+    row
+  }
+)
+
+#' Plot contour around cells based on their density
+#'
+#' @inheritParams ggrepel::geom_label_repel
+#' @inheritParams geom_cell_contour
+geom_cell_contour_label <- function(
+  mapping = NULL,
+  data = construct_get_cell_info(),
+  relative_bandwidth = 0.2,
+  min.segment.length = Inf,
+  ...,
+  show.legend = NA
+) {
+  mapping <- aesIntersect(mapping, aes_(x=~x, y=~y))
+
+  layer(
+    data = data,
+    mapping = mapping,
+    stat = StatCellContourLabel,
+    geom = GeomCellContourLabel,
+    position = "identity",
+    show.legend = show.legend,
+    inherit.aes = FALSE,
+    params = lst(
+      na.rm = FALSE,
+      relative_bandwidth = relative_bandwidth,
+      min.segment.length = min.segment.length,
+      ...
+    )
+  )
+}

R/geom_expression.R

@@ -1,7 +1,7 @@
-GeomExpressionRaster <- ggproto(
-  "GeomExpressionRaster",
-  GeomRaster,
-  default_aes = aesIntersect(GeomRaster$default_aes, aes(color = "grey80"))
+GeomExpressionTile <- ggproto(
+  "GeomExpressionTile",
+  GeomTile,
+  default_aes = aesIntersect(GeomTile$default_aes, aes(color = "grey80"))
 )
 
 #' Plotting expression
@@ -13,11 +13,12 @@ GeomExpressionRaster <- ggproto(
 #' @rdname geom_cell
 #'
 #' @export
-geom_expression_raster <- function(
+geom_expression_tile <- function(
   mapping = NULL,
   ...,
   show.legend = NA,
-  data = construct_get_expression_info()
+  rescale = dynutils::scale_quantile,
+  data = construct_get_expression_info(rescale = rescale)
 ) {
   assign("mapping", mapping, envir = environment(data)) # place the mapping in the data environment
 
@@ -26,7 +27,7 @@ geom_expression_raster <- function(
     data = data,
     mapping = mapping,
     stat = StatIdentity,
-    geom = GeomExpressionRaster,
+    geom = GeomExpressionTile,
     position = "identity",
     show.legend = show.legend,
     inherit.aes = FALSE,
@@ -37,18 +38,34 @@ geom_expression_raster <- function(
   )
 }
 
-construct_get_expression_info <- function() {
+construct_get_expression_info <- function(rescale) {
   function(data) {
     feature_info <- attr(data, "data")$feature_info
     cell_info <- attr(data, "data")$cell_info
 
-    expression <- get_expression(attr(data, "data")$dataset)[cell_info$cell_id, feature_info$feature_id]
+    # fetch data
+    expression <- get_expression(attr(data, "data")$dataset)
 
-    expression_info <- reshape2::melt(as.matrix(expression), varnames = c("cell_id", "feature_id"), value.name = "expression") %>%
-      mutate(cell_id = as.character(cell_id), feature_id = as.character(feature_id))
+    # ensure exact subset and ordering
+    expression <- expression[cell_info$cell_id, feature_info$feature_id, drop = FALSE]
 
+    # rescale expression
+    expression <- rescale(expression)
+
+    # reshape into data frame
+    expression_info <- reshape2::melt(
+      as.matrix(expression),
+      varnames = c("cell_id", "feature_id"),
+      value.name = "expression"
+    ) %>%
+      mutate(
+        cell_id = as.character(.data$cell_id),
+        feature_id = as.character(.data$feature_id)
+      )
+
+    # add info data frames
     expression_info %>%
       left_join(feature_info, "feature_id") %>%
-      left_join(cell_info %>% select(-y), "cell_id")
+      left_join(cell_info %>% select(-.data$y), "cell_id")
   }
 }

R/geom_feature_.R

@@ -0,0 +1,45 @@
+GeomFeatureText <- ggproto(
+  "GeomFeatureText",
+  GeomText,
+  default_aes = aesIntersect(aes(hjust = 1), GeomText$default_aes)
+)
+
+#' Plotting feature text
+#'
+#' @param mapping Set of aesthetic mappings created by aes().
+#' @param data A function created by [get_feature_info_constructor()].
+#' @param show.legend Whether to show a legend for this geom
+#'
+#' @rdname geom_feature
+#'
+#' @export
+geom_feature_text <- function(
+  mapping = NULL,
+  data = construct_get_feature_info(),
+  ...,
+  show.legend = TRUE
+) {
+  mapping <- aesIntersect(mapping, aes_(x=~0, y=~y, label=~feature_id))
+  layer(
+    data = data,
+    mapping = mapping,
+    stat = StatIdentity,
+    geom = GeomFeatureText,
+    position = "identity",
+    show.legend = show.legend,
+    inherit.aes = FALSE,
+    params = lst(
+      na.rm = FALSE,
+      ...
+    )
+  )
+}
+
+
+
+construct_get_feature_info <- function() {
+  function(data) {
+    attr(data, "data")$feature_info
+  }
+}
+

R/geom_milestone_.R

@@ -1,13 +1,13 @@
 GeomMilestoneLabel <- ggproto(
   "GeomMilestoneLabel",
   GeomLabel,
-  default_aes = aesIntersect(aes(fill = "#111111CC", fontface = "bold"), GeomLabel$default_aes)
+  default_aes = aesIntersect(aes(color = "white", fill = "#111111CC", fontface = "bold"), GeomLabel$default_aes)
 )
 
 #' Plotting milestones
 #'
 #' @param mapping Set of aesthetic mappings created by aes().
-#' @param data A function created by [get_cell_info_constructor()].
+#' @param data A function created by [construct_get_cell_info()].
 #' @param show.legend Whether to show a legend for this geom
 #'
 #' @rdname geom_milestone
@@ -17,7 +17,7 @@ geom_milestone_label <- function(
   mapping = NULL,
   data = get_milestone_info,
   ...,
-  show.legend = TRUE
+  show.legend = FALSE
 ) {
   mapping <- aesIntersect(mapping, aes_(x=~x, y=~y, label=~label))
   layer(
@@ -35,19 +35,50 @@ geom_milestone_label <- function(
   )
 }
 
+
+GeomMilestonePoint <- ggproto(
+  "GeomMilestoneLabel",
+  GeomPoint,
+  draw_panel = function(self, data, panel_params, coord, shadow, ...) {
+    original_draw_panel <- GeomPoint$draw_panel
+
+    # draw point ---------
+    grob_point <- original_draw_panel(data = data, coord = coord, panel_params = panel_params, ...)
+
+    # draw shadows --------
+    assert_that(is.list(shadow) || isFALSE(shadow), msg = "shadow should be a list created by shadow_defaults() or FALSE")
+    if (is.list(shadow)) {
+      grob_point_shadow <- original_draw_panel(
+        data = data %>% mutate(colour = shadow$color, size = size + shadow$size), panel_params = panel_params, coord = coord, ...)
+    } else {
+      grob_point_shadow <- grid::grob()
+    }
+
+    # combine grobs
+    grid::gList(
+      grob_point_shadow,
+      grob_point
+    )
+  },
+  # this function is just here so that shadow becomes a parameter
+  draw_group = function(data, panel_params, coord, shadow) {
+  }
+)
+
 #' @export
 geom_milestone_point <- function(
   mapping = NULL,
   data = get_milestone_info,
   position = "identity",
   show.legend = NA,
   size = 10,
+  shadow = shadow_defaults(),
   ...
 ) {
-  mapping <- aesIntersect(mapping, aes_(x=~x, y=~y))
-  layer(data = data, mapping = mapping, stat = StatIdentity, geom = GeomPoint,
+  mapping <- aesIntersect(mapping, aes_(x=~x, y=~y, color=~milestone_id))
+  layer(data = data, mapping = mapping, stat = StatIdentity, geom = GeomMilestonePoint,
         position = position, show.legend = show.legend, inherit.aes = FALSE,
-        params = lst(na.rm = FALSE, size = size, ...)
+        params = lst(na.rm = FALSE, size = size, shadow = shadow, ...)
   )
 }
 

R/geom_trajectory_segments.R

@@ -2,7 +2,7 @@ GeomTrajectorySegments <- ggproto(
   "GeomTrajectorySegments",
   GeomPath,
   default_aes = aesIntersect(aes(linejoin = "mitre", lineend = "square"), GeomPath$default_aes),
-  draw_panel = function(self, data, panel_params, coord, arrow = NULL, arrow_size = 1, shadow = dynplot::shadow(),  ...) {
+  draw_panel = function(self, data, panel_params, coord, arrow = NULL, arrow_size = 1, shadow = shadow_defaults(),  ...) {
     original_draw_panel <- GeomPath$draw_panel
 
     # draw path ---------
@@ -11,15 +11,17 @@ GeomTrajectorySegments <- ggproto(
     # draw arrows -------
     # select the rows at draw_arrow and draw_arrow - 1
     data_arrows <- filter(data, draw_arrow | lead(draw_arrow)) %>%
-      mutate(group = ceiling(row_number() / 2)) %>%
-      mutate(size = size * arrow_size)
+      mutate(group = ceiling(row_number() / 2))
+
+    # the size of the arrow is scaled with the arrow size and segment size parameters (the first one of the latter)
+    arrow$length <- arrow$length * arrow_size * data_arrows$size[1]
 
     if (nrow(data_arrows) > 1) {
       grob_arrows <- original_draw_panel(data = data_arrows, panel_params = panel_params, coord = coord, arrow = arrow, lineend = "butt", linejoin = "mitre", ...)
     } else {grob_arrows <- grid::grob()}
 
     # draw shadows ------
-    assert_that(is.list(shadow) || isFALSE(shadow), msg = "shadow should be a list created by shadow() or FALSE")
+    assert_that(is.list(shadow) || isFALSE(shadow), msg = "shadow should be a list created by shadow_defaults() or FALSE")
     if (is.list(shadow)) {
       grob_path_shadow <- original_draw_panel(
         data = data %>% mutate(colour = shadow$color, size = size + shadow$size), panel_params = panel_params, coord = coord, arrow = NULL, lineend = data$lineend[[1]], linejoin = data$linejoin[[1]], ...)
@@ -37,7 +39,8 @@ GeomTrajectorySegments <- ggproto(
       grob_path
     )
   },
-  draw_group = function(data, panel_params, coord, arrow = NULL, arrow_size = 1, shadow = "black") {
+  # this function is just here so that shadow, arrow and arrow_size becomes a parameter
+  draw_group = function(data, panel_params, coord, arrow = NULL, arrow_size = 1, shadow) {
   },
   draw_key = function(data, params, size) {
     data$linetype[is.na(data$linetype)] <- 0
@@ -58,7 +61,7 @@ GeomTrajectorySegments <- ggproto(
 #' @inheritParams ggplot2::geom_segment
 #' @param position_arrow Where to place the arrows within the segments. Typically these are functions created by [position_trajectory_arrows_middle()] or [position_trajectory_arrows_boundaries()].
 #' @param arrow_size The size of the arrow relative to the line size.
-#' @param shadow Shadow specification as created by [shadow()]
+#' @param shadow Shadow specification as created by [shadow_defaults()]
 #' @param data A function created by [construct_get_segment_info()].
 #'
 #' @export
@@ -67,7 +70,7 @@ geom_trajectory_segments <- function(
   position_arrow = position_trajectory_arrows_middle(),
   arrow = ggplot2::arrow(length = ggplot2::unit(0.5, "cm"), type = "closed"),
   arrow_size = 1,
-  shadow = if("colour" %in% names(mapping)) {dynplot::shadow()} else {FALSE},
+  shadow = if("colour" %in% names(mapping)) shadow_defaults() else FALSE,
   ...,
   data = construct_get_segment_info(position_arrow),
   show.legend = NA
@@ -81,7 +84,7 @@ geom_trajectory_segments <- function(
     position = "identity",
     show.legend = show.legend,
     inherit.aes = FALSE,
-    params = lst(arrow, arrow_size, shadow = shadow, ...)
+    params = lst(arrow, arrow_size, shadow, ...)
   )
 }
 
@@ -93,6 +96,10 @@ get_edge_info <- function(data) {
 }
 construct_get_segment_info <- function(position_arrow = position_trajectory_arrows_middle()) {
   function(data) {
+    if (is.null(attr(data, "data")$segment_info)) {
+      stop("No segment information in the data, did you add a trajectory?")
+    }
+
     segment_info <- attr(data, "data")$segment_info
     segment_info <- segment_info %>% position_arrow()
 
@@ -107,14 +114,15 @@ construct_get_segment_info <- function(position_arrow = position_trajectory_arro
 # used to position the arrows nicely
 calculate_trajectory_segment_length <- function(data) {
   data %>%
-    group_by(from, to) %>%
+    group_by(.data$from, .data$to) %>%
     mutate(
       length = c(
         0,
-        sqrt((head(y, -1) - tail(y, -1))**2 + (head(x, -1) - tail(x, -1))**2)
+        sqrt((head(.data$y, -1) - tail(.data$y, -1))**2 + (head(.data$x, -1) - tail(.data$x, -1))**2)
       ),
-      cumlength = cumsum(length)
-    )
+      cumlength = cumsum(.data$length)
+    ) %>%
+    ungroup()
 }
 
 
@@ -126,7 +134,7 @@ position_trajectory_arrows_middle <- function() {
       calculate_trajectory_segment_length() %>%
       group_by(from, to) %>%
       mutate(
-        draw_arrow = (row_number() == which.min(abs(cumlength - last(cumlength)/2))) & row_number() > 1
+        draw_arrow = (row_number() == which.min(abs(cumlength - last(cumlength)/2))) & (row_number() > 1) & (directed)
       ) %>%
       ungroup()
   }
@@ -143,7 +151,7 @@ position_trajectory_arrows_boundaries <- function(quantile = 0.3) {
       mutate(
         draw_arrow_start = (row_number() == which.min(abs(cumlength - last(cumlength) * quantile))) & row_number() > 1,
         draw_arrow_end = (row_number() == which.min(abs(cumlength - last(cumlength) * (1-quantile)))) & row_number() > 1,
-        draw_arrow = draw_arrow_start | draw_arrow_end
+        draw_arrow = (draw_arrow_start | draw_arrow_end) & (row_number() > 1) & (directed)
       ) %>%
       ungroup()
   }
@@ -152,9 +160,6 @@ position_trajectory_arrows_boundaries <- function(quantile = 0.3) {
 
 
 #' @export
-shadow <- function(size = 1, color = "#222222") {
-  lst(
-    size,
-    color
-  )
+shadow_defaults <- function(size = 1, color = "#222222") {
+  lst(size, color)
 }

R/geom_velocity_.R → R/geom_velocity_arrow.R

@@ -39,7 +39,7 @@ geom_velocity_arrow <- function(
   data = construct_get_velocity_info(stat),
   show.legend = NA
 ) {
-  mapping <- aesIntersect(mapping, aes_(x=~x, y=~y, xend=~x_projected, yend=~y_projected, length=~length))
+  mapping <- aesIntersect(mapping, aes_(x=~x, y=~y, xend=~x_future, yend=~y_future, length=~length))
   layer(
     data = data,
     mapping = mapping,
@@ -61,7 +61,7 @@ construct_get_velocity_info <- function(position) {
   get_velocity_info <- function(data) {
     cell_positions <- attr(data, "data")$cell_info
     assert_that(
-      c("x", "y", "x_projected", "y_projected") %all_in% colnames(cell_positions),
+      c("x", "y", "x_future", "y_future") %all_in% colnames(cell_positions),
       msg = "This layout does not contain information on velocity"
     )
 
@@ -71,28 +71,36 @@ construct_get_velocity_info <- function(position) {
 }
 
 
-#' @param cell_positions Dataframe contains at least x, y, x_projected and y_projected
+#' @rdname embed_arrows
 embed_arrows_cells <- function(cell_positions) {
   cell_positions %>%
     mutate(
-      length = sqrt((y_projected - y)**2 + (x_projected - x)**2),
+      length = sqrt((y_future - y)**2 + (x_future - x)**2),
     )
 }
 
+#' @rdname embed_arrows
+embed_arrows_random <- function(cell_positions, n_cells = 100) {
+  embed_arrows_cells(cell_positions) %>%
+    sample_n(min(nrow(cell_positions), n_cells))
+}
 
 
-#' @param cell_positions The dimensionality reduction which contains at least x, y, x_projected and y_projected
+#' Embed velocity arrows in a grid
+#' @param cell_positions The dimensionality reduction which contains at least x, y, x_future and y_future
 #' @param grid_n Number of rows and columns in the grid
-#' @param grid_sd Standard deviation for smoothing arrows
+#' @param grid_bandwidth Bandwidth relative to the spacing of the grid
+#'
+#' @rdname embed_arrows
 embed_arrows_grid <- function(
   cell_positions,
   grid_n = c(15, 15),
-  grid_sd = NULL,
+  grid_bandwidth = 1/3,
   max_arrow_length = NULL,
   filter = rlang::quo(mass > max(mass) * 0.1)
 ) {
   assert_that(is.data.frame(cell_positions))
-  assert_that(c("x", "y", "x_projected", "y_projected") %all_in% colnames(cell_positions))
+  assert_that(c("x", "y", "x_future", "y_future") %all_in% colnames(cell_positions))
 
   if (length(grid_n) == 1) {
     grid_n <- c(grid_n, grid_n)
@@ -104,31 +112,31 @@ embed_arrows_grid <- function(
   grid_h <- grid_n[2]
 
   # calculate grid points
-  range_x <- range(unlist(cell_positions[, c("x", "x_projected")]))
-  range_y <- range(unlist(cell_positions[, c("y", "y_projected")]))
+  range_x <- range(unlist(cell_positions[, c("x", "x_future")]), na.rm = TRUE)
+  range_y <- range(unlist(cell_positions[, c("y", "y_future")]), na.rm = TRUE)
   grid_x <- seq(range_x[1],range_x[2],length.out=grid_w)
   grid_y <- seq(range_y[1],range_y[2],length.out=grid_h)
 
   diff_x <- grid_x[2] - grid_x[1]
   diff_y <- grid_y[2] - grid_y[1]
 
-  if(is.null(grid_sd)) {
-    grid_sd <- sqrt((diff_x)^2 + (diff_y)^2)/3
-  }
+  grid_sd <- sqrt((diff_x)^2 + (diff_y)^2) * grid_bandwidth
   if(is.null(max_arrow_length)) {
     max_arrow_length <- min(c(diff_x, diff_y))
   }
 
   cell_positions_difference <- tibble(
-    x = cell_positions$x_projected - cell_positions$x,
-    y = cell_positions$y_projected - cell_positions$y
+    x = cell_positions$x_future - cell_positions$x,
+    y = cell_positions$y_future - cell_positions$y
   )
 
   # calculate for each gaussian the smoothed arrow using a gaussian kernel
   garrows <- map_dfr(grid_x, function(x) {
     # cell distances and weights to each grid point
     cd <- sqrt(outer(cell_positions$y,-grid_y,'+')^2 + (x-cell_positions$x)^2)
-    cw <- dnorm(cd,sd=grid_sd)
+    # cw <- dexp(cd, 100)
+    # cw <- dnorm(cd,sd=grid_sd)
+    cw <- cd < grid_sd
 
     # calculate the actual arrow
     gw <- Matrix::colSums(cw)
@@ -157,13 +165,14 @@ embed_arrows_grid <- function(
       length = length * norm,
       x_difference = x_difference * norm,
       y_difference = y_difference * norm,
-      x_projected = x + x_difference,
-      y_projected = y + y_difference,
+      x_future = x + x_difference,
+      y_future = y + y_difference,
     )
 
   garrows
 }
 
+
 #' @export
 stat_velocity_cells <- dynutils::inherit_default_params(
   list(embed_arrows_cells),
@@ -186,7 +195,7 @@ stat_velocity_grid <- dynutils::inherit_default_params(
         embed_arrows_grid(
           attr(data, "data")$cell_info,
           grid_n = grid_n,
-          grid_sd = grid_sd,
+          grid_bandwidth = grid_bandwidth,
           max_arrow_length = max_arrow_length,
           filter = filter
         )
@@ -198,17 +207,13 @@ formals(stat_velocity_grid) <- formals(embed_arrows_grid)[2:length(formals(embed
 
 #' @export
 stat_velocity_random <- dynutils::inherit_default_params(
-  list(embed_arrows_cells),
-  function(sample_n = 100, ...) {
+  list(embed_arrows_random),
+  function(...) {
     list(
       data = function(data) {
-        embedding <- embed_arrows_cells(
-          attr(data, "data")$cell_info
-        )
-        embedding %>%
-          sample_n(min(nrow(embedding), !!sample_n))
+        embed_arrows_random(attr(data, "data")$cell_info, n_cells = n_cells)
       }
     )
   }
 )
-formals(stat_velocity_random) <- formals(embed_arrows_cells)[2:length(formals(embed_arrows_cells))]
+formals(stat_velocity_random) <- formals(embed_arrows_random)[2:length(formals(embed_arrows_random))]

R/geom_velocity_stream.R

@@ -0,0 +1,157 @@
+GeomVelocityStream <- ggproto(
+  "GeomVelocityStream",
+  GeomPath,
+  default_aes = aesIntersect(GeomPath$default_aes, aes(color = "black", linejoin = "mitre", lineend = "butt", length = 1)),
+  draw_panel = function(self, data, panel_params, coord, arrow = NULL, arrow_size = 1, shadow = shadow_defaults(),  ...) {
+    original_draw_panel <- GeomPath$draw_panel
+
+    # draw path ---------
+    grob_path <- original_draw_panel(data = data, panel_params = panel_params, coord = coord, arrow = NULL, lineend = data$lineend[[1]], linejoin = data$linejoin[[1]], ...)
+
+    # draw arrows -------
+    # select the two rows closest to percentage 0.5
+    data_arrows <- data %>%
+      group_by(group) %>%
+      mutate(
+        draw_arrow = (row_number() == which.min(abs(percentage - 0.5)))
+      ) %>%
+      filter(
+        draw_arrow | lag(draw_arrow)
+      ) %>%
+      ungroup()
+
+    # the size of the arrow is scaled with the arrow size and segment size parameters (the first one of the latter)
+    arrow$length <- arrow$length * arrow_size * data_arrows$size[1]
+
+    if (nrow(data_arrows) > 1) {
+      grob_arrows <- original_draw_panel(data = data_arrows, panel_params = panel_params, coord = coord, arrow = arrow, lineend = "butt", linejoin = "mitre", ...)
+    } else {grob_arrows <- grid::grob()}
+
+    # combine grobs
+    grid::gList(
+      grob_arrows,
+      grob_path
+    )
+  },
+
+  draw_group = function(data, panel_params, coord, arrow = NULL, arrow_size = 1, shadow = "black") {
+  },
+  required_aes = union(GeomPath$required_aes, c("percentage"))
+)
+
+#' Plotting velocity
+#'
+#' @inheritParams ggplot2::geom_segment
+#' @param stat Where to place the arrows, such as for every cell ([stat_velocity_cells()]) or using a grid ([stat_velocity_grid()])
+#' @param data A function created by [construct_get_velocity_info()]
+#'
+#'
+#' @rdname geom_velocity
+#'
+#' @export
+geom_velocity_stream <- function(
+  mapping = NULL,
+  stat = stat_velocity_stream(),
+  arrow = ggplot2::arrow(length = ggplot2::unit(0.5, "cm"), type = "closed"),
+  arrow_size = 1,
+  ...,
+  data = construct_get_velocity_info(stat),
+  show.legend = NA
+) {
+  mapping <- aesIntersect(mapping, aes_(x=~x, y=~y, group=~line, percentage=~percentage))
+
+  layer(
+    data = data,
+    mapping = mapping,
+    stat = StatIdentity,
+    geom = GeomVelocityStream,
+    position = "identity",
+    show.legend = show.legend,
+    inherit.aes = FALSE,
+    params = lst(
+      arrow,
+      arrow_size,
+      na.rm = FALSE,
+      ...
+    )
+  )
+}
+
+#' @rdname embed_arrows
+embed_arrows_stream <- function(
+  cell_positions,
+  grid_bandwidth = 1/3,
+  filter = rlang::quo(mass > max(mass) * 0.1)
+) {
+  matplotlib <- reticulate::import("matplotlib")
+
+  ax <- matplotlib$pyplot$axes()
+
+  grid_arrows <- embed_arrows_grid(cell_positions, filter = TRUE, grid_n = 20)
+  grid_arrows <- grid_arrows %>%
+    mutate(
+      filter = rlang::eval_tidy(filter, data = .)
+    ) %>%
+    mutate(
+      x_difference = ifelse(filter, x_difference, 0),
+      y_difference = ifelse(filter, y_difference, 0)
+    )
+
+  streamplot <- matplotlib$axes$mstream$streamplot(
+    ax,
+    reticulate::np_array(unique(grid_arrows$x)),
+    reticulate::np_array(unique(grid_arrows$y)),
+    reticulate::np_array(matrix(grid_arrows$x_difference, nrow = length(unique(grid_arrows$x)), byrow = F)),
+    reticulate::np_array(matrix(grid_arrows$y_difference, nrow = length(unique(grid_arrows$x)), byrow = F))
+  )
+  vertices <- streamplot$lines$get_segments()
+
+  streamplot_data <- vertices %>% map_dfr(function(matrix) {
+    tibble(
+      x = matrix[1],
+      y = matrix[3],
+      x_future = matrix[2],
+      y_future = matrix[4]
+    )
+  })
+
+  # group into lines and remove duplicate coordinates
+  streamplot_data <- streamplot_data %>%
+    mutate(diff = (lead(x) - x_future) + (lead(y) - y_future)) %>%
+    mutate(line = lag(cumsum(diff != 0), default = 0)) %>%
+    mutate(order = row_number()) %>%
+    group_by(x, y) %>% # remove duplicate x and y coordinates
+    slice(1) %>%
+    ungroup() %>%
+    arrange(order)
+
+  # calculate other statistics for each line
+  streamplot_data <- streamplot_data %>%
+    mutate(length = sqrt((x_future - x) ** 2 + (y_future - y) ** 2)) %>%
+    group_by(line) %>%
+    mutate(
+      cumlength = cumsum(length),
+      percentage = cumlength / max(cumlength)
+    )
+
+  # streamplot_data <- streamplot_data %>% filter(line == 25)
+
+  streamplot_data
+}
+
+#' @export
+stat_velocity_stream <- dynutils::inherit_default_params(
+  list(embed_arrows_stream),
+  function(...) {
+    list(
+      data = function(data) {
+        embed_arrows_stream(
+          attr(data, "data")$cell_info,
+          grid_bandwidth = grid_bandwidth,
+          filter = filter
+        )
+      }
+    )
+  }
+)
+formals(stat_velocity_stream) <- formals(embed_arrows_stream)[2:length(formals(embed_arrows_stream))]

R/helper_calculate_segments_from_edges.R

@@ -0,0 +1,23 @@
+# partitions edges into smaller segments
+calculate_segments_from_edges <- function(edge_positions, n_segments_per_edge = 100) {
+  n_segments_per_edge <- 100
+  segments <- pmap(edge_positions, function(from, to, comp_1_from, comp_2_from, comp_1_to, comp_2_to, ...) {
+    segment_progressions <- tibble(
+      from = from,
+      to = to,
+      percentage = seq(0, 1, length.out = n_segments_per_edge),
+      point_id = paste0(from, "_", to, "_", seq_len(n_segments_per_edge))
+    )
+    segment_positions <- tibble(
+      comp_1 = seq(comp_1_from, comp_1_to, length.out = n_segments_per_edge),
+      comp_2 = seq(comp_2_from, comp_2_to, length.out = n_segments_per_edge),
+      point_id = paste0(from, "_", to, "_", seq_len(n_segments_per_edge))
+    )
+    lst(segment_progressions, segment_positions)
+  })
+
+  lst(
+    segment_progressions = map_dfr(segments, "segment_progressions"),
+    segment_positions = map_dfr(segments, "segment_positions")
+  )
+}

R/helper_linearise_cells.R

@@ -195,51 +195,74 @@ calculate_connections <- function(linearised) {
   # get all connections that are necessary
   # direct connections are those that are reachable without up
   connections <- crossing(
-    linearised$edge_positions %>% select(from, comp_1_from),
-    linearised$edge_positions %>% select(to, comp_1_to)
+    linearised$edge_positions %>% select(.data$from, .data$comp_1_from),
+    linearised$edge_positions %>% select(.data$to, .data$comp_1_to)
   ) %>% filter(
-    from == to,
-    comp_1_from != comp_1_to
+    .data$from == .data$to,
+    .data$comp_1_from != .data$comp_1_to
   ) %>% mutate(
-    comp_1_diff = abs(comp_1_to-comp_1_from)
+    comp_1_diff = abs(.data$comp_1_to - .data$comp_1_from)
   ) %>%
-    arrange(comp_1_diff) %>%
-    select(milestone_id = from, comp_1_from, comp_1_to, comp_1_diff) %>%
+    arrange(.data$comp_1_diff) %>%
+    select(
+      milestone_id = .data$from,
+      .data$comp_1_from,
+      .data$comp_1_to,
+      .data$comp_1_diff
+    ) %>%
     mutate(
       level = NA,
-      direct = near(comp_1_diff, linearised$margin)
+      direct = near(.data$comp_1_diff, linearised$margin)
     )
 
   for (i in seq_len(nrow(connections))) {
     connection <- connections %>% extract_row_to_list(i)
 
     overlapping_connections <- connections %>%
       filter(
-        dplyr::row_number() < i,
-        pmax(comp_1_from, comp_1_to) > min(connection$comp_1_from, connection$comp_1_to),
-        pmin(comp_1_from, comp_1_to) < max(connection$comp_1_from, connection$comp_1_to)
+        row_number() < i,
+        pmax(.data$comp_1_from, .data$comp_1_to) > min(connection$comp_1_from, connection$comp_1_to),
+        pmin(.data$comp_1_from, .data$comp_1_to) < max(connection$comp_1_from, connection$comp_1_to)
       )
 
     if (nrow(overlapping_connections)) {
-      connections$level[i] <- max(overlapping_connections$level) + 1
+      connections$level[[i]] <- max(overlapping_connections$level) + 1
     } else {
-      if (connections$direct[i]) {
-        connections$level[i] <- 0
+      if (connections$direct[[i]]) {
+        connections$level[[i]] <- 0
       } else {
-        connections$level[i] <- 1
+        connections$level[[i]] <- 1
       }
     }
   }
 
   # calculate connection positions
-  connections_direct <- connections %>% filter(direct)
-  connections_indirect <- connections %>% filter(!direct)
+  connections_direct <- connections %>% filter(.data$direct)
+  connections_indirect <- connections %>% filter(!.data$direct)
 
   connection_positions <- bind_rows(
-    connections_direct %>% mutate(connection_ix = 1, comp_2_from = 0, comp_2_to = 0),
-    connections_indirect %>% mutate(comp_1_to = comp_1_from, comp_2_from = 0, comp_2_to = level, connection_ix = 1),
-    connections_indirect %>% mutate(comp_2_from = level, comp_2_to = level, connection_ix = 2),
-    connections_indirect %>% mutate(comp_1_from = comp_1_to, comp_2_from = level, comp_2_to = 0, connection_ix = 3)
+    connections_direct %>% mutate(
+      connection_ix = 1,
+      comp_2_from = 0,
+      comp_2_to = 0
+    ),
+    connections_indirect %>% mutate(
+      comp_1_to = .data$comp_1_from,
+      comp_2_from = 0,
+      comp_2_to = .data$level,
+      connection_ix = 1
+    ),
+    connections_indirect %>% mutate(
+      comp_2_from = .data$level,
+      comp_2_to = .data$level,
+      connection_ix = 2
+    ),
+    connections_indirect %>% mutate(
+      comp_1_from = .data$comp_1_to,
+      comp_2_from = .data$level,
+      comp_2_to = 0,
+      connection_ix = 3
+    )
   )
 
   connection_positions

R/helper_new_scale.R

@@ -1,6 +1,16 @@
+#' Add a new scale
 #' @export
-
+#' @rdname new_scale
 new_scale <- ggnewscale::new_scale
 
 #' @export
-new_scale_fillcolour <- function() {list(new_scale("colour"), new_scale("fill"))}
+#' @rdname new_scale
+new_scale_color <- ggnewscale::new_scale_color
+
+#' @export
+#' @rdname new_scale
+new_scale_colour <- ggnewscale::new_scale_colour
+
+#' @export
+#' @rdname new_scale
+new_scale_fill <- ggnewscale::new_scale_fill

R/helper_rename_dimred_xy.R

@@ -0,0 +1,5 @@
+rename_dimred_xy <- function(df) {
+  colnames(df) <- gsub("^comp_1", "x", colnames(df))
+  colnames(df) <- gsub("^comp_2", "y", colnames(df))
+  df
+}

R/layout_heatmap.R

@@ -16,6 +16,8 @@ layout_heatmap <- function(
   )
 }
 
+#' @importFrom stats kmeans
+#' @importFrom Matrix t
 #' @export
 get_features <- function(dataset, n_modules = 5, num_features = 100) {
   feature_importances <- dynfeature::calculate_overall_feature_importance(dataset)
@@ -36,7 +38,7 @@ get_features <- function(dataset, n_modules = 5, num_features = 100) {
   feature_modules
 }
 
-
+#' @importFrom stats cor
 #' @export
 layout_modules <- function(dataset, feature_modules, cell_layout, margin = 0.02) {
   margin <- margin * nrow(feature_modules)

R/package.R

@@ -1,4 +1,6 @@
-#' Plot all the trajectories
+#' Modular trajectory plotting
+#'
+#' Finetune visualisations of trajectories using ggplot2's grammar of graphics principles.
 #'
 #' @import dplyr
 #' @import tidyr
@@ -12,7 +14,9 @@
 #' @import purrr
 #' @importFrom assertthat assert_that
 #' @importFrom magrittr %<>% %$% set_rownames set_colnames
+#' @importFrom utils head tail
+#' @importFrom rlang %|% quo_text quo eval_tidy
 #'
 #' @docType package
-#' @name dynplot
+#' @name dynplot2
 NULL

R/palettes.R

@@ -1,6 +1,7 @@
 #' @importFrom RColorBrewer brewer.pal
 #' @importFrom grDevices rainbow
 #' @importFrom shades hue
+#' @importFrom utils head
 milestone_palette_list <- list(
   cubeHelix = function(n) cubeHelix(n = n),
   Set3 = function(n) {
@@ -15,7 +16,8 @@ milestone_palette_list <- list(
       # milestone_palette_list$cubeHelix(n)
       all_colors <- grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)]
       all_colors <- all_colors[order(shades::hue(all_colors))][-c(1:2)] # sort and remove white/black
-      all_colors[seq(0, length(all_colors), length(all_colors)/(n+1)) %>% ceiling() %>% head(-1)]
+      ix <- ceiling(seq(0, length(all_colors), length(all_colors)/(n+1)))
+      all_colors[head(ix, -1)]
     }
   }
 )
@@ -24,6 +26,8 @@ milestone_palette_list <- list(
 #' @param n The number of colours to be in the palette.
 #'
 #' @rdname get_milestone_palette_names
+#'
+#' @export
 milestone_palette <- function(n, name = "auto") {
   milestone_palette_list[[name]](n)
 }

R/scale_expression.R

@@ -1,20 +1,113 @@
 # TODO: add rescaling, e.g. quantile
 
-ScaleExpressionFillColour <- ggproto(
-  "ScaleExpressionFillColour",
-  scale_colour_distiller(type = "div", palette = "RdBu"),
-  aesthetics = c("fill", "colour"),
-  map = function(self, x, limits = self$get_limits()) {
-    self$super()$map(x, limits = limits)
-  },
-  oob = scales::squish
+common_colorbar_legend <- guide_colorbar(
+  title.position = "top",
+  title.hjust = 0.5,
+
+  frame.colour = "#333333",
+
+  ticks = FALSE
+)
+
+
+# from https://carto.com/carto-colors/
+#' @importFrom stringr str_split
+ScaleExpressionColour <- ggproto(
+  "ScaleExpressionColour",
+  scale_colour_gradientn(colours = stringr::str_split("#fcde9c,#faa476,#f0746e,#e34f6f,#dc3977,#b9257a,#7c1d6f", ",")[[1]]),
+  # scale_colour_distiller(type = "seq", palette = "OrRd", direction = 1),
+  aesthetics = c("colour"),
+  oob = scales::squish,
+  guide = common_colorbar_legend
+)
+
+
+#' Scale expression values
+#' @export
+#' @rdname scale_expression
+scale_expression_colour <- scale_expression_color <- function(name = "Expression", ...){
+  ggproto(
+    "ScaleExpressionColour",
+    ScaleExpressionColour,
+    name = name,
+    ...
+  )
+}
+
+#' @export
+#' @rdname scale_expression
+scale_expression_color <- scale_expression_colour
+
+#' @importFrom stringr str_split
+ScaleExpressionFill <- ggproto(
+  "ScaleExpressionFill",
+  scale_fill_gradientn(colours = stringr::str_split("#fcde9c,#faa476,#f0746e,#e34f6f,#dc3977,#b9257a,#7c1d6f", ",")[[1]]),
+  # scale_colour_distiller(type = "seq", palette = "OrRd", direction = 1),
+  aesthetics = c("fill"),
+  oob = scales::squish,
+  guide = common_colorbar_legend
 )
 
+
 #' @export
-scale_expression_fillcolour <- function(...){
+#' @rdname scale_expression
+scale_expression_fill <- function(name = "Expression", ...){
   ggproto(
-    "",
-    ScaleExpressionFillColour,
+    "ScaleExpressionFil",
+    ScaleExpressionFill,
+    name = name,
+    ...
+  )
+}
+
+
+
+#' @export
+scale_velocity_fill <- function(name = "Velocity", guide = common_colorbar_legend, ...) {
+  center_limits <- function(limits) {
+    if (abs(limits[[1]]) > abs(limits[[2]])) {
+      limits[[2]] <- abs(limits[[1]])
+    } else if (abs(limits[[2]]) > abs(limits[[1]])) {
+      limits[[1]] <- -abs(limits[[2]])
+    }
+    limits
+  }
+
+  scale <- scale_color_distiller(
+    palette = "RdBu",
+    aesthetics = c("fill"),
+    limits = center_limits,
+    name = name,
+    guide = guide,
+    ...
+  )
+}
+
+#' @export
+scale_velocity_color <- function(name = "Velocity", guide = common_colorbar_legend, ...) {
+  center_limits <- function(limits) {
+    if (abs(limits[[1]]) > abs(limits[[2]])) {
+      limits[[2]] <- abs(limits[[1]])
+    } else if (abs(limits[[2]]) > abs(limits[[1]])) {
+      limits[[1]] <- -abs(limits[[2]])
+    }
+    limits
+  }
+
+  scale <- scale_color_distiller(
+    palette = "RdBu",
+    aesthetics = c("color"),
+    limits = center_limits,
+    name = name,
+    guide = guide,
     ...
   )
 }
+
+
+
+
+ScaleContinuousCentered <- ggproto(
+  "ScaleContinuousCentered",
+  ScaleContinuous
+)

R/scale_milestone.R

@@ -22,12 +22,7 @@ ScaleMilestoneFillColour <- ggproto(
     }
 
     self$milestone_ids <- milestone_ids
-    if (length(milestone_ids) > 0) {
-      self$milestone_colors <- milestone_palette(length(milestone_ids)) %>%
-        set_names(milestone_ids) %>%
-        col2rgb() %>%
-        t()
-    }
+    self$milestone_colors <- define_milestone_colors(self$milestone_colors, self$milestone_ids)
   },
   train = function(self, x) {
     print("train")
@@ -47,12 +42,11 @@ ScaleMilestoneFillColour <- ggproto(
     }
 
     y <- map_chr(x, color_milestone_percentages, milestone_colors = self$milestone_colors)
+
     y
   },
   get_breaks = function(self) {
-    map(self$milestone_ids, function(milestone_id) {
-      tibble(milestone_id = milestone_id, percentage = 1)
-    })
+    milestone_percentage_breaks(self$milestone_ids)
   },
   get_labels = function(self, breaks) {
     self$milestone_ids
@@ -61,12 +55,66 @@ ScaleMilestoneFillColour <- ggproto(
   values = NULL
 )
 
+
+ScaleMilestoneFill <- ggproto(
+  "ScaleMilestoneFill",
+  ScaleMilestoneFillColour,
+  aesthetics = c("fill")
+)
+
+ScaleMilestoneColor <- ggproto(
+  "ScaleMilestoneColor",
+  ScaleMilestoneFillColour,
+  aesthetics = c("colour")
+)
+
+#' Milestone scales
+#' @name scale_milestones
+
+#' @rdname scale_milestones
+#' @export
+scale_milestones_fill <- function(name = "Milestone", milestone_colors = NULL) {
+  ggproto(NULL, ScaleMilestoneFill, name = name, milestone_colors = milestone_colors)
+}
+
+#' @rdname scale_milestones
+#' @export
+scale_milestones_color <- function(name = "Milestone", milestone_colors = NULL) {
+  ggproto(NULL, ScaleMilestoneColor, name = name, milestone_colors = milestone_colors)
+}
+
+#' @rdname scale_milestones
+#' @export
+scale_milestones_colour <- scale_milestones_color
+
+
+#' Helper functions for coloring milestones
+#' @rdname helpers_milestone_coloring
+#' @importFrom grDevices col2rgb
 #' @export
-scale_milestones_fillcolour <- function(name = "Milestone") {
-  ggproto(NULL, ScaleMilestoneFillColour, name = name)
+define_milestone_colors <- function(milestone_colors, milestone_ids) {
+  if (length(milestone_ids) > 0 && is.null(milestone_colors)) {
+    milestone_colors <- milestone_palette(length(milestone_ids)) %>%
+      set_names(milestone_ids) %>%
+      col2rgb() %>%
+      t()
+  } else if(length(milestone_ids) > 0 && is.character(milestone_colors)) {
+    # convert from color to rgb matrix
+    milestone_colors <- milestone_colors %>%
+      col2rgb() %>%
+      t()
+  }
+
+  milestone_colors
 }
 
+#' @rdname helpers_milestone_coloring
+#' @param milestone_percentages A tibble of milestone percentages of a particular cell
+#' @param milestone_colors The matrix linking milestones to RGB, as created by define_milestone_colors
+#' @export
 color_milestone_percentages <- function(milestone_percentages, milestone_colors) {
+  assert_that(!is.null(milestone_colors))
+
   mix_colors <- function(milid, milpct) {
     color_rgb <- apply(milestone_colors[milid,,drop = FALSE], 2, function(x) sum(x * milpct))
     color_rgb[color_rgb < 0] <- 0
@@ -76,3 +124,11 @@ color_milestone_percentages <- function(milestone_percentages, milestone_colors)
 
   mix_colors(as.character(milestone_percentages$milestone_id), milestone_percentages$percentage)
 }
+
+#' @rdname helpers_milestone_coloring
+#' @export
+milestone_percentage_breaks <- function(milestone_ids) {
+  map(milestone_ids, function(milestone_id) {
+    tibble(milestone_id = milestone_id, percentage = 1)
+  })
+}

R/scale_x_.R

@@ -0,0 +1,15 @@
+# TODO: train so that text does not get clipped
+
+ScaleXHeatmap <- ggproto(
+  "ScaleXHeatmap",
+  scale_x_continuous(expand = c(0.5, 0.5))
+)
+
+#' @export
+scale_x_heatmap <- function(...){
+  ggproto(
+    "",
+    ScaleXHeatmap,
+    ...
+  )
+}

R/selectors.R

@@ -1,7 +1,57 @@
 #' @export
-select_expression <- function(feature_id) {
+#' @importFrom stringr str_glue
+select_feature_expression <- function(feature_id, d, expression_source = "expression", scale = dynutils::scale_quantile) {
+  assert_that(
+    !missing(d) && class(d) == "rlang_data_pronoun",
+    msg = str_glue("The second argument of select_feature_expression should be {crayon::italic('.data')} or {crayon::italic('d = .data')}")
+  )
+
+  # to get the correct expression, we need two additional data objects:
+  # - the dataset from which we can get the expression
+  # - the order of the cell ids (which may be duplicated, i.e. in the case of facetting)
+  data <- get("data", envir = d$d[[1]])
+  cell_ids <- d$cell_id
+
+  # make sure only one feature_id is requested, or the feature_ids are coming from facet_data
+  assert_that(length(feature_id) == 1 || (length(feature_id) %% nrow(data)) == 0)
+
   assert_that(!is.null(data))
-  assert_that(feature_id %in% attr(data, "data")$dataset$feature_info$feature_id)
+  assert_that(all(feature_id %in% attr(data, "data")$dataset$feature_info$feature_id))
+
+  expression <- get_expression(attr(data, "data")$dataset, expression_source)[, unique(feature_id), drop=F]
+
+  # scale expression
+  assert_that(is.function(scale))
+  expression <- scale(expression)
+
+  if (length(feature_id) == 1) {
+    feature_id <- rep(feature_id, length(cell_ids))
+  }
+
+  expression_molten <- expression %>%
+    as.matrix() %>%
+    reshape2::melt(varnames = c("cell_id", "feature_id"), value.name = "expression") %>%
+    mutate(cell_id = as.character(cell_id), feature_id = as.character(feature_id))
+
+  # make sure feature_id is a multiple of cell_ids
+  assert_that((length(feature_id) %% length(cell_ids)) == 0)
+
+  # do the left join to get the required expression
+  expression_df <- tibble(cell_id = rep(cell_ids, length(feature_id)/length(cell_ids)), feature_id = feature_id) %>%
+    left_join(expression_molten, c("cell_id", "feature_id"))
+
+  # create output with attributes
+  out <- expression_df$expression
+  out
+}
+
+#' @importFrom stringr str_glue
+#' @export
+select_feature_velocity <- function(feature_id, d) {
+  assert_that(
+    !missing(d) && class(d) == "rlang_data_pronoun",
+    msg = str_glue("The second argument of select_feature_velocity should be {crayon::italic('.data')} or {crayon::italic('d = .data')}")
+  )
 
-  get_expression(attr(data, "data")$dataset, "expression")[attr(data, "data")$cell_info$cell_id, feature_id]
+  select_feature_expression(feature_id, d, "velocity_vector")
 }

R/themes.R

@@ -16,7 +16,7 @@ theme_clean <- function() {
 #' @keywords plot_helpers
 #'
 #' @export
-theme_graph <- function() {
+theme_dynplot <- function() {
   theme_void()
 }
 
@@ -27,5 +27,5 @@ theme_graph <- function() {
 #'
 #' @export
 empty_plot <- function() {
-  ggplot(tibble(x = character())) + geom_point(aes(x, x)) + theme_graph()
+  ggplot(tibble(x = character())) + geom_point(aes(x, x)) + theme_dynplot()
 }

R/utils.R

@@ -1,6 +1,6 @@
 # copied from package rje because it was causing problems
-cubeHelix <- function (n, start = 0.5, r = -1.5, hue = 1, gamma = 1)
-{
+#' @importFrom grDevices rgb
+cubeHelix <- function (n, start = 0.5, r = -1.5, hue = 1, gamma = 1) {
   M = matrix(c(-0.14861, -0.29227, 1.97294, 1.78277, -0.90649, 0), ncol = 2)
   lambda = seq(0, 1, length.out = n)
   l = rep(lambda^gamma, each = 3)

README.Rmd

@@ -1,21 +1,23 @@
 ---
 output: github_document
+editor_options: 
+  chunk_output_type: console
 ---
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
 ```{r setup1, include=FALSE}
-knitr::opts_chunk$set(fig.path=".readme_files/", warning=FALSE, message=FALSE, error=FALSE, echo = TRUE)
+knitr::opts_chunk$set(fig.path = "man/figures/README_", warning=FALSE, message=FALSE, error=FALSE, echo = TRUE)
 ```
 
-<a href="https://travis-ci.org/dynverse/dynplot"><img src="https://travis-ci.org/dynverse/dynplot.svg" align="left"></a> <a href="https://codecov.io/gh/dynverse/dynplot"> <img src="https://codecov.io/gh/dynverse/dynplot/branch/master/graph/badge.svg" align="left" /></a>  [**ℹ️ Tutorials**](https://dynverse.org)  &nbsp;  &nbsp; [**ℹ️ Reference documentation**](https://dynverse.org/reference/dynplot) <br><img src="man/figures/logo.png" align="right" />
+<a href="https://travis-ci.org/dynverse/dynplo2"><img src="https://travis-ci.org/dynverse/dynplo2.svg" align="left"></a> <a href="https://codecov.io/gh/dynverse/dynplo2"> <img src="https://codecov.io/gh/dynverse/dynplo2/branch/master/graph/badge.svg" align="left" /></a>  [**ℹ️ Tutorials**](https://dynverse.org)  &nbsp;  &nbsp; [**ℹ️ Reference documentation**](https://dynverse.org/reference/dynplo2) <br><img src="man/figures/logo.png" align="right" />
 
 # Common visualisation tools for single-cell trajectories
 
 ## Installation
 
 ```r
-devtools::install_github("dynverse/dynplot")
+devtools::install_github("dynverse/dynplot2")
 ```
 
 On linux, udunits2 has to be installed:
@@ -29,98 +31,117 @@ The package provides different ways to plot both the topology and cellular prope
 
 ```{r, echo=FALSE}
 library(dplyr)
-library(dynplot)
+library(dynplot2)
 library(dynutils)
 library(dynwrap)
 library(tidyverse)
 
 library(patchwork)
 
+set.seed(1)
+
 data(example_bifurcating)
-trajectory <- example_bifurcating
-
-trajectory <- trajectory %>% add_root()
-grouping <- trajectory$prior_information$groups_id
-groups <- tibble(dynwrap::group_onto_nearest_milestones(trajectory)) %>% 
-  mutate(color=dynplot:::milestone_palette_list$auto(n()))
-features_oi <- apply(as.matrix(trajectory$counts), 2, sd) %>% sort() %>% names() %>% tail(10)
-feature_oi <- features_oi[[10]]
+dataset <- example_bifurcating %>% add_root() %>% add_dimred(dimred = dyndimred::dimred_landmark_mds)
 ```
 
+```{r}
+dynplot_dimred(dataset) +
+  geom_cell_point(aes(colour = select_feature_expression("G1", d = .data))) +
+  scale_expression_color() + # a scale has to be given here, otherwise error
+  new_scale_color() +
+  geom_trajectory_segments(aes(colour = milestone_percentages), size = 2) +
+  geom_milestone_label(aes(fill = milestone_id)) +
+  scale_milestones_color() +
+  scale_milestones_fill()
+```
 
-```{r cells, fig.width=15, fig.height=10, echo=FALSE}
-empty <- ggplot() + theme_void()
-question <- empty + geom_text(aes(0, 0, label="????"))
-unavailable <- empty + geom_text(aes(0, 0, label="----"))
-
-title <- list(
-  dynplot:::empty_plot(), 
-  # empty + geom_label(aes(0, 0, label="Topology")),
-  empty + geom_label(aes(0, 0, label="Ordering")),
-  empty + geom_label(aes(0, 0, label="Grouping/clustering")),
-  empty + geom_label(aes(0, 0, label="Expression of \n a single gene")), 
-  empty + geom_label(aes(0, 0, label="Pseudotime"))
-)
-
-dendro <- list(
-  empty + geom_label(aes(0, 0, label="Dendrogram")),
-  # plot_dendro(trajectory),
-  plot_dendro(trajectory, "milestone"),
-  plot_dendro(trajectory, grouping=grouping, groups=groups),
-  plot_dendro(trajectory, feature_oi=feature_oi),
-  plot_dendro(trajectory, "pseudotime")
-)
-
-onedim <- list(
-    empty + geom_label(aes(0, 0, label="Onedim")),
-    # plot_onedim(trajectory),
-    plot_onedim(trajectory, "milestone"),
-    plot_onedim(trajectory, grouping=grouping, groups=groups),
-    plot_onedim(trajectory, feature_oi=feature_oi),
-    plot_onedim(trajectory, "pseudotime")
-  )
-
-graph <- list(
-  empty + geom_label(aes(0, 0, label="Graph")),
-  # plot_graph(trajectory),
-  plot_graph(trajectory, "milestone"),
-  plot_graph(trajectory, grouping=grouping, groups=groups),
-  plot_graph(trajectory, feature_oi=feature_oi),
-  plot_graph(trajectory, "pseudotime")
-)
-
-dimred <- list(
-  empty + geom_label(aes(0, 0, label="Dimensionality reduction")),
-  # plot_dimred(trajectory),
-  plot_dimred(trajectory, "milestone"),
-  plot_dimred(trajectory, grouping=grouping, groups=groups),
-  plot_dimred(trajectory, feature_oi=feature_oi),
-  plot_dimred(trajectory, "pseudotime")
-)
-
-wrap_plots(
-  c(title, dendro, onedim, graph, dimred),
-  nrow=5,
-  byrow=TRUE
-) & theme(legend.position = "none")
+```{r}
+# dynplot_dimred(dataset) +
+#   geom_cell_point(color = "grey80") +
+#   geom_trajectory_segments(aes(colour = milestone_percentages), size = 2) +
+#   geom_milestone_label(aes(fill = milestone_id)) +
+#   scale_milestones_fill() +
+#   scale_milestones_colour() +
+#   geom_velocity_arrow(stat = stat_velocity_grid(grid_n = 20))
 ```
 
+```{r}
+dynplot_graph(dataset) +
+  geom_trajectory_divergence() +
+  geom_trajectory_segments(size = 2, color = "#333333", arrow_size = 0.5) +
+  geom_cell_point(aes(colour = milestone_percentages)) +
+  geom_milestone_label(aes(fill = milestone_id)) +
+  scale_milestones_fill() +
+  scale_milestones_color()
+```
+
+```{r}
+dynplot_onedim(dataset, equal_cell_width = TRUE, margin = 0.05) +
+  geom_trajectory_segments(aes(color = milestone_percentages), size = 5, color = "#333333") +
+  geom_trajectory_connection() +
+  scale_milestones_colour() +
+  new_scale_colour() +
+  geom_cell_point(aes(colour = select_feature_expression("G2", d = .data))) +
+  scale_expression_colour() +
+  new_scale_colour() +
+  geom_milestone_label(aes(y = -0.1, fill = milestone_id, hjust = as.integer(type == "end"))) +
+  scale_milestones_fill()
+```
+
+```{r}
+dynplot_dendro(dataset) +
+  geom_trajectory_segments(aes(color = milestone_percentages), size = 5, color = "#333333") +
+  geom_trajectory_connection() +
+  scale_milestones_colour() +
+  new_scale_colour() +
+  geom_cell_point(aes(colour = select_feature_expression("G1", d = .data))) +
+  scale_expression_colour() +
+  new_scale_colour() +
+  # geom_cell_point(aes(colour = milestone_percentages)) +
+  geom_milestone_label(aes(fill = milestone_id)) +
+  scale_milestones_fill()
+```
+
+```{r}
+dynplot_dendro(dataset) +
+  geom_trajectory_segments(aes(color = milestone_percentages), size = 5, color = "#333333") +
+  geom_trajectory_connection() +
+  scale_milestones_colour() +
+  new_scale_colour() +
+  geom_cell_point(aes(y = y + dynutils::scale_minmax(select_feature_expression("G2", d=.data)) * 0.5 - 0.25, colour = select_feature_expression("G2", d = .data))) +
+  scale_expression_colour() +
+  new_scale_colour() +
+  # geom_cell_point(aes(colour = milestone_percentages)) +
+  geom_milestone_label(aes(fill = milestone_id)) +
+  scale_milestones_fill()
+```
 
-And to plot the expression and feature importances of many genes along the trajectory
 
-```{r heatmap, echo=FALSE, fig.width=10, fig.height=6}
-cell_feature_importance <- dynfeature::calculate_cell_feature_importance(trajectory %>% dynwrap::add_waypoints())
-plot_heatmap(trajectory, grouping = trajectory$prior_information$grouping_assignment)
+```{r}
+# cell_layout <- layout_onedim(dataset)
+# feature_modules <- get_features(dataset)
+# feature_layout <- layout_modules(dataset, feature_modules = feature_modules, cell_layout = cell_layout)
+# layout <- layout_heatmap(dataset, feature_layout = feature_layout)
+# 
+# dynplot(dataset, layout = layout) +
+#   geom_trajectory_segments(aes(color = milestone_percentages)) +
+#   geom_trajectory_connection() +
+#   geom_milestone_label(aes(fill = milestone_id, hjust = as.integer(type == "end"))) +
+#   scale_milestones_colour() +
+#   new_scale_fill() +
+#   geom_expression_raster() +
+#   scale_expression_fill() +
+#   new_scale_fill() +
+#   geom_tile(aes(x = x, y = 1))
 ```
 
 ## Latest changes
 
-Check out `news(package = "dynwrap")` or [NEWS.md](inst/NEWS.md) for a full list of changes.
+Check out `news(package = "dynplot2")` or [NEWS.md](NEWS.md) for a full list of changes.
 
-<!-- This section gets automatically generated from inst/NEWS.md, and also generates inst/NEWS -->
+<!-- This section gets automatically generated from NEWS.md -->
 
 ```{r news, echo=FALSE, results="asis"}
-dynutils::update_news()
 cat(dynutils::recent_news())
 ```
 

README.md

@@ -1,19 +1,19 @@
 
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
-<a href="https://travis-ci.org/dynverse/dynplot"><img src="https://travis-ci.org/dynverse/dynplot.svg" align="left"></a>
-<a href="https://codecov.io/gh/dynverse/dynplot">
-<img src="https://codecov.io/gh/dynverse/dynplot/branch/master/graph/badge.svg" align="left" /></a>
+<a href="https://travis-ci.org/dynverse/dynplo2"><img src="https://travis-ci.org/dynverse/dynplo2.svg" align="left"></a>
+<a href="https://codecov.io/gh/dynverse/dynplo2">
+<img src="https://codecov.io/gh/dynverse/dynplo2/branch/master/graph/badge.svg" align="left" /></a>
 [**ℹ️ Tutorials**](https://dynverse.org)     [**ℹ️ Reference
-documentation**](https://dynverse.org/reference/dynplot)
+documentation**](https://dynverse.org/reference/dynplo2)
 <br><img src="man/figures/logo.png" align="right" />
 
 # Common visualisation tools for single-cell trajectories
 
 ## Installation
 
 ``` r
-devtools::install_github("dynverse/dynplot")
+devtools::install_github("dynverse/dynplot2")
 ```
 
 On linux, udunits2 has to be installed:
@@ -26,20 +26,117 @@ On linux, udunits2 has to be installed:
 The package provides different ways to plot both the topology and
 cellular properties of a trajectory:
 
-![](.readme_files/cells-1.png)<!-- -->
+``` r
+dynplot_dimred(dataset) +
+  geom_cell_point(aes(colour = select_feature_expression("G1", d = .data))) +
+  scale_expression_color() + # a scale has to be given here, otherwise error
+  new_scale_color() +
+  geom_trajectory_segments(aes(colour = milestone_percentages), size = 2) +
+  geom_milestone_label(aes(fill = milestone_id)) +
+  scale_milestones_color() +
+  scale_milestones_fill()
+```
 
-And to plot the expression and feature importances of many genes along
-the trajectory
+![](man/figures/README_unnamed-chunk-2-1.png)<!-- -->
 
-![](.readme_files/heatmap-1.png)<!-- -->
+``` r
+# dynplot_dimred(dataset) +
+#   geom_cell_point(color = "grey80") +
+#   geom_trajectory_segments(aes(colour = milestone_percentages), size = 2) +
+#   geom_milestone_label(aes(fill = milestone_id)) +
+#   scale_milestones_fill() +
+#   scale_milestones_colour() +
+#   geom_velocity_arrow(stat = stat_velocity_grid(grid_n = 20))
+```
+
+``` r
+dynplot_graph(dataset) +
+  geom_trajectory_divergence() +
+  geom_trajectory_segments(size = 2, color = "#333333", arrow_size = 0.5) +
+  geom_cell_point(aes(colour = milestone_percentages)) +
+  geom_milestone_label(aes(fill = milestone_id)) +
+  scale_milestones_fill() +
+  scale_milestones_color()
+```
+
+![](man/figures/README_unnamed-chunk-4-1.png)<!-- -->
+
+``` r
+dynplot_onedim(dataset, equal_cell_width = TRUE, margin = 0.05) +
+  geom_trajectory_segments(aes(color = milestone_percentages), size = 5, color = "#333333") +
+  geom_trajectory_connection() +
+  scale_milestones_colour() +
+  new_scale_colour() +
+  geom_cell_point(aes(colour = select_feature_expression("G2", d = .data))) +
+  scale_expression_colour() +
+  new_scale_colour() +
+  geom_milestone_label(aes(y = -0.1, fill = milestone_id, hjust = as.integer(type == "end"))) +
+  scale_milestones_fill()
+```
+
+![](man/figures/README_unnamed-chunk-5-1.png)<!-- -->
+
+``` r
+dynplot_dendro(dataset) +
+  geom_trajectory_segments(aes(color = milestone_percentages), size = 5, color = "#333333") +
+  geom_trajectory_connection() +
+  scale_milestones_colour() +
+  new_scale_colour() +
+  geom_cell_point(aes(colour = select_feature_expression("G1", d = .data))) +
+  scale_expression_colour() +
+  new_scale_colour() +
+  # geom_cell_point(aes(colour = milestone_percentages)) +
+  geom_milestone_label(aes(fill = milestone_id)) +
+  scale_milestones_fill()
+```
+
+![](man/figures/README_unnamed-chunk-6-1.png)<!-- -->
+
+``` r
+dynplot_dendro(dataset) +
+  geom_trajectory_segments(aes(color = milestone_percentages), size = 5, color = "#333333") +
+  geom_trajectory_connection() +
+  scale_milestones_colour() +
+  new_scale_colour() +
+  geom_cell_point(aes(y = y + dynutils::scale_minmax(select_feature_expression("G2", d=.data)) * 0.5 - 0.25, colour = select_feature_expression("G2", d = .data))) +
+  scale_expression_colour() +
+  new_scale_colour() +
+  # geom_cell_point(aes(colour = milestone_percentages)) +
+  geom_milestone_label(aes(fill = milestone_id)) +
+  scale_milestones_fill()
+```
+
+![](man/figures/README_unnamed-chunk-7-1.png)<!-- -->
+
+``` r
+# cell_layout <- layout_onedim(dataset)
+# feature_modules <- get_features(dataset)
+# feature_layout <- layout_modules(dataset, feature_modules = feature_modules, cell_layout = cell_layout)
+# layout <- layout_heatmap(dataset, feature_layout = feature_layout)
+# 
+# dynplot(dataset, layout = layout) +
+#   geom_trajectory_segments(aes(color = milestone_percentages)) +
+#   geom_trajectory_connection() +
+#   geom_milestone_label(aes(fill = milestone_id, hjust = as.integer(type == "end"))) +
+#   scale_milestones_fillcolour() +
+#   new_scale_fillcolour() +
+#   geom_expression_raster() +
+#   scale_expression_fillcolour() +
+#   new_scale_fillcolour() +
+#   geom_tile(aes(x = x, y = 1))
+```
 
 ## Latest changes
 
-Check out `news(package = "dynwrap")` or [NEWS.md](inst/NEWS.md) for a
-full list of
-changes.
+Check out `news(package = "dynplot2")` or [NEWS.md](NEWS.md) for a full
+list of changes.
+
+<!-- This section gets automatically generated from NEWS.md -->
+
+### Recent changes in dynplot2 1.0.0 (unreleased)
 
-<!-- This section gets automatically generated from inst/NEWS.md, and also generates inst/NEWS -->
+  - Spin-off from dynplot, dynplot2 provides a more modular interface to
+    visualize trajectories
 
 ### Recent changes in dynplot 1.0.0 (28-03-2019)