### matchingVariograms Example Setup

Source: https://jef.works/STcompare/reference/matchingVariograms.html

Example setup for preparing data and parameters before calling the matchingVariograms function.

```R
data(quakes)

#remove duplicated positions
quakes_data <- quakes[!duplicated(cbind(quakes$lat, quakes$long)),]

lat <- quakes_data$lat
long <- quakes_data$long
X <- quakes_data$depth
ids <- 1:length(X)

#number of permutations
B <- 3

#sequence of deltas
delta <- seq(0.1,0.9,0.1)
```

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

### Install STcompare from GitHub

Source: https://jef.works/STcompare/articles/getting-started-with-STcompare.html

Install the STcompare R package from the JEFworks-Lab GitHub repository using remotes::install_github(). Requires the remotes package.

```R
require(remotes)
remotes::install_github('JEFworks-Lab/STcompare')
```

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

### Install STcompare with remotes from GitHub

Source: https://jef.works/STcompare/articles/Install.html

Use the 'remotes' package to install STcompare from GitHub. This is the currently recommended method. The 'build_vignettes' argument is set to FALSE.

```R
require(remotes)
remotes::install_github('JEFworks-Lab/STcompare', build_vignettes = FALSE)
```

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

### spatialSimilarity Usage Example

Source: https://jef.works/STcompare/reference/spatialSimilarity.html

Example demonstrating how to rasterize gene expression data and compute spatial similarity between two experiments.

```R
data(speKidney)
##### Rasterize to get pixels at matched spatial locations #####
rastKidney <- SEraster::rasterizeGeneExpression(speKidney,
               assay_name = 'counts', resolution = 0.2, fun = "mean",
               BPPARAM = BiocParallel::MulticoreParam(), square = FALSE)

s <- spatialSimilarity(list(rastKidney$A, rastKidney$C))
```

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

### Define example genes and shared pixels

Source: https://jef.works/STcompare/articles/brain-MERFISH-10x-visium.html

Select example genes and identify the intersection of spatial coordinates between the MERFISH and Visium datasets to ensure one-to-one comparisons.

```R
geneSig <- svgSigPos[1]
geneNotSig <- nonsvgNotSig[1]

#get the shared pixels between the two samples to show only the pixels that are used for the one-to-one comparisons
sharedPixels <- intersect(rownames(SpatialExperiment::spatialCoords(rast$MERFISH)),
                           rownames(SpatialExperiment::spatialCoords(rast$Visium)))
```

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

### Example: rasterize, compute similarity, and plot pixel classes

Source: https://jef.works/STcompare/reference/pixelClass.html

Complete example workflow. First rasterizes gene expression from speKidney using SEraster::rasterizeGeneExpression, then computes spatial similarity with spatialSimilarity(), and finally calls pixelClass to generate the plot. The output shows a coordinate system message.

```r
data(speKidney)
##### Rasterize to get pixels at matched spatial locations #####
rastKidney <- SEraster::rasterizeGeneExpression(speKidney,
                assay_name = 'counts', resolution = 0.2, fun = "mean",
                BPPARAM = BiocParallel::MulticoreParam(), square = FALSE)
s <- spatialSimilarity(list(rastKidney$A, rastKidney$B))
pixelClass(s, "Gene")
```

```output
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
```

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

### Example: viladomatCorrelation with quakes data

Source: https://jef.works/STcompare/reference/viladomatCorrelation.html

Example using the `quakes` dataset. Removes duplicated positions, prepares the data matrix, sets delta sequence, maxDistPrctile, and nPermutations, then calls `viladomatCorrelation`.

```r
data(quakes)

#remove duplicated positions
quakes_data <- quakes[!duplicated(cbind(quakes$lat, quakes$long)),]

data <- cbind(quakes_data$depth,
              quakes_data$mag,
              quakes_data$lat,
              quakes_data$long)

#sequence of deltas
delta <- seq(0.1,0.9,0.1)

# maximum distance for the variogram set at the 25% percentile of
# the distribution of pairs of distances:
maxDistPrctile <- 0.25

#number of permutations
nPermutations <- 3

resultsPermuteX <- viladomatCorrelation(data, delta, maxDistPrctile, nPermutations)
```

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

### Load STcompare and required R packages

Source: https://jef.works/STcompare/articles/brain-MERFISH-10x-visium.html

R setup block that loads STcompare and supporting libraries for the tutorial, including rhdf5, Matrix, ggplot2, dplyr, patchwork, SEraster, rjson, and scatterbar.

```R
library(STcompare)
library(rhdf5)
library(Matrix)
library(ggplot2)
library(dplyr)
library(patchwork)
library(BiocGenerics)
library(SEraster)
library(rjson)
library(scatterbar)
devtools::load_all()
#> ℹ Loading STcompare
#> Warning in fun(libname, pkgname): no display name and no $DISPLAY environment
#> variable
```

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

### Example: Calculating and Plotting Gene Expression Correlation

Source: https://jef.works/STcompare/reference/plotCorrelationGeneExp.html

This example demonstrates rasterizing gene expression data, calculating spatial correlation coefficients, and generating scatterplots for two different gene expression comparisons.

```R
data(speKidney)

##### Rasterize to get pixels at matched spatial locations #####
rastKidney <- SEraster::rasterizeGeneExpression(speKidney,
               assay_name = 'counts', resolution = 0.2, fun = "mean",
               BPPARAM = BiocParallel::MulticoreParam(), square = FALSE)

##### Use STcompare to calculate Pearson's correlation coefficient #####
rastGexpListAB <- list(A = rastKidney$A, B = rastKidney$B)
rastGexpListAC <- list(A = rastKidney$A, C = rastKidney$C)

negCorrelation <- spatialCorrelationGeneExp(rastGexpListAB)
#> 1: Gene
posCorrelation <- spatialCorrelationGeneExp(rastGexpListAC)
#> 1: Gene

negCorrelation
#>      correlationCoef   pValueNaive pValuePermuteX pValuePermuteY
#> Gene      -0.9472813 5.652003e-136              0              0
#>      deltaStarMedianX deltaStarMedianY   deltaStarX   deltaStarY
#> Gene              0.2              0.2 0.2, 0.2.... 0.2, 0.2....
#>      nullCorrelationsX nullCorrelationsY
#> Gene      -0.16074....      -0.08458....
posCorrelation
#>      correlationCoef   pValueNaive pValuePermuteX pValuePermuteY
#> Gene       0.9431531 1.409195e-133              0              0
#>      deltaStarMedianX deltaStarMedianY   deltaStarX   deltaStarY
#> Gene              0.2              0.2 0.2, 0.4.... 0.2, 0.2....
#>      nullCorrelationsX nullCorrelationsY
#> Gene      -0.00632....      0.052321....

expAB <- plotCorrelationGeneExp(rastGexpListAB, negCorrelation, "Gene")
expAC <- plotCorrelationGeneExp(rastGexpListAC, posCorrelation, "Gene")

expAB
#> Ignoring unknown labels:
#> • fill : "Data"

expAC
#> Ignoring unknown labels:
#> • fill : "Data"
```

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

### Load required libraries for STcompare

Source: https://jef.works/STcompare/articles/getting-started-with-STcompare.html

Loads STcompare and its dependencies (SpatialExperiment, SEraster, ggplot2, patchwork, viridis, reshape2) needed for the examples on this page.

```R
library(STcompare)
library(SpatialExperiment)
library(SEraster)
library(ggplot2)
library(patchwork)
library(viridis)
library(reshape2)
```

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

### Calculate spatial correlation for earthquake data

Source: https://jef.works/STcompare/reference/spatialCorrelation.html

This example demonstrates how to prepare the quakes dataset by removing duplicated positions and calculating the spatial correlation between depth and magnitude.

```R
data(quakes)

#remove duplicated positions
quakes_data <- quakes[!duplicated(cbind(quakes$lat, quakes$long)),]

cor <- spatialCorrelation(X = quakes_data$depth,
                          Y = quakes_data$mag,
                          pos = cbind(quakes_data$lat, quakes_data$long))
cor
#>     correlationCoef  pValueNaive pValuePermuteX pValuePermuteY deltaStarMedianX
#> cor       -0.230266 1.779041e-13           0.03           0.13              0.8
#>     deltaStarMedianY   deltaStarX   deltaStarY nullCorrelationsX
#> cor              0.5 0.9, 0.8.... 0.3, 0.7....      -0.15504....
#>     nullCorrelationsY
#> cor      -0.10144....
```

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

### Calculate spatial correlation within a sample

Source: https://jef.works/STcompare/reference/spatialCorrelationGeneExpWithinSample.html

This example demonstrates how to rasterize gene expression data and then calculate spatial correlation. Note that the provided example code results in a method dispatch error due to the input format.

```R
data(speKidney)

##### Rasterize to get pixels at matched spatial locations #####
rastKidney <- SEraster::rasterizeGeneExpression(speKidney,
               assay_name = 'counts', resolution = 0.2, fun = "mean",
               BPPARAM = BiocParallel::MulticoreParam(), square = FALSE)

sc_within_sample <- spatialCorrelationGeneExpWithinSample(
                      input = rastKidney,
                      assayName = "A"
                      )
#> Error in (function (classes, fdef, mtable) {    methods <- .findInheritedMethods(classes, fdef, mtable)    if (length(methods) == 1L)         return(methods[[1L]])    else if (length(methods) == 0L) {        cnames <- paste0("\"", vapply(classes, as.character,             ""), "\"", collapse = ", ")        stop(gettextf("unable to find an inherited method for function %s for signature %s",             sQuote(fdef@generic), sQuote(cnames)), domain = NA)    }    else stop("Internal error in finding inherited methods; didn't return a unique method",         domain = NA)})(list("list"), new("standardGeneric", .Data = function (x,     ...) standardGeneric("spatialCoords"), generic = structure("spatialCoords", package = "SpatialExperiment"),     package = "SpatialExperiment", group = list(), valueClass = character(0),     signature = "x", default = NULL, skeleton = (function (x,         ...)     stop(gettextf("invalid call in method dispatch to '%s' (no default method)",         "spatialCoords"), domain = NA))(x, ...)), <environment>): unable to find an inherited method for function ‘spatialCoords’ for signature ‘"list"’
sc_within_sample
#> Error: object 'sc_within_sample' not found
```

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

### savePlots usage example

Source: https://jef.works/STcompare/reference/savePlots.html

Complete example: rasterizes gene expression data, computes spatial similarity, then calls savePlots to generate and save plots. Includes the expected console output with coordinate system warnings and a missing values warning.

```r
data(speKidney)
#' ##### Rasterize to get pixels at matched spatial locations #####
rastKidney <- SEraster::rasterizeGeneExpression(speKidney,
                assay_name = 'counts', resolution = 0.2, fun = "mean",
                BPPARAM = BiocParallel::MulticoreParam(), square = FALSE)
s <- spatialSimilarity(list(rastKidney$A, rastKidney$B))
plts <- savePlots("Gene", s, rastKidney)
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
plts
#> $Gene
#> Warning: Removed 15 rows containing missing values or values outside the scale range
#> (`geom_point()`).

#> 
```

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

### Calculate spatial correlation with permutations

Source: https://jef.works/STcompare/reference/spatialCorrelationGeneExpIterPermutations.html

This example demonstrates rasterizing gene expression data from a SpatialExperiment object and then calculating spatial correlation using the spatialCorrelationGeneExpIterPermutations function.

```R
data(speKidney)
if (FALSE) { # \dontrun{
rastKidney <- SEraster::rasterizeGeneExpression(
  speKidney,
  assay_name = "counts",
  resolution = 0.2,
  fun = "mean",
  BPPARAM = BiocParallel::MulticoreParam(),
  square = FALSE
)

rastGexpListAB <- list(A = rastKidney$A, B = rastKidney$B)

corr <- spatialCorrelationGeneExpIterPermutations(
  rastGexpListAB,
  nPermutations = c(100, 1000),
  nThreads = 5
)

negCorrelation
} # }
```

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

### Select example genes for analysis

Source: https://jef.works/STcompare/articles/acute-kidney-injury-10x-visium-rasterized.html

Selects genes with the lowest similarity for positively correlated SVG genes and the smallest difference in percent dissimilarity for negatively correlated SVG genes.

```R
# for the significantly positively correlated svg genes, get gene with lowest similarity (if tied break tie by highest correlation coefficient)
gene1 <- sim_corr_2 %>%
  dplyr::filter(gene %in% svgSigPos) %>% # filter to only significantly positively correlated genes
  dplyr::arrange(desc(correlation)) %>% # if there are ties in similarity, the gene with the highest correlation coefficient will be at the top
  dplyr::arrange(similar) %>% # arrange in ascending order of similarity so that the gene with the lowest similarity is at the top, 
  dplyr::pull(gene) %>%
  head(1)

# for the significantly negatively correlated svg genes, get gene with smallest difference in percent dissimilarity
gene2 <- sim_corr_2 %>%
  dplyr::filter(gene %in% svgSigNeg) %>% # filter to only significantly negatively correlated genes
  dplyr::mutate(diff = abs(percentDissimilarityX - percentDissimilarityY)) %>% # create a new column for the absolute difference between percentDissimilarityX and percentDissimilarityY
  dplyr::arrange(diff) %>% # arrange in ascending order of the difference column so that the gene with the smallest difference is at the top
  dplyr::pull(gene) %>%
  head(1)
```

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

### matchingVariograms

Source: https://jef.works/STcompare/reference/matchingVariograms.html

Creates a permutation of spatial dataset X by starting with a randomly shuffled version of X and smoothing it until the variogram of this permutation matches the variogram of the unmodified X.

```APIDOC
## matchingVariograms

### Description
Function to create a permutation of spatial dataset X by starting with a randomly shuffled version of X and then smoothing it until the variogram of this permutation matches the variogram of the unmodified X.

### Usage
```
matchingVariograms(
  X.randomized,
  long,
  lat,
  delta,
  target_variog,
  prctile,
  ids,
  i,
  seed = 0
)
```

### Arguments
- **X.randomized** (numeric vector) - Required - Random permutation of the values of X across its x,y-coordinates.
- **long** (numeric vector) - Required - y-coordinates of points in dataset X.
- **lat** (numeric vector) - Required - x-coordinates of points in dataset X.
- **delta** (numeric vector) - Required - Given a point in dataset X, the percentage of neighbors which should be within the smoothing kernel. This can be a list or a single numeric. Input to `locfit::lp` as `nn`.
- **target_variog** (variogram) - Required - Variogram of dataset X that will be used as target when doing the matching. See geoR::variog.
- **prctile** (numeric) - Required - The max distance that should be used for evaluating the fit between the variograms from the permuted data and the target variogram. Can be determined from calculating all distances between the points and then calculating a specified quantile. Input to `geoR::variog` as `max.dist`.
- **ids** (integer) - Required - Indices from 1 to the number of points in the dataset X.
- **i** (integer) - Required - The number indicating the ith permutation.
- **seed** (integer) - Optional - Seed for the random number generator. Default `0`.

### Value
The output is returned as a list.
- **residus** (numeric vector) - Sum of squares of the residuals between the target variogram and the permutation variogram for each delta.
- **delta.star.id** (integer) - The index of the delta that minimizes the residual sum of squares.
- **hat.X.delta.star** (numeric vector) - Permutation generated from optimal delta, whose variogram best matches the variogram of X.

### References
Adapted from: Viladomat, Júlia et al. "Assessing the significance of global and local correlations under spatial autocorrelation: a nonparametric approach." Biometrics vol. 70,2 (2014): 409-18. doi:10.1111/biom.12139

### Examples
```
data(quakes)

#remove duplicated positions
quakes_data <- quakes[!duplicated(cbind(quakes$lat, quakes$long)),]

lat <- quakes_data$lat
long <- quakes_data$long
X <- quakes_data$depth
ids <- 1:length(X)

#number of permutations
B <- 3

#sequence of deltas
delta <- seq(0.1,0.9,0.1)
```
```

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

### Create SpatialExperiment objects for Visium and MERFISH

Source: https://jef.works/STcompare/articles/brain-MERFISH-10x-visium.html

Initializes SpatialExperiment objects for both Visium and MERFISH datasets, including one-hot encoding for MERFISH cluster labels.

```R
Visium_SE_ct <- SpatialExperiment::SpatialExperiment(
  assays = list(celltypes = t(ctVisium)),
  spatialCoords = as.matrix(SpatialExperiment::spatialCoords(Visium_SE))

)
# make new spatial experiment for MERFISH with one-hot encoded cluster labels as the assay
one_hot <- model.matrix(~ x - 1, data = data.frame( x = ctMERFISH[intersect(colnames(MERFISH_SE), names(ctMERFISH))]))
colnames(one_hot) <- gsub("x", "", colnames(one_hot))

MERFISH_SE_ct <- SpatialExperiment::SpatialExperiment(
  assays = list(celltypes = t(one_hot)),
  spatialCoords = as.matrix(SpatialExperiment::spatialCoords(MERFISH_SE)[intersect(colnames(MERFISH_SE), names(ctMERFISH)), ])
)
```

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

### Download and load Visium cell type proportions

Source: https://jef.works/STcompare/articles/brain-MERFISH-10x-visium.html

Downloads deconvolved cell type proportions for Visium spots from Zenodo and formats the column names to match the cell type matrix.

```R
#load ctVisium: matrix of deconvolved cell type proportions for Visium spots, rows are the spot barcodes and columns are cell type names
zenodo_url <- "https://zenodo.org/records/19582556/files/STalign_Visium_cell_type_annotations.csv.gz?download=1"
temp <- tempfile(fileext = ".csv.gz")
download.file(zenodo_url, destfile = temp, mode = "wb")

ctVisium <- read.csv(gzfile(temp), row.names = 1, stringsAsFactors = FALSE)
colnames(ctVisium) <- rownames(cmat)
ctVisium[1:5, 1:5]
#>                         CT-1 CT-2 CT-3      CT-4      CT-5
#> AAACAGAGCGACTCCT-1 0.0000000    0    0 0.0000000 0.2784436
#> AAACCCGAACGAAATC-1 0.0535788    0    0 0.2691073 0.0000000
#> AAACCGGGTAGGTACC-1 0.0000000    0    0 0.0000000 0.0000000
#> AAACCGTTCGTCCAGG-1 0.1393211    0    0 0.0000000 0.0000000
#> AAACGAAGAACATACC-1 0.0000000    0    0 0.0000000 0.1962521
```

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

### Download and Load Visium Data

Source: https://jef.works/STcompare/articles/brain-MERFISH-10x-visium.html

Downloads the Visium FFPE mouse brain dataset from 10X Genomics, extracts the tarball, and reads the count matrix, barcodes, and feature names into R.

```R
# The Visium can be found on 10X Genomics:https://www.10xgenomics.com/datasets/adult-mouse-brain-ffpe-1-standard-1-3-0 

# Visium counts  ###############################################################

url <- "https://cf.10xgenomics.com/samples/spatial-exp/1.3.0/Visium_FFPE_Mouse_Brain/Visium_FFPE_Mouse_Brain_filtered_feature_bc_matrix.tar.gz"

temp_tar <- tempfile(fileext = ".tar.gz")
download.file(url, destfile = temp_tar, mode = "wb")

temp_dir <- tempdir()
untar(temp_tar, exdir = temp_dir)

matrix_path <- file.path(temp_dir, "filtered_feature_bc_matrix", "matrix.mtx.gz")
barcodes_path <- file.path(temp_dir, "filtered_feature_bc_matrix", "barcodes.tsv.gz")
features_path <- file.path(temp_dir, "filtered_feature_bc_matrix", "features.tsv.gz")

Visium.gexp.raw <- Matrix::readMM(matrix_path)
Visium.barcodes <- read.csv(barcodes_path, header=FALSE)
Visium.genenames <- read.csv(features_path, sep='\t', header=FALSE)


dim(Visium.gexp.raw)
#> [1] 19465  2264
length(Visium.barcodes[,1])
#> [1] 2264
length(Visium.genenames[,2])
#> [1] 19465
rownames(Visium.gexp.raw) <- Visium.genenames[,2]
colnames(Visium.gexp.raw) <- Visium.barcodes[,1]
Visium.gexp.raw[1:5, 1:5]
#> 5 x 5 sparse Matrix of class "dgTMatrix"
#>        AAACAGAGCGACTCCT-1 AAACCCGAACGAAATC-1 AAACCGGGTAGGTACC-1
#> Xkr4                    1                  .                  1
#> Rp1                     .                  .                  .
#> Sox17                   .                  .                  .
#> Lypla1                  1                  .                  .
#> Tcea1                   1                  .                  1
#>        AAACCGTTCGTCCAGG-1 AAACGAAGAACATACC-1
#> Xkr4                    .                  .
#> Rp1                     .                  .
#> Sox17                   .                  .
#> Lypla1                  .                  .
#> Tcea1                   1                  1
```

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

### Create SpatialExperiment Objects

Source: https://jef.works/STcompare/articles/brain-MERFISH-10x-visium.html

Constructs SpatialExperiment objects for MERFISH and Visium datasets, including counts and normalized assays.

```R
MERFISH_SE <- SpatialExperiment::SpatialExperiment(
  assays = list(counts = MERFISH.gexp[genes.have, rownames(MERFISH.pos.good)],
                libnorm = MERINGUE::normalizeCounts(MERFISH.gexp[genes.have, rownames(MERFISH.pos.good)], log = FALSE)),
  spatialCoords = as.matrix(MERFISH.pos.good),
  
)
#> Converting to sparse matrix ...
#> Normalizing matrix with 40697 cells and 466 genes.
#> normFactor not provided. Normalizing by library size.
#> Using depthScale 1e+06

Visium_SE <- SpatialExperiment::SpatialExperiment(
  assays = list(counts = Visium.gexp.raw[genes.have,],
                libnorm = MERINGUE::normalizeCounts(Visium.gexp.raw[genes.have,], log=FALSE)),
  spatialCoords = as.matrix(Visium.pos.good)

)
#> Converting to sparse matrix ...
#> Normalizing matrix with 2264 cells and 466 genes.
#> normFactor not provided. Normalizing by library size.
#> Using depthScale 1e+06
```

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

### Load and inspect simRanPatternRasts dataset

Source: https://jef.works/STcompare/reference/simRanPatternRasts.html

Loads the simRanPatternRasts dataset and displays the first element, a SpatialExperiment object. The output shows the object structure with 1 row and 280 columns, an assay named 'pixelval', and spatial coordinates. The subsequent attempt to access the assay fails because the 'assays' function is not found, likely due to missing package loading.

```r
data(simRanPatternRasts)
simRanPatternRasts[[1]]
#> class: SpatialExperiment 
#> dim: 1 280 
#> metadata(0):
#> assays(1): pixelval
#> rownames(1): 1
#> rowData names(0):
#> colnames(280): pixel151 pixel152 ... pixel588 pixel589
#> colData names(6): num_cell cellID_list ... geometry sample_id
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> spatialCoords names(2) : x y
#> imgData names(0):
assays(simRanPatternRasts[[1]])$pixelval[1, 1:5]
#> Error in assays(simRanPatternRasts[[1]]): could not find function "assays"
```

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

### Visualize spatial gene expression and correlation

Source: https://jef.works/STcompare/articles/acute-kidney-injury-10x-visium-rasterized.html

Plots spatial gene expression patterns for shared pixels and generates correlation plots for the selected example genes.

```R
#get the shared pixels between the two samples to show only the pixels that are used for the one-to-one comparisons
sharedPixels <- intersect(rownames(SpatialExperiment::spatialCoords(rast$AKI_ctrl)),
                           rownames(SpatialExperiment::spatialCoords(rast$AKI_aki)))

#use SEraster::plotRaster to plot the spatial pattern of gene expression for each gene in each sample for only the shared pixels
gene1_plt_ctrl  <- SEraster::plotRaster(rast$AKI_ctrl[gene1, sharedPixels], assay_name = "CPM", name = "CPM expression")
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
gene1_plt_aki <- SEraster::plotRaster(rast$AKI_aki[gene1, sharedPixels], assay_name = "CPM", name = "CPM expression")
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.

gene2_plt_ctrl  <- SEraster::plotRaster(rast$AKI_ctrl[gene2, sharedPixels], assay_name = "CPM", name = "CPM expression")
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.
gene2_plt_aki <- SEraster::plotRaster(rast$AKI_aki[gene2, sharedPixels], assay_name = "CPM", name = "CPM expression")
#> Coordinate system already present.
#> ℹ Adding new coordinate system, which will replace the existing one.


# use the plotCorrelationGeneExp function to generate the correlation plots 

gene1_corr <- plotCorrelationGeneExp(
  speList = input, 
  spatialCorrelation = kidneyCorrelation, 
  geneName = gene1, 
  assayName = "CPM"
  )
gene2_corr <- plotCorrelationGeneExp(
  speList = input, 
  spatialCorrelation = kidneyCorrelation, 
  geneName = gene2, 
  assayName = "CPM"
  )

# use patchwork to put three plots together for each gene
gene1_plt_ctrl | gene1_plt_aki | gene1_corr
#> Ignoring unknown labels:
#> • fill : "Data"
```

```R
gene2_plt_ctrl | gene2_plt_aki | gene2_corr
#> Ignoring unknown labels:
#> • fill : "Data"
```

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

### Download and process control counts

Source: https://jef.works/STcompare/articles/acute-kidney-injury-10x-visium-rasterized.html

Downloads control count data from Zenodo, reads the HDF5 matrix, and converts it to a dgCMatrix.

```R
# Control counts ###############################################################

zenodo_url <- "https://zenodo.org/records/19074288/files/NL3_filtered_feature_bc_matrix.h5?download=1"
temp_h5 <- tempfile(fileext = ".h5")
download.file(zenodo_url, temp_h5, mode = "wb")

# Read the CSC-encoded matrix components from the HDF5 file
file_barcodes <- as.character(rhdf5::h5read(temp_h5, "matrix/barcodes"))
barcodes_matrix <- rhdf5::h5read(temp_h5, "matrix") # stores the sparse matrix components 

# convert CSC encoding into an R dgCMatrix
control_counts <- Matrix::sparseMatrix(
  dims = barcodes_matrix$shape,
  i = as.numeric(barcodes_matrix$indices),
  p = as.numeric(barcodes_matrix$indptr),
  x = as.numeric(barcodes_matrix$data),
  index1 = FALSE
)
colnames(control_counts) <- file_barcodes
rownames(control_counts) <- barcodes_matrix[["features"]]$name

unlink(temp_h5)
control_counts[1:5, 1:6]
```

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

### Create SpatialExperiment objects

Source: https://jef.works/STcompare/articles/acute-kidney-injury-10x-visium-rasterized.html

Initializes SpatialExperiment objects for both control and AKI datasets using their respective counts matrices and aligned spatial coordinates.

```R
AKI_ctrl_SE <- SpatialExperiment::SpatialExperiment(
  assays = list(counts = control_counts),
  spatialCoords = as.matrix(ctrl_pos_rot[,1:2]),
)

AKI_aki_SE <- SpatialExperiment::SpatialExperiment(
  assays = list(counts = aki_counts),
  spatialCoords = as.matrix(aki_pos_aligned[,1:2])
)
```

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

### Example: linearRegression on rasterized kidney data

Source: https://jef.works/STcompare/reference/linearRegression.html

Rasterizes gene expression from speKidney, computes spatial similarity, and generates the regression plot for 'Gene'. The output includes a warning about removed rows.

```r
data(speKidney)
##### Rasterize to get pixels at matched spatial locations #####
rastKidney <- SEraster::rasterizeGeneExpression(speKidney,
               assay_name = 'counts', resolution = 0.2, fun = "mean",
               BPPARAM = BiocParallel::MulticoreParam(), square = FALSE)

s <- spatialSimilarity(list(rastKidney$A, rastKidney$C))
linearRegression(s, "Gene")
#> Warning: Removed 14 rows containing missing values or values outside the scale range
#> (`geom_point()`).
```

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

### Load required libraries for STcompare analysis

Source: https://jef.works/STcompare/articles/acute-kidney-injury-10x-visium-rasterized.html

Loads STcompare, SpatialExperiment, MERINGUE, rhdf5, ggplot2, dplyr, patchwork, Matrix, BiocGenerics, and devtools. The output shows package loading messages and masking warnings.

```R
library(STcompare)
library(SpatialExperiment)
#> Loading required package: SingleCellExperiment
#> Loading required package: SummarizedExperiment
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#> 
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#> 
#>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#>     colWeightedMeans, colWeightedMedians, colWeightedSds,
#>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#>     rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#>     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
#>     table, tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#> 
#>     findMatches
#> The following objects are masked from 'package:base':
#> 
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: GenomeInfoDb
#> Loading required package: Biobase
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#> 
#>     rowMedians
#> The following objects are masked from 'package:matrixStats':
#> 
#>     anyMissing, rowMedians
library(MERINGUE)
library(rhdf5)
library(ggplot2)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:Biobase':
#> 
#>     combine
#> The following objects are masked from 'package:GenomicRanges':
#> 
#>     intersect, setdiff, union
#> The following object is masked from 'package:GenomeInfoDb':
#> 
#>     intersect
#> The following objects are masked from 'package:IRanges':
#> 
#>     collapse, desc, intersect, setdiff, slice, union
#> The following objects are masked from 'package:S4Vectors':
#> 
#>     first, intersect, rename, setdiff, setequal, union
#> The following objects are masked from 'package:BiocGenerics':
#> 
#>     combine, intersect, setdiff, union
#> The following object is masked from 'package:matrixStats':
#> 
#>     count
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(patchwork)
library(Matrix)
#> 
#> Attaching package: 'Matrix'
#> The following object is masked from 'package:S4Vectors':
#> 
#>     expand
library(BiocGenerics)
devtools::load_all()
#> ℹ Loading STcompare
#> Warning in fun(libname, pkgname): no display name and no $DISPLAY environment
#> variable
```

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

### R example: spatialCorrelationGeneExp with speKidney

Source: https://jef.works/STcompare/reference/spatialCorrelationGeneExp.html

Rasterizes gene expression from the speKidney dataset using SEraster::rasterizeGeneExpression, then computes Pearson correlations with spatialCorrelationGeneExp for two gene pairs. The output shows correlation coefficients, p-values, and delta star summaries.

```R
data(speKidney)

##### Rasterize to get pixels at matched spatial locations #####
rastKidney <- SEraster::rasterizeGeneExpression(speKidney,
               assay_name = 'counts', resolution = 0.2, fun = "mean",
               BPPARAM = BiocParallel::MulticoreParam(), square = FALSE)

##### Use STcompare to calculate Pearson's correlation coefficient #####
rastGexpListAB <- list(A = rastKidney$A, B = rastKidney$B)
rastGexpListAC <- list(A = rastKidney$A, C = rastKidney$C)

negCorrelation <- spatialCorrelationGeneExp(rastGexpListAB, nThreads = 5)
#> 1: Gene
posCorrelation <- spatialCorrelationGeneExp(rastGexpListAC, nThreads = 5)
#> 1: Gene

negCorrelation
#>      correlationCoef   pValueNaive pValuePermuteX pValuePermuteY
#> Gene      -0.9472813 5.652003e-136              0              0
#>      deltaStarMedianX deltaStarMedianY   deltaStarX   deltaStarY
#> Gene              0.2              0.2 0.2, 0.2.... 0.2, 0.2....
#>      nullCorrelationsX nullCorrelationsY
#> Gene      -0.16074....      -0.08458....
posCorrelation
#>      correlationCoef   pValueNaive pValuePermuteX pValuePermuteY
#> Gene       0.9431531 1.409195e-133              0              0
#>      deltaStarMedianX deltaStarMedianY   deltaStarX   deltaStarY
#> Gene              0.2              0.2 0.2, 0.4.... 0.2, 0.2....
#>      nullCorrelationsX nullCorrelationsY
#> Gene      -0.00632....      0.052321....
```

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

### Download and inspect 10x Visium HDF5 data

Source: https://jef.works/STcompare/articles/acute-kidney-injury-10x-visium-rasterized.html

Downloads the AKI dataset from Zenodo and lists the internal structure of the HDF5 file using rhdf5.

```R
# The Visium mouse kidney datasets can be found Zenodo: https://doi.org/10.5281/zenodo.17676991
# IL3 - ischemia acute kidney injury dataset 
# NL3 - control dataset 

# AKI counts ###################################################################
zenodo_url <- "https://zenodo.org/records/19074288/files/IL3_filtered_feature_bc_matrix.h5?download=1"
temp_h5 <- tempfile(fileext = ".h5")
download.file(zenodo_url, temp_h5, mode = "wb")

# this is a list of the internal groups and datasets stored in the HDF5 file 
# we are interested in the "matrix/barcodes" dataset 
rhdf5::h5ls(temp_h5)
#>               group          name       otype  dclass      dim
#> 0                 /        matrix   H5I_GROUP                 
#> 1           /matrix      barcodes H5I_DATASET  STRING     2970
#> 2           /matrix          data H5I_DATASET INTEGER 10857645
#> 3           /matrix      features   H5I_GROUP                 
#> 4  /matrix/features _all_tag_keys H5I_DATASET  STRING        1
#> 5  /matrix/features  feature_type H5I_DATASET  STRING    32285
#> 6  /matrix/features        genome H5I_DATASET  STRING    32285
#> 7  /matrix/features            id H5I_DATASET  STRING    32285
#> 8  /matrix/features          name H5I_DATASET  STRING    32285
#> 9           /matrix       indices H5I_DATASET INTEGER 10857645
#> 10          /matrix        indptr H5I_DATASET INTEGER     2971
#> 11          /matrix         shape H5I_DATASET INTEGER        2
```

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

### Download and Parse Visium Spatial Data

Source: https://jef.works/STcompare/articles/brain-MERFISH-10x-visium.html

Downloads the Visium FFPE Mouse Brain dataset, extracts the spatial files, and reads the tissue positions into a data frame.

```R
url <- "https://cf.10xgenomics.com/samples/spatial-exp/1.3.0/Visium_FFPE_Mouse_Brain/Visium_FFPE_Mouse_Brain_spatial.tar.gz"

temp_tar <- tempfile(fileext = ".tar.gz")
download.file(url, destfile = temp_tar, mode = "wb")

temp_dir <- tempdir()
untar(temp_tar, exdir = temp_dir)

pos_path <- file.path(temp_dir, "spatial", "tissue_positions_list.csv")
scale_path <- file.path(temp_dir, "spatial", "scalefactors_json.json")

Visium.pos.info <- read.csv(pos_path, header = FALSE, row.names = 1)
head(Visium.pos.info)
#>                    V2 V3 V4   V5    V6
#> ACGCCTGACACGCGCT-1  0  0  0 5247 21461
#> TACCGATCCAACACTT-1  0  1  1 5391 21208
#> ATTAAAGCGGACGAGC-1  0  0  2 5537 21460
#> GATAAGGGACGATTAG-1  0  1  3 5682 21207
#> GTGCAAATCACCAATA-1  0  0  4 5827 21460
#> TGTTGGCTGGCGGAAG-1  0  1  5 5972 21207
```

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

### Load tissue positions from Zenodo

Source: https://jef.works/STcompare/articles/acute-kidney-injury-10x-visium-rasterized.html

Downloads and reads the tissue position CSV files for AKI and control samples into R data frames.

```R
zenodo_url <- "https://zenodo.org/records/19074288/files/IL3_tissue_positions.csv?download=1"
aki_pos <- read.csv(zenodo_url, header = TRUE, stringsAsFactors = FALSE)

zenodo_url <- "https://zenodo.org/records/19074288/files/NL3_tissue_positions.csv?download=1"
ctrl_pos <- read.csv(zenodo_url, header = TRUE, stringsAsFactors = FALSE)

head(ctrl_pos)
#>              barcode array_row array_col
#> 1 AAACAAGTATCTCCCA-1       102        50
#> 2 AAACAGAGCGACTCCT-1        94        14
#> 3 AAACAGCTTTCAGAAG-1         9        43
#> 4 AAACAGGGTCTATATT-1        13        47
#> 5 AAACATTTCCCGGATT-1        97        61
#> 6 AAACCCGAACGAAATC-1       115        45
head(aki_pos)
#>              barcode array_row array_col
#> 1 AAACAAGTATCTCCCA-1       102        50
#> 2 AAACAGAGCGACTCCT-1        94        14
#> 3 AAACAGCTTTCAGAAG-1         9        43
#> 4 AAACAGGGTCTATATT-1        13        47
#> 5 AAACATTTCCCGGATT-1        97        61
#> 6 AAACCCGAACGAAATC-1       115        45
```

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

### Load MERFISH cell type annotations

Source: https://jef.works/STcompare/articles/brain-MERFISH-10x-visium.html

Downloads the cell type annotations for MERFISH cells.

```R
#load ctMERFISH: named vector of cell type annotations for MERFISH cells, names are cell ids
zenodo_url <- "https://zenodo.org/records/19582556/files/STalign_S2R3_cell_type_annotations.csv.gz?download=1"
temp <- tempfile(fileext = ".csv.gz")
download.file(zenodo_url, destfile = temp, mode = "wb")
```