### Install spacelink from GitHub with devtools

Source: https://hanbyul-lee.github.io/spacelink/index.html

Installs the development version of spacelink from GitHub. Requires R ≥ 4.0.0 and a supported OS (Windows, macOS, Linux). Uncomment the devtools installation line if devtools is not already installed.

```r
# install.packages("devtools")
devtools::install_github("hanbyul-lee/spacelink")
```

--------------------------------

### Load Example Dataset

Source: https://hanbyul-lee.github.io/spacelink/articles/AD_data_analysis.html

Load control and late-stage AD Visium samples and inspect the dimensions of the count matrices.

```R
control_spe <- readRDS("50500136_1.rds")
control_counts <- as.matrix(counts(control_spe))
control_spatial_coords <- as.data.frame(spatialCoords(control_spe))

AD_spe <- readRDS("18414513_1.rds")
AD_counts <- as.matrix(counts(AD_spe))
AD_spatial_coords <- as.data.frame(spatialCoords(AD_spe))

dim(control_counts)
```

```R
## [1] 36601  1663
```

```R
dim(AD_counts)
```

```R
## [1] 36601  1273
```

--------------------------------

### Example: Running spacelink on simulated data

Source: https://hanbyul-lee.github.io/spacelink/reference/spacelink.html

Complete example that generates synthetic spatial expression data and runs spacelink. The output shows ESV, pval, and padj for each gene. Note the \donttest{} wrapper indicates the example is not run during package checks.

```r
# \donttest{
library(spacelink)
set.seed(123)
n_spots <- 100
n_genes <- 10

coords <- expand.grid(
  x = seq(0, 10, length.out = sqrt(n_spots)),
  y = seq(0, 10, length.out = sqrt(n_spots))
)

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 / 2); chol_K <- chol(K)
for (i in 1:(n_genes / 2)) {
  expr_data[i, ] <- rnorm(n_spots) + i * rnorm(n_spots) %*% chol_K
}
for (i in (n_genes / 2 + 1):n_genes) {
  expr_data[i, ] <- rnorm(n_spots)
}

results <- spacelink(normalized_counts = expr_data, spatial_coords = coords)
print(results[, c("ESV", "pval", "padj")])
#>                ESV         pval         padj
#> Gene_1  0.36240416 1.316941e-02 2.633882e-02
#> Gene_2  0.67450493 1.641875e-05 5.472917e-05
#> Gene_3  0.63469294 2.190094e-04 5.475235e-04
#> Gene_4  0.67450493 1.096747e-06 5.483735e-06
#> Gene_5  0.67450493 2.671183e-08 2.671183e-07
#> Gene_6  0.00000000 7.483181e-01 8.645015e-01
#> Gene_7  0.00000000 8.050144e-01 8.645015e-01
#> Gene_8  0.06886313 3.269677e-01 5.449462e-01
#> Gene_9  0.00000000 8.351624e-01 8.645015e-01
#> Gene_10 0.00000000 8.645015e-01 8.645015e-01
# }
```

--------------------------------

### Install remaining spacelink dependencies from CRAN

Source: https://hanbyul-lee.github.io/spacelink/index.html

Installs the other dependency packages from CRAN. If automatic installation fails, install these packages manually.

```r
install.packages(c("Rcpp", "RcppML", "pracma", "RANN", "fields", "Matrix"))
```

--------------------------------

### Load and inspect Visium example dataset

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_visium_data.html

Load the built-in Visium human DLPFC dataset and inspect the gene expression counts, spatial coordinates, and cell type proportions.

```R
data(Visium_human_DLPFC)
counts <- Visium_human_DLPFC$counts
spatial_coords <- Visium_human_DLPFC$spatial_coords
cell_type_proportions <- Visium_human_DLPFC$cell_type_proportions
head(counts[,1:10])
```

```R
## 6 x 10 sparse Matrix of class "dgTMatrix"
##   [[ suppressing 10 column names ‘AAACAAGTATCTCCCA-1’, ‘AAACAATCTACTAGCA-1’, ‘AAACACCAATAACTGC-1’ ... ]]                               
## MIR1302-2HG . . . . . . . . . .
## FAM138A     . . . . . . . . . .
## OR4F5       . . . . . . . . . .
## AL627309.1  . . . . . . . . . .
## AL627309.3  . . . . . . . . . .
## AL627309.2  . . . . . . . . . .
```

```R
head(spatial_coords)
```

```R
##                    pxl_col_in_fullres pxl_row_in_fullres
## AAACAAGTATCTCCCA-1               9791               8468
## AAACAATCTACTAGCA-1               5769               2807
## AAACACCAATAACTGC-1               4068               9505
## AAACAGAGCGACTCCT-1               9271               4151
## AAACAGCTTTCAGAAG-1               3393               7583
## AAACAGGGTCTATATT-1               3665               8064
```

```R
head(cell_type_proportions)
```

```R
##                     endothelial       L2_3           L4           L6 oligodendrocyte          OPC       astro         Inh    microglia
## AAACAAGTATCTCCCA-1 2.331552e-02 0.86383709 9.941496e-05 8.994572e-05     0.008841064 8.994572e-05 0.037094727 0.060440671 0.0061916188
## AAACAATCTACTAGCA-1 7.952809e-02 0.55623941 1.284478e-04 1.413909e-02     0.001048564 1.284478e-04 0.083634323 0.264896732 0.0002568957
## AAACACCAATAACTGC-1 4.144100e-05 0.07909515 4.144100e-05 8.593492e-02     0.826629114 7.799136e-05 0.006276291 0.000124323 0.0017793364
## AAACAGAGCGACTCCT-1 4.419690e-04 0.70660998 1.437067e-04 1.437067e-04     0.012711896 1.437067e-04 0.119882025 0.133047245 0.0268757663
## AAACAGCTTTCAGAAG-1 1.398113e-02 0.52170271 8.882893e-05 2.874930e-01     0.050319200 8.882893e-05 0.026019353 0.096453571 0.0038533732
## AAACAGGGTCTATATT-1 8.872097e-05 0.37248441 8.872097e-05 3.500929e-01     0.206785555 8.872097e-05 0.029330741 0.040862758 0.0001774419
```

--------------------------------

### Example: running spacelink_ctSVG on simulated data

Source: https://hanbyul-lee.github.io/spacelink/reference/spacelink_ctSVG.html

Complete example that generates simulated spatial expression data and runs spacelink_ctSVG on a focal cell type. The output shows p-values and ESV scores for each gene.

```R
# \donttest{
library(spacelink)
set.seed(123)
n_spots    <- 100
n_genes    <- 10
n_celltypes <- 2

coords <- expand.grid(
  x = seq(0, 10, length.out = sqrt(n_spots)),
  y = seq(0, 10, length.out = sqrt(n_spots))
)

cell_type_props <- matrix(nrow = n_spots, ncol = n_celltypes)
cell_type_props[1:(n_spots / 2), 1]            <- runif(n_spots / 2, 0, 0.5)
cell_type_props[(n_spots / 2 + 1):n_spots, 1]  <- runif(n_spots / 2, 0.5, 1)
cell_type_props[, 2] <- 1 - cell_type_props[, 1]
colnames(cell_type_props) <- paste0("Cell_type_", 1:n_celltypes)

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)
Sigma <- chol(K) %*% diag(cell_type_props[, "Cell_type_1"])
for (i in 1:(n_genes / 2))
  expr_data[i, ] <- rnorm(n_spots) + i * rnorm(n_spots) %*% Sigma
Sigma <- chol(K) %*% diag(cell_type_props[, "Cell_type_2"])
for (i in (n_genes / 2 + 1):n_genes)
  expr_data[i, ] <- rnorm(n_spots) + rnorm(n_spots) %*% Sigma

results <- spacelink_ctSVG(
  normalized_counts     = expr_data,
  spatial_coords        = coords,
  cell_type_proportions = cell_type_props,
  focal_cell_type       = "Cell_type_1"
)
print(results[, c("pval", "ESV")])
#>                 pval          ESV
#> Gene_1  8.204439e-02 4.167043e-01
#> Gene_2  5.121033e-03 6.063774e-01
#> Gene_3  5.702071e-04 6.964893e-01
#> Gene_4  5.175750e-06 6.582678e-01
#> Gene_5  1.856661e-04 6.579634e-01
#> Gene_6  8.726797e-01 3.665291e-06
#> Gene_7  6.295618e-01 2.876897e-03
#> Gene_8  8.989732e-01 8.545377e-06
#> Gene_9  7.814450e-01 3.790261e-06
#> Gene_10 7.961633e-01 3.368537e-06
# }
```

--------------------------------

### Source File License Header Template

Source: https://hanbyul-lee.github.io/spacelink/LICENSE.html

Attach this notice to the start of each source file to state the exclusion of warranty and provide a pointer to the full license text.

```text
<one line to give the program's name and a brief idea of what it does.>
Copyright (C) <year>  <name of author>

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <http://www.gnu.org/licenses/>.
```

--------------------------------

### Raw count matrix expr_data (first 10 genes)

Source: https://hanbyul-lee.github.io/spacelink/articles/spacelink_vignette.html

Example output of the raw count matrix `expr_data` for the first 10 genes and 10 spots. This is the input data that must be log-transformed and normalized before running `spacelink`.

```R
##         [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## Gene_1     1    0    5    1    0    4    1    0    0     0
## Gene_2     0    0    0    1    0   29    5    9    2     0
## Gene_3     1    1    0    0    5    0    1    6    1     4
## Gene_4     2    1    1    1    8    0    3    0    2     0
## Gene_5     0    2    6    1    1   27    9    7    1     3
## Gene_6     0    0    0    0    0    0    6    3    0     0
## Gene_7     1    4    2    1    2    1    2    5    1     0
## Gene_8     1    1    2    3    0    4    2    1    9     3
## Gene_9     0    0    6    3    7    0    0    4    1     1
## Gene_10    2    0    2    0    1    0    1    0    3     0
```

--------------------------------

### Run global Spacelink on selected genes

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_visium_data.html

Runs spacelink on six example genes to identify spatially variable genes. Prints results ordered by ESV_adj descending, showing time, pval, padj, ESV, and ESV_adj.

```R
gene_list <- c("CNP","EIF5B","HOPX","KRT19","PHB","SHISA5")
global_results <- spacelink(
  normalized_counts = normalized_counts[gene_list,],
  spatial_coords = spatial_coords
)
print(global_results[order(global_results$ESV_adj, decreasing=T),c("time","pval","padj","ESV","ESV_adj")])
```

--------------------------------

### Install gaston dependency from CRAN

Source: https://hanbyul-lee.github.io/spacelink/index.html

Installs the gaston dependency package from CRAN. This package may need to be installed manually if automatic installation fails.

```r
install.packages("gaston")
```

--------------------------------

### Output of spacelink_ctSVG for oligodendrocytes

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_cosmx_data.html

Displays the resulting data frame with columns for time, pval, padj, ESV, and ESV_adj for each gene. The values are based on the 6 example genes and will differ from a full-gene-set run.

```R
##           time         pval         padj          ESV     ESV_adj
## GFAP   138.568 0.000000e+00 0.000000e+00 1.812070e-01 0.181206965
## MBP    139.349 0.000000e+00 0.000000e+00 1.069256e-01 0.106925570
## CNP    154.706 0.000000e+00 0.000000e+00 5.807032e-02 0.058070319
## HOPX   139.915 8.518875e-72 1.277831e-71 5.427712e-03 0.005427712
## ANKRD6 148.591 5.112956e-01 5.112956e-01 1.752922e-07 0.000000000
## PSG4   140.014 2.226631e-01 2.671957e-01 2.500886e-03 0.000000000
```

--------------------------------

### Prepare spatial coordinates with distance categories

Source: https://hanbyul-lee.github.io/spacelink/articles/AD_data_analysis.html

Prepares spatial coordinates data frames for both control and AD samples by adding expression, layer, and distance columns, and categorizing distance into 'near' (<=150) or 'far' (>150) relative to neuritic plaques. This setup is required before plotting the spatial gene expression.

```R
colnames(control_spatial_coords) = c("x","y")
control_spatial_coords$expression = control_counts['ENSG00000067225',]
control_spatial_coords$layer = control_spe$SpatialCluster
control_spatial_coords$distance = control_spe$MinDistNP
control_spatial_coords <- control_spatial_coords %>%
  mutate(distance_category = case_when(
    distance <= 150 ~ "near",
    distance > 150  ~ "far",
    TRUE            ~ NA_character_
  ))

colnames(AD_spatial_coords) = c("x","y")
AD_spatial_coords$expression = AD_counts['ENSG00000067225',]
AD_spatial_coords$layer = AD_spe$SpatialCluster
AD_spatial_coords$distance = AD_spe$MinDistNP
AD_spatial_coords <- AD_spatial_coords %>%
  mutate(distance_category = case_when(
    distance <= 150 ~ "near",
    distance > 150  ~ "far",
    TRUE            ~ NA_character_
  ))
```

--------------------------------

### Sparse Matrix Output of Normalized Counts

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_cosmx_data.html

Example output from a dgCMatrix sparse matrix, showing normalized counts for genes (rows) across cells (columns). The output is truncated with suppressed column names and shows values like 0.6931472 and 1.098612.

```R
## 6 x 10 sparse Matrix of class "dgCMatrix"
##   [[ suppressing 10 column names ‘c_1_100_131’, ‘c_1_100_156’, ‘c_1_100_240’ ... ]]
##                                                                                 
## A1BG  .         . .         .         .         .         . . .         1.098612
## A2M   .         . 0.6931472 .         .         .         . . .         .       
## AAAS  .         . .         .         0.6931472 .         . . .         .       
## AAK1  0.6931472 . .         0.6931472 0.6931472 .         . . .         .       
## AAMP  0.6931472 . .         0.6931472 .         .         . . 0.6931472 .       
## AARS1 0.6931472 . .         .         .         0.6931472 . . .         .  
```

--------------------------------

### Interactive Terminal Startup Notice Template

Source: https://hanbyul-lee.github.io/spacelink/LICENSE.html

Use this notice for programs that perform terminal interaction to inform users about the lack of warranty and redistribution rights upon startup.

```text
<program>  Copyright (C) <year>  <name of author>
This program comes with ABSOLUTELY NO WARRANTY; for details type 'show w'.
This is free software, and you are welcome to redistribute it
under certain conditions; type 'show c' for details.
```

--------------------------------

### Load Spacelink Library

Source: https://hanbyul-lee.github.io/spacelink/articles/quick_start.html

Initial step to load the package into the R environment.

```R
library(spacelink)
```

--------------------------------

### Benchmark Spacelink Runtime and Memory

Source: https://hanbyul-lee.github.io/spacelink/articles/runtime_memory_usage.html

This script defines a benchmarking function to measure elapsed time and peak RAM usage for Spacelink and Spacelink-lite across varying spot counts. It uses the peakRAM package and repeats experiments 10 times to calculate median values.

```R
library(spacelink)
library(peakRAM)
library(dplyr)
library(ggplot2)

set.seed(123)
scales  <- c(small = 1000, medium = 5000, large = 100000)
coords_list <- list()
coords_list[['small']]  <- cbind(kronecker(rep(1,25), 1:40),  kronecker(1:25,  rep(1,40)))
coords_list[['medium']] <- cbind(kronecker(rep(1,50), 1:100), kronecker(1:50,  rep(1,100)))
coords_list[['large']]  <- cbind(kronecker(rep(1,250),1:400), kronecker(1:250, rep(1,400)))

n_iter <- 10

run_benchmark <- function(scale_name, method_name, lite_flag, n_iter) {
  n_spots <- scales[[scale_name]]
  coords  <- coords_list[[scale_name]]
  
  times   <- numeric(n_iter)
  peaks   <- numeric(n_iter)
  
  for (i in seq_len(n_iter)) {
    norm_data <- matrix(rnorm(n_spots), nrow = 1, ncol = n_spots)
    
    pr <- peakRAM(
      spacelink(
        normalized_counts = norm_data,
        spatial_coords    = coords,
        lite              = lite_flag
      )
    )
    
    times[i] <- pr$Elapsed_Time_sec
    peaks[i] <- pr$Peak_RAM_Used
  }
  
  data.frame(
    scale       = scale_name,
    n_spots     = n_spots,
    method      = method_name,
    time_s      = median(times),
    peak_mem = median(peaks)
  )
}

benchmark_results <- lapply(names(scales), function(scale_name) {
  cat("Benchmarking scale:", scale_name, "\n")
  
  if (scale_name == "large") {
    run_benchmark(scale_name, "spacelink_lite", lite_flag = TRUE,  n_iter = n_iter)
  } else {
    bind_rows(
      run_benchmark(scale_name, "spacelink_full", lite_flag = FALSE, n_iter = n_iter),
      run_benchmark(scale_name, "spacelink_lite", lite_flag = TRUE,  n_iter = n_iter)
    )
  }
})

bm_df <- bind_rows(benchmark_results) |>
  mutate(
    scale  = factor(scale, levels = names(scales)),
    method = recode_factor(method,
                           "spacelink_full" = "Spacelink",
                           "spacelink_lite" = "Spacelink-lite")
  )
```

--------------------------------

### Load Datasets

Source: https://hanbyul-lee.github.io/spacelink/articles/disease_informativeness_evaluation.html

Load the four required CSV files into the R environment and display the first few rows of the spacelink_results table.

```R
spacelink_results <- read.csv("data/spacelink_results.csv")
pops_scores <- read.csv("data/pops_scores.csv",row.names=1)
pops_trait_category <- read.csv("data/pops_trait_category.csv")
confounders <- read.csv("data/confounders.csv")
head(spacelink_results)
```

--------------------------------

### Load Required R Packages

Source: https://hanbyul-lee.github.io/spacelink/articles/AD_data_analysis.html

Load the necessary libraries for spatial data analysis, visualization, and Spacelink functionality.

```R
library(SpatialExperiment) 
library(Seurat)
library(dplyr)
library(ggplot2)
library(viridis)
library(patchwork)
library(spacelink) 
```

--------------------------------

### spacelink function usage

Source: https://hanbyul-lee.github.io/spacelink/reference/spacelink.html

Shows the function signature and default arguments for spacelink. Use this to understand the required and optional parameters.

```r
spacelink(
  normalized_counts,
  spatial_coords,
  covariates = NULL,
  lite = FALSE,
  grid_size = NULL,
  kernel = "Exponential",
  n_lengthscales = 5,
  n_workers = 1
)
```

--------------------------------

### Run spacelink_ctSVG cell-type-specific analysis

Source: https://hanbyul-lee.github.io/spacelink/articles/overview.html

Runs cell-type-specific SVG hypothesis testing with spacelink_ctSVG. Requires normalized_counts, spatial_coords, cell_type_proportions (spots x cell-types), and focal_cell_type matching a column name. lite = FALSE uses exact computation; n_workers = 1 runs sequentially.

```R
cell_type_results <- spacelink_ctSVG(
  normalized_counts = normalized_expression_matrix,
  spatial_coords = coordinates,
  cell_type_proportions = cell_type_proportions,
  focal_cell_type = cell_type_name,
  lite = FALSE,
  n_workers = 1
)
```

--------------------------------

### Run Global Spacelink with lite=TRUE and grid_size=0.1

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_cosmx_data.html

R code to run Spacelink at the global tissue level on a subset of six genes. Uses lite=TRUE for grid-based spatial binning to reduce memory and computation time, and sets grid_size=0.1 to achieve a meaningful reduction in effective spatial units for densely packed CosMx cells. The results include padj and ESV_adj columns for significance and effect size.

```R
library(spacelink)

gene_list <- c('ANKRD6','CNP','GFAP','HOPX','MBP','PSG4')
global_results <- spacelink(
  normalized_counts=normalized_counts[gene_list,], 
  spatial_coords=spatial_coords, 
  lite=TRUE, 
  grid_size=0.1)
print(global_results[order(global_results$ESV_adj, decreasing=T),c("time","pval","padj","ESV","ESV_adj")])
```

--------------------------------

### Visualize Benchmark Results

Source: https://hanbyul-lee.github.io/spacelink/articles/runtime_memory_usage.html

Generates log-scale plots for runtime and peak memory usage using ggplot2 based on the benchmark data frame.

```R
# --- Plot: Runtime ---
p_time <- ggplot(bm_df, aes(x = n_spots, y = time_s, color = method, group = method)) +
  geom_line(linewidth = 1, na.rm = TRUE) +
  geom_point(size = 3) +
  scale_x_log10(
    breaks = scales,
    labels = paste0(names(scales), "\n(n=", format(scales, big.mark = ","), ")")
  ) +
  scale_y_log10() +
  labs(
    title = "Spacelink Runtime by Scale",
    x     = "Number of Spots",
    y     = "Median Runtime (seconds, log scale)",
    color = "Method"
  ) +
  theme_bw(base_size = 13)

# --- Plot: Peak Memory ---
p_mem <- ggplot(bm_df, aes(x = n_spots, y = peak_mem / (1024^2), color = method, group = method)) +
  geom_line(linewidth = 1, na.rm = TRUE) +
  geom_point(size = 3) +
  scale_x_log10(
    breaks = scales,
    labels = paste0(names(scales), "\n(n=", format(scales, big.mark = ","), ")")
  ) +
  scale_y_log10() +
  labs(
    title = "Spacelink Peak Memory by Scale",
    x     = "Number of Spots",
    y     = "Peak Memory Used (MB, log scale)",
    color = "Method"
  ) +
  theme_bw(base_size = 13)

print(p_time)
print(p_mem)
```

--------------------------------

### Run spacelink global analysis

Source: https://hanbyul-lee.github.io/spacelink/articles/overview.html

Runs global SVG hypothesis testing with spacelink. Requires normalized_counts (genes x spots) and spatial_coords (spots x 2). lite = FALSE uses exact computation; n_workers = 1 runs sequentially. The head() call orders results by ESV_adj descending and shows pval, padj, ESV, and ESV_adj columns.

```R
global_results <- spacelink(
  normalized_counts = normalized_expression_matrix,
  spatial_coords = coordinates,
  lite = FALSE,
  n_workers = 1
)

# View results
head(global_results[order(global_results$ESV_adj, decreasing=T), c("pval","padj","ESV","ESV_adj")])
```

--------------------------------

### Inspect pops_trait_category Output

Source: https://hanbyul-lee.github.io/spacelink/articles/disease_informativeness_evaluation.html

The output of the head(pops_trait_category) command.

```text
##         domain                               var                         var_clean
## 1 Neurological            PASS_ADHD_Demontis2018             ADHD (Demontis, 2018)
## 2 Neurological                    PASS_Alzheimer                         Alzheimer
## 3 Neurological        PASS_Alzheimers_Jansen2019         Alzheimers (Jansen, 2019)
## 4 Neurological                     PASS_Anorexia                          Anorexia
## 5 Neurological                       PASS_Autism                            Autism
## 6 Neurological PASS_BipolarDisorder_Ruderfer2018 Bipolar Disorder (Ruderfer, 2018)
```

--------------------------------

### spacelink_ctSVG function usage

Source: https://hanbyul-lee.github.io/spacelink/reference/spacelink_ctSVG.html

Shows the function signature and all arguments for spacelink_ctSVG. Use this to call the function with required and optional parameters.

```R
spacelink_ctSVG(
  normalized_counts,
  spatial_coords,
  cell_type_proportions,
  focal_cell_type,
  covariates = NULL,
  lite = FALSE,
  grid_size = NULL,
  kernel = "Exponential",
  n_lengthscales = 5,
  c1 = 0.1,
  c2 = 0.1,
  n_workers = 1
)
```

--------------------------------

### Simulate Gaussian process data and compare length-scale errors (R)

Source: https://hanbyul-lee.github.io/spacelink/articles/lengthscale_estimation.html

Simulates normal expression data from a Gaussian process with a single exponential covariance kernel and compares length-scale estimation errors between nnSVG and Spacelink. The true length scale is set to 5 times the minimum distance, and errors are computed as absolute differences divided by the minimum distance.

```R
library(MASS)
library(nnSVG)
library(spacelink)

data("Visium_human_DLPFC")
spatial_coords <- Visium_human_DLPFC$spatial_coords
range_all <- max(apply(spatial_coords, 2, function(col) diff(range(col))))
spatial_coords <- apply(spatial_coords, 2, function(col) (col - min(col)) / range_all)

D <- as.matrix(dist(spatial_coords)); N <- nrow(spatial_coords)
minD <- min(D[D!=0])

phi <- 1/(5*minD)
snr <- 0.25
tau <- 0.5

K <- exp(-phi*D); chol_K <- chol(K)

set.seed(123)
norm_expr <- matrix(1, 50, N) + sqrt(tau)*matrix(rnorm(50*N), 50, N) + sqrt(tau*snr)*mvrnorm(50, rep(0,N), K)

nnsvg_result <- nnSVG(norm_expr, spatial_coords)
spacelink_result <- spacelink(norm_expr, spatial_coords)

nnsvg_lengthscale_error <- abs(1/(nnsvg_result[,"phi"]) - 1/phi)/minD
spacelink_error <- abs(1/apply(spacelink_result, 1, function(x)x[paste0("phi",which.max(x[grepl("sigma",names(x))]))]) - 1/phi)/minD
```

--------------------------------

### Simulate spatial gene expression data in R

Source: https://hanbyul-lee.github.io/spacelink/articles/spacelink_vignette.html

Simulates a spatial transcriptomics dataset with 10 genes across 500 spots. The first six genes are spatially variable with varying length scales and signal-to-noise ratios, while the last four are non-SVGs. The final expression matrix is Poisson-distributed.

```R
set.seed(123)
n_spots <- 500
n_genes <- 10

# Obtain spatial coordinates
coords <- expand.grid(
  x = seq(0, 10, length.out = 25),
  y = seq(0, 10, length.out = 20)
)

expr_data <- matrix(nrow = n_genes, ncol = n_spots)
rownames(expr_data) <- paste0("Gene_", 1:n_genes)
D <- as.matrix(dist(coords))

# Simulate Gene_1, Gene_2, Gene_3
lengthscale = 0.5
snr = c(0.5, 1, 2)
K <- exp(-D/lengthscale); chol_K <- chol(K)
for (i in 1:3) {
  expr_data[i, ] <- rnorm(n_spots) + snr[i] * rnorm(n_spots) %*% chol_K
}

# Simulate Gene_4, Gene_5, Gene_6
lengthscale = 2
snr = c(0.5, 1, 2)
K <- exp(-D/lengthscale); chol_K <- chol(K)
for (i in 1:3) {
  expr_data[i+3, ] <- rnorm(n_spots) + snr[i] * rnorm(n_spots) %*% chol_K
}

# Simulate Gene_7, Gene_8, Gene_9, Gene_10
for (i in 7:10) {
  expr_data[i, ] <- rnorm(n_spots)
}

expr_data <- apply(exp(expr_data), 1:2, function(x)rpois(1,x))
expr_data[,1:10]
```

--------------------------------

### Run Spacelink Analysis

Source: https://hanbyul-lee.github.io/spacelink/articles/spacelink_vignette.html

Executes the spacelink function with normalized counts and spatial coordinates, then displays the adjusted p-values and ESV scores.

```R
library(spacelink)

results <- spacelink(normalized_counts = norm_data, spatial_coords = coords)
results[, c("padj", "ESV")]
```

--------------------------------

### Generate Mock Spatial Data

Source: https://hanbyul-lee.github.io/spacelink/articles/quick_start.html

Creates a synthetic dataset with spatial coordinates, cell-type proportions, and gene expression counts for testing.

```R
# Mock data: 5 SVGs (all cell-type-1 specific ctSVGs) and 5 non-SVGs
set.seed(123)
spatial_coords <- expand.grid(x = seq(0, 10, length.out = 20), y = seq(0, 10, length.out = 20))
cell_type_proportions <- data.frame("Cell_type_1"=c(runif(200,0,0.5),runif(200,0.5,1)),"Cell_type_2"=rep(NA,400)); cell_type_proportions$Cell_type_2 = 1 - cell_type_proportions$Cell_type_1
K <- exp(-as.matrix(dist(spatial_coords))); Sigma <- 2 * chol(K) %*% diag(cell_type_proportions$Cell_type_1); counts <- apply(exp(rbind(matrix(rnorm(2000),5,400) %*% Sigma, matrix(rnorm(2000),5,400))), 1:2, function(x)rpois(1,x)); rownames(counts) <- paste0("Gene_", 1:10)
normalized_counts <- log1p(counts)/log(2)
```

--------------------------------

### Normalize counts with SCTransform

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_visium_data.html

Creates a Seurat object and applies SCTransform with vst.flavor = 'v2' to normalize counts. Extracts the normalized data matrix and displays the first 10 columns.

```R
seurat_obj <- CreateSeuratObject(counts = counts)
seurat_norm = SCTransform(seurat_obj, vst.flavor = "v2", verbose = FALSE)
normalized_counts <- seurat_norm@assays$SCT$data
head(normalized_counts[,1:10])
```

--------------------------------

### Inspect Regression Results

Source: https://hanbyul-lee.github.io/spacelink/articles/disease_informativeness_evaluation.html

Displays the first few rows of the resulting tidy data frame containing the regression results.

```R
head(spacelink_reg_res)
```

--------------------------------

### Inspect spacelink_results Output

Source: https://hanbyul-lee.github.io/spacelink/articles/disease_informativeness_evaluation.html

The output of the head(spacelink_results) command showing the first six rows of the dataset.

```text
##   dataset cell_type  gene          pval          padj        ESV    ESV_adj
## 1  Cortex    Global  A1BG  0.000000e+00  0.000000e+00 0.99800827 0.99800827
## 2  Cortex    Global   A2M 6.187195e-268 2.259640e-267 0.08264960 0.08264960
## 3  Cortex    Global  AAAS 3.417640e-167 9.833155e-167 0.08194720 0.08194720
## 4  Cortex    Global  AAK1  0.000000e+00  0.000000e+00 0.68435610 0.68435610
## 5  Cortex    Global  AAMP  1.571541e-38  2.331317e-38 0.03354861 0.03354861
## 6  Cortex    Global AARS1  0.000000e+00  0.000000e+00 0.25700653 0.25700653
```

--------------------------------

### Simulate ctSVG data with Gaussian process model

Source: https://hanbyul-lee.github.io/spacelink/articles/spacelink_ctSVG_vignette.html

Simulates gene expression data for two cell types with cell-type-specific spatially variable genes (ctSVGs) and non-SVGs. Sets a seed, defines spatial coordinates, cell type proportions, and generates expression counts via a Gaussian process model with Poisson noise. The code creates 10 genes: six cell-type-1 ctSVGs with varying length scales and signal-to-noise ratios, two cell-type-2 ctSVGs, and two non-SVGs.

```R
set.seed(123)
n_spots    <- 500
n_genes    <- 10
n_celltypes <- 2

# Obtain spatial coordinates
coords <- expand.grid(
  x = seq(0, 10, length.out = 25),
  y = seq(0, 10, length.out = 20)
)

# Obtain cell type proportions
cell_type_props <- matrix(nrow = n_spots, ncol = n_celltypes)
colnames(cell_type_props) <- paste0("Cell_type_", 1:n_celltypes)
cell_type_props[1:(n_spots/2), 'Cell_type_1'] <- runif(n_spots / 2, 0, 0.5)
cell_type_props[(n_spots/2+1):n_spots, 'Cell_type_1'] <- runif(n_spots / 2, 0.5, 1)
cell_type_props[, 'Cell_type_2'] <- 1 - cell_type_props[, 'Cell_type_1']

expr_data <- matrix(nrow = n_genes, ncol = n_spots)
rownames(expr_data) <- paste0("Gene_", 1:n_genes)
D <- as.matrix(dist(coords))

# Simulate Gene_1, Gene_2, Gene_3
lengthscale = 0.5
snr = c(1, 2, 5)
K <- exp(-D/lengthscale)
Sigma <- chol(K) %*% diag(cell_type_props[, "Cell_type_1"])
for (i in 1:3) {
  expr_data[i, ] <- rnorm(n_spots) + snr[i] * rnorm(n_spots) %*% Sigma
}

# Simulate Gene_4, Gene_5, Gene_6
lengthscale = 2
snr = c(1, 2, 5)
K <- exp(-D/lengthscale)
Sigma <- chol(K) %*% diag(cell_type_props[, "Cell_type_1"])
for (i in 1:3) {
  expr_data[i+3, ] <- rnorm(n_spots) + snr[i] * rnorm(n_spots) %*% Sigma
}

# Simulate Gene_7, Gene_8
lengthscale = 1
snr = c(1, 5)
K <- exp(-D/lengthscale)
Sigma <- chol(K) %*% diag(cell_type_props[, "Cell_type_2"])
for (i in 1:2) {
  expr_data[i+6, ] <- rnorm(n_spots) + snr[i] * rnorm(n_spots) %*% Sigma
}

# Simulate Gene_9, Gene_10
for (i in 9:10) {
  expr_data[i, ] <- rnorm(n_spots)
}

expr_data <- apply(exp(expr_data), 1:2, function(x)rpois(1,x))
expr_data[,1:10]
```

--------------------------------

### Output of normalized counts head

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_visium_data.html

Shows the first 6 rows and 10 columns of the normalized counts matrix, a sparse matrix of class dgCMatrix.

```R
## 6 x 10 sparse Matrix of class "dgCMatrix"
##   [[ suppressing 10 column names ‘AAACAAGTATCTCCCA-1’, ‘AAACAATCTACTAGCA-1’, ‘AAACACCAATAACTGC-1’ ... ]]                                                    ##                             
## HES4  0.6931472 .        . .         .         .         . .         .         .
## ISG15 .         .        . .         0.6931472 .         . 0.6931472 .         .
## AGRN  0.6931472 1.386294 . .         0.6931472 .         . 0.6931472 .         .
## SDF4  0.6931472 .        . 0.6931472 .         .         . .         0.6931472 .
## ACAP3 0.6931472 .        . .         .         .         . .         0.6931472 .
## DVL1  .         .        . .         .         0.6931472 . .         .         .
```

--------------------------------

### Visualize oligodendrocyte proportions and gene expression with ESV comparison

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_cosmx_data.html

Plots the spatial distribution of oligodendrocytes using a magma color scale, then for the genes with highest and lowest ct-ESV, plots their expression with a log1p-transformed gradient. Titles include global ESV and ct-ESV for comparison.

```R
gene_names <- rownames(cell_type_results)[c(which.max(cell_type_results$ESV_adj),
                                            which.min(cell_type_results$ESV))]

# Visualize oligodendrocyte proportions
df <- as.data.frame(cbind(spatial_coords, expr = cell_type_info[,'oligodendrocyte']))
colnames(df) <- c('X0', 'X1', 'prop')
print(ggplot(df %>% arrange(prop), aes(x = X0, y = X1, color = prop)) +
        geom_point(size = 0.1) +
        ggtitle("Oligodendrocyte") + theme_void() +
        scale_color_viridis(option="magma", direction=-1, limits=c(0,1)) +
        theme(plot.title = element_text(face = "italic", size = 20), legend.position="none"))

# Visualize gene expression (compare ESV vs ct-ESV)
for(gene_name in gene_names){
  df <- as.data.frame(cbind(spatial_coords, expr = counts[gene_name,]))
  colnames(df) <- c('X0', 'X1', 'expr')
  title <- paste0(gene_name,
                ", ESV:", round(global_results[gene_name,'ESV'],4),
                ", ct-ESV:", round(cell_type_results[gene_name,'ESV'],4))
  print(ggplot(df, aes(x = X0, y = X1, color = expr)) +
        geom_point(size = 0.1) +
        ggtitle(title) + theme_void() + 
        scale_color_gradient(low = "grey97", high = "red", trans = "log1p", name = "counts") +
        theme(plot.title = element_text(face = "italic", size = 20)))
}
```

--------------------------------

### Inspect pops_scores

Source: https://hanbyul-lee.github.io/spacelink/articles/disease_informativeness_evaluation.html

Inspect the first five columns of the pops_scores matrix.

```R
head(pops_scores[,1:5])
```

--------------------------------

### Inspect pops_scores Output

Source: https://hanbyul-lee.github.io/spacelink/articles/disease_informativeness_evaluation.html

The output of the head(pops_scores[,1:5]) command.

```text
##         PASS_ADHD_Demontis2018 PASS_Alzheimer PASS_Alzheimers_Jansen2019 PASS_Anorexia PASS_Autism
## OR4F5                  0.00000        0.48576                    0.00000       0.80940     0.23031
## SAMD11                 0.77787        0.00000                    0.14254       0.00000     0.00000
## NOC2L                  0.36135        0.00000                    0.69540       0.00000     0.00000
## KLHL17                 0.00000        0.00000                    0.00000       0.00000     0.00000
## PLEKHN1                0.00000        0.37651                    0.24570       0.00000     0.00000
## PERM1                  0.00000        0.35478                    0.00000       0.29606     0.00000
```

--------------------------------

### Run cell-type-specific Spacelink

Source: https://hanbyul-lee.github.io/spacelink/articles/illustration_on_visium_data.html

Runs spacelink_ctSVG with oligodendrocyte as the focal cell type on the same six genes. Prints results ordered by ESV_adj descending.

```R
cell_type_results <- spacelink_ctSVG(
  normalized_counts = normalized_counts[gene_list,],
  spatial_coords = spatial_coords,
  cell_type_proportions = cell_type_proportions,
  focal_cell_type = "oligodendrocyte"
)

print(cell_type_results[order(cell_type_results$ESV_adj, decreasing=T),])
```

--------------------------------

### Normalize count data

Source: https://hanbyul-lee.github.io/spacelink/articles/spacelink_ctSVG_vignette.html

Apply log normalization to the expression data using log1p and log(2). This step is required before running the spacelink_ctSVG method.

```R
norm_data <- log1p(expr_data)/log(2)
norm_data[,1:10]
```