#' Performs spatial tests on spatial cytometry data.
#'
#' @param cells A SegmentedCells or data frame that contains at least the
#' variables x and y, giving the location coordinates of each cell, and cellType.
#' @param condition Vector of conditions to be tested corresponding to each image
#' if cells is a data frame.
#' @param subject Vector of subject IDs corresponding to each image if cells is
#' a data frame.
#' @param covariates Vector of covariate names that should be included in the
#' mixed effects model as fixed effects.
#' @param from vector of cell types which you would like to compare to the to vector
#' @param to vector of cell types which you would like to compare to the from vector
#' @param dist The distance at which the statistic is obtained.
#' @param integrate Should the statistic be the integral from 0 to dist, or the
#' value of the L curve at dist.
#' @param nsim Number of simulations to perform. If empty, the p-value from lmerTest is used.
#' @param verbose logical indicating whether to output messages.
#' @param weights logical indicating whether to include weights based on cell counts.
#' @param weightsByPair logical indicating whether weights should be calculated for each cell type pair.
#' @param weightFactor numeric that controls the convexity of the weight function.
#' @param window Should the window around the regions be 'square', 'convex' or 'concave'.
#' @param window.length A tuning parameter for controlling the level of concavity
#' when estimating concave windows.
#' @param BPPARAM A BiocParallelParam object.
#' @param sigma A numeric variable used for scaling when fitting inhomogeneous L-curves.
#' @param minLambda Minimum value for density for scaling when fitting inhomogeneous L-curves.
#' @param Rs A vector of the radii that the measures of association should be calculated.
#' @param fast A logical describing whether to use a fast approximation of the
#' inhomogeneous L-curves.
#' @param edgeCorrect A logical indicating whether to perform edge correction.
#' @param includeZeroCells A logical indicating whether to include cells with
#' zero counts in the pairwise association calculation.
#' @param imageID The image ID if using SingleCellExperiment.
#' @param cellType The cell type if using SingleCellExperiment.
#' @param spatialCoords The spatial coordinates if using a SingleCellExperiment.
#' @param ... Other options to pass to bootstrap.
#'
#' @return Data frame of p-values.
#' @export
#'
#' @examples
#' data("diabetesData")
#'
#' # Test with random effect for patient on only one pairwise combination of cell types.
#' spicy(diabetesData, condition = "stage", subject = "case",
#' from = "Tc", to = "Th")
#'
#' # Test all pairwise combination of cell types without random effect of patient.
#' #spicyTest <- spicy(diabetesData, condition = "stage", subject = "case")
#'
#' # Test all pairwise combination of cell types with random effect of patient.
#' #spicy(diabetesData, condition = "condition", subject = "subject")
#'
#' # Test all pairwise combination of cell types with random effect of patient using
#' # a bootstrap to calculate significance.
#' #spicy(diabetesData, condition = "stage", subject = "case", nsim = 10000)
#'
#' @aliases
#' spicy
#' spicy,spicy-method
#' @importFrom mgcv gam ti
#' @importFrom BiocParallel SerialParam
#' @importFrom scam scam
spicy <- function(cells,
condition = NULL,
subject = NULL,
covariates = NULL,
from = NULL,
to = NULL,
dist = NULL,
integrate = TRUE,
nsim = NULL,
verbose = TRUE,
weights = TRUE,
weightsByPair = FALSE,
weightFactor = 1,
window = "convex",
window.length = NULL,
BPPARAM = BiocParallel::SerialParam(),
sigma = NULL,
Rs = NULL,
minLambda = 0.05,
fast = TRUE,
edgeCorrect = TRUE,
includeZeroCells = FALSE,
imageID = "imageID",
cellType = "cellType",
spatialCoords = c("x", "y"),
...) {
if(is(cells, "SingleCellExperiment")|is(cells,"SpatialExperiment"))
cells <- extractSpicyInfo(cells,
imageID = imageID,
cellType = cellType,
spatialCoords = spatialCoords,
condition = condition,
subject = subject,
covariates = covariates)
if (!is(cells, "SegmentedCells")) {
stop('cells needs to be a SegmentedCells object')
}
if(is.null(from)|is.null(to)){
if (is.null(from))
from <- as.character(unique(cellType(cells)))
if (is.null(to))
to <- as.character(unique(cellType(cells)))
# from <- as.character(unique(from))
# to <- as.character(unique(to))
m1 <- rep(from, times = length(to))
m2 <- rep(to, each = length(from))
labels <- paste(m1, m2, sep = "__")
}else{
m1 <- from
m2 <- to
labels <- paste(m1, m2, sep = "__")
if(any(duplicated(labels)))stop("There are duplicated from-to pairs")
}
if (any((!to %in% cellType(cells)) |
(!from %in% cellType(cells))))
stop("to and from need to be cell type in your SegmentedCells")
nCells <- table(imageID(cells), cellType(cells))
## Find pairwise associations
if(fast){
pairwiseAssoc <- getPairwise(cells,
Rs = Rs,
sigma = sigma,
window = window,
window.length = window.length,
minLambda = minLambda,
from = from,
to = to,
fast = fast,
edgeCorrect = edgeCorrect,
includeZeroCells = includeZeroCells,
BPPARAM = BPPARAM)
pairwiseAssoc <- as.data.frame(pairwiseAssoc)
pairwiseAssoc <- pairwiseAssoc[labels]
}else{
MoreArgs1 <- list(cells = cells, dist = dist, window = window, window.length = window.length, fast = FALSE)
if (verbose)
message("Calculating pairwise spatial associations")
pairwiseAssoc <- mapply(getPairwise,
from = m1,
to = m2,
MoreArgs = MoreArgs1)
colnames(pairwiseAssoc) <- labels
}
weightFunction <- getWeightFunction(pairwiseAssoc, nCells, m1, m2, BPPARAM, weights, weightsByPair, weightFactor)
pairwiseAssoc <- as.list(data.frame(pairwiseAssoc))
names(pairwiseAssoc) <- labels
## Linear model
if (is.null(subject) & !is.null(condition)) {
if (verbose)
message("Testing for spatial differences across conditions")
MoreArgs2 <-
list(
cells = cells,
condition = condition,
covariates = covariates,
cellCounts = table(imageID(cells), cellType(cells)),
pheno = as.data.frame(imagePheno(cells))
)
linearModels <- mapply(
spatialLM,
spatAssoc = pairwiseAssoc,
from = m1,
to = m2,
weightFunction = weightFunction,
MoreArgs = MoreArgs2,
SIMPLIFY = FALSE
)
df <- cleanLM(linearModels, nsim, BPPARAM = BPPARAM)
}
## Mixed effects model
if ((!is.null(subject)) & !is.null(condition)) {
if (verbose)
message(
"Testing for spatial differences across conditions accounting for multiple images per subject"
)
MoreArgs2 <-
list(
cells = cells,
subject = subject,
condition = condition,
covariates = covariates,
cellCounts = table(imageID(cells), cellType(cells)),
pheno = as.data.frame(imagePheno(cells))
)
mixed.lmer <- mapply(
spatialMEM,
spatAssoc = pairwiseAssoc,
from = m1,
to = m2,
weightFunction = weightFunction,
MoreArgs = MoreArgs2,
SIMPLIFY = FALSE
)
mixed.lmer <- lapply(mixed.lmer, function(x){
if(is(x,"lmerModLmerTest")){
if(x@devcomp$cmp["REML"]==-Inf){return(NA)}}
x
})
df <- cleanMEM(mixed.lmer, nsim, BPPARAM = BPPARAM)
}
df$pairwiseAssoc <- pairwiseAssoc
df$comparisons <- data.frame(from = m1,
to = m2,
labels = labels)
df$weights <- weightFunction
df$nCells <- nCells
df <- methods::new('SpicyResults', df)
df
}
#' @importFrom dplyr bind_rows
cleanLM <- function(linearModels, nsim, BPPARAM) {
if (length(nsim) > 0) {
boot <- BiocParallel::bplapply(linearModels, spatialLMBootstrap, nsim = nsim, BPPARAM = BPPARAM)
#p <- do.call(rbind, p)
tBoot <- lapply(boot, function(coef) {
if(length(coef)>1){
coef <- as.data.frame(t(coef))
coef <-
split(coef, c("coefficient", "se", "statistic", "p.value"))
}else{
coef <- list(coefficient = NA, se = NA, statistic = NA, p.value = NA)
}
coef
})
df <- apply(do.call(rbind, tBoot), 2, function(x)
do.call(rbind, x))
df <- lapply(df, function(x) {
rownames(x) <- names(linearModels)
x
})
} else {
tLm <- lapply(linearModels, function(LM) {
if(is(LM,"lm")){
coef <- as.data.frame(t(summary(LM)$coef))
coef <-
split(coef, c("coefficient", "se", "statistic", "p.value"))
}else{
n <- data.frame(NA)
colnames(n) <- "(Intercept)"
coef <- list(coefficient = n, se = n, statistic = n, p.value = n)
}
coef
})
df <- do.call("rbind", tLm)
df <- suppressWarnings(apply(df, 2, function(x){
dplyr::bind_rows(x)
}))
df <- lapply(df, function(x) {
rownames(x) <- names(linearModels)
x
})
}
df
}
#' @importFrom dplyr bind_rows
cleanMEM <- function(mixed.lmer, nsim, BPPARAM) {
if (length(nsim) > 0) {
boot <- BiocParallel::bplapply(mixed.lmer, spatialMEMBootstrap, nsim = nsim, BPPARAM = BPPARAM)
#boot <- lapply(mixed.lmer, spatialMEMBootstrap, nsim = nsim)
#p <- do.call(rbind, p)
tBoot <- lapply(boot, function(coef) {
if(length(coef)>1){
coef <- as.data.frame(t(coef))
coef <-
split(coef, c("coefficient", "se", "statistic", "p.value"))
}else{
n <- data.frame(NA)
colnames(n) <- "(Intercept)"
coef <- list(coefficient = n, df = n, p.value = n, se = n, statistic = n)
}
coef
})
df <- apply(do.call(rbind, tBoot), 2, function(x)
do.call(rbind, x))
df <- lapply(df, function(x) {
rownames(x) <- names(mixed.lmer)
x
})
} else {
tLmer <- lapply(mixed.lmer, function(lmer) {
if(is(lmer,"lmerMod")){
coef <- as.data.frame(t(summary(lmer)$coef))
coef <-
split(coef,
c("coefficient", "se", "df", "statistic", "p.value"))
}else{
n <- data.frame(NA)
colnames(n) <- "(Intercept)"
coef <- list(coefficient = n, df = n, p.value = n, se = n, statistic = n)
}
coef
})
df <- do.call("rbind", tLmer)
df <- suppressWarnings(apply(df, 2, function(x){
dplyr::bind_rows(x)
}))
df <- lapply(df, function(x) {
rownames(x) <- names(mixed.lmer)
x
})
}
df
}
#' Get statistic from pairwise L curve of a single image.
#'
#' @param cells A SegmentedCells or data frame that contains at least the
#' variables x and y, giving the location coordinates of each cell, and cellType.
#' @param from The 'from' cellType for generating the L curve.
#' @param to The 'to' cellType for generating the L curve.
#' @param dist The distance at which the statistic is obtained.
#' @param window Should the window around the regions be 'square', 'convex' or 'concave'.
#' @param window.length A tuning parameter for controlling the level of concavity
#' when estimating concave windows.
#' @param Rs A vector of the radii that the measures of association should be calculated.
#' @param sigma A numeric variable used for scaling when fitting inhomogeneous L-curves.
#' @param minLambda Minimum value for density for scaling when fitting inhomogeneous L-curves.
#' @param fast A logical describing whether to use a fast approximation of the
#' inhomogeneous L-curves.
#' @param edgeCorrect A logical indicating whether to perform edge correction.
#' @param includeZeroCells A logical indicating whether to include cells with
#' zero counts in the pairwise association calculation.
#' @param BPPARAM A BiocParallelParam object.
#' @param imageID The imageID if using a SingleCellExperiment or SpatialExperiment.
#' @param cellType The cellType if using a SingleCellExperiment or SpatialExperiment.
#' @param spatialCoords The spatialCoords if using a SingleCellExperiment or SpatialExperiment.
#'
#' @return Statistic from pairwise L curve of a single image.
#'
#'
#' @examples
#' data("diabetesData")
#' pairAssoc <- getPairwise(diabetesData[1,])
#' @export
#' @importFrom BiocParallel bplapply
getPairwise <- function(cells, from = NULL, to = NULL, dist = NULL, window = "convex", window.length = NULL, Rs = c(20, 50, 100), sigma = NULL, minLambda = 0.05, fast = TRUE, edgeCorrect = TRUE, includeZeroCells = TRUE, BPPARAM=BiocParallel::SerialParam(),
imageID = "imageID", cellType = "cellType", spatialCoords = c("x", "y")) {
cells2 <- prepCellSummary(cells, spatialCoords, cellType, imageID, bind = FALSE)
if(is.null(from))from <- levels(cells2$cellType)
if(is.null(to))to <- levels(cells2$cellType)
if(fast){
pairwiseVals <- BiocParallel::bplapply(cells2,
inhomLPair,
Rs = Rs,
sigma = sigma,
window = window,
window.length = window.length,
minLambda = minLambda,
from = from,
to = to,
edgeCorrect = edgeCorrect,
includeZeroCells = includeZeroCells,
BPPARAM=BPPARAM)
return(do.call("rbind",pairwiseVals ))
}else{
pairwiseVals <- lapply(cells2,
getStat,
from = from,
to = to,
dist = dist, window, window.length)
return(unlist(pairwiseVals))
}
}
#' Get proportions from a SegmentedCells, SingleCellExperiment, SpatialExperiment or data.frame.
#'
#' @param cells SegmentedCells, SingleCellExperiment, SpatialExperiment or data.frame
#' @param feature The feature of interest
#' @param imageID The imageID's
#'
#' @return Proportions
#'
#'
#' @examples
#' data("diabetesData")
#' prop <- getProp(diabetesData)
#' @export
#' @importFrom SummarizedExperiment colData
getProp <- function(cells, feature = "cellType", imageID = "imageID") {
if (is.data.frame(cells)) {
df <- cells[,c(imageID, feature)]
}
if (is(cells, "SingleCellExperiment")|is(cells, "SpatialExperiment")) {
df <- as.data.frame(SummarizedExperiment::colData(cells))[,c(imageID, feature)]
}
if (is(cells, "SegmentedCells")) {
cellSummary <- cellSummary(cells, bind = TRUE)
df <- as.data.frame(cellSummary[,c(imageID, feature)])
}
tab <- table(df[,imageID], df[,feature])
tab <- sweep(tab,1,rowSums(tab), "/")
as.data.frame.matrix(tab)
}
#' @importFrom spatstat.explore Lcross
getStat <- function(cells, from, to, dist, window, window.length) {
pppCell <- pppGenerate(cells, window, window.length)
L <- tryCatch({
spatstat.explore::Lcross(pppCell,
from = from,
to = to,
correction = "best")
}, error = function(e) {
})
if (!is(L,"fv")) {
return(NA)
}
if (is.null(dist))
dist <- max(L$r)
theo = L$theo[L$r <= dist]
iso = L$iso[L$r <= dist]
mean(iso - theo)
}
# Performs bootstrapping to estimate p-value.
#' @importFrom lme4 fixef
#' @importFrom lmerTest lmer
#' @importFrom lme4 lmerControl
#' @importFrom stats formula weights sd
spatialMEMBootstrap <- function(mixed.lmer, nsim = 19) {
functionToReplicate <- function(x) {
toGet <- sample(nrow(x@frame), replace = TRUE)
spatialData <- x@frame[toGet,]
spatialData$weights <- spatialData$`(weights)`
mixed.lmer1 <- suppressWarnings(suppressMessages(tryCatch({lmerTest::lmer(stats::formula(x),
data = spatialData,
weights = weights,
control = lme4::lmerControl(calc.derivs = FALSE))},
error = function(e){})))
if (!is(mixed.lmer1,"lmerMod")) {
return(rep(NA,ncol(coef(mixed.lmer)[[1]])))
}
summary(mixed.lmer1)$coef[, "t value"]
}
if(!is(mixed.lmer, "lmerMod"))return(NA)
stats <- replicate(nsim, functionToReplicate(x = mixed.lmer))
stats <- t(stats)
fe <- lme4::fixef(mixed.lmer)
# s1 <- sweep(stats,2,sign(fe),"*")
# pval <- colMeans(s1 < 0, na.rm = TRUE)*2 #+ colMeans(sweep(s1,2,2*abs(summary(mixed.lmer)$coef[, "t value"]),">"))
pval <- mapply(pmin,colMeans(stats < 0,na.rm = TRUE), colMeans(stats > 0,na.rm = TRUE), MoreArgs = list(na.rm = TRUE))*2
df <-
data.frame(
coefficient = fe,
se = apply(stats, 2, stats::sd, na.rm = TRUE),
statistic = fe/apply(stats, 2, stats::sd, na.rm = TRUE),
p.value = pval
)
df
}
#' @importFrom stats p.adjust
.show_SpicyResults <- function(df) {
pval <- as.data.frame(df$p.value)
cond <- colnames(pval)[grep('condition', colnames(pval))]
message(df$test)
message("Number of cell type pairs: ", nrow(pval), "\n")
message("Number of differentially localised cell type pairs: \n")
if (nrow(pval) == 1)
print(sum(pval[cond] < 0.05, na.rm = TRUE))
if (nrow(pval) > 1)
print(colSums(apply(pval[cond], 2, p.adjust, 'fdr') < 0.05, na.rm = TRUE))
}
setMethod("show", methods::signature(object = "SpicyResults"), function(object) {
.show_SpicyResults(object)
})
#' @importFrom lmerTest lmer
#' @importFrom stats predict weights
#' @importFrom methods is
spatialMEM <-
function(spatAssoc,
from,
to,
cells,
subject,
condition,
covariates,
weightFunction,
cellCounts,
pheno) {
spatAssoc[is.na(spatAssoc)] <- 0
count1 <- cellCounts[, from]
count2 <- cellCounts[, to]
# filter <- !is.na(spatAssoc)
#
# if (sum(filter) < 3)
# return(NA)
spatialData <-
data.frame(spatAssoc,
condition = pheno[, condition],
subject = pheno[, subject],
pheno[covariates])
formula <- 'spatAssoc ~ condition + (1|subject)'
if (!is.null(covariates))
formula <-
paste('spatAssoc ~ condition + (1|subject)',
paste(covariates, collapse = '+'),
sep = "+")
spatialData$weights = weightFunction
mixed.lmer <- suppressWarnings(suppressMessages(tryCatch({lmerTest::lmer(stats::formula(formula),
data = spatialData,
weights = weights)},
error = function(e) {
})))
if (!is(mixed.lmer,"lmerMod")) {
return(NA)
}
mixed.lmer
}
#' @importFrom stats predict lm
spatialLM <-
function(spatAssoc,
from,
to,
cells,
condition,
covariates,
weightFunction,
cellCounts,
pheno) {
count1 <- cellCounts[, from]
count2 <- cellCounts[, to]
#print(pheno)
spatialData <-
data.frame(spatAssoc, condition = pheno[, condition], pheno[,covariates])
# spatialData <- spatialData[filter, ]
# count1 <- count1[filter]
# count2 <- count2[filter]
#print(spatialData)
formula <- 'spatAssoc ~ condition'
if (!is.null(covariates))
formula <-
paste('spatAssoc ~ condition',
paste(covariates, collapse = '+'),
sep = "+")
lm1 <- tryCatch({stats::lm(stats::formula(formula),
data = spatialData,
weights = weightFunction)},
error = function(e) {
})
if (!is(lm1,"lm")) {
return(NA)
}
lm1
}
#' @importFrom stats sd coef
spatialLMBootstrap <- function(linearModels, nsim=19) {
functionToReplicate <- function(x) {
toGet <- sample(nrow(x$model), replace = TRUE)
spatAssocBoot <- x$model$spatAssoc[toGet]
conditionBoot <- x$model$condition[toGet]
weightsBoot <- x$weights[toGet]
spatialDataBoot <- data.frame(spatAssoc = spatAssocBoot,
condition = conditionBoot)
lm1 <- tryCatch({stats::lm(spatAssoc ~ condition,
data = spatialDataBoot,
weights = weightsBoot)},
error = function(e) {
})
if (!is(lm1,"lm")) {
return(NA)
}
lm1$coefficients[2]
}
if(!is(linearModels, "lm"))return(NA)
stats <- replicate(nsim, functionToReplicate(x = linearModels))
df <- coef(summary(linearModels))
pval <- pmin(mean(stats < 0,na.rm = TRUE), mean(stats > 0, na.rm = TRUE), na.rm = TRUE) * 2
df[2,4] <- pval
df
}
###########################
#
# Generate distances
#
###########################
#' @importFrom spatstat.geom owin convexhull ppp
#' @importFrom concaveman concaveman
makeWindow <-
function(data,
window = "square",
window.length = NULL) {
data = data.frame(data)
ow <-
spatstat.geom::owin(xrange = range(data$x), yrange = range(data$y))
if (window == "convex") {
p <- spatstat.geom::ppp(data$x, data$y, ow)
ow <- spatstat.geom::convexhull(p)
}
if (window == "concave") {
message("Concave windows are temperamental. Try choosing values of window.length > and < 1 if you have problems.")
if(is.null(window.length)){
window.length <- (max(data$x) - min(data$x))/20
}else{
window.length <- (max(data$x) - min(data$x))/20 * window.length
}
dist <- (max(data$x) - min(data$x)) / (length(data$x))
bigDat <-
do.call("rbind", lapply(as.list(as.data.frame(t(data[, c("x", "y")]))), function(x)
cbind(
x[1] + c(0, 1, 0,-1,-1, 0, 1,-1, 1) * dist,
x[2] + c(0, 1, 1, 1,-1,-1,-1, 0, 0) * dist
)))
ch <-
concaveman::concaveman(bigDat,
length_threshold = window.length,
concavity = 1)
poly <- as.data.frame(ch[nrow(ch):1, ])
colnames(poly) <- c("x", "y")
ow <-
spatstat.geom::owin(
xrange = range(poly$x),
yrange = range(poly$y),
poly = poly
)
}
ow
}
#' @importFrom spatstat.geom ppp
pppGenerate <- function(cells, window, window.length) {
ow <- makeWindow(cells, window, window.length)
pppCell <- spatstat.geom::ppp(
cells$x,
cells$y,
window = ow,
marks = cells$cellType
)
pppCell
}
#' @importFrom spatstat.explore density.ppp
#' @importFrom spatstat.geom closepairs nearest.valid.pixel area ppp
#' @importFrom tidyr pivot_longer
#' @importFrom dplyr left_join
inhomLPair <-
function (data,
Rs = c(20, 50, 100),
sigma = NULL,
window = "convex",
window.length = NULL,
minLambda = 0.05,
from = NULL,
to = NULL,
edgeCorrect = TRUE,
includeZeroCells = TRUE) {
ow <- makeWindow(data, window, window.length)
X <-
spatstat.geom::ppp(
x = data$x,
y = data$y,
window = ow,
marks = data$cellType
)
if (is.null(Rs))
Rs = c(20, 50, 100)
# if (is.null(sigma))
# sigma = 100000
maxR <- min(ow$xrange[2]- ow$xrange[1], ow$yrange[2]- ow$yrange[1])/2.01
Rs <- unique(pmin(c(0, sort(Rs)),maxR))
if(!is.null(sigma)){
den <- spatstat.explore::density.ppp(X, sigma = sigma)
den <- den / mean(den)
den$v <- pmax(den$v, minLambda)
}
if(is.null(from)) from <- levels(data$cellType)
if(is.null(to)) to <- levels(data$cellType)
use <- data$cellType %in% c(from, to)
data <- data[use,]
X <- X[use,]
p <- spatstat.geom::closepairs(X, max(Rs), what = "ijd", distinct = FALSE)
n <- X$n
p$j <- data$cellID[p$j]
p$i <- data$cellID[p$i]
cT <- data$cellType
names(cT) <- data$cellID
p$d <- cut(p$d, Rs, labels = Rs[-1], include.lowest = TRUE)
# inhom density
p$wt <- rep(1,length(p$d))
if(!is.null(sigma)){
np <- spatstat.geom::nearest.valid.pixel(X$x, X$y, den)
w <- den$v[cbind(np$row, np$col)]
names(w) <- data$cellID
p$wt <- 1/w[p$j]*mean(w)
rm(np)
}
lam <- table(data$cellType)/spatstat.geom::area(X)
# p$wt <- as.numeric(p$wt/lam[cT[p$j]])
p$cellTypeJ <- cT[p$j]
p$cellTypeI <- cT[p$i]
p$i <- factor(p$i, levels = data$cellID)
if(edgeCorrect){
edge <- sapply(Rs[-1],function(x)borderEdge(X,x), simplify = FALSE)
edge <- do.call("cbind",edge)
edge <- as.data.frame(edge)
colnames(edge) <- Rs[-1]
edge$i <- data$cellID
edge <- tidyr::pivot_longer(edge,-i,"d")
p <- dplyr::left_join(as.data.frame(p), edge, c("i", "d"))
}else{
p <- as.data.frame(p)
p$value <- 1
}
p$d <- factor(p$d, levels = Rs[-1])
p <- p[p$i != p$j, ]
use <- p$cellTypeI %in% from & p$cellTypeJ %in% to
p <- p[use,]
r <- inhomL(p, lam, X, Rs)
wt <- r$wt
names(wt) <- paste(r$cellTypeI, r$cellTypeJ, sep = "__")
m1 <- rep(from, times = length(to))
m2 <- rep(to, each = length(from))
labels <- paste(m1, m2, sep = "__")
assoc <- rep(-sum(Rs), length(labels))
names(assoc) <- labels
if(!includeZeroCells)assoc[!(m1%in%X$marks&m2%in%X$marks)] = NA
assoc[names(wt)] <- wt
names(assoc) <- labels
assoc
}
#' @importFrom data.table as.data.table setkey CJ .SD ":="
#' @importFrom spatstat.geom area
inhomL <-
function (p, lam, X, Rs) {
r <- data.table::as.data.table(p)
r$wt <- r$wt/r$value/as.numeric(lam[r$cellTypeJ])/as.numeric(lam[r$cellTypeI])/spatstat.geom::area(X)
r <- r[,j:=NULL]
r <- r[,value:=NULL]
r <- r[,i:=NULL]
data.table::setkey(r, d, cellTypeI, cellTypeJ)
r <- r[data.table::CJ(d, cellTypeI, cellTypeJ, unique = TRUE)
][, lapply(.SD, sum), by = .(d, cellTypeI, cellTypeJ)
][is.na(wt), wt := 0]
r <- r[, wt := cumsum(wt), by = list(cellTypeI, cellTypeJ)]
r <- r[, list(wt=sum(sqrt(wt/pi))), by=.(cellTypeI, cellTypeJ)]
r$wt <- r$wt - sum(Rs)
r <- as.data.frame(r)
r
}
#' @importFrom spatstat.geom union.owin border inside.owin solapply intersect.owin area discs
borderEdge <- function(X, maxD){
W <-X$window
bW <- spatstat.geom::union.owin(spatstat.geom::border(W,maxD, outside = FALSE),
spatstat.geom::border(W,2, outside = TRUE))
inB <- spatstat.geom::inside.owin(X$x, X$y, bW)
e <- rep(1, X$n)
if(any(inB)){
circs <-spatstat.geom:: discs(X[inB], maxD, separate = TRUE)
circs <- spatstat.geom::solapply(circs, spatstat.geom::intersect.owin, X$window)
areas <- unlist(lapply(circs, spatstat.geom::area))/(pi*maxD^2)
e[inB] <- areas
}
e
}
#' @importFrom scam scam
#' @importFrom stats quantile
calcWeights <- function(rS, M1, M2, nCells, weightFactor){
count1 <- as.vector(nCells[, M1])
count2 <- as.vector(nCells[, M2])
rS <- as.vector(rS)
toWeight <- !is.na(rS)
resSqToWeight <- rS[toWeight]
count1ToWeight <- count1[toWeight]
count2ToWeight <- count2[toWeight]
if(length(count1ToWeight)<=20){
warning("A cell type pair is seen less than 20 times, not using weights for this pair.")
#weightFunction <- lm(log10(resSqToWeight)~ 1)
weightFunction <- rep(1, nrow(pairwiseAssoc))
return(weightFunction)
}
weightFunction <- scam::scam(log10(resSqToWeight+1)~ s(log10(count1ToWeight+1),bs="mpd")+s(log10(count2ToWeight+1), bs="mpd"))#, optimizer = "nlm.fd")
z1 <- suppressWarnings(stats::predict(weightFunction, data.frame(count1ToWeight = as.numeric(count1),
count2ToWeight = as.numeric(count2))))
#w <- 1 / sqrt(z1 - min(z1) + 1)
#z1[which(as.numeric(count1)==0|as.numeric(count2)==0)] <- NA
w <- 1/pmax(z1, stats::quantile(z1[z1>0.1], 0.01, na.rm = TRUE))
w <- w / mean(w, na.rm = TRUE)
w^weightFactor
}
#' @importFrom BiocParallel bpmapply
getWeightFunction <- function(pairwiseAssoc, nCells, m1, m2, BPPARAM, weights, weightsByPair, weightFactor){
if(!weights){
# weightFunction <- sapply(colnames(pairwiseAssoc), function(x){NULL}, simplify = FALSE, USE.NAMES = TRUE)
# return(weightFunction)
weightFunction <- rep(1, nrow(pairwiseAssoc)*ncol(pairwiseAssoc))
pair <- rep(colnames(pairwiseAssoc), each = nrow(pairwiseAssoc))
weightFunction <- split(weightFunction, pair)
return(weightFunction)
}
resSq <- apply(pairwiseAssoc, 2, function(x){
if(stats::sd(x, na.rm = TRUE)>0){
return((x - mean(x, na.rm = TRUE))^2)
}else{
return(rep(NA,length(x)))
}
}
)
if(weightsByPair){
weightFunction <- BiocParallel::bpmapply(calcWeights, rS = as.list(as.data.frame(resSq)), M1 = m1, M2 = m2, BPPARAM=BPPARAM, MoreArgs = list(nCells = nCells, weightFactor), SIMPLIFY = FALSE)
}else{
weightFunction <- calcWeights(m1, m2, rS = resSq, nCells, weightFactor)
pair <- rep(colnames(pairwiseAssoc), each = nrow(pairwiseAssoc))
weightFunction <- split(weightFunction, pair)
# weightFunction <- sapply(colnames(pairwiseAssoc), function(x, weightFunction){weightFunction}, simplify = FALSE, USE.NAMES = TRUE, weightFunction = weightFunction)
}
weightFunction[colnames(pairwiseAssoc)]
}
#' @importFrom SummarizedExperiment colData
#' @importFrom SpatialExperiment spatialCoords
prepCellSummary <- function(cells, spatialCoords, cellType, imageID, bind = FALSE){
if (is.data.frame(cells)) {
cells <- SegmentedCells(cells,
spatialCoords = spatialCoords,
cellTypeString = cellType,
imageIDString = imageID)
cellSummary <- cellSummary(cells, bind = bind)
}
if (class(cells) == "SingleCellExperiment") {
cells <- colData(cells)
cells <- SegmentedCells(cells,
spatialCoords = spatialCoords,
cellTypeString = cellType,
imageIDString = imageID)
cellSummary <- cellSummary(cells, bind = bind)
}
if (class(cells) == "SpatialExperiment") {
cells <- cbind(colData(cells), spatialCoords(cells))
cells <- SegmentedCells(cells,
spatialCoords = spatialCoords,
cellTypeString = cellType,
imageIDString = imageID)
cellSummary <- cellSummary(cells, bind = bind)
}
if (is(cells, "SegmentedCells")) {
cellSummary <- cellSummary(cells, bind = bind)
}
cellSummary
}
extractSpicyInfo <- function(cells,
imageID = imageID,
cellType = cellType,
spatialCoords = spatialCoords,
condition = condition,
subject = subject,
covariates = covariates){
extra <- c(condition, subject, covariates)
if (class(cells) == "SingleCellExperiment") {
cells <- colData(cells)
colnames(cells)[colnames(cells)%in%extra] <- paste0("phenotype_",colnames(cells)[colnames(cells)%in%extra])
cells <- SegmentedCells(cells,
spatialCoords = spatialCoords,
cellTypeString = cellType,
imageIDString = imageID)
}
if (class(cells) == "SpatialExperiment") {
cells <- cbind(colData(cells), spatialCoords(cells))
colnames(cells)[colnames(cells)%in%extra] <- paste0("phenotype_",colnames(cells)[colnames(cells)%in%extra])
cells <- SegmentedCells(cells,
spatialCoords = spatialCoords,
cellTypeString = cellType,
imageIDString = imageID)
}
cells
}
#' Perform a simple wilcoxon-rank-sum test or t-test on the columns of a data fram
#'
#' @param df A data.frame or SingleCellExperiment, SpatialExperiment
#' @param condition The condition of interest
#' @param type The type of test, "wilcox" or "ttest"
#' @param feature Can be used to calculate the proportions of this feature for each image
#' @param imageID The imageID's if presenting a SingleCellExperiment
#'
#' @return Proportions
#'
#'
#' @examples
#'
#' # Test for an association with long-duration diabetes
#' # This is clearly ignoring the repeated measures...
#' data("diabetesData")
#' props <- getProp(diabetesData)
#' condition <- imagePheno(diabetesData)$stage
#' names(condition) <- imagePheno(diabetesData)$imageID
#' condition <- condition[condition%in% c("Long-duration", "Onset")]
#' test <- colTest(props[names(condition), ], condition)
#' @export
#' @importFrom SummarizedExperiment colData
#' @importFrom SingleCellExperiment SingleCellExperiment
#' @importFrom stats wilcox.test t.test
#' @importFrom S4Vectors as.data.frame
colTest <- function(df, condition, type = "ttest", feature = NULL, imageID = "imageID"){
if(is(df, "SingleCellExperiment")|is(df, "SpatialExperiment")){
if(is.null(feature))stop("'feature' is still null")
x <- df@colData
x <- x[,c(imageID,condition)]
x <- unique(x)
condition <- x[[condition]]
names(condition) <- x[[imageID]]
df <- getProp(df, imageID = imageID, feature = feature)
condition <- condition[rownames(df)]
}
test <- apply(df, 2, function(x){
if(type == "wilcox")test <- stats::wilcox.test(x ~ condition)
if(type == "ttest")test <- stats::t.test(x ~ condition)
signif(c(test$estimate, tval <- test$statistic, pval = test$p.value),2)
})
test <- as.data.frame(t(test))
test$adjPval <- signif(stats::p.adjust(test$pval, "fdr"),2)
test$cluster <- rownames(test)
test <- test[order(test$pval),]
test
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.