Add timeline component and plotting functions

NAMESPACE

@@ -11,12 +11,16 @@ export(deconvSingle)
 export(findHVG)
 export(findPure)
 export(getCenters)
+export(getClusterOrder)
 export(getEstimates)
 export(getMajorMatrix)
 export(getPval)
 export(getTheta)
+export(getTimeline)
 export(heatmapKseq)
 export(plotContTable)
+export(plotGeneTimeline)
+export(plotMultipleGeneTimeline)
 export(predictTheta)
 export(runDescend)
 export(scaleRowSums)

R/SOUP.R

@@ -307,3 +307,5 @@ getSimilarity <- function(expr, type="log") {
   return (A) 
 }
 
+
+

R/SOUPtime.R

@@ -0,0 +1,65 @@
+
+#' Estimate Developmental Timeline
+#' 
+#' Estimate the developmental timeline for cells from SOUP soft membership.
+#' 
+#' @param membership An n.cell-by-K soft membership matrix
+#' @param centers A K-by-n.gene cluster center matrix
+#' @param k.start Which cluster is the starting point? 
+#'     One, and only one of \code{k.start} and \code{k.end} must be specified.
+#' @param k.end Which cluster is the ending point?
+#'     One, and only one of \code{k.start} and \code{k.end} must be specified.
+#'     
+#' @return A vector of timepoints for all cells.
+#' @export
+getTimeline <- function(membership, centers, 
+                        k.start=NULL, k.end=NULL) {
+
+  K = nrow(centers)
+  if (ncol(membership) != K) {
+    stop("Error: Membership must have the same ncol as centers.\n")
+  }
+
+  ## re-order clusters
+  cor.center = cor(t(centers))
+  if (!is.null(k.start)) {
+    k.order = getClusterOrder(cor.center, k.start)
+  } else if (!is.null(k.end)) {
+    ## start from end, then reverse
+    k.order = getClusterOrder(cor.center, k.end)
+    k.order = k.order[c(K:1)]
+  }
+
+  ## timeline
+  timeline = rowSums(membership[, k.order] %*% diag(c(1:K)))
+  names(timeline) <- row.names(membership)
+
+  return(timeline)
+}
+
+#' Order Clusters
+#' 
+#' Re-order cluster labels, starting from \code{k.start}, 
+#' and consecutively append the next one with the highest correlation in cluster centers.
+#' 
+#' @param cor.center The K-by-K correlation matrix among cluster centers
+#' @param k.start The cluster to start with
+#' 
+#' @return A list with length K that re-orders the clusters. 
+#' @export
+getClusterOrder <- function(cor.center, ## K-by-K correlation matrix of centers
+                            k.start=1 ## which cluster to start with
+) {
+  K = ncol(cor.center)
+  ## order clusters
+  colnames(cor.center) = c(1:K)
+  k.order = c(k.start)
+  k.cur = k.start
+  while (length(k.order) < K) {
+    k.next = names(which.max(cor.center[k.cur, -k.order]))
+    k.next = as.numeric(k.next)
+    k.order = c(k.order, k.next)
+    k.cur = k.next
+  }
+  return (k.order)
+}

R/utils_plot.R

@@ -166,4 +166,66 @@ plotContTable <- function(est_label, true_label, short.names=NULL, xlab="Referen
     )
 
   return(g)
-}
+}
+
+#' Plot Gene Expression along Trajectory
+#' 
+#' @param expr The cell-by-gene expression matrix to be visualized
+#' @param marker.gene The gene to be visualized; must be one of the columns of \code{expr}
+#' @param timeline Timepoints for all cells, a vector
+#' @param x.title (optional) Label of x-axis
+#' @param y.title (optional) Label of y-axis
+#' @param title (optional) Plot title
+#' 
+#' @return A ggplot object
+#' @export
+plotGeneTimeline <- function(expr, marker.gene, timeline, 
+                             x.title="SOUP trajectory", y.title="Expression", title="") {
+  df = data.frame(Expression=expr[, marker.gene],
+                  Trajectory=timeline)
+  g = ggplot(df, aes(x=Trajectory, y=Expression)) +
+    geom_point(aes(color=Trajectory)) +
+    labs(x=x.title, y=y.title, title=title) +
+    scale_colour_gradient(low="#56B1F7", high="#132B43") +
+    theme_bw() +
+    theme(panel.background = element_blank(),
+          panel.border = element_blank(),
+          panel.grid.minor = element_blank(),
+          panel.grid.major = element_blank(),
+          axis.line = element_line(colour = "black"),
+          axis.title.x=element_text(size=12),
+          axis.title.y=element_text(size=12),
+          legend.title=element_text(size=12),
+          legend.text=element_text(size=12),
+          plot.title = element_text(size=12, face="bold", hjust = 0.5))
+  return(g)
+}
+
+#' Plot Multiple Gene Expressions along Trajectory
+#' 
+#' Expression levels of multiple genes along trajectory, arranged on a grid of plots.
+#' 
+#' @param expr The cell-by-gene expression matrix to be visualized
+#' @param genelist A list of genes to be visualized; must be among the columns of \code{expr}
+#' @param timeline Timepoints for all cells, a vector
+#' @param x.title (optional) Label of x-axis
+#' @param y.title (optional) Label of y-axis
+#' @param nrow (optional) number of rows in the plot grid.
+#' @param ncol (optional) number of cols in the plot grid.
+#' 
+#' @return A ggplot object
+#' @export
+plotMultipleGeneTimeline <- function(expr, genelist, timeline,
+                                     x.title="SOUP trajectory", y.title="Expression", 
+                                     nrow=NULL, ncol=NULL) {
+  g.list = list()
+  for (gene in genelist) {
+    g.list = c(g.list, 
+               list(plotGeneTimeline(expr=expr, marker.gene=gene, timeline=timeline, 
+                                x.title=x.title, y.title=y.title, title=gene)))
+  }
+
+  g.multi = ggarrange(plotlist=g.list, nrow=nrow, ncol=ncol,
+                      common.legend=TRUE, legend="right")
+  return(g.multi)
+}