@@ -153,14 +153,8 @@ runTSCAN <- function(inSCE,

     sce <- runScranSNN(inSCE, useReducedDim = useReducedDim, seed = seed)

     cluster <- colData(sce)$scranSNN_cluster

     rm(sce)

-  } else if (length(cluster) == 1 & is.character(cluster)) {

-    if (!cluster %in% names(colData(inSCE))) {

-      stop("Specified cluster not found in `colData(inSCE)`.")

-    }

-    cluster <- colData(inSCE)[[cluster]]

-  } else if (length(cluster) != ncol(inSCE)) {

-    stop("`cluster` should either has length equal to `ncol(inSCE)` or be a ",

-         "colData variable.")

+  } else {

+    cluster <- .manageCellVar(inSCE, var = cluster)

   }

   if (length(unique(cluster)) == 1) {

     stop("Only one cluster found. Unable to calculate.")

@@ -443,6 +443,7 @@ plotTSCANPseudotimeHeatmap <- function(inSCE,

 #' Use the viridis color scale to match with `plotTSCANResult` embedding color

 #' @param x numeric vector of pseudotime

 #' @return function object returned by circlize::colorRamp2

+#' @noRd

 .viridisPseudoTimeColor <- function(x) {

   a <- max(x, na.rm = TRUE)

   b <- min(x, na.rm = TRUE)

@@ -862,6 +863,7 @@ plotTSCANDimReduceFeatures <- function(

 #' @param useCluster Vector of clusters, that are going to be the vertices for

 #' edge finding.

 #' @return data.frame, subset rows of line.data

+#' @noRd

 .getClustersLineData <- function(line.data, useCluster) {

   idx <- vapply(useCluster, function(x){

     grepl(paste0("^", x, "--"), line.data$edge) |

@@ -442,6 +442,7 @@ plotTSCANPseudotimeHeatmap <- function(inSCE,

 #' Basing on pseudotime is presented within range from 0 to 100

 #' Use the viridis color scale to match with `plotTSCANResult` embedding color

 #' @param x numeric vector of pseudotime

+#' @return function object returned by circlize::colorRamp2

 .viridisPseudoTimeColor <- function(x) {

   a <- max(x, na.rm = TRUE)

   b <- min(x, na.rm = TRUE)

@@ -652,10 +652,7 @@ plotTSCANClusterPseudo <- function(inSCE, useCluster, useReducedDim = "UMAP",

   by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = clusters)

   line.data <- TSCAN::reportEdges(by.cluster, mst = results$mst,

                                   clusters = NULL, use.dimred = useReducedDim)

-  line.data.sub <- line.data[grepl(paste0("^", useCluster, "--"),

-                                   line.data$edge) |

-                               grepl(paste0("--", useCluster, "$"),

-                                     line.data$edge),]

+  line.data.sub <- .getClustersLineData(line.data, useCluster)

   # Get Branch pseudotime

   tscan.pseudo <- TSCAN::orderCells(results$maptscan, mst = results$mst,

@@ -694,6 +691,9 @@ plotTSCANClusterPseudo <- function(inSCE, useCluster, useReducedDim = "UMAP",

 #' @param inSCE Input \linkS4class{SingleCellExperiment} object.

 #' @param useCluster Choose a cluster used for identifying DEG with

 #' \code{\link{runTSCANClusterDEAnalysis}}. Required.

+#' @param pathIndex Specifies one of the branching paths from \code{useCluster}

+#' and plot the top DEGs on this path. Ususally presented by the terminal

+#' cluster of a path. By default \code{NULL} plot top DEGs of all paths.

 #' @param useReducedDim A single character for the matrix of 2D embedding.

 #' Should exist in \code{reducedDims} slot. Default \code{"UMAP"}.

 #' @param topN Integer. Use top N genes identified. Default \code{9}.

@@ -802,55 +802,50 @@ plotTSCANDimReduceFeatures <- function(

   combinePlot = c("all", "none"))

 {

   results <- getTSCANResults(inSCE, analysisName = "Pseudotime")

-  if (!is.null(useCluster) && length(useCluster) != 1)

-    stop("`useCluster`: can only use one cluster.")

   combinePlot <- match.arg(combinePlot)

   clusters <- colData(inSCE)$TSCAN_clusters

   by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = clusters)

   line.data <- TSCAN::reportEdges(by.cluster, mst = results$mst,

                                   clusters = NULL, use.dimred = useReducedDim)

   if (!is.null(useCluster)) {

-    line.data <- line.data[grepl(paste0("^", useCluster, "--"),

-                                 line.data$edge) |

-                             grepl(paste0("--", useCluster, "$"),

-                                   line.data$edge),]

+    line.data <- .getClustersLineData(line.data, useCluster)

     inSCE <- inSCE[, clusters %in% useCluster]

   }

   if (by == "rownames") by <- NULL

   plotList <- list()

-  if (!is.na(features)) {

-    for (f in features) {

-      g <- plotSCEDimReduceFeatures(inSCE, feature = f,

-                                    useAssay = useAssay,

-                                    featureLocation = by,dim1 = 1, dim2 = 2,

-                                    featureDisplay = featureDisplay,

-                                    reducedDimName = useReducedDim,

-                                    title = f) +

-        ggplot2::geom_line(data = line.data,

-                           ggplot2::aes_string(x = colnames(line.data)[2],

-                                               y = colnames(line.data)[3],

-                                               group = "edge"),

-                           inherit.aes = FALSE)

-      # Text labeling code credit: scater

-      reddim <- as.data.frame(SingleCellExperiment::reducedDim(inSCE,

-                                                               useReducedDim))

-      text_out <- scater::retrieveCellInfo(inSCE, "TSCAN_clusters",

-                                           search = "colData")

-      by_text_x <- vapply(split(reddim[,1], text_out$val),

-                          stats::median, FUN.VALUE = 0)

-      by_text_y <- vapply(split(reddim[,2], text_out$val),

-                          stats::median, FUN.VALUE = 0)

-      g <- g + ggrepel::geom_text_repel(

-        data = data.frame(x = by_text_x,

-                          y = by_text_y,

-                          label = names(by_text_x)),

-        mapping = ggplot2::aes_string(x = "x",

-                                      y = "y",

-                                      label = "label"),

-        inherit.aes = FALSE

-      )

-      plotList[[f]] <- g

-    }

+  features <- stats::na.omit(features)

+  for (f in features) {

+    # Should not enter the loop if features is length zero after NA omit

+    g <- plotSCEDimReduceFeatures(inSCE, feature = f,

+                                  useAssay = useAssay,

+                                  featureLocation = by,dim1 = 1, dim2 = 2,

+                                  featureDisplay = featureDisplay,

+                                  reducedDimName = useReducedDim,

+                                  title = f) +

+      ggplot2::geom_line(data = line.data,

+                         ggplot2::aes_string(x = colnames(line.data)[2],

+                                             y = colnames(line.data)[3],

+                                             group = "edge"),

+                         inherit.aes = FALSE)

+    # Text labeling code credit: scater

+    reddim <- as.data.frame(SingleCellExperiment::reducedDim(inSCE,

+                                                             useReducedDim))

+    text_out <- scater::retrieveCellInfo(inSCE, "TSCAN_clusters",

+                                         search = "colData")

+    by_text_x <- vapply(split(reddim[,1], text_out$val),

+                        stats::median, FUN.VALUE = 0)

+    by_text_y <- vapply(split(reddim[,2], text_out$val),

+                        stats::median, FUN.VALUE = 0)

+    g <- g + ggrepel::geom_text_repel(

+      data = data.frame(x = by_text_x,

+                        y = by_text_y,

+                        label = names(by_text_x)),

+      mapping = ggplot2::aes_string(x = "x",

+                                    y = "y",

+                                    label = "label"),

+      inherit.aes = FALSE, na.rm = TRUE,

+    )

+    plotList[[f]] <- g

   }

   if (combinePlot == "all") {

     if (length(plotList) > 0)

@@ -858,3 +853,22 @@ plotTSCANDimReduceFeatures <- function(

   }

   return(plotList)

 }

+

+#' Extract edges that connects to the clusters of interest

+#' @param line.data data.frame of three columns. First column should be named

+#' "edge", the other two are coordinates of one of the vertices that this

+#' edge connects.

+#' @param useCluster Vector of clusters, that are going to be the vertices for

+#' edge finding.

+#' @return data.frame, subset rows of line.data

+.getClustersLineData <- function(line.data, useCluster) {

+  idx <- vapply(useCluster, function(x){

+    grepl(paste0("^", x, "--"), line.data$edge) |

+      grepl(paste0("--", x, "$"), line.data$edge)

+  }, logical(nrow(line.data)))

+  # `idx` returned was c("matrix", "array") class. Each column is a logical

+  # vector for edge selection of each cluster.

+  # Next, do "or" for each row.

+  idx <- rowSums(idx) > 0

+  line.data[idx,]

+}

@@ -21,7 +21,7 @@

 #' \code{listTSCANTerminalNodes}.

 #'

 #' When \code{analysisName = "ClusterDEAnalysis"}, returns the list result from

-#' \code{\link{runTSCANBranchDEG}}. Here \code{pathName} needs to match

+#' \code{\link{runTSCANClusterDEAnalysis}}. Here \code{pathName} needs to match

 #' with the \code{useCluster} argument when running the algorithm.

 #' @export

 #' @rdname getTSCANResults

@@ -50,10 +50,11 @@ setMethod("getTSCANResults", signature(x = "SingleCellExperiment"),

     results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN$Pseudotime

   else {

     if (is.null(pathName)) {

-      stop("Have to specify `pathName` for TSCAN DE analyses")

+      results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN[[analysisName]]

+    } else {

+      pathName <- as.character(pathName)

+      results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN[[analysisName]][[pathName]]

     }

-    pathName <- as.character(pathName)

-    results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN[[analysisName]][[pathName]]

   }

   return(results)

 })

@@ -161,7 +162,9 @@ runTSCAN <- function(inSCE,

     stop("`cluster` should either has length equal to `ncol(inSCE)` or be a ",

          "colData variable.")

   }

-

+  if (length(unique(cluster)) == 1) {

+    stop("Only one cluster found. Unable to calculate.")

+  }

   message(date(), " ... Running TSCAN to estimate pseudotime")

   inSCE <- scran::computeSumFactors(inSCE, clusters = cluster)

   by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = cluster)

@@ -205,12 +208,12 @@ runTSCAN <- function(inSCE,

   edgeList <- mst

   for (i in seq_along(unique(colData(inSCE)$TSCAN_clusters))){

     degree <- sum(edgeList[i] != 0)

-    if (degree > 2) branchClusters <- c(branchClusters, i)

+    if (degree > 1) branchClusters <- c(branchClusters, i)

   }

   ## Save these results in a list and then make S4 accessor that passes the

   ## entire list

-  result <- list(pseudo = tscan.pseudo, mst = mst, clusters = cluster,

+  result <- list(pseudo = tscan.pseudo, mst = mst,

                  maptscan = map.tscan,

                  pathClusters = pathClusters,

                  branchClusters = branchClusters,

@@ -246,7 +249,8 @@ runTSCAN <- function(inSCE,

 plotTSCANResults <- function(inSCE, useReducedDim = "UMAP") {

   results <- getTSCANResults(inSCE, analysisName = "Pseudotime")

-  by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = results$clusters)

+  clusters <- colData(inSCE)$TSCAN_clusters

+  by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = clusters)

   line.data <- TSCAN::reportEdges(by.cluster, mst = results$mst,

                                   clusters = NULL, use.dimred = useReducedDim)

@@ -543,12 +547,12 @@ plotTSCANPseudotimeGenes <- function(inSCE,

 #' data("mouseBrainSubsetSCE", package = "singleCellTK")

 #' mouseBrainSubsetSCE <- runTSCAN(inSCE = mouseBrainSubsetSCE,

 #'                                 useReducedDim = "PCA_logcounts")

-#' mouseBrainSubsetSCE <- runTSCANBranchDEG(inSCE = mouseBrainSubsetSCE,

+#' mouseBrainSubsetSCE <- runTSCANClusterDEAnalysis(inSCE = mouseBrainSubsetSCE,

 #'                                          useCluster = 1)

-runTSCANBranchDEG <- function(inSCE,

-                              useCluster,

-                              useAssay = "logcounts",

-                              fdrThreshold = 0.05){

+runTSCANClusterDEAnalysis <- function(inSCE,

+                                      useCluster,

+                                      useAssay = "logcounts",

+                                      fdrThreshold = 0.05) {

   results <- getTSCANResults(inSCE, analysisName = "Pseudotime")

   message(date(), " ... Finding DEG between TSCAN branches")

   tscan.pseudo <- TSCAN::orderCells(results$maptscan, mst = results$mst,

@@ -637,14 +641,15 @@ runTSCANBranchDEG <- function(inSCE,

 #' data("mouseBrainSubsetSCE", package = "singleCellTK")

 #' mouseBrainSubsetSCE <- runTSCAN(inSCE = mouseBrainSubsetSCE,

 #'                                 useReducedDim = "PCA_logcounts")

-#' plotTSCANBranchPseudo(mouseBrainSubsetSCE, useCluster = 1,

-#'                       useReducedDim = "TSNE_logcounts")

-plotTSCANBranchPseudo <- function(inSCE, useCluster, useReducedDim = "UMAP",

+#' plotTSCANClusterPseudo(mouseBrainSubsetSCE, useCluster = 1,

+#'                        useReducedDim = "TSNE_logcounts")

+plotTSCANClusterPseudo <- function(inSCE, useCluster, useReducedDim = "UMAP",

                                    combinePlot = c("all", "none")){

   # Get the plotting info of edges

   results <- getTSCANResults(inSCE, analysisName = "Pseudotime")

   combinePlot <- match.arg(combinePlot)

-  by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = results$clusters)

+  clusters <- colData(inSCE)$TSCAN_clusters

+  by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = clusters)

   line.data <- TSCAN::reportEdges(by.cluster, mst = results$mst,

                                   clusters = NULL, use.dimred = useReducedDim)

   line.data.sub <- line.data[grepl(paste0("^", useCluster, "--"),

@@ -681,20 +686,20 @@ plotTSCANBranchPseudo <- function(inSCE, useCluster, useReducedDim = "UMAP",

 ###  plot gene of interest in branch cluster

 ###################################################

-#' @title Plot features identified by \code{\link{runTSCANBranchDEG}} on cell

-#' 2D embedding with MST overlaid

+#' @title Plot features identified by \code{\link{runTSCANClusterDEAnalysis}} on

+#' cell 2D embedding with MST overlaid

 #' @description A wrapper function which plot the top features expression

-#' identified by \code{\link{runTSCANBranchDEG}} on the 2D embedding of the

-#' cells cluster used in the analysis. The related MST edges are overlaid.

+#' identified by \code{\link{runTSCANClusterDEAnalysis}} on the 2D embedding of

+#' the cells cluster used in the analysis. The related MST edges are overlaid.

 #' @param inSCE Input \linkS4class{SingleCellExperiment} object.

 #' @param useCluster Choose a cluster used for identifying DEG with

-#' \code{\link{runTSCANBranchDEG}}. Required.

+#' \code{\link{runTSCANClusterDEAnalysis}}. Required.

 #' @param useReducedDim A single character for the matrix of 2D embedding.

 #' Should exist in \code{reducedDims} slot. Default \code{"UMAP"}.

 #' @param topN Integer. Use top N genes identified. Default \code{9}.

 #' @param useAssay A single character for the feature expression matrix. Should

 #' exist in \code{assayNames(inSCE)}. Default \code{NULL} for using the one used

-#' in \code{\link{runTSCANBranchDEG}}.

+#' in \code{\link{runTSCANClusterDEAnalysis}}.

 #' @param featureDisplay Specify the feature ID type to display. Users can set

 #' default value with \code{\link{setSCTKDisplayRow}}. \code{NULL} or

 #' \code{"rownames"} specifies the rownames of \code{inSCE}. Other character

@@ -703,8 +708,8 @@ plotTSCANBranchPseudo <- function(inSCE, useCluster, useReducedDim = "UMAP",

 #' will combine plots of each feature into a single \code{.ggplot} object,

 #' while \code{"none"} will output a list of plots. Default \code{"all"}.

 #' @return A \code{.ggplot} object of cell scatter plot, colored by the

-#' expression of a gene identified by \code{\link{runTSCANBranchDEG}}, with the

-#' layer of trajectory.

+#' expression of a gene identified by \code{\link{runTSCANClusterDEAnalysis}},

+#' with the layer of trajectory.

 #' @export

 #' @author Yichen Wang

 #' @importFrom S4Vectors metadata

@@ -712,13 +717,14 @@ plotTSCANBranchPseudo <- function(inSCE, useCluster, useReducedDim = "UMAP",

 #' data("mouseBrainSubsetSCE", package = "singleCellTK")

 #' mouseBrainSubsetSCE <- runTSCAN(inSCE = mouseBrainSubsetSCE,

 #'                                 useReducedDim = "PCA_logcounts")

-#' mouseBrainSubsetSCE <- runTSCANBranchDEG(inSCE = mouseBrainSubsetSCE,

-#'                                          useCluster = 1)

-#' plotTSCANBranchDEG(mouseBrainSubsetSCE, useCluster = 1,

-#'                    useReducedDim = "TSNE_logcounts")

-plotTSCANBranchDEG <- function(

+#' mouseBrainSubsetSCE <- runTSCANClusterDEAnalysis(inSCE = mouseBrainSubsetSCE,

+#'                                                  useCluster = 1)

+#' plotTSCANClusterDEG(mouseBrainSubsetSCE, useCluster = 1,

+#'                     useReducedDim = "TSNE_logcounts")

+plotTSCANClusterDEG <- function(

   inSCE,

   useCluster,

+  pathIndex = NULL,

   useReducedDim = "UMAP",

   topN = 9,

   useAssay = NULL,

@@ -731,12 +737,17 @@ plotTSCANBranchDEG <- function(

                                pathName = useCluster)

   if (is.null(DEResults)) {

     stop("Branch DE result of specified cluster not found. Run ",

-         "`runTSCANBranchDEG()` with the same `useCluster` first.")

+         "`runTSCANClusterDEAnalysis()` with the same `useCluster` first.")

   }

   combinePlot <- match.arg(combinePlot)

-  deg <- do.call(rbind, DEResults$DEgenes)

-  deg <- deg[order(deg$FDR),]

-  features <- rownames(deg)[seq(min(topN, nrow(deg)))]

+  if (is.null(pathIndex)) {

+    deg <- do.call(rbind, DEResults$DEgenes)

+    deg <- deg[order(deg$FDR),]

+  } else {

+    pathIndex <- as.character(pathIndex)

+    deg <- DEResults$DEgenes[[pathIndex]]

+  }

+  features <- rownames(deg)[seq(min(topN, nrow(deg), na.rm = TRUE))]

   if (is.null(useAssay)) useAssay <- DEResults$useAssay

   plotTSCANDimReduceFeatures(inSCE = inSCE, features = features,

                              useReducedDim = useReducedDim,

@@ -794,7 +805,8 @@ plotTSCANDimReduceFeatures <- function(

   if (!is.null(useCluster) && length(useCluster) != 1)

     stop("`useCluster`: can only use one cluster.")

   combinePlot <- match.arg(combinePlot)

-  by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = results$clusters)

+  clusters <- colData(inSCE)$TSCAN_clusters

+  by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = clusters)

   line.data <- TSCAN::reportEdges(by.cluster, mst = results$mst,

                                   clusters = NULL, use.dimred = useReducedDim)

   if (!is.null(useCluster)) {

@@ -802,44 +814,46 @@ plotTSCANDimReduceFeatures <- function(

                                  line.data$edge) |

                              grepl(paste0("--", useCluster, "$"),

                                    line.data$edge),]

-    inSCE <- inSCE[, colData(inSCE)$TSCAN_clusters %in% useCluster]

+    inSCE <- inSCE[, clusters %in% useCluster]

   }

   if (by == "rownames") by <- NULL

-  if (!is.null(featureDisplay) &&

-      featureDisplay == "rownames")

-    featureDisplay <- NULL

   plotList <- list()

-  for (f in features) {

-    g <- plotSCEDimReduceFeatures(inSCE, feature = f,

-                                  featureLocation = by,

-                                  featureDisplay = featureDisplay,

-                                  reducedDimName = useReducedDim,

-                                  title = f) +

-      ggplot2::geom_line(data = line.data,

-                         ggplot2::aes_string(x = colnames(line.data)[2],

-                                             y = colnames(line.data)[3]),

-                         inherit.aes = FALSE)

-    # Text labeling code credit: scater

-    reddim <- as.data.frame(SingleCellExperiment::reducedDim(inSCE,

-                                                             useReducedDim))

-    text_out <- scater::retrieveCellInfo(inSCE, "TSCAN_clusters",

-                                         search = "colData")

-    by_text_x <- vapply(split(reddim[,1], text_out$val),

-                        stats::median, FUN.VALUE = 0)

-    by_text_y <- vapply(split(reddim[,2], text_out$val),

-                        stats::median, FUN.VALUE = 0)

-    g <- g + ggrepel::geom_text_repel(

-      data = data.frame(x = by_text_x,

-                        y = by_text_y,

-                        label = names(by_text_x)),

-      mapping = ggplot2::aes_string(x = "x",

-                                    y = "y",

-                                    label = "label"),

-      inherit.aes = FALSE

-    )

-    plotList[[f]] <- g

+  if (!is.na(features)) {

+    for (f in features) {

+      g <- plotSCEDimReduceFeatures(inSCE, feature = f,

+                                    useAssay = useAssay,

+                                    featureLocation = by,dim1 = 1, dim2 = 2,

+                                    featureDisplay = featureDisplay,

+                                    reducedDimName = useReducedDim,

+                                    title = f) +

+        ggplot2::geom_line(data = line.data,

+                           ggplot2::aes_string(x = colnames(line.data)[2],

+                                               y = colnames(line.data)[3],

+                                               group = "edge"),

+                           inherit.aes = FALSE)

+      # Text labeling code credit: scater

+      reddim <- as.data.frame(SingleCellExperiment::reducedDim(inSCE,

+                                                               useReducedDim))

+      text_out <- scater::retrieveCellInfo(inSCE, "TSCAN_clusters",

+                                           search = "colData")

+      by_text_x <- vapply(split(reddim[,1], text_out$val),

+                          stats::median, FUN.VALUE = 0)

+      by_text_y <- vapply(split(reddim[,2], text_out$val),

+                          stats::median, FUN.VALUE = 0)

+      g <- g + ggrepel::geom_text_repel(

+        data = data.frame(x = by_text_x,

+                          y = by_text_y,

+                          label = names(by_text_x)),

+        mapping = ggplot2::aes_string(x = "x",

+                                      y = "y",

+                                      label = "label"),

+        inherit.aes = FALSE

+      )

+      plotList[[f]] <- g

+    }

   }

   if (combinePlot == "all") {

+    if (length(plotList) > 0)

     plotList <- cowplot::plot_grid(plotlist = plotList)

   }

   return(plotList)

@@ -1,61 +1,82 @@

 ###################################################

 ###  Setting accessor functions

 ###################################################

-#' @title getTSCANResults accessor function 

+#' @title getTSCANResults accessor function

+#' @description SCTK allows user to access all TSCAN related results with

+#' \code{"getTSCANResults"}. See details.

 #' @param x Input \linkS4class{SingleCellExperiment} object.

-#' @param analysisName Algorithm name implemented

-#' @param pathName Sub folder name within the \code{analysisName}

-#' @param value Value to be stored within the \code{pathName} or 

+#' @param analysisName Algorithm name implemented, should be one of

+#' \code{"Pseudotime"}, \code{"DEG"}, or \code{"ClusterDEAnalysis"}.

+#' @param pathName Sub folder name within the \code{analysisName}. See details.

+#' @param value Value to be stored within the \code{pathName} or

 #' \code{analysisName}

+#' @details

+#' When \code{analysisName = "Pseudotime"}, returns the list result from

+#' \code{\link{runTSCAN}}, including the MST structure.

+#'

+#' When \code{analysisName = "DEG"}, returns the list result from

+#' \code{\link{runTSCANDEG}}, including \code{DataFrame}s containing genes that

+#' increase/decrease along each the pseudotime paths. \code{pathName} indicates

+#' the path index, the available options of which can be listed by

+#' \code{listTSCANTerminalNodes}.

+#'

+#' When \code{analysisName = "ClusterDEAnalysis"}, returns the list result from

+#' \code{\link{runTSCANBranchDEG}}. Here \code{pathName} needs to match

+#' with the \code{useCluster} argument when running the algorithm.

 #' @export

 #' @rdname getTSCANResults

 #' @return Get or set TSCAN results

+#' @examples

+#' data("mouseBrainSubsetSCE", package = "singleCellTK")

+#' mouseBrainSubsetSCE <- runTSCAN(inSCE = mouseBrainSubsetSCE,

+#'                                 useReducedDim = "PCA_logcounts")

+#' results <- getTSCANResults(mouseBrainSubsetSCE, "Pseudotime")

 setGeneric("getTSCANResults", signature = "x",

-           function(x, analysisName=NULL, pathName = NULL) {

+           function(x, analysisName = NULL, pathName = NULL) {

              standardGeneric("getTSCANResults")

            }

 )

 #' @export

 #' @rdname getTSCANResults

-#' 

-setMethod("getTSCANResults", signature(x = "SingleCellExperiment"), 

-          function(x, analysisName=NULL, pathName = NULL){

+#'

+setMethod("getTSCANResults", signature(x = "SingleCellExperiment"),

+          function(x, analysisName = NULL, pathName = NULL){

   result.names <- listTSCANResults(x)

   if(!analysisName %in% result.names) {

-    stop("The analysis was not found in the results for tool 'Trajectory'")  

-  }

-  if (is.null(pathName)){

-    results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN[[analysisName]]

-    

+    stop("The analysis was not found for TSCAN results")

   }

-  else{

+  if (analysisName == "Pseudotime")

+    results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN$Pseudotime

+  else {

+    if (is.null(pathName)) {

+      stop("Have to specify `pathName` for TSCAN DE analyses")

+    }

+    pathName <- as.character(pathName)

     results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN[[analysisName]][[pathName]]

-    

   }

   return(results)

 })

 #' @export

 #' @rdname getTSCANResults

-#' 

-setGeneric("getTSCANResults<-", 

-           function(x, analysisName, pathName = NULL, value) 

+#'

+setGeneric("getTSCANResults<-",

+           function(x, analysisName, pathName = NULL, value)

              standardGeneric("getTSCANResults<-")

            )

 #' @export

 #' @rdname getTSCANResults

-setReplaceMethod("getTSCANResults", signature(x = "SingleCellExperiment"), 

+setReplaceMethod("getTSCANResults", signature(x = "SingleCellExperiment"),

                  function(x, analysisName, pathName = NULL, value) {

-  

-  if (is.null(pathName)){

-    S4Vectors::metadata(x)$sctk$Traj$TSCAN[[analysisName]] <- value

-    

-  }

-  else{

+  if (analysisName == "Pseudotime")

+    S4Vectors::metadata(x)$sctk$Traj$TSCAN$Pseudotime <- value

+  else {

+    if (is.null(pathName))

+      stop("Have to specify `pathName` for TSCAN DE analyses")

+    pathName <- as.character(pathName)

     S4Vectors::metadata(x)$sctk$Traj$TSCAN[[analysisName]][[pathName]] <- value

-    

   }

   return(x)

 })

@@ -70,810 +91,756 @@ setGeneric("listTSCANResults", signature = "x",

 #' @rdname getTSCANResults

 setMethod("listTSCANResults", "SingleCellExperiment", function(x){

   all.results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN

-  if(is.null(all.results) || length(names(all.results)) == 0) {

-    stop("No results from 'Trajectory' are found. Please run `runTSCAN` first.") 

+  if (is.null(all.results) ||

+      !"Pseudotime" %in% names(all.results)) {

+    stop("No TSCAN result found. Please run `runTSCAN` first.")

   }

   return(names(all.results))

 })

 #' @export

 #' @rdname getTSCANResults

-#' 

 setGeneric("listTSCANTerminalNodes", signature = "x",

-           function(x, analysisName=NULL, pathName = NULL) {

+           function(x) {

              standardGeneric("listTSCANTerminalNodes")

            }

 )

 #' @export

 #' @rdname getTSCANResults

-#' 

-setMethod("listTSCANTerminalNodes", signature(x = "SingleCellExperiment"), 

-          function(x, analysisName=NULL, pathName = NULL){

-            

-  result.names <- listTSCANResults(x)

-  if(!analysisName %in% result.names) {

-    stop("The analysis was not found in the results for tool 'Trajectory'")  

-  }

-  

-  results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN[[analysisName]][["pathClusters"]]

+setMethod("listTSCANTerminalNodes", signature(x = "SingleCellExperiment"),

+          function(x){

+  if (is.null(S4Vectors::metadata(x)$sctk$Traj$TSCAN$Pseudotime))

+    stop("No TSCAN result found. Please run `runTSCAN()` first.")

+  results <- S4Vectors::metadata(x)$sctk$Traj$TSCAN$Pseudotime$pathClusters

   return(names(results))

-          

   })

 ###################################################

 ###  STEP 1:: creating cluster and MST

 ###################################################

-#' @title Run runTSCAN function to obtain pseudotime values for cells

-#' @description Wrapper for obtaining a pseudotime ordering of the cells by 

-#' projecting them onto the MST

-#'

+#' @title Run TSCAN to obtain pseudotime values for cells

+#' @description Wrapper for obtaining a pseudotime ordering of the cells by

+#' projecting them onto the minimum spanning tree (MST)

 #' @param inSCE Input \linkS4class{SingleCellExperiment} object.

-#' @param cluster Grouping for each cell in \code{inSCE}. A user may input a 

-#' vector equal length to the number of the samples in \code{inSCE}, or can be 

-#' retrieved from the \code{colData} slot. Default \code{NULL}.

-#' @param seed An integer. Set the seed for random process that happens only in 

-#' "random" generation. Default \code{12345}.

-#' @param useReducedDim Character. Saved dimension reduction name in 

-#' \code{inSCE} object. Required. Used for specifying which low-dimension 

-#' representation to perform the clustering algorithm and building nearest 

-#' neighbor graph on. Default \code{"PCA"}

-#'

-#' @return A \linkS4class{SingleCellExperiment} object with pseudotime ordering 

-#' of the cells along the paths

-#' @export 

+#' @param useReducedDim Character. A low-dimension representation in

+#' \code{reducedDims}, will be used for both clustering if \code{cluster} not

+#' specified and MST construction. Default \code{"PCA"}.

+#' @param cluster Grouping for each cell in \code{inSCE}. A vector with equal

+#' length to the number of the cells in \code{inSCE}, or a single character for

+#' retriving \code{colData} variable. Default \code{NULL}, will run

+#' \code{runScranSNN} to obtain.

+#' @param seed An integer. Random seed for clustering if \code{cluster} is not

+#' specified. Default \code{12345}.

+#' @return The input \code{inSCE} object with pseudotime ordering of the cells

+#' along the paths and the cluster label stored in \code{colData}, and other

+#' unstructured information in \code{metadata}.

+#' @export

 #' @author Nida Pervaiz

-#' @importFrom utils head

 #' @examples

-#' data("scExample", package = "singleCellTK")

-#' sce <- subsetSCECols(sce, colData = "type != 'EmptyDroplet'")

-#' rowData(sce)$Symbol <- rowData(sce)$feature_name

-#' rownames(sce) <- rowData(sce)$Symbol

-#' sce <- runNormalization(sce, 

-#'                         normalizationMethod = "LogNormalize", 

-#'                         useAssay = "counts", 

-#'                         outAssayName = "logcounts")

-#' sce <- runDimReduce(inSCE = sce, 

-#'                     method = "scaterPCA", 

-#'                     useAssay = "logcounts", 

-#'                     reducedDimName = "PCA")

-#' sce <- runDimReduce(inSCE = sce, 

-#'                     method = "rTSNE", 

-#'                     useReducedDim = "PCA", 

-#'                     reducedDimName = "TSNE")

-#' sce <- runTSCAN (inSCE = sce, 

-#'                  useReducedDim = "PCA", 

-#'                  seed = NULL)

-runTSCAN <- function(inSCE, 

-                     useReducedDim, 

-                     cluster = NULL, 

-                     seed = 12345) { 

-  

-  if (is.null(seed)) {

-    if(is.null(cluster)){

-      inSCE <- runScranSNN(inSCE, useReducedDim = useReducedDim)

-      cluster <- colData(inSCE)$"scranSNN_cluster"

+#' data("mouseBrainSubsetSCE", package = "singleCellTK")

+#' mouseBrainSubsetSCE <- runTSCAN(inSCE = mouseBrainSubsetSCE,

+#'                                 useReducedDim = "PCA_logcounts")

+runTSCAN <- function(inSCE,

+                     useReducedDim = "PCA",

+                     cluster = NULL,

+                     seed = 12345) {

+  if (is.null(cluster)) {

+    # DON'T RETURN TO `inSCE`

+    # It really overwrites existing cluster labels

+    sce <- runScranSNN(inSCE, useReducedDim = useReducedDim, seed = seed)

+    cluster <- colData(sce)$scranSNN_cluster

+    rm(sce)

+  } else if (length(cluster) == 1 & is.character(cluster)) {

+    if (!cluster %in% names(colData(inSCE))) {

+      stop("Specified cluster not found in `colData(inSCE)`.")

     }

+    cluster <- colData(inSCE)[[cluster]]

+  } else if (length(cluster) != ncol(inSCE)) {

+    stop("`cluster` should either has length equal to `ncol(inSCE)` or be a ",

+         "colData variable.")

   }

-  else{

-    withr::with_seed(   

-      seed, 

-      if(is.null(cluster)){

-        inSCE <- runScranSNN(inSCE, useReducedDim = useReducedDim)

-        cluster <- colData(inSCE)$"scranSNN_cluster"

-      })

-  }

+

+  message(date(), " ... Running TSCAN to estimate pseudotime")

   inSCE <- scran::computeSumFactors(inSCE, clusters = cluster)

-  by.cluster <- scuttle::aggregateAcrossCells(inSCE , ids = cluster)

-  centroids <- SingleCellExperiment::reducedDim(by.cluster, useReducedDim)         

-  mst <- TSCAN::createClusterMST(centroids, clusters = NULL)      

-  

-  # Map each cell to the closest edge on the MST, reporting also the distance to 

+  by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = cluster)

+  centroids <- SingleCellExperiment::reducedDim(by.cluster, useReducedDim)

+  mst <- TSCAN::createClusterMST(centroids, clusters = NULL)

+

+  # Map each cell to the closest edge on the MST, reporting also the distance to

   # the corresponding vertices.

-  map.tscan <- TSCAN::mapCellsToEdges(inSCE , mst = mst, 

-                                      use.dimred = useReducedDim , 

+  map.tscan <- TSCAN::mapCellsToEdges(inSCE , mst = mst,

+                                      use.dimred = useReducedDim,

                                       clusters = cluster)

-  

-  # Compute a pseudotime for each cell lying on each path through the MST from a 

+

+  # Compute a pseudotime for each cell lying on each path through the MST from a

   # given starting node.

   tscan.pseudo <- TSCAN::orderCells(map.tscan, mst)

-  common.pseudo <- TrajectoryUtils::averagePseudotime(tscan.pseudo) 

-  

-  colData(inSCE)$"TSCAN_clusters" <- cluster

-  colData(inSCE)$"TSCAN_pseudotime" <- common.pseudo 

-  

+  common.pseudo <- TrajectoryUtils::averagePseudotime(tscan.pseudo)

+

+  colData(inSCE)$TSCAN_clusters <- cluster

+  colData(inSCE)$TSCAN_pseudotime <- common.pseudo

+

   pathClusters <- list()

   pathList <- data.frame()

   for(i in seq(ncol(tscan.pseudo))){

     pathIndex = as.character(colnames(tscan.pseudo)[i])

-    inSCE[[paste0("Path_",pathIndex,"_pseudotime") ]] <- TrajectoryUtils::pathStat(tscan.pseudo)[,pathIndex]

-    nonNA <- inSCE[,!is.na(inSCE[[paste0("Path_",pathIndex,"_pseudotime") ]])]

-    pathClusters[[pathIndex]] <- (sort(unique(colData(nonNA)$TSCAN_clusters)))

-    

+    colDataVarName <- paste0("Path_",pathIndex,"_pseudotime")

+    inSCE[[colDataVarName]] <- TrajectoryUtils::pathStat(tscan.pseudo)[,pathIndex]

+    nonNA.idx <- !is.na(colData(inSCE)[[colDataVarName]])

+    pathClusters[[pathIndex]] <- unique(colData(inSCE)$TSCAN_clusters[nonNA.idx])

     pathList[i,1] <- pathIndex

-    pathList[i,2] <- toString (pathClusters[[pathIndex]])

-    message("Cluster involved in path index ",pathIndex," are: ", pathClusters[i])

+    pathList[i,2] <- toString(pathClusters[[pathIndex]])

+    message(date(), " ...   Clusters involved in path index ", pathIndex,

+            " are: ", paste(pathClusters[[i]], collapse = ", "))

   }

-  

+

   maxWidth <- max(stringr::str_length(pathList[, 1]))

-  choiceList <- pathList[, 1]

-  names(choiceList) <- paste0(stringr::str_pad(pathList[, 1], width=maxWidth, 

-                                               side="right"), 

-                              "|", pathList[, 2])

-  

-  branchClustersList <- c()

+  choiceList <- paste0(stringr::str_pad(pathList[, 1], width=maxWidth,

+                                        side="right"),

+                       "|", pathList[, 2])

+

+  branchClusters <- c()

   edgeList <- mst

   for (i in seq_along(unique(colData(inSCE)$TSCAN_clusters))){

-    clustEdges <- sum(edgeList[i] != 0)

-    if(clustEdges > 1){

-      branchClustersList <- sort(BiocGenerics::append(i,branchClustersList))

-    }

-    

+    degree <- sum(edgeList[i] != 0)

+    if (degree > 2) branchClusters <- c(branchClusters, i)

   }

-  

-  ## Save these results in a list and then make S4 accessor that passes the 

+

+  ## Save these results in a list and then make S4 accessor that passes the

   ## entire list

-  result <- list(pseudo = tscan.pseudo, mst = mst, clusters = cluster, 

-                 by.cluster = by.cluster, maptscan = map.tscan, 

-                 pathClusters = pathClusters, 

-                 branchClustersList = branchClustersList, 

-                 pathIndexList = names(choiceList))

-  getTSCANResults(inSCE, analysisName = "Pseudotime") <- result 

-  

-  message("Number of estimated paths is ", ncol(tscan.pseudo) )

-  

+  result <- list(pseudo = tscan.pseudo, mst = mst, clusters = cluster,

+                 maptscan = map.tscan,

+                 pathClusters = pathClusters,

+                 branchClusters = branchClusters,

+                 pathIndexList = choiceList)

+  getTSCANResults(inSCE, analysisName = "Pseudotime") <- result

+

+  message(date(), " ...   Number of estimated paths is ", ncol(tscan.pseudo))

+

   return (inSCE)

-}  

+}

 ###################################################

 ###  plot pseudotime values

 ###################################################

-#' @title Plot MST pseudotime values for cells

-#' @description A wrapper function which visualizes outputs from the 

-#' \code{\link{runTSCAN}} function. Plots the pseudotime ordering of the cells 

-#' by projecting them onto the MST

-#'

+#' @title Plot MST pseudotime values on cell 2D embedding

+#' @description A wrapper function which visualizes outputs from the

+#' \code{\link{runTSCAN}} function. Plots the pseudotime ordering of the cells

+#' and project them onto the MST.

 #' @param inSCE Input \linkS4class{SingleCellExperiment} object.

-#' @param useReducedDim Saved dimension reduction name in \code{inSCE} object. 

+#' @param useReducedDim Saved dimension reduction name in \code{inSCE} object.

 #' Required.

-#'

-#' @return A plot with the pseudotime ordering of the cells by projecting them 

-#' onto the MST. 

-#' @export 

+#' @return A \code{.ggplot} object with the pseudotime ordering of the cells

+#' colored on a cell 2D embedding, and the MST path drawn on it.

+#' @export

 #' @author Nida Pervaiz

-#' @importFrom utils head

 #' @examples

-#' data("scExample", package = "singleCellTK")

-#' sce <- subsetSCECols(sce, colData = "type != 'EmptyDroplet'")

-#' rowData(sce)$Symbol <- rowData(sce)$feature_name

-#' rownames(sce) <- rowData(sce)$Symbol

-#' sce <- runNormalization(sce, 

-#'                         normalizationMethod = "LogNormalize", 

-#'                         useAssay = "counts", 

-#'                         outAssayName = "logcounts")

-#' sce <- runDimReduce(inSCE = sce, 

-#'                     method = "scaterPCA", 

-#'                     useAssay = "logcounts", 

-#'                     reducedDimName = "PCA")

-#' sce <- runDimReduce(inSCE = sce, 

-#'                     method = "rTSNE", 

-#'                     useReducedDim = "PCA", 

-#'                     reducedDimName = "TSNE")

-#' sce <- runTSCAN (inSCE = sce, 

-#'                  useReducedDim = "PCA", 

-#'                  seed = NULL)

-#' plotTSCANResults(inSCE = sce, 

-#'                  useReducedDim = "TSNE")

-

-plotTSCANResults <- function(inSCE, 

-                             useReducedDim){

-  

+#' data("mouseBrainSubsetSCE", package = "singleCellTK")

+#' mouseBrainSubsetSCE <- runTSCAN(inSCE = mouseBrainSubsetSCE,

+#'                                 useReducedDim = "PCA_logcounts")

+#' plotTSCANResults(inSCE = mouseBrainSubsetSCE,

+#'                  useReducedDim = "TSNE_logcounts")

+plotTSCANResults <- function(inSCE, useReducedDim = "UMAP") {

+

   results <- getTSCANResults(inSCE, analysisName = "Pseudotime")

-  line.data <- TSCAN::reportEdges(results$by.cluster, mst = results$mst, 

-                                  clusters = NULL, use.dimred = useReducedDim) 

-  

-  a <- b <- c <- NULL

-  a = unlist(line.data[,2])

-  b = unlist(line.data[,3])

-  c = unlist(line.data[,1])

-  

-  scater::plotReducedDim(inSCE, dimred = useReducedDim, 

-                         colour_by = I(colData(inSCE)$"TSCAN_pseudotime"), 

-                         text_by = "TSCAN_clusters", text_colour = "black") + 

-    ggplot2::geom_path(data = line.data, ggplot2::aes(x = a, y = b, group = c)) 

-}

+  by.cluster <- scuttle::aggregateAcrossCells(inSCE, ids = results$clusters)

+  line.data <- TSCAN::reportEdges(by.cluster, mst = results$mst,

+                                  clusters = NULL, use.dimred = useReducedDim)

+  scater::plotReducedDim(inSCE, dimred = useReducedDim,