## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
    eval = TRUE,
    collapse = TRUE,
    comment = "#>",
    out.width = "100%",
    dev = "png",
    dpi = 50,
    fig.height = 4.2,
    fig.width = 5.6
)

## ----install_pkg, eval=FALSE--------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# 
# BiocManager::install("SpNeigh")

## ----loadPackage--------------------------------------------------------------
library(SpNeigh)
library(ggplot2)

## -----------------------------------------------------------------------------
coords <- readRDS(system.file("extdata", "MouseBrainCoords.rds",
    package = "SpNeigh"
))
head(coords)

## -----------------------------------------------------------------------------
rownames(coords) <- coords$cell

## -----------------------------------------------------------------------------
logNorm_expr <- readRDS(system.file("extdata", "LogNormExpr.rds",
    package = "SpNeigh"
))
class(logNorm_expr)

## -----------------------------------------------------------------------------
new_cell_types <- c(
    "0" = "Meninges",
    "1" = "Cerebral_cortex",
    "2" = "White_matter",
    "3" = "Cerebral_cortex",
    "4" = "Cerebral_cortex",
    "5" = "Thalamus",
    "6" = "Cerebral_cortex",
    "7" = "Hippocampus",
    "8" = "Hippocampus",
    "9" = "Hippocampus",
    "10" = "White_matter",
    "11" = "Cerebral_cortex"
)

## -----------------------------------------------------------------------------
coords_sub <- subset(coords, cluster %in% c("0", "2"))
coords_sub <- as.matrix(coords_sub[, c("x", "y")])

## -----------------------------------------------------------------------------
metadata_sub <- subset(coords[, c("cell", "cluster")], cluster %in% c("0", "2"))

## -----------------------------------------------------------------------------
spe <- SpatialExperiment::SpatialExperiment(
    assay = list("logcounts" = logNorm_expr),
    colData = metadata_sub,
    spatialCoords = coords_sub
)
spe

## -----------------------------------------------------------------------------
spe_coords <- extractCoords(spe)

## -----------------------------------------------------------------------------
spe_logcounts <- SingleCellExperiment::logcounts(spe)

## -----------------------------------------------------------------------------
library(Seurat)

## -----------------------------------------------------------------------------
seu_sp <- CreateSeuratObject(
    assay = "Spatial",
    counts = logNorm_expr,
    meta.data = metadata_sub
)
seu_sp

## -----------------------------------------------------------------------------
LayerData(seu_sp, assay = "Spatial", layer = "data") <- logNorm_expr

## -----------------------------------------------------------------------------
cents <- CreateCentroids(coords_sub[, c("x", "y")])

fov <- CreateFOV(
    coords = list("centroids" = cents),
    type = c("centroids"),
    assay = "Spatial"
)

seu_sp[["fov"]] <- fov

## -----------------------------------------------------------------------------
seu_sp$seurat_clusters <- seu_sp$cluster

## -----------------------------------------------------------------------------
seu_sp_coords <- extractCoords(seu_sp)

## -----------------------------------------------------------------------------
seu_sp_log_expr <- GetAssayData(seu_sp)

## ----PlotBoundary_spe, fig.width=6.6, fig.height=4.2--------------------------
plotBoundary(spe, one_cluster = "2")

## ----PlotBoundary_seu, fig.width=6.6, fig.height=4.2--------------------------
plotBoundary(seu_sp, one_cluster = "2", eps = 120)

## -----------------------------------------------------------------------------
bon_points_c2 <- getBoundary(
    data = coords,
    one_cluster = 2,
    eps = 120,
    minPts = 10
)
table(bon_points_c2$region_id)

## ----PlotWithBoundary_C2, fig.width=6.6, fig.height=4.2-----------------------
plotBoundary(coords, boundary = bon_points_c2)

## -----------------------------------------------------------------------------
# plotBoundary(coords) + addBoundary(boundary = bon_points_c2)

## ----PlotRegion_C2,fig.height = 2.8,  fig.width = 5.6-------------------------
bon_polys_c2 <- buildBoundaryPoly(bon_points_c2)
plotRegion(bon_polys_c2)

## -----------------------------------------------------------------------------
ring_regions <- getRingRegion(boundary = bon_points_c2)

## ----PlotRing_C2,fig.height = 2.8,  fig.width = 5.6---------------------------
plotRegion(ring_regions)

## -----------------------------------------------------------------------------
cells_ring <- getCellsInside(data = coords, boundary = ring_regions)
cells_ring

## ----PlotCellsInsideRing_C2,fig.height = 3.6,  fig.width = 5.6----------------
plotCellsInside(cells_ring, point_size = 0.2)

## -----------------------------------------------------------------------------
stats_ring <- statsCellsInside(cells_ring)

## ----PlotStatsBar_Proportion_C2_Ring, fig.height = 4.2,  fig.width = 6.6------
plotStatsBar(stats_ring, stat_column = "proportion")

## ----PlotStats_Donut_C2_Ring, fig.height = 3.5,  fig.width = 5.6--------------
plotStatsPie(stats_ring, plot_donut = TRUE)

## -----------------------------------------------------------------------------
interaction_matrix <- computeSpatialInteractionMatrix(
    subset(coords, cell %in% cells_ring$cell)
)
interaction_matrix

## ----Heatmap_InteractionMatrix_C2_Ring, fig.height = 4.2,  fig.width = 5.6----
plotInteractionMatrix(interaction_matrix)

## ----PlotCellsInside_SplitByCluster_C2_Ring, fig.width=8.4, fig.height=4.2----
plotCellsInside(cells_ring) +
    facet_wrap(~cluster) +
    Seurat::RotatedAxis() +
    addBoundaryPoly(ring_regions, linewidth_boundary = 0.1)

## -----------------------------------------------------------------------------
cells_inside <- getCellsInside(data = coords, boundary = bon_points_c2)
cells_inside

## -----------------------------------------------------------------------------
cells_in <- subset(cells_inside, region_id == 1 & cluster == 2)[["cell"]]
cells_out <- subset(cells_ring, region_id == 1 & cluster == 2)[["cell"]]

## -----------------------------------------------------------------------------
tab <- runLimmaDE(
    exp_mat = logNorm_expr,
    min.pct = 0.25,
    cells_reference = cells_in,
    cells_target = cells_out
)

## -----------------------------------------------------------------------------
head(tab[, c("gene", "logFC", "adj.P.Val", "pct.reference", "pct.target")])

## -----------------------------------------------------------------------------
set.seed(123)

## ----Expression_Slc17a7_c2, fig.width = 8, fig.height = 2.5-------------------
plotExpression(
    data = coords[colnames(logNorm_expr), ],
    exp_mat = logNorm_expr,
    genes = tab$gene[1:2],
    sub_plot = TRUE,
    one_cluster = 2,
    shuffle = TRUE,
    point_size = 0.1,
    angle_x_label = 45
)

## ----Expression_Slc17a7_InOutBoundary,fig.height = 2.8, fig.width = 5.6-------
plotExpression(
    data = coords[colnames(logNorm_expr), ],
    exp_mat = logNorm_expr,
    genes = "Slc17a7",
    sub_plot = TRUE,
    shuffle = TRUE,
    sub_cells = c(cells_in, cells_out),
    point_size = 0.3
) +
    addBoundaryPoly(bon_polys_c2[1, ], linewidth_boundary = 0.5)

## ----PlotExpression_Sox10,fig.height = 2.8,  fig.width = 5.6------------------
plotExpression(
    data = coords[colnames(logNorm_expr), ],
    exp_mat = logNorm_expr,
    sub_cells = c(cells_in, cells_out),
    sub_plot = TRUE,
    genes = "Sox10",
    shuffle = TRUE,
    point_size = 0.3
) +
    addBoundaryPoly(bon_polys_c2[1, ], linewidth_boundary = 0.5)

## -----------------------------------------------------------------------------
cells_c0 <- subset(coords, cluster == 0)[, "cell"]

## ----Plot_NoBoundary_C0, fig.width=8*0.7+1, fig.height=6*0.7------------------
plotBoundary(coords[cells_c0, ])

## -----------------------------------------------------------------------------
bon_points_c0 <- getBoundary(data = coords, one_cluster = 0)
bon_polys_c0 <- buildBoundaryPoly(bon_points_c0)

## ----PlotBoundary_C0, fig.width=6.6, fig.height=4.2---------------------------
plotEdge(boundary_poly = bon_polys_c0, linewidth_boundary = 0.6)

## -----------------------------------------------------------------------------
weights_bon <- computeBoundaryWeights(
    data = coords,
    cell_ids = cells_c0,
    boundary = bon_points_c0
)

## -----------------------------------------------------------------------------
weights_bon[1:3]

## -----------------------------------------------------------------------------
weights_cen <- computeCentroidWeights(data = coords, cell_ids = cells_c0)

## ----fig.width=6.6, fig.height=4.2--------------------------------------------
plotWeights(data = coords, weights = weights_bon, point_size = 0.8) +
    labs(title = "Boundary weights")

## -----------------------------------------------------------------------------
tab_sp <- runSpatialDE(
    exp_mat = logNorm_expr[, cells_c0],
    spatial_distance = weights_bon,
    cell_ids = cells_c0
)

## -----------------------------------------------------------------------------
head(tab_sp)

## ----Expression_top_genes_boundary, fig.width = 8, fig.height = 2.5-----------
plotExpression(
    data = coords[colnames(logNorm_expr), ],
    exp_mat = logNorm_expr,
    genes = tab_sp$gene[1:2],
    sub_plot = TRUE,
    one_cluster = 0,
    shuffle = TRUE,
    point_size = 0.1,
    angle_x_label = 45
)

## ----Expression_along_boundary_weights_c0, fig.width=6.4, fig.height=4--------
plotSpatialExpression(
    exp_mat = logNorm_expr[, cells_c0],
    spatial_distance = weights_bon,
    scale_method = "minmax",
    genes = tab_sp$gene[1:10],
    label_x = "Boundary weights"
)

## -----------------------------------------------------------------------------
SEI_bon <- computeSpatialEnrichmentIndex(
    exp_mat = logNorm_expr[, cells_c0],
    weights = weights_bon
)
head(SEI_bon)

## ----Expression_top_genes_SEI_boundary, fig.width = 8, fig.height = 2.5-------
plotExpression(
    data = coords[colnames(logNorm_expr), ],
    exp_mat = logNorm_expr,
    genes = SEI_bon$gene[1:2],
    sub_plot = TRUE,
    one_cluster = 0,
    shuffle = TRUE,
    point_size = 0.1,
    angle_x_label = 45
)

## -----------------------------------------------------------------------------
sessionInfo()

