MERINGUE

https://jef.works/MERINGUE/
Admin
MERINGUE is an R package for characterizing spatial gene expression heterogeneity and cell-cell...
Tokens:15,653
Snippets:138
Trust Score:9
Update:1 month ago
Show doc for...
Tokens:
Raw
### Install MERINGUE via remotes

Source: https://jef.works/MERINGUE

Recommended installation method for the MERINGUE package using the remotes R package.

```R
# install.packages(remotes)
require(remotes)
remotes::install_github('JEFworks-Lab/MERingue', build_vignettes = TRUE)
```

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

### Visualize gene expression in space

Source: https://jef.works/MERINGUE/spatial_clustering

Plots the expression of two example genes (first and 50th) in spatial coordinates to show that they mark either only the blue group or both the red and green groups.

```R
par(mfrow=c(1,2))
plotEmbedding(pos, colors=scale(pcs[,1])[,1],
              main='Cell-Types Gene Expression', xlab='Spatial X', ylab='Spatial Y')
plotEmbedding(pos, colors=scale(pcs[,50])[,1],
              main='Cell-Types Gene Expression', xlab='Spatial X', ylab='Spatial Y')
```

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

### Load MERINGUE and BCL Dataset

Source: https://jef.works/MERINGUE/BCL_analysis

Load the MERINGUE library and the BCL dataset, then inspect the structure of the positions and counts.

```R
suppressMessages(library(MERINGUE))

data(BCL) # object contains counts and positions for all 4 sections. Will only need sections 1-3
head(BCL$pos)
BCL$counts[1:5,1:5]
```

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

### Sample Genes for Demonstration

Source: https://jef.works/MERINGUE/BCL_analysis

Subsets the expression matrix to 1000 genes for demonstration purposes.

```R
# Sample 1000 genes for demonstrative purposes only
set.seed(0)
test <- sample(rownames(mat), 1000)
mat <- mat[test,]
```

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

### Load MERINGUE package

Source: https://jef.works/MERINGUE/simulation

Loads the MERINGUE package while suppressing startup messages. Required before using functions like getSpatialNeighbors, spatialCrossCor, and spatialCrossCorTest.

```R
suppressMessages(library(MERINGUE))
```

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

### Simulate three spatially distinct cell groups

Source: https://jef.works/MERINGUE/spatial_clustering

Loads the MERINGUE library and sets up plotting parameters. Simulates 300 cells in three spatially distinct groups using rnorm and assigns group labels.

```R
library(MERINGUE)
par(mfrow=c(1,1), mar=rep(4,4))

## simulate 3 spatially but 2 transcriptionally distinct groups
N <- 300
## Three spatially distinct groups
pos1 <- cbind(rnorm(N/3), rnorm(N/3))
pos2 <- cbind(rnorm(N/3, 10), rnorm(N/3))
pos3 <- cbind(rnorm(N/3, 10), rnorm(N/3, 10))
pos <- rbind(rbind(pos1, pos2), pos3)
group <- c(rep(1, N/3), rep(2, N/3), rep(3, N/3))
names(group) <- rownames(pos) <- paste0('cell', 1:N)
plotEmbedding(pos, groups=group, main='Cell-Types in Space', xlab='Spatial X', ylab='Spatial Y')
```

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

### Visualize 3D network

Source: https://jef.works/MERINGUE/drosophila_3D_analysis

Render the 3D neighbor network using RGL and save the output to a file.

```R
## plot
require(rgl)
MERINGUE:::plotNetwork3D(pos, N, size=1)
## save plot
rgl.snapshot(paste0(dir(), 'neighbors_3D-flat.png'))
```

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

### Visualize marker genes as a heatmap

Source: https://jef.works/MERINGUE/mOB_analysis

Prepares gene groups and column/row colors to visualize marker genes in a heatmap.

```R
dg.genes <- unlist(dg.sig)
ggroup <- unlist(lapply(1:length(dg.sig), function(i) { 
  rep(names(dg.sig)[i], length(dg.sig[[i]]))
}))
names(ggroup) <- dg.genes
ggroup <- factor(ggroup)

# Plot
ccol <- rainbow(length(levels(annot)))[annot]
names(ccol) <- names(annot) # column colors
gcol <- rainbow(length(levels(ggroup)), v=0.5)[ggroup]
names(gcol) <- names(ggroup) # row colors

m <- as.matrix(mat[dg.genes, names(sort(annot))])
m <- winsorize(t(scale(t(m))))
heatmap(m, scale="none", 
          Colv=NA, Rowv=NA, labRow=NA, labCol=NA,
          ColSideColors=ccol[colnames(m)],
          RowSideColors=gcol[rownames(m)],
          col=colorRampPalette(c('blue', 'white', 'red'))(100)
)
```

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

### Load mOB data and extract positions and counts

Source: https://jef.works/MERINGUE/mOB_analysis

Loads the MERINGUE package and the mOB dataset, extracting the spot positions (pos) and counts matrix (cd) for downstream analysis.

```R
suppressMessages(library(MERINGUE))
data(mOB)
pos <- mOB$pos
cd <- mOB$counts
```

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

### Simulate Gene1 and Gene2 expression gradients

Source: https://jef.works/MERINGUE/simulation

Simulates expression for Gene1 and Gene2 as jittered versions of a sorted random gradient, then plots both on the simulated spatial layout. Uses set.seed(0) and set.seed(1) for reproducibility.

```R
par(mfrow=c(1,2), mar=rep(5,4))
gexp0 <- sort(abs(rnorm(N)))

set.seed(0)
gexp1 <- jitter(gexp0, amount = 0.5)
names(gexp1) <- rownames(pos)
plotEmbedding(posw, col=gexp1,  
              main='Gene1')

set.seed(1)
gexp2 <- jitter(gexp0, amount = 0.5)
names(gexp2) <- rownames(pos)
plotEmbedding(posw, col=gexp2,  
              main='Gene2')
```

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

### Simulate transcriptionally identical groups with gene expression matrix

Source: https://jef.works/MERINGUE/spatial_clustering

Creates a gene expression matrix where two of the three spatial groups are transcriptionally identical. 30 genes are upregulated in the blue group and 60 genes are upregulated in the green and red groups. Visualizes the expression matrix with a heatmap.

```R
M <- 30
## But two are transcriptionally identical
pcs12 <- matrix(rnorm(N*2/3*M), N*2/3, M)
pcs3 <- matrix(rnorm(N*1/3*M, 10), N*1/3, M)
pcs <- rbind(pcs12, pcs3)
pcs <- cbind(pcs, abs(10-pcs))
colnames(pcs) <- paste0('gene', 1:ncol(pcs))
rownames(pcs) <- rownames(pos)
heatmap(pcs, scale="none", Rowv = NA, Colv=NA,
        RowSideColors=MERINGUE:::fac2col(group),
        col=colorRampPalette(c('blue', 'white', 'red'))(100))
```

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

### Perform spatially informed clustering and visualize results

Source: https://jef.works/MERINGUE/spatial_clustering

Runs getSpatiallyInformedClusters with k=50, alpha=0, beta=0, displays the cluster table, and plots the clusters in spatial and UMAP embeddings. The output shows three clusters of 100 cells each.

```R
com2 <- getSpatiallyInformedClusters(pcs, W=W, k=50, alpha=0, beta=0)
table(com2)
plotEmbedding(pos, groups=com2, main='Spatially Aware Transcriptional Clusters', xlab='Spatial X', ylab='Spatial Y')
```

```R
plotEmbedding(emb, groups=com2, main='Spatially Aware Transcriptional Clusters', xlab='UMAP X', ylab='UMAP Y')
```

```R
## com2
##   1   2   3 
## 100 100 100
```

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

### Filter and normalize counts with cleanCounts and normalizeCounts

Source: https://jef.works/MERINGUE/mOB_analysis

Filters out poor spots and genes (fewer than 100 counts) using cleanCounts, then normalizes to counts per million (CPM) with normalizeCounts. The output shows the filtering results and normalization parameters.

```R
# Remove poor datasets and genes
counts <- cleanCounts(counts = cd, 
                      min.reads = 100, 
                      min.lib.size = 100, 
                      plot=TRUE,
                      verbose=TRUE)
```

```R
pos <- pos[colnames(counts),]

# CPM normalize
mat <- normalizeCounts(counts = counts, 
                       log=FALSE,
                       verbose=TRUE)
```

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

### Simulate and plot Gene5 and Gene6 spatial gradients

Source: https://jef.works/MERINGUE/simulation

Simulates two genes with spatial gradients (Gene5 left-to-right, Gene6 top-to-bottom) using jitter on coordinates, then plots their embeddings. Uses set.seed(0) and set.seed(1) for reproducibility.

```R
par(mfrow=c(1,2), mar=rep(2,4))
set.seed(0)
gexp5 <- jitter(pos[,1], amount = 0.5)
names(gexp5) <- rownames(pos)
plotEmbedding(posw, col=gexp5,  
              main='Gene5')
```

```R
set.seed(1)
gexp6 <- jitter(pos[,2], amount = 0.5)
names(gexp6) <- rownames(pos)
plotEmbedding(posw, col=gexp6,  
              main='Gene6')
```

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

### Compare spatially aware and unaware clustering and run differential expression

Source: https://jef.works/MERINGUE/spatial_clustering

Compares the spatially aware clusters (com2) with the unaware clusters (com) using table, then runs differential expression on the cells in clusters 1 and 2. The output shows no significant differentially expressed genes (0 in each cluster).

```R
## what's different between our spatially aware and unaware clustering
print(table(com, com2))
## try to identify subtly differentially expressed genes
vi <- com2 %in% c(1,2)
cols <- factor(com2[vi])
diffgexp <- getDifferentialGenes(t(pcs[vi,]), cols)
```

```R
## Warning in getDifferentialGenes(t(pcs[vi, ]), cols): NaNs produced
```

```R
sigdiffgexp <- sapply(diffgexp, function(x) sum(x$p.adj < 0.05))
sigdiffgexp
```

```R
##    com2
## com   1   2   3
##   1 100 100   0
##   2   0   0 100
## [1] "Running differential expression with 2 clusters ... "
## [1] "Summarizing results ... "
## 1 2 
## 0 0
```

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

### Simulate Gene3 expression pattern

Source: https://jef.works/MERINGUE/simulation

Simulates expression for Gene3 based on a subset of cells (even x-coordinates) along a gradient, then plots it on the spatial layout. This is the first part of Scenario 2, which contrasts with the general cross-correlation result.

```R
par(mfrow=c(1,2), mar=rep(5,4))
num <- pos[,1]

vi <- (num %% 2) == 0
exp3 <- as.numeric(vi)
exp3 <- exp3 * sort(abs(rnorm(length(exp3)))+1)
names(exp3) <- rownames(pos)
plotEmbedding(posw, col=exp3,  
              main='Gene3')
```

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

### Perform standard transcriptional clustering with getClusters

Source: https://jef.works/MERINGUE/spatial_clustering

Uses MERINGUE's getClusters function with k=50 to identify transcriptionally distinct clusters, then visualizes the clusters in both spatial and UMAP coordinates.

```R
## First, with standard spatially-unaware cluster detection
par(mfrow=c(1,2))
com <- getClusters(pcs, k=50)
plotEmbedding(pos, groups=com, main='Transcriptional Clusters', xlab='Spatial X', ylab='Spatial Y')
plotEmbedding(emb, groups=com, main='Transcriptional Clusters', xlab='UMAP X', ylab='UMAP Y')
```

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

### Simulate cells in space

Source: https://jef.works/MERINGUE/simulation

Simulates a 15x15 grid of cell positions, adds jitter and warping, then plots the resulting spatial layout with plotEmbedding. The variable posw is used in later spatial analyses.

```R
# 15x15 grid of cells
N <- 15^2
pos <- t(combn(c(1:sqrt(N), rev(1:sqrt(N))), 2))
pos <- unique(pos)
rownames(pos) <- paste0('cell', 1:N)
colnames(pos) <- c('x', 'y')
# jitter
posj <- jitter(pos, amount = 0.5)
# induce warping
posw <- 1.1^posj
# plot
par(mfrow=c(1,1), mar=rep(5,4))
plotEmbedding(posw,  main='Simulated Cells in Space')
```

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

### Load Drosophila dataset

Source: https://jef.works/MERINGUE/drosophila_3D_analysis

Load the drosophila dataset and extract spatial positions and gene expression matrices.

```R
data("drosophila")
pos <- drosophila$pos
gexp <- drosophila$mat
```

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

### Visualize gene expression and sLISA plots

Source: https://jef.works/MERINGUE/mOB_analysis

Identifies genes upregulated in the Granule Cell Layer and spatially aggregated in Spatial Pattern 1, then plots their expression and signed LISA values.

```R
# look at genes 
gs <- intersect(spatgexp[["Spatial Pattern 2"]], dg.sig[['1: Granular Cell Layer']])
# order by degree of spatial clustering
gs <- gs[order(I[gs,]$observed, decreasing=TRUE)]

# plot
par(mfrow=c(2,2), mar=rep(2,4))
invisible(lapply(gs[1:2], function(g) {
  gexp <- winsorize(scale(mat[g,])[,1])
  plotEmbedding(pos, colors=gexp, cex=2, main=g, verbose=FALSE)
  slisa <- winsorize(signedLisa(gexp, w))
  plotEmbedding(pos, colors=slisa, gradientPalette=colorRampPalette(c('darkgreen', 'white', 'darkorange'))(100), cex=2)
}))
```

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

### Visualize pattern association heatmap

Source: https://jef.works/MERINGUE/drosophila_3D_analysis

Generate a heatmap showing the association between gene spatial patterns.

```R
# Look at pattern association
gcol <- rainbow(length(levels(ggroup$groups)), v=0.5)[ggroup$groups]
names(gcol) <- names(ggroup$groups)
heatmap(scc[ggroup$hc$labels, ggroup$hc$labels], scale='none', 
        Colv=as.dendrogram(ggroup$hc), 
        Rowv=as.dendrogram(ggroup$hc), 
        labCol=NA,
        RowSideColors=gcol[ggroup$hc$labels],
        ColSideColors=gcol[ggroup$hc$labels],
        col=colorRampPalette(c('white', 'black'))(100)
)
```

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

### Parse Tissue Sections

Source: https://jef.works/MERINGUE/BCL_analysis

Extract position and count data for the first three sections of the BCL dataset.

```R
## parse out the 3 sections
BCL1.pos <- BCL$pos[BCL$pos$slice==1,][,c("x", "y")]
BCL2.pos <- BCL$pos[BCL$pos$slice==2,][,c("x", "y")]
BCL3.pos <- BCL$pos[BCL$pos$slice==3,][,c("x", "y")]
BCL1.counts <- BCL$counts[,BCL$pos$slice==1]
BCL2.counts <- BCL$counts[,BCL$pos$slice==2]
BCL3.counts <- BCL$counts[,BCL$pos$slice==3]
```

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

### Create spatial adjacency weight matrix

Source: https://jef.works/MERINGUE/mOB_analysis

Generates a binary adjacency weight matrix based on neighbor relationships within a specified distance.

```R
# Get neighbor-relationships
w <- getSpatialNeighbors(pos, filterDist = 2.5)
plotNetwork(pos, w)
```

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

### Compute spatial neighbors and plot network

Source: https://jef.works/MERINGUE/spatial_clustering

Computes the spatial neighbor graph W using getSpatialNeighbors with filterDist = 2, then plots the network with plotNetwork. This is the first step in integrating spatial information into clustering.

```R
W <- getSpatialNeighbors(pos, filterDist = 2)
plotNetwork(pos, W)
```

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

### Plot spatially clustered genes as heatmap

Source: https://jef.works/MERINGUE/mOB_analysis

Visualizes gene expression across voxels, ordered by cell-type and spatial pattern.

```R
# Plot as heatmap
sp.genes <- unlist(lapply(levels(ggroup$groups), function(x) {
  names(ggroup$groups)[ggroup$groups==x]
}))
ccol <- rainbow(length(levels(annot)))[annot]
names(ccol) <- names(annot)

m <- as.matrix(mat[sp.genes,names(sort(annot))])
m <- winsorize(t(scale(t(m))))
heatmap(m, scale="none", 
          Colv=NA, Rowv=NA, labRow=NA, labCol=NA,
          ColSideColors=ccol[colnames(m)],
          RowSideColors=gcol[rownames(m)],
          col=colorRampPalette(c('blue', 'white', 'red'))(100)
)
```

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

### Execute Individual Section Analysis

Source: https://jef.works/MERINGUE/BCL_analysis

Applies the helper function to each section to identify spatially autocorrelated genes.

```R
# Analyze each section using helper function
par(mfrow=c(1,3), mar=rep(2,4))
L1 <- helper(posList[[1]], mat[, rownames(posList[[1]])])

## Number of significantly autocorrelated genes: 92

## ...driven by > 12.7 cells: 48

L2 <- helper(posList[[2]], mat[, rownames(posList[[2]])])

## Number of significantly autocorrelated genes: 87

## ...driven by > 12.5 cells: 44

L3 <- helper(posList[[3]], mat[, rownames(posList[[3]])])

## Number of significantly autocorrelated genes: 67

## ...driven by > 13.15 cells: 29
```

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

### Compare and Plot Consistent Genes

Source: https://jef.works/MERINGUE/BCL_analysis

Identifies genes significant in multiple sections and consistent across sections, then plots their expression.

```R
gdups <- unlist(list('L1'=L1$sig.genes, 
                     'L2'=L2$sig.genes, 
                     'L3'=L3$sig.genes
))
## genes that are significantly spatially variable in multiple sections
gdups <- gdups[duplicated(gdups)]
## genes that are significantly spatially variable across sections
gcross <- cross$sig.genes
## both 
gall <- intersect(gcross, gdups)
## order by Moran's I
gall <- gall[order(L1$I[gall,]$observed, decreasing=TRUE)]

# Plot few consistent genes
invisible(lapply(gall[1:2], function(g1) {
  par(mfrow=c(1,3), mar=rep(2,4))
  plotEmbedding(posList[[1]], 
                colors=winsorize(scale(mat[g1, rownames(posList[[1]])])[,1]), 
                main=paste0(g1, '\n section 1'),
                cex=2, verbose=FALSE)
  plotEmbedding(posList[[2]], 
                colors=winsorize(scale(mat[g1, rownames(posList[[2]])])[,1]), 
                main=paste0(g1, '\n section 2'), 
                cex=2, verbose=FALSE)
  plotEmbedding(posList[[3]], 
                colors=winsorize(scale(mat[g1, rownames(posList[[3]])])[,1]), 
                main=paste0(g1, '\n section 3'), 
                cex=2, verbose=FALSE)
}))
```

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

### Visualize pair-rule gene embeddings

Source: https://jef.works/MERINGUE/drosophila_3D_analysis

Plot spatial embeddings for a specific subset of pair-rule genes.

```R
par(mfrow=c(2,4), mar=rep(2,4))
## pair-rule genes
sub <- c('eve', 'prd', 'odd', 'run', 'ftz', 'trn', 'h', 'dpn')
invisible(lapply(sub, function(g) {
  plotEmbedding(pos[,c(1,3)], 
                col=gexp[g,], 
                main=g)
}))
```

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

### Construct spatial neighbor network

Source: https://jef.works/MERINGUE/drosophila_3D_analysis

Generate a neighbor adjacency matrix based on 3D spatial coordinates with a distance filter.

```R
N <- getSpatialNeighbors(pos, filterDist = 10, verbose=TRUE)
```

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

### Define cell types and plot them in space

Source: https://jef.works/MERINGUE/simulation

Defines cell-type A as cells expressing Gene3 and cell-type B as cells expressing Gene4, then plots the cell types on the embedding.

```R
ctA <- names(gexp3)[gexp3>0]
ctB <- names(gexp4)[gexp4>0]
```

```R
# double check mutually exclusive
print(intersect(ctA, ctB))
```

```R
# plot
par(mfrow=c(1,1), mar=rep(5,4))
cellType <- factor(rownames(posw) %in% ctA)
levels(cellType) <- c('ctB', 'ctA')
names(cellType) <- rownames(posw)
plotEmbedding(posw, groups=cellType, 
              show.legend=TRUE,  
              main='Cell Types in Space')
```

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

### Spatially-unaware analysis: PCA, tSNE, clustering, and annotation

Source: https://jef.works/MERINGUE/mOB_analysis

Performs dimensionality reduction (PCA, tSNE), graph-based clustering, and manual annotation of clusters with cell types. Includes plotting the embeddings.

```R
# Dimensionality reduction by PCA on log10 CPM expression values
pcs.info <- prcomp(t(log10(as.matrix(mat)+1)), center=TRUE)
nPcs <- 5
pcs <- pcs.info$x[,1:nPcs]

# 2D embedding by tSNE
emb <- Rtsne::Rtsne(pcs,
             is_distance=FALSE,
             perplexity=30,
             num_threads=1,
             verbose=FALSE)$Y
rownames(emb) <- rownames(pcs)

# Graph-based cluster detection
k <- 30
com <- getClusters(pcs, k, weight=TRUE)

# Manually annotate identified clusters with cell-types
annot <- as.character(com); names(annot) <- names(com)
annot[com==4] <- '1: Granule Cell Layer'
annot[com==1] <- '2: Mitral Cell Layer'
annot[com==3] <- '3: Outer Plexiform Layer'
annot[com==2] <- '4: Glomerular Layer'
annot[com==5] <- '5: Olfactory Nerve Layer'
annot <- as.factor(annot)

# Plot
par(mfrow=c(1,2), mar=rep(1,4))
plotEmbedding(emb, groups=annot, 
              show.legend=TRUE, xlab=NA, ylab=NA,
              verbose=FALSE)
plotEmbedding(pos, groups=annot, 
              cex=1, xlab=NA, ylab=NA,
              verbose=FALSE)
```

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

### Define Helper for Section Analysis

Source: https://jef.works/MERINGUE/BCL_analysis

Defines a helper function to calculate spatial neighbors, plot the network, and filter for significant spatially aggregated genes.

```R
helper <- function(pos, mat) {
  w <- getSpatialNeighbors(pos, filterDist = 3)
  plotNetwork(pos, w)
  # get spatially clustered genes
  I <- getSpatialPatterns(mat, w)
  # filter for significant hits driven by more than 10% of cells
  results.filter <- filterSpatialPatterns(mat = mat,
                                          I = I,
                                          w = w,
                                          adjustPv = TRUE,
                                          alpha = 0.05,
                                          minPercentCells = 0.05,
                                          verbose = TRUE)
  # return results
  list(I=I, sig.genes=results.filter)
}
```

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

### Visualize gene expression embeddings

Source: https://jef.works/MERINGUE/drosophila_3D_analysis

Plot gene expression across different 2D projections of the 3D spatial data.

```R
par(mfrow=c(1,3), mar=rep(2,4))
g <- 'apt'
plotEmbedding(pos[,c(1,2)], col=gexp[g,], main='X-Y')
plotEmbedding(pos[,c(3,2)], col=gexp[g,], main='Z-Y')
plotEmbedding(pos[,c(1,3)], col=gexp[g,], main='X-Z')
```

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

### Identify primary spatial patterns

Source: https://jef.works/MERINGUE/mOB_analysis

Uses groupSigSpatialPatterns to cluster genes based on spatial cross-correlation. Requires pre-defined pos, mat, and scc objects.

```R
par(mfrow=c(2,2), mar=rep(2,4))
ggroup <- groupSigSpatialPatterns(pos = pos, 
                                  mat = as.matrix(mat[results.filter,]), 
                                  scc = scc, 
                                  power = 1, 
                                  hclustMethod = 'ward.D', 
                                  deepSplit = 2,
                                  zlim=c(-1.5,1.5))
```

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

### Filter and Normalize Counts

Source: https://jef.works/MERINGUE/BCL_analysis

Filter genes and spots with low counts, then normalize the resulting matrix.

```R
cc <- cleanCounts(counts, min.reads = 100, min.lib.size = 100, plot=TRUE)

mat <- normalizeCounts(cc, log=FALSE, verbose=TRUE)

posList[[1]] <- posList[[1]][intersect(rownames(posList[[1]]), colnames(mat)),]
posList[[2]] <- posList[[2]][intersect(rownames(posList[[2]]), colnames(mat)),]
posList[[3]] <- posList[[3]][intersect(rownames(posList[[3]]), colnames(mat)),]
```

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

### Plot scatterplot and compute Pearson correlation for Gene5 and Gene6

Source: https://jef.works/MERINGUE/simulation

Creates a scatterplot of Gene5 versus Gene6 with a linear regression line, then computes Pearson's product-moment correlation. The test output shows a non-significant correlation (cor = 0.005302307, p-value = 0.937).

```R
# Plot
par(mfrow=c(1,1), mar=rep(2,4))
plot(gexp5, gexp6, main='Scatterplot of\nGene5 versus Gene6')
abline(lm(gexp5~gexp6))
```

```R
# Compute cross correlation
cor.test(gexp5, gexp6)
```

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

### Plot Section Embeddings

Source: https://jef.works/MERINGUE/BCL_analysis

Visualizes the spatial embedding for three sections using plotEmbedding.

```R
# Plot
par(mfrow=c(1,3), mar=rep(2,4))
plotEmbedding(posList[[1]], groups=section, main='section 1', cex=2)
plotEmbedding(posList[[2]], groups=section, main='section 2', cex=2)
plotEmbedding(posList[[3]], groups=section, main='section 3', cex=2)
```

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

### Visualize spatial cross-correlation matrix

Source: https://jef.works/MERINGUE/mOB_analysis

Generates a heatmap of the spatial cross-correlation matrix for the identified gene groups.

```R
# Look at pattern association
gcol <- rainbow(length(levels(ggroup$groups)), v=0.5)[ggroup$groups]
names(gcol) <- names(ggroup$groups)
heatmap(scc[ggroup$hc$labels, ggroup$hc$labels], scale='none', 
        Colv=as.dendrogram(ggroup$hc), 
        Rowv=as.dendrogram(ggroup$hc), 
        labRow=NA, labCol=NA,
        ColSideColors=gcol[ggroup$hc$labels],
        RowSideColors=gcol[ggroup$hc$labels],
        col=colorRampPalette(c('black', 'white'))(100)
)
```

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

### Analyze intra-cell-type spatial heterogeneity

Source: https://jef.works/MERINGUE/mOB_analysis

Restricts analysis to voxels in the Granular Cell Layer to identify genes exhibiting spatial clustering within that specific cell-type.

```R
# Restrict to just voxels corresponding to the Granular Cell Layer.
sub <- names(annot)[annot == '1: Granule Cell Layer']
w.sub <- getSpatialNeighbors(pos[sub,], filterDist = 2.5)
plotNetwork(pos[sub,], w.sub)
```

```R
gs <- dg.sig[['1: Granule Cell Layer']]
I.sub <- getSpatialPatterns(mat[gs, sub], w.sub)
results.filter.sub <- filterSpatialPatterns(mat = mat[, sub],
                                        I = I.sub,
                                        w = w.sub,
                                        adjustPv = TRUE,
                                        alpha = 0.05,
                                        minPercentCells = 0.05,
                                        verbose = TRUE)
```

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

### Perform Multi-section Analysis

Source: https://jef.works/MERINGUE/BCL_analysis

Identifies genes with consistent spatial patterns across sections using cross-layer neighbors.

```R
# K-mutual nearest neighbors across sections
cw <- getCrossLayerNeighbors(posList, k=3)
# Look for genes exhibiting spatial autocorrelation across layers
I <- getSpatialPatterns(mat, cw)
results.filter <- filterSpatialPatterns(mat = mat,
                                        I = I,
                                        w = cw,
                                        adjustPv = TRUE,
                                        alpha = 0.05,
                                        minPercentCells = 0.05/4,
                                        verbose = TRUE)

## Number of significantly autocorrelated genes: 85

## ...driven by > 9.5875 cells: 79

cross <- list(I=I, sig.genes=results.filter)
```

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

### Inter-cell-type spatial cross-correlation

Source: https://jef.works/MERINGUE/simulation

Computes the inter-cell-type weight matrix and then the inter-cell-type spatial cross-correlation statistic and its test.

```R
par(mfrow=c(1,1), mar=rep(5,4))
weightIc <- getInterCellTypeWeight(ctA, ctB, 
                                   weight, posw, 
                                   plot=TRUE, 
                                   main='Adjacency Weight Matrix\nBetween Cell-Types')
```

```R
spatialCrossCor(gexp3, gexp4, weightIc)
interCellTypeSpatialCrossCor(gexp3, gexp4, ctA, ctB, weightIc)
```

```R
spatialCrossCorTest(gexp3, gexp4, weightIc)
```

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

### Compare spatial patterns with differential expression

Source: https://jef.works/MERINGUE/mOB_analysis

Calculates hypergeometric test significance for overlaps between spatial patterns and differentially expressed genes, then visualizes the results as a heatmap.

```R
# Compare two different types of identifying genes
diffgexp <- dg.sig
spatgexp <- lapply(levels(ggroup$groups), function(x) {
  names(ggroup$groups)[ggroup$groups==x]
})
names(spatgexp) <- paste0('Spatial Pattern ', levels(ggroup$groups))

# Assess significance of overlap by hypergeometric test
sigoverlap <- do.call(rbind, lapply(1:length(spatgexp), function(i) {
  so <- unlist(lapply(1:length(diffgexp), function(j) {
    
    #x = # of genes in common between two groups.
    #n = # of genes in group 1.
    #D = # of genes in group 2.
    #N = total genes
    x <- length(intersect(spatgexp[[i]], diffgexp[[j]])) ## shared
    n <- length(spatgexp[[i]])
    D <- length(diffgexp[[j]])
    N <- nrow(counts) ## total
    
    phyper(x, D, N-D, n, lower.tail=FALSE)
  }))
  names(so) <- names(diffgexp)
  return(so)
}))
rownames(sigoverlap) <- 1:length(spatgexp)

# Visualize as heatmap
pvo <- sigoverlap
rownames(pvo) <- paste0('spatial: ', rownames(sigoverlap))
colnames(pvo) <- colnames(sigoverlap)
pvo[pvo < 1e-6] <- 1e-6 # prevent Infs
pvo <- -log10(pvo)
# order for diagonal
matrix.sort <- function(matrix) {
  row.max <- apply(matrix,1,which.max)
  if(all(table(row.max) != 1)) stop("Ties cannot be resolved")
  matrix[names(sort(row.max)),]
}
pvo <- matrix.sort(pvo)
heatmap(pvo, 
        col=colorRampPalette(c('white', 'black'))(100), 
        scale="none", Rowv=NA, Colv=NA, margins = c(25,15))
```

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

### Combine Data Across Sections

Source: https://jef.works/MERINGUE/BCL_analysis

Merge counts and positions from multiple sections into a unified matrix and list.

```R
# Get common set of genes for the sections
genes.have <- Reduce(intersect, list(
  rownames(BCL1.counts),
  rownames(BCL2.counts),
  rownames(BCL3.counts)
))

# Combine into large counts matrix. Counts of genes for all the sections being assessed
counts <- cbind(
  BCL1.counts[genes.have,],
  BCL2.counts[genes.have,],
  BCL3.counts[genes.have,]
)

# section factor
section <-  c(
  rep('L1', ncol(BCL1.counts)),
  rep('L2', ncol(BCL2.counts)),
  rep('L3', ncol(BCL3.counts))
)
names(section) <- colnames(counts)

# List of positions
posList <- list(
  BCL1.pos[colnames(BCL1.counts),],
  BCL2.pos[colnames(BCL2.counts),],
  BCL3.pos[colnames(BCL3.counts),]
)
```

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

### Scatterplot and general cross-correlation of Gene3 vs Gene4

Source: https://jef.works/MERINGUE/simulation

Plots a scatterplot of Gene3 versus Gene4 and computes the Pearson correlation, which is significantly negative.

```R
# Plot
par(mfrow=c(1,1), mar=rep(5,4))
plot(gexp3, gexp4, 
     main='Scatterplot of\nGene3 versus Gene4')
abline(lm(gexp4~gexp3))
```

```R
# Compute cross correlation
cor.test(gexp3, gexp4)
```

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

### Analyze spatial patterns

Source: https://jef.works/MERINGUE/drosophila_3D_analysis

Identify significantly spatially variable genes and filter the results.

```R
## analyze
results <- getSpatialPatterns(gexp, N)
filter <- filterSpatialPatterns(mat = gexp, 
                                I = results, 
                                w = N,
                                alpha = 0.05,
                                details = TRUE, 
                                minPercentCells = 0.05)

print(head(filter))
```

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

### Group genes by spatial patterns

Source: https://jef.works/MERINGUE/drosophila_3D_analysis

Perform spatial cross-correlation analysis and group genes based on their spatial patterns.

```R
scc <- spatialCrossCorMatrix(gexp[sub,], N)
method = 'ward.D'
par(mfrow=c(1,2), mar=rep(2,4))
ggroup <- groupSigSpatialPatterns(pos[,c(1,3)], gexp, scc,
                                    hclustMethod = method,
                                    deepSplit = 2,
                                    binSize = 50,
                                    power = 1)
```

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

### Perform UMAP dimensionality reduction on transcriptome

Source: https://jef.works/MERINGUE/spatial_clustering

Applies UMAP to the gene expression matrix with min_dist=0.5 and visualizes the resulting embedding colored by the expression of the first and 50th genes.

```R
######### See if we can integrate spatial information in cluster detection
par(mfrow=c(1,2))
emb <- uwot::umap(pcs, min_dist = 0.5)
rownames(emb) <- rownames(pcs)
plotEmbedding(emb, col=scale(pcs[,1])[,1], 
              main='Cell-Types Gene Expression', xlab='UMAP X', ylab='UMAP Y')
plotEmbedding(emb, col=scale(pcs[,50])[,1],
              main='Cell-Types Gene Expression', xlab='UMAP X', ylab='UMAP Y')
```

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

### Identify and plot inconsistent genes

Source: https://jef.works/MERINGUE/BCL_analysis

Calculates the set difference between duplicate and cross-section genes, then visualizes the spatial embedding of the top inconsistent gene across three sections.

```R
gdiff <- setdiff(gdups, gcross)
gdiff <- gdiff[order(L1$I[gdiff,]$observed, decreasing=TRUE)]

# Plot inconsistent gene
invisible(lapply(gdiff[1], function(g1) {
  par(mfrow=c(1,3), mar=rep(2,4))
  plotEmbedding(posList[[1]], 
                colors=winsorize(scale(mat[g1, rownames(posList[[1]])])[,1]), 
                main=paste0(g1, '\n section 1'),
                cex=2, verbose=FALSE)
  plotEmbedding(posList[[2]], 
                colors=winsorize(scale(mat[g1, rownames(posList[[2]])])[,1]), 
                main=paste0(g1, '\n section 2'), 
                cex=2, verbose=FALSE)
  plotEmbedding(posList[[3]], 
                colors=winsorize(scale(mat[g1, rownames(posList[[3]])])[,1]), 
                main=paste0(g1, '\n section 3'), 
                cex=2, verbose=FALSE)
}))
```

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

### Plot Gene1 versus Gene2 scatterplot

Source: https://jef.works/MERINGUE/simulation

Plots Gene1 versus Gene2 expression as a scatterplot with a linear regression line. Shows the positive relationship between the two genes.

```R
# Plot
par(mfrow=c(1,1), mar=rep(5,4))
plot(gexp1, gexp2, 
     main='Scatterplot of\nGene1 versus Gene2')
abline(lm(gexp2~gexp1))
```

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

### Spatial cross-correlation of Gene3 vs Gene4

Source: https://jef.works/MERINGUE/simulation

Computes the spatial cross-correlation statistic and its test, showing a positive spatial correlation despite the negative general correlation.

```R
# Compute spatial cross correlation
spatialCrossCor(gexp3, gexp4, weight)
par(mfrow=c(1,1), mar=rep(5,4))
spatialCrossCorTest(gexp3, gexp4, weight, 
                    plot=TRUE)
```

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

### Compute spatial cross-correlation for Gene1 and Gene2

Source: https://jef.works/MERINGUE/simulation

Computes spatial neighbors with getSpatialNeighbors using filterDist = 1, plots the network, then computes the spatial cross-correlation statistic and its test. The output shows spatialCrossCor = 0.8257857 and a p-value of 0.000999001, indicating a significant positive spatial cross-correlation.

```R
weight <- getSpatialNeighbors(posw, filterDist = 1)
plotNetwork(posw, weight)
```

```R
spatialCrossCor(gexp1, gexp2, weight)
par(mfrow=c(1,1), mar=rep(5,4))
spatialCrossCorTest(gexp1, gexp2, weight, 
                    plot=TRUE)
```

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

### Filter spatial patterns

Source: https://jef.works/MERINGUE/mOB_analysis

Applies multiple-testing correction and filters genes based on p-value and minimum percentage of probe spots.

```R
results.filter <- filterSpatialPatterns(mat = mat[test,],
                                        I = I,
                                        w = w,
                                        adjustPv = TRUE,
                                        alpha = 0.05,
                                        minPercentCells = 0.05,
                                        verbose = TRUE)
```

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

### Plot Gene4 expression

Source: https://jef.works/MERINGUE/simulation

Creates a numeric vector for Gene4 expression based on odd/even positions and plots it on the embedding.

```R
vi <- (num %% 2) == 1
gexp4 <- as.numeric(vi)
gexp4 <- gexp4 * sort(abs(rnorm(length(gexp4)))+1)
names(gexp4) <- rownames(pos)
plotEmbedding(posw, col=gexp4,  
              main='Gene4')
```

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

### Compute and test spatial cross-correlation for Gene5 and Gene6

Source: https://jef.works/MERINGUE/simulation

Computes the spatial cross-correlation and performs a permutation test with plotting. The output shows a spatial cross-correlation value of 0.006021392 and a p-value of 0.9320679, indicating no significant spatial cross-correlation.

```R
# Compute spatial cross correlation
spatialCrossCor(gexp5, gexp6, weight)
par(mfrow=c(1,1), mar=rep(2,4))
spatialCrossCorTest(gexp5, gexp6, weight, 
                    plot=TRUE)
```

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

### Identify differentially expressed genes

Source: https://jef.works/MERINGUE/mOB_analysis

Samples 2000 genes and identifies significantly upregulated genes per cluster using the Wilcox test.

```R
set.seed(0)
test <- sample(rownames(mat), 2000)

# Identify significantly differentially upregulated genes
# in each identified cluster by Wilcox test
dg <- getDifferentialGenes(as.matrix(mat[test,]), annot)
dg.sig <- lapply(dg, function(x) {
  x <- x[x$p.adj < 0.05,]
  x <- na.omit(x)
  x <- x[x$highest,]
  rownames(x)
})
print(lapply(dg.sig, length))
```

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

### Compute general cross-correlation with cor.test

Source: https://jef.works/MERINGUE/simulation

Computes Pearson's product-moment correlation between Gene1 and Gene2. The output shows a significant positive correlation (cor = 0.847839, p-value < 2.2e-16).

```R
# Compute cross correlation
cor.test(gexp1, gexp2)
```

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

### Identify spatially auto-correlated genes

Source: https://jef.works/MERINGUE/mOB_analysis

Calculates spatial patterns for the sampled genes using the adjacency weight matrix.

```R
# Identify sigificantly spatially auto-correlated genes
start_time <- Sys.time()
I <- getSpatialPatterns(mat[test,], w)
end_time <- Sys.time()
print(end_time - start_time)
```
Tokens:
Raw