Merge pull request #51 from YosefLab/staging

Merge Staging to Master

Author: deto

.travis.yml

@@ -5,7 +5,9 @@ cache: packages
 
 r:
     - bioc-release
-    - bioc-devel
+    - devel
+    - oldrel
+    - 3.4
 
 sudo: false
 

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: VISION
 Title: Functional interpretation of single cell RNA-seq latent manifolds
-Version: 1.0.1
+Version: 1.1.0
 Authors@R: c(person("Matt", "Jones", email = "mattjones315@gmail.com", role = c("aut", "cre")),
 			 person("David", "Detomaso", email = "david.detomaso@berkeley.edu", role = c("aut", "cre")),
 			 person("Tal", "Ashuach", email = "tal_ashuach@berkeley.edu", role = c("aut")),
@@ -13,7 +13,6 @@ Imports:
 	fastICA,
 	igraph,
 	jsonlite,
-	jug,
 	loe,
 	logging,
     parallel,
@@ -22,6 +21,7 @@ Imports:
 	matrixStats,
 	mclust,
 	methods,
+    plumber,
 	Rcpp,
 	rsvd,
 	stats,

NAMESPACE

@@ -25,7 +25,6 @@ exportMethods(saveAndViewResults)
 exportMethods(viewResults)
 import(Matrix)
 import(Rcpp)
-import(jug)
 import(loe)
 import(logging)
 import(methods)
@@ -65,6 +64,8 @@ importFrom(mclust,Mclust)
 importFrom(mclust,mclustBIC)
 importFrom(parallel,mclapply)
 importFrom(pbmcapply,pbmclapply)
+importFrom(plumber,forward)
+importFrom(plumber,plumber)
 importFrom(rsvd,rsvd)
 importFrom(stats,chisq.test)
 importFrom(stats,kmeans)

NEWS.md

@@ -1,3 +1,11 @@
+# VISION 1.1.0
+
+* Added differential expression testing to the output report
+    * Thanks to @Yanay1 (Yanay Rosen)
+* Added gene-signature importance calculation
+* More flexible projection plotting
+* Added UMAP as a projection method
+
 # VISION 1.0.1
 
 * Bugfixes related to caching in the output report

R/AllClasses.R

@@ -98,9 +98,11 @@ Vision <- setClass("Vision",
         projection_methods = "character",
         Projections = "list",
         TrajectoryProjections = "list", # list of TrajectoryProjection
+        SigGeneImportance = "list",
         pools = "list",
         inputProjections = "list",
         name = "character",
+        num_neighbors = "numericORNULL",
         latentSpace = "matrix",
         latentTrajectory = "Trajectory",
         version = "numeric",
@@ -126,9 +128,11 @@ Vision <- setClass("Vision",
         projection_methods = character(),
         Projections = list(),
         TrajectoryProjections = list(),
+        SigGeneImportance = list(),
         pools = list(),
         inputProjections = list(),
         name = "",
+        num_neighbors = NULL, 
         latentSpace = matrix(NA, 1, 1),
         latentTrajectory = NULL,
         version = 1.1,

R/AnalysisFunctions.R

@@ -13,8 +13,9 @@ clusterCells <- function(object) {
     n_workers <- getOption("mc.cores")
     n_workers <- if (is.null(n_workers)) 2 else n_workers
 
+    K <- min(object@num_neighbors, 30)
     kn <- ball_tree_knn(res,
-                        min(round(sqrt(nrow(res))), 30),
+                        K,
                         n_workers)
 
     cl <- louvainCluster(kn, res)
@@ -67,7 +68,8 @@ poolCells <- function(object,
                                               cellsPerPartition = object@cellsPerPartition,
                                               filterInput = object@projection_genes,
                                               filterThreshold = object@threshold,
-                                              latentSpace = object@latentSpace)
+                                              latentSpace = object@latentSpace, 
+                                              K = object@num_neighbors)
 
         object@pools <- pools
     }
@@ -232,6 +234,82 @@ calcSignatureScores <- function(object,
     return(object)
 }
 
+
+#' calculate gene-signature importance
+#'
+#' For each signature, the contribution of each gene to the signature score
+#' is evaluated by calculating the covariance between signature scores and expression
+#' The correlation of genes with a negative sign in the signature are inverted.
+#'
+#' @importFrom pbmcapply pbmclapply
+#' @importFrom matrixStats colSds
+#' @importFrom matrixStats rowSds
+#'
+#' @param object the VISION object
+#' @return the VISION object, with SigGeneImportance slot populated
+evalSigGeneImportance <- function(object){
+
+    message("Evaluating signature-gene importance...\n")
+
+    sigScores <- object@sigScores
+
+    if (length(object@sigData) == 0) {
+        object@SigGeneImportance <- list()
+        return(object)
+    }
+
+    if (length(sigScores) <= 1){
+        stop(
+            sprintf("Signature scores have not yet been computed.  `calcSignatureScores` must be run before running `evalSigGeneImportance`")
+            )
+    }
+
+    normExpr <- getNormalizedCopy(object@exprData, object@sig_norm_method)
+
+    # Center each column of sigScores first
+
+    mu <- colMeans(sigScores)
+
+    sigScores <- t(sigScores)
+    sigScores <- (sigScores - mu)
+    sigScores <- t(sigScores)
+
+    # Center each row of normExpr
+    mu <- rowMeans(normExpr)
+
+    normExpr <- (normExpr - mu)
+
+    # Compute Covariances
+    sigData <- object@sigData
+
+    sigGene <- function(signame) {
+        sigdata <- sigData[[signame]]
+
+        genes <- sigdata@sigDict
+
+        sigvals <- sigScores[, signame]
+
+        geneIndices <- match(names(genes), rownames(normExpr))
+
+        corr <- sigGeneInner(sigvals, normExpr, geneIndices)
+
+        names(corr) <- names(genes)
+
+        corr <- corr * genes
+
+        return(corr)
+    }
+
+    sigs <- colnames(sigScores)
+    res <- pbmclapply(setNames(sigs, sigs), sigGene)
+
+    object@SigGeneImportance <- res
+
+
+    return(object)
+}
+
+
 #' Computes the latent space of the expression matrix using PCA
 #'
 #' @param object the VISION object for which compute the latent space
@@ -302,13 +380,24 @@ generateProjections <- function(object) {
   projections <- generateProjectionsInner(object@exprData,
                                      object@latentSpace,
                                      projection_genes = object@projection_genes,
-                                     projection_methods = object@projection_methods)
+                                     projection_methods = object@projection_methods, 
+                                     K = object@num_neighbors)
 
   # Add inputProjections
   for (proj in names(object@inputProjections)){
       projections[[proj]] <- object@inputProjections[[proj]]
   }
 
+  # Make sure all projections have column names
+  n <- names(projections)
+  projections <- lapply(setNames(n, n), function(pname){
+      proj <- projections[[pname]]
+      if (is.null(colnames(proj))){
+          colnames(proj) <- paste0(pname, "-", seq_len(ncol(proj)))
+      }
+      return(proj)
+  })
+
   object@Projections <- projections
 
   message("")
@@ -337,7 +426,8 @@ analyzeLocalCorrelations <- function(object, signatureBackground = NULL) {
                                 object@latentSpace,
                                 object@sigScores,
                                 object@metaData,
-                                signatureBackground)
+                                signatureBackground,
+                                object@num_neighbors)
 
   message("Clustering signatures...\n")
   sigClusters <- clusterSignatures(object@sigScores,

R/Microclusters.R

@@ -18,6 +18,7 @@
 #' If greater than 1, this specifies the number of cells in which a gene must be detected
 #' for it to be used when computing PCA. If less than 1, this instead specifies the proportion of cells needed
 #' @param latentSpace (Optional) Latent space to be used instead of PCA numeric matrix cells x components
+#' @param K Number of neighbors to use for finding pools. 
 #' @importFrom Matrix tcrossprod
 #' @importFrom Matrix rowMeans
 #' @return pooled cells - named list of vectors - cells in each supercell
@@ -26,7 +27,7 @@ applyMicroClustering <- function(
                          exprData, cellsPerPartition=10,
                          filterInput = "fano",
                          filterThreshold = round(ncol(exprData) * 0.05),
-                         latentSpace = NULL) {
+                         latentSpace = NULL, K=round(sqrt(ncol(exprData)))) {
 
     if (is.data.frame(exprData)){
         exprData <- data.matrix(exprData)
@@ -92,7 +93,7 @@ applyMicroClustering <- function(
     n_workers <- if (is.null(n_workers)) 2 else n_workers
 
     kn <- ball_tree_knn(res,
-                        min(round(sqrt(nrow(res))), 30),
+                        min(K, 30),
                         n_workers)
 
     cl <- louvainCluster(kn, res)

R/Projections.R

@@ -25,17 +25,17 @@ registerMethods <- function(projection_methods, numCells) {
 #' @param latentSpace numeric matrix cells x components
 #' @param projection_genes character vector of gene names to use for projections
 #' @param projection_methods character vector of projection methods to use
+#' @param K Number of neighbors to use in projections.
 #' @return list of Projection objects
-generateProjectionsInner <- function(expr, latentSpace, projection_genes=NULL, projection_methods = NULL) {
+generateProjectionsInner <- function(expr, latentSpace, projection_genes=NULL, projection_methods = NULL, K = round(sqrt(ncol(expr)))) {
 
     NUM_CELLS <- nrow(latentSpace)
     methodList <- registerMethods(projection_methods, NUM_CELLS)
 
     projections <- list()
 
-    projections[["PCA: 1,2"]] <- latentSpace[, c(1, 2)]
-    projections[["PCA: 1,3"]] <- latentSpace[, c(1, 3)]
-    projections[["PCA: 2,3"]] <- latentSpace[, c(2, 3)]
+    projections[["Latent Space"]] <- latentSpace
+
     N <- length(methodList)
     for (i in seq_len(N)){
         method <- names(methodList)[i]
@@ -55,6 +55,9 @@ generateProjectionsInner <- function(expr, latentSpace, projection_genes=NULL, p
 
             res <- methodList[[method]](exprData)
             projections[[method]] <- res
+        } else if (method == "UMAP") {
+            res <- methodList[[method]](t(latentSpace), K)
+            projections[[method]] <- res
         } else { ## run on reduced data
             res <- methodList[[method]](t(latentSpace))
             projections[[method]] <- res
@@ -248,10 +251,11 @@ applySpectralEmbedding <- function(exprData) {
 #' Performs UMAP on data
 #'
 #' @param exprData Expression data, NUM_GENES x NUM_SAMPLES
-#'
+#' @param K Number of neighbors to use in UMAP projection.
 #' @return Reduced data NUM_SAMPLES x NUM_COMPONENTS
-applyUMAP <- function(exprData) {
+applyUMAP <- function(exprData, K) {
 
+	print(requireNamespace('uwot'))
     if (!requireNamespace("uwot", quietly = TRUE)){
         stop("Package \"uwot\" needed to run UMAP.  Please install it using:\n\n   devtools::install_github(\"jlmelville/uwot\")\n\n",
             call. = FALSE)
@@ -261,7 +265,7 @@ applyUMAP <- function(exprData) {
     n_workers <- getOption("mc.cores")
     n_workers <- if (is.null(n_workers)) 2 else n_workers
 	res <- uwot::umap(
-        ndataT, n_neighbors = round(sqrt(ncol(exprData))),
+        ndataT, n_neighbors = K,
         n_threads = n_workers, ret_nn = T
     )
 	res <- res$embedding
@@ -389,7 +393,7 @@ clipBottom <- function(x, mi) {
 #' @importFrom Matrix sparseMatrix
 #' @importFrom matrixStats rowMaxs
 #' @param object matrix to use for KNN
-#' @param K number of neighbors to compute this for
+#' @param K Number of neighbors to consider.
 #' @return a list of two items:
 #'          indices: matrix, cells X neighbors
 #'              Each row specifies indices of nearest neighbors
@@ -398,6 +402,7 @@ clipBottom <- function(x, mi) {
 setMethod("computeKNNWeights", signature(object = "matrix"),
     function(object, K = round(sqrt(nrow(object)))) {
 
+        print(K)
         n_workers <- getOption("mc.cores")
         n_workers <- if (is.null(n_workers)) 2 else n_workers
 
@@ -406,6 +411,7 @@ setMethod("computeKNNWeights", signature(object = "matrix"),
         d <- k[[2]]
 
         sigma <- rowMaxs(d)
+        sigma[sigma == 0] <- 1.0  # Can happen if all neighbors at same point as cell
         sparse_weights <- exp(-1 * (d * d) / sigma ^ 2)
 
         # Normalize row sums = 1

R/RcppExports.R

@@ -41,3 +41,11 @@ multMat <- function(X, Y) {
     .Call('_VISION_multMat', PACKAGE = 'VISION', X, Y)
 }
 
+sigGeneInner <- function(sigScores, exp, geneIndices) {
+    .Call('_VISION_sigGeneInner', PACKAGE = 'VISION', sigScores, exp, geneIndices)
+}
+
+wilcox_subset <- function(vals, indicesA, indicesB) {
+    .Call('_VISION_wilcox_subset', PACKAGE = 'VISION', vals, indicesA, indicesB)
+}
+