Spatial Transcriptomics and Cell Type Decomposition

library(Seurat)
library(dplyr)
library(ggplot2)
library(patchwork)

dcis <- Load10X_Spatial(data.dir = here::here("data/dcis"),
                filename = "filtered_feature_bc_matrix.h5")
dcis_norm <- SCTransform(dcis, assay = "Spatial", verbose = FALSE)
dcis_obj <- RunPCA(dcis_norm, assay = "SCT", verbose = FALSE)
# compute K nearest neighbors (KNN)
dcis_obj <- FindNeighbors(dcis_obj, reduction = "pca", dims = 1:30)
# Leiden algorithm for community detection
dcis_obj <- FindClusters(dcis_obj, verbose = FALSE)
# PCA result is the default UMAP input, use dimensions 1:30 as input features
dcis_obj <- RunUMAP(dcis_obj, reduction = "pca", dims = 1:30)

annotation_data_dcis <- read.csv(here::here("data/dcis/Pathology.csv"), header = TRUE, row.names = 1)
dcis_obj <- AddMetaData(dcis_obj, metadata = annotation_data_dcis)
dcisPlot <- DimPlot(dcis_obj, group.by = "Pathology", label = TRUE, pt.size = 0.5) + NoLegend()
annotation_data_ids <- read.csv(here::here("data/ids/Pathology.csv"), header = TRUE, row.names = 1)
ids_obj <- AddMetaData(ids_obj, metadata = annotation_data_ids)
idsPlot <- DimPlot(ids_obj, group.by = "Pathology", label = TRUE, pt.size = 0.5) + NoLegend()
dcisPlot + idsPlot

merged_obj <- merge(x = dcis, y = ids, add.cell.ids = c("DCIS", "IDS"))
merged_obj <- NormalizeData(merged_obj)
merged_obj <- FindVariableFeatures(merged_obj)
merged_obj <- ScaleData(merged_obj)
merged_obj <- RunPCA(merged_obj)
merged_obj <- IntegrateLayers(object = merged_obj,
                              method = RPCAIntegration,
                              orig.reduction = "pca",
                              new.reduction = "integrated.rpca",
                              verbose = FALSE)
obj <- FindNeighbors(merged_obj, reduction = "integrated.rpca", dims = 1:30)
obj <- FindClusters(obj, resolution = 0.8, cluster.name = "rpca_clusters")

obj <- RunUMAP(obj, reduction = "integrated.rpca", dims = 1:30, reduction.name = "umap.rpca")
obj <- BuildClusterTree(obj)
obj <- JoinLayers(obj)

DimPlot(
  obj,
  reduction = "umap.rpca",
  group.by = "Pathology",
  combine = FALSE, label.size = 2
)

DimPlot(
  obj,
  reduction = "umap.rpca",
  label = TRUE,
  group.by = "rpca_clusters",
  combine = FALSE, label.size = 2
)

SpatialDimPlot(obj, group.by = 'Pathology')

SpatialDimPlot(obj, group.by = 'seurat_clusters')

PlotClusterTree(obj)

cluster_id <- 13
degs <- FindMarkers(obj, ident.1 = cluster_id, min.pct = 0.25, logfc.threshold = 0.25)

SpatialFeaturePlot(object = obj, 
                   features = rownames(degs)[1:2], 
                   alpha = c(0.1, 1), ncol = 2)