R/pseudotime.R

In sbrn3/disscat: Calculates Distances Between Single Cell RNA Seq Categories

###############################################
# Pseudotime
###############################################

#' Returns a slingshot based on the plsda of a seurat object
#'
#' @param seurat_object seurat_object containing counts and metadata
#' @param reduction Dimensionality reduction data to use
#' @param dims Dimensions to use
#' @param group.by Categories that it will draw trajectories through. If there
#' is no preference then it is possible to select "recluster" where it will
#' cluster cells and generate groups using Mclust.
#' @param output Output object. Can choose "sce"(Single Cell Experiment S4),
#' "slingshot" (Slingshot S4) or "seurat"(Seurat S4)
#'
#' @import mclust
#'
#' @return a sce which contains the slingshot data
#' @export
RunSlingshot <- function(seurat_object,
                         reduction,
                         dims = 1:2,
                         group.by = "recluster",
                         output = "seurat",
                         thresh = 0.001) {
    # Extract counts from seurat object. Only works for integrated data.
    count.data <-as.matrix(seurat_object@assays[[seurat_object@active.assay]]@data)

    # Create a single cell experiment object from count data
    sce.object <-
      SingleCellExperiment::SingleCellExperiment(assays = base::list(counts = count.data))

    # Add reduced dimensions to single cell experiment
    rd1 <- as.matrix(seurat_object@reductions[[reduction]]@cell.embeddings[,dims])

    SingleCellExperiment::reducedDims(sce.object) <- S4Vectors::SimpleList(reduced_data = rd1)

    # Cluster cells (Seurat requires clustering)
    if (group.by == "recluster") {
      SummarizedExperiment::colData(sce.object)[["label"]] <-
        mclust::Mclust(rd1)$classification
    } else {
      SummarizedExperiment::colData(sce.object)[["label"]] <-
        as.factor(as.matrix(seurat_object[[group.by]]))
    }
    # Run slingshot
    sce.slingshot <- slingshot::slingshot(sce.object,
                               clusterLabels = "label",
                               reducedDim = 'reduced_data',
                               thresh = thresh)

    # Output
    if (output == "seurat") {
        # Add cell data
        seurat_object <- AddMetaDataMatrix(seurat_object,
                                           t(as.matrix(SummarizedExperiment::colData(sce.slingshot))))

        location <- paste("slingshot", reduction, sep = "_")
        seurat_object@tools[[location]] <-
          slingshot::SlingshotDataSet(sce.slingshot)
        return(seurat_object)
    } else if (output == "sce") {
        return(sce.slingshot)
    } else if (output == "slingshot") {
        return(slingshot::SlingshotDataSet(sce.slingshot))
    } else {
        stop("Select a valid output type: 'seurat', 'sce', 'slingshot'")
    }
}

#' Return slingshot object stored in a seurat object
#'
#' @param seurat_object A seurat object
#' @param reduction Dimensionality reduction method used for slingshot
#'
#' @return A SingleCellExperiment with Slingshot data
GetSlingshot <- function(seurat_object, reduction) {
    location <- paste("slingshot", reduction, sep = "_")
    return(seurat_object@tools[[location]])
}

#' Plots slingshot curves stored in a seurat object.
#'
#' @param seurat_object A seurat object
#' @param reduction Dimensionality reduction method used for slingshot
#' @param group.by Groups used for slingshot analysis
#' @param lineages Draw curve linear connections between discreet groups
#'
#' @return Scatter with slingshot lineage(s). Coloured by groups.
#' @export
#'
PlotSlingshot <- function(seurat_object,
                          reduction,
                          group.by = NULL,
                          lineages = FALSE) {
    numeric.factor <- as.numeric(as.factor(as.matrix(seurat_object[[group.by]])))
    if (length(unique(as.matrix(seurat_object[[group.by]]))) >= 10) {
        colors <- grDevices::colorRampPalette(RColorBrewer::brewer.pal(11,'Spectral')[-6])(100)
        plotcol <- colors[cut(numeric.factor, breaks=100)]
    } else {
        different <- length(unique(as.matrix(seurat_object[[group.by]])))

        # Setup colours for plotting
        plotcol <- RColorBrewer::brewer.pal(different,'Set1')[numeric.factor]

    }

    plot.data <- as.data.frame(seurat_object@reductions[[reduction]]@cell.embeddings[,1:2])
    if (!lineages) {
        p <- slingshot::plot(plot.data, col = plotcol, pch=16, asp = 0.8) +
          slingshot::lines(GetSlingshot(seurat_object, reduction), lwd=2, col='black')
    } else {
        p <- slingshot::plot(plot.data, col = plotcol, pch=16, asp = 0.8) +
          slingshot::lines(GetSlingshot(seurat_object, reduction),
                           lwd=2,
                           type = 'lineages',
                           col='black')
    }

    return(p)
}

#' Create density plot from seurat object
#'
#' @description
#' Creates a density plot of the a continuous variable of a seurat object
#' and is able to split it by a category. There is also an option to show the
#' means of each group.
#'
#' @param seurat_object A seurat object
#' @param x X axis category
#' @param group.by Category to split the density plot by
#' @param title Title of plot
#' @param subtitle Subtitle of plot
#' @param show.mean Whether to show the category means
#'
#' @export
#' @return ggplot density plot
DenPlot <- function(seurat_object,
                    x,
                    group.by,
                    title = NULL,
                    subtitle = NULL,
                    show.mean = FALSE) {
    plot.data <- seurat_object@meta.data[c(x, group.by)]
    colnames(plot.data) <- c("dimension", "group.by")
    plot.data$group.by <- as.factor(plot.data$group.by)
    p <- ggplot2::ggplot(data = plot.data, ggplot2::aes(x = plot.data$dimension,
                                                 fill = plot.data$group.by,
                                                 color = group.by)) +
      ggplot2::geom_density(alpha = 0.2) +
      ggplot2::scale_fill_discrete(name = group.by) +
      ggplot2::scale_color_discrete(name = group.by) +
      ggplot2::labs(x = x, title = title, subtitle = subtitle)
    if (show.mean) {
        mu <- plyr::ddply(plot.data, "group.by", plyr::summarise,
                    grp.mean=mean(stats::na.omit(dimension)))
        p <- p + ggplot2::geom_vline(data=mu,
                                     ggplot2::aes(xintercept=mu$grp.mean,
                                                  color=group.by),
                                     linetype="dashed")
    }
    return(p)
}

#' Transfer Slingshot pseudotime dimension to other seurat objects
#'
#' @description
#' If you would like to cast a naive dataset into the same pseudotime as some
#' other dataset. This does require that you have performed the same
#' dimensionality reduction on both datasets.
#'
#' @param reference Seurat object
#' @param query Seurat object
#' @param reduction Reduction method. Default is "plsda"
#' @param dims Dimensions used in slingshot pseudotime calculation.
#'
#' @return Query dataset but with pseudotime values
#' @export
#'
TransferSlingshot <- function(reference,
                              query,
                              reduction = "plsda",
                              dims = 1:10) {
  # Extract data
  reference.slingshot <- reference@tools$slingshot_plsda
  query.embeddings <- Seurat::Embeddings(query, reduction = "plsda")

  predicted.slingshot <- slingshot::predict(reference.slingshot,
                                 query.embeddings[,dims])

  for (i in 1:length(predicted.slingshot@curves)) {
    pseudotime_values <- predicted.slingshot@curves[[i]]$lambda
    location <- paste("slingPseudotime", i, sep = "_")
    query[[location]] <- pseudotime_values
  }

  query@tools[["predicted_slingshot"]] <- predicted.slingshot
  return(query)
}


#' Merges Seurat Objects
#'
#' @description
#' Merges seurat objects but also includes dimensionality reductions
#'
#' @param seurat_1 Seurat Object
#' @param seurat_2 Seurat Object
#' @param reduction Reductions to include
#'
#' @return Merged Seurat object
#' @export
MergeSeurat <- function(seurat_1, seurat_2, reduction = c("pca", "umap", "plsda")) {
  seurat_combined <- merge(seurat_1, seurat_2)
  for (i in 1:length(reduction)) {
    # does the reduction exist in both objects?
    if (is.null(seurat_1@reductions[[reduction[[i]]]])) {
      stop(paste("Reduction", reduction[[i]], "does not exist in the first object"))
    }
    if (is.null(seurat_2@reductions[[reduction[[i]]]])) {
      stop(paste("Reduction", reduction[[i]], "does not exist in the second object"))
    }

    # seurat 1 embeddings
    x_1 <- as.matrix(seurat_1@reductions[[reduction[[i]]]]@cell.embeddings)
    # seurat 2 embeddings
    x_2 <- as.matrix(seurat_2@reductions[[reduction[[i]]]]@cell.embeddings)

    # Are the number of columns the same in both objects
    if (ncol(x_1) != ncol(x_2)) {
      stop(paste("The number of dimensions for each experiment in",
                 reduction[[i]], "must be the same"))
    }

    # merge the embeddings (rbind)
    x <- rbind(x_1, x_2)
    # Input into object
    dim_reduce_object <-   Seurat::CreateDimReducObject(embeddings = x,
                                                        key = paste(reduction[[i]], "_", sep= ""),
                                                        assay = Seurat::DefaultAssay(seurat_combined))
    seurat_combined@reductions[[reduction[[i]]]] <- dim_reduce_object
  }
  return(seurat_combined)
}

#' Calculates the pseudotime of a object with multiple time points
#'
#' @description
#' Generally meant to enable you to find the pseudotime of a single cell type
#' within a tissue when you have multiple time points.
#'
#' @param seurat.object Seurat object of interest
#' @param subset.data Name of metadata dataset to subset the data by
#' @param subset.value Name of metadata value to subset the data by
#' @param group.data Name of metadata data that desribes time points
#' @param group.begin Starting time point
#' @param group.end Ending time point
#' @param comp number of plsda components to use
#' @param near.zero.var see mixomics::plsda
#'
#' @return List containing the dataset with original pslda called 'extremes'
#' the final dataset with pseudotime values calculated.
#' @export
RunPseudotime <- function(seurat.object,
                       subset.data = NULL,
                       subset.value= NULL,
                       group.data,
                       group.begin,
                       group.end,
                       comp = 2,
                       near.zero.var = TRUE,
                       verbose = FALSE,
                       thresh = 0.001) {
  output.list <- vector(mode = "list")

  if (!is.null(subset.data)) {
    if (verbose == TRUE) {
      message(paste("Subsetting", subset.data, "to", subset.value, sep = ""))
    }
    # Subset to cell type
    cell.type <- as.vector(seurat.object[[subset.data]] == subset.value)
    seurat.object <- seurat.object[,cell.type]
    output.list[["name"]] <- subset.data
  }

  # Subset to extreme groups
  if (verbose == TRUE) {
    message(paste("Subsetting", group.data, "to", group.begin, "and", group.end, sep = " "))
  }
  extreme.groups <- as.vector(seurat.object[[group.data]] == group.begin | seurat.object[[group.data]] == group.end)
  seurat.object.extreme <- seurat.object[,extreme.groups]

  # PLSDA on extreme groups
  if (verbose == TRUE) {
    message(paste("Running plsda on", group.begin, "and", group.end, sep = " "))
  }
  seurat.object.extreme <- RunPLSDA(seurat.object.extreme, group = group.data,
                                    comp = comp,
                                    near.zero.var = near.zero.var)

  # Cast all groups into plsda
  if (verbose == TRUE) {
    message("Transfering PLSDA variates to all data")
  }
  seurat.object <- TransferPLSDA(seurat.object.extreme, seurat.object)

  # Run Slingshot on all groups
  if (verbose == TRUE) {
    message("Running slingshot pseudotime ordering")
  }
  seurat.object <- RunSlingshot(seurat.object,
                                reduction = "plsda",
                                dims = 1:comp,
                                group.by = group.data,
                                thresh = thresh)

  output.list[["extreme"]]   <- seurat.object.extreme
  output.list[["pseudotime"]] <- seurat.object
  return(output.list)
}

#' Plot some of the pseudotime results
#'
#' @param pseudotime.list Results from RunPseudotime
#' @param group.by Group to colour by
#' @param reductions Dimensionality reductions to plot
#'
#' @return List of plots
#' @export
VisualisePseudotime <- function(pseudotime.list,
                                group.by,
                                reductions = c("pca", "plsda")) {
  seurat.object <- pseudotime.list[["pseudotime"]]
  seurat.object.extreme <- pseudotime.list[["extreme"]]
  output.list <- vector(mode = "list")
  for (i in 1:length(reductions)) {
    name <- paste(reductions[i], group.by, sep = ".")


    output.list[[paste(name, "all", sep = ".")]] <- Seurat::DimPlot(seurat.object, reduction = reductions[i], group.by = group.by)
    output.list[[paste(name, "extreme")]] <- Seurat::DimPlot(seurat.object, reduction = reductions[i], group.by = group.by, label = TRUE)
  }
}