Add function for spatial autocorrelation for SpatialExperiment object

NAMESPACE

@@ -10,6 +10,7 @@ export(gsvaParam)
 export(igsva)
 export(plageParam)
 export(readGMT)
+export(spatCor)
 export(ssgseaParam)
 export(zscoreParam)
 exportClasses(GsvaExprData)
@@ -24,6 +25,7 @@ exportMethods(filterGeneSets)
 exportMethods(geneSetSizes)
 exportMethods(geneSets)
 exportMethods(gsva)
+exportMethods(spatCor)
 import(methods)
 importClassesFrom(Biobase,ExpressionSet)
 importClassesFrom(DelayedArray,DelayedArray)
@@ -62,13 +64,16 @@ importFrom(S4Vectors,SimpleList)
 importFrom(SingleCellExperiment,SingleCellExperiment)
 importFrom(SpatialExperiment,SpatialExperiment)
 importFrom(SummarizedExperiment,SummarizedExperiment)
+importFrom(SummarizedExperiment,assay)
 importFrom(graphics,plot)
 importFrom(methods,new)
 importFrom(parallel,splitIndices)
 importFrom(sparseMatrixStats,colRanks)
 importFrom(sparseMatrixStats,rowRanges)
+importFrom(stats,dist)
 importFrom(stats,ecdf)
 importFrom(stats,na.omit)
+importFrom(stats,pnorm)
 importFrom(stats,rnorm)
 importFrom(stats,rpois)
 importFrom(stats,sd)

R/AllGenerics.R

@@ -49,3 +49,12 @@ setGeneric("gsvaAnnotation", function(object) standardGeneric("gsvaAnnotation"))
 
 setGeneric("gsvaAssayNames", function(object) standardGeneric("gsvaAssayNames"))
 
+## spatial methods
+
+#' @export
+setGeneric("spatCor",
+           function(spe, ...) standardGeneric("spatCor"))
+
+
+
+

R/spatial.R

@@ -0,0 +1,105 @@
+#' @importFrom stats dist pnorm
+#' @importFrom SummarizedExperiment assay
+
+### ----- Method for computing Spatial Autocorrelation for SpatialExperiment object -----
+
+#' @title Compute Spatial Autocorrelation for SpatialExperiment objects
+#' 
+#' @description Computes spatial autocorrelation using Moran's I statistic
+#' for a \code{SpatialExperiment} object, using an inverse squared distance weight matrix as default,
+#' or an inverse distance weight matrix as an alternative. It also tests for spatial autocorrelation
+#' assuming normality.
+#' 
+#' @param spe An object of \code{SpatialExperiment} class.
+#' @param alternative A character string specifying the alternative hypothesis tested against the null hypothesis of no spatial autocorrelation;
+#'  must be one of "two.sided", "less", or "greater", or any unambiguous abbreviation of these.
+#' @param na.rm A logical indicating whether missing values should be removed.
+#' @param squared A logical indicating whether the inverse distance weight matrix should be squared or not.
+#' 
+#' @return A \code{data.frame} with the same row names as the original \code{SpatialExperiment} object.
+#' Columns include the observed Moran's I statistic, the expected Moran's I statistic under no spatial autocorrelation, the expected
+#' standard deviation under no spatial autocorrelation, and the p-value of the test.
+#' 
+#' @aliases spatCor spatCor,SpatialExperiment-method
+#' @name spatCor
+#' @rdname spatCor
+#' @exportMethod spatCor
+#' @export
+
+setMethod("spatCor", signature("SpatialExperiment"),
+          function(spe, na.rm = FALSE, alternative = "two.sided", squared = TRUE) {
+            weight_list <- .spe_dist_weight_matrix(spe, squared)
+            logc <- assay(spe)
+
+            spe_Moran <- list()
+            rowns <- rownames(spe)
+            spe_Moran <- lapply(rowns, function(x){
+              .internal_moran(logc[rowns == x, ], weight_list, na.rm = na.rm, alternative = alternative)
+            })
+            names(spe_Moran) <- rownames(spe)
+            df_res <- do.call(rbind, lapply(spe_Moran, as.data.frame))
+            return(df_res)
+          })
+
+.internal_moran <- function (x, weight_list, na.rm = FALSE, alternative = "two.sided") 
+{
+  weight <- weight_list$weight
+  n <- length(x)
+  ei <- -1/(n - 1)
+  nas <- is.na(x)
+  if (any(nas)) {
+    if (na.rm) {
+      x <- x[!nas]
+      n <- length(x)
+      weight <- weight[!nas, !nas]
+    }
+    else {
+      warning("'x' has missing values: maybe you wanted to set na.rm = TRUE?")
+      return(list(observed = NA, expected = ei, sd = NA, 
+                  p.value = NA))
+    }
+  }
+  m <- mean(x, na.rm = na.rm)
+  y <- x - m
+  cv <- sum(weight * y %o% y, na.rm = na.rm)
+  v <- sum(y^2, na.rm = na.rm)
+  s <- weight_list$s
+  obs <- (n/s) * (cv/v)
+  k <- (sum(y^4)/n)/(v/n)^2
+  s.sq <- weight_list$s.sq
+  S1 <- weight_list$S1
+  S2 <- weight_list$S2
+  sdi <- sqrt((n * ((n^2 - 3 * n + 3) * S1 - n * S2 + 3 * s.sq) - 
+                 k * (n * (n - 1) * S1 - 2 * n * S2 + 6 * s.sq))/((n - 
+                                                                     1) * (n - 2) * (n - 3) * s.sq) - 1/((n - 1)^2))
+  alternative <- match.arg(alternative, c("two.sided", "less", 
+                                          "greater"))
+  pv <- pnorm(obs, mean = ei, sd = sdi)
+  if (alternative == "two.sided") 
+    pv <- if (obs <= ei) 
+      2 * pv
+  else 2 * (1 - pv)
+  if (alternative == "greater") 
+    pv <- 1 - pv
+  list(observed = obs, expected = ei, sd = sdi, p.value = pv)
+}
+
+
+
+.spe_dist_weight_matrix<-function(spe, squared = TRUE){
+  xy <- spatialCoords(spe)
+  weight<-as.matrix(dist(xy))
+  if(squared == TRUE)
+    weight=1/weight^2
+  else weight = 1/weight
+  diag(weight)<-0
+  s <- sum(weight)
+  ROWSUM <- rowSums(weight)
+  ROWSUM[ROWSUM == 0] <- 1
+  return(list(weight=weight,
+              ROWSUM = ROWSUM,
+              s = s,
+              S1 = 0.5 * sum((2*weight)^2),
+              S2 = sum((ROWSUM + colSums(weight))^2),
+              s.sq = s^2))
+}