Performs global spatial gene variability analysis across multiple length scales using kernel-based methods.
Usage
spacelink_global(
normalized_counts,
spatial_coords,
covariates = NULL,
lengthscales = NULL,
n_lengthscales = 5,
M = 1,
n_workers = 1
)Arguments
- normalized_counts
Normalized count matrix. Rows and columns indicate genes and spots, respectively.
- spatial_coords
Spatial coordinates matrix. Rows indicate spots.
- covariates
Spatial features (without intercept). Rows indicate spots.
- lengthscales
If NULL, length-scales are calculated using the two-step approach.
- n_lengthscales
Number of length-scales (i.e., number of kernels). Default is 5.
- M
Controls the minimum length-scale. Minimum length-scale is set to be the minimum distance multiplied by M. Default is 1.
- n_workers
Number of workers for parallelization. Default is 1.
Value
A data frame containing:
- tau.sq
Nugget variance component
- sigma.sq1, sigma.sq2, ...
Spatial variance components for each kernel
- raw_ESV
Effective Spatial Variability score
- pval
Combined p-value from score tests
- padj
Benjamini-Hochberg adjusted p-values
- ESV
ESV adjusted for multiple testing (0 if padj > 0.05)
Examples
# Set up simulated data
library(spacelink)
set.seed(123)
n_spots <- 100
n_genes <- 10
# Generate spatial coordinates (2D grid)
coords <- expand.grid(
x = seq(0, 10, length.out = sqrt(n_spots)),
y = seq(0, 10, length.out = sqrt(n_spots))
)
# Simulate (normalized) expression data with spatial autocorrelation
expr_data <- matrix(nrow=n_genes, ncol=n_spots)
rownames(expr_data) <- paste0("Gene_", 1:n_genes)
D <- as.matrix(dist(coords))
K <- exp(-D); chol_K <- chol(K)
for(i in 1:(n_genes/2)){
expr_data[i,] <- matrix(rnorm(n_spots, 0, 1), nrow=1, ncol=n_spots)
expr_data[i,] <- expr_data[i,] + i*matrix(rnorm(n_spots, 0, 1), 1, n_spots)%*%chol_K
}
# Simulate (normalized) expression data without spatial autocorrelation
for(i in (n_genes/2+1):n_genes){
expr_data[i,] <- i*matrix(rnorm(n_spots, 0, 1), nrow=1, ncol=n_spots)
}
global_results <- spacelink_global(normalized_counts = expr_data, spatial_coords = coords)
print(global_results[, c("raw_ESV", "pval", "ESV", "padj")])
#> raw_ESV pval ESV padj
#> Gene_1 0.15061194 0.158662165 0.0000000 0.31732433
#> Gene_2 0.36958214 0.013996958 0.3695821 0.04665653
#> Gene_3 0.24545260 0.063870746 0.0000000 0.15967687
#> Gene_4 0.42172756 0.007834405 0.4217276 0.03917202
#> Gene_5 0.43943594 0.004087093 0.4394359 0.03917202
#> Gene_6 0.00000000 0.748318133 0.0000000 0.86450152
#> Gene_7 0.00000000 0.805014384 0.0000000 0.86450152
#> Gene_8 0.06886313 0.326967744 0.0000000 0.54494624
#> Gene_9 0.00000000 0.835162432 0.0000000 0.86450152
#> Gene_10 0.00000000 0.864501520 0.0000000 0.86450152