library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)
library(hdf5r)
library(Rfast2)

# load data with spatial coordinates and low-res image
brain <- Load10X_Spatial("spatial-mouse-sagittal-anterior/",filename="V1_Mouse_Brain_Sagittal_Anterior_Section_2_filtered_feature_bc_matrix.h5")

# basic control plot
plot1 <- VlnPlot(brain, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")
wrap_plots(plot1, plot2)

# normalization
brain <- SCTransform(brain, assay = "Spatial", verbose = FALSE)

# display the expression of two chosen genes
SpatialFeaturePlot(brain, features = c("Hpca", "Ttr")) #  Hpca = hippocampus marker, Ttr = choroid plexus
plot <- SpatialFeaturePlot(brain, features = c("Ttr")) + theme(legend.text = element_text(size = 0),
                                                               legend.title = element_text(size = 20), legend.key.size = unit(1, "cm"))
plot

# cluster the Visium spots to define areas or regions
brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)
brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30)
brain <- FindClusters(brain, verbose = FALSE)
brain <- RunUMAP(brain, reduction = "pca", dims = 1:30)
p1 <- DimPlot(brain, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)
p1 + p2

# monochrome display of chosen regions
SpatialDimPlot(brain, cells.highlight = CellsByIdentities(object = brain, idents = c(2, 1, 4, 3, 5, 8)),
               facet.highlight = TRUE, ncol = 3)

# find genes that display differential expression between regions 3 and 4, plot top 3
de_markers <- FindMarkers(brain, ident.1 = 3, ident.2 = 4)
SpatialFeaturePlot(object = brain, features = rownames(de_markers)[1:3], alpha=c(1,1),ncol = 3)

# an example of gene with same expression in regions 3 and 4
SpatialFeaturePlot(object = brain, features = "Camk2n1")

# find all the variable genes
brain <- FindSpatiallyVariableFeatures(brain, assay = "SCT", features = VariableFeatures(brain)[1:1000],
                                       selection.method = "moransi")
top.features <- head(SpatiallyVariableFeatures(brain, selection.method = "moransi"), 6)
SpatialFeaturePlot(brain, features = top.features, ncol = 3, alpha = c(0.1, 1))

# load a cortex single-cell atlas
allen_reference <- readRDS("allen_cortex.rds")
allen_reference <- SCTransform(allen_reference, ncells = 3000, verbose = FALSE) %>%
  RunPCA(verbose = FALSE) %>%
  RunUMAP(dims = 1:30)
DimPlot(allen_reference, group.by = "subclass", label = TRUE)

# score the brain regions
anchors <- FindTransferAnchors(reference = allen_reference, query = brain, normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors, refdata = allen_reference$subclass, prediction.assay = TRUE,
                                  weight.reduction = brain[["pca"]], dims = 1:30)
brain[["predictions"]] <- predictions.assay
DefaultAssay(brain) <- "predictions"
SpatialFeaturePlot(brain, features = c("L2/3 IT", "L4", "Astro", "L5 PT"), pt.size.factor = 1.6, ncol = 2, crop = TRUE)