library(Seurat)
library(SPARK)

so <- readRDS("01_spatial_qc.rds")

# Step 2.1: Dimension reduction and clustering
so <- RunPCA(so, assay="SCT", verbose=FALSE)
so <- FindNeighbors(so, reduction="pca", dims=1:30)
so <- FindClusters(so, verbose=FALSE, resolution={resolution})
so <- RunUMAP(so, reduction="pca", dims=1:30)

# Step 2.2: Spatial cluster visualization
SpatialDimPlot(so, label=TRUE, label.size=3, pt.size.factor=1.5)
ggsave("02_spatial_clusters.pdf", width=8, height=6)

# Step 2.3: Find SVGs with SPARK
raw_count <- GetAssayData(so, layer="counts")
info <- cbind.data.frame(x=so@images[[Images(so)[1]]]@coordinates$col,
                         y=so@images[[Images(so)[1]]]@coordinates$row)
rownames(info) <- colnames(raw_count)

spark <- CreateSPARKObject(counts=raw_count, location=info, percentage=0.1, min_total_counts=100)
spark <- spark.vc(spark, covariates=NULL, lib_size=NULL, num_core=4, verbose=FALSE)
spark <- spark.test(spark, check_positive=TRUE, verbose=FALSE)
write.csv(spark@res_mtest, "02_SVG_results.csv")

saveRDS(so, "02_spatial_clustered.rds")

library(RCTD)
library(Seurat)

so <- readRDS("02_spatial_clustered.rds")
sc_ref <- readRDS("{sc_reference}")

# Step 3.1: Prepare RCTD reference
counts <- GetAssayData(sc_ref, layer="counts")
cell_types <- sc_ref$cell_type
names(cell_types) <- colnames(sc_ref)
nUMI <- sc_ref$nCount_RNA
names(nUMI) <- colnames(sc_ref)

reference <- Reference(counts, cell_types, nUMI)

# Step 3.2: Prepare spatial data
spatial_counts <- GetAssayData(so, layer="counts")
spatial_coords <- so@images[[Images(so)[1]]]@coordinates[,c("col","row")]
colnames(spatial_coords) <- c("x","y")

puck <- SpatialRNA(spatial_coords, spatial_counts, colSums(spatial_counts))

# Step 3.3: Run RCTD
my_rctd <- create.RCTD(puck, reference, max_cores=4)
my_rctd <- run.RCTD(my_rctd, doublet_mode="doublet")

# Step 3.4: Extract results
results <- my_rctd@results
weights <- as.data.frame(results$weights)
weights_norm <- weights / rowSums(weights)
so@meta.data <- cbind(so@meta.data, weights_norm)

# Step 3.5: Visualize deconvolution
SpatialFeaturePlot(so, features=colnames(weights_norm)[1:4], pt.size.factor=1.5, ncol=2)
ggsave("03_deconvolution.pdf", width=12, height=8)

saveRDS(so, "03_spatial_niche.rds")

library(Seurat)
library(ggplot2)

so <- readRDS("03_spatial_niche.rds")

# Step 4.1: Integrate metabolite data (MALDI-MSI)
maldi_data <- read.csv("{maldi_csv}", row.names=1)
so[["MALDI"]] <- CreateAssayObject(counts=maldi_data[, colnames(so)])
so <- NormalizeData(so, assay="MALDI")

# Step 4.2: Weighted Nearest Neighbor analysis
so <- FindMultiModalNeighbors(so, reduction.list=list("pca","MALDI.pca"),
                              dims.list=list(1:30, 1:10), modality.weight.name="MALDI.weight")
so <- RunUMAP(so, nn.name="weighted.nn", reduction.name="wnn.umap", reduction.key="wnnUMAP_")

# Step 4.3: Visualize
DimPlot(so, reduction="wnn.umap", label=TRUE)
ggsave("04_multimodal_umap.pdf", width=8, height=6)

# Step 4.4: Correlation between gene and metabolite
DefaultAssay(so) <- "SCT"
gene_expr <- FetchData(so, vars="Hmgcs2")
met_expr <- FetchData(so, vars="beta_hydroxybutyrate")
cor_result <- cor.test(gene_expr$Hmgcs2, met_expr$beta_hydroxybutyrate)
print(paste("Hmgcs2 vs beta-HB correlation:", round(cor_result$estimate, 3)))

saveRDS(so, "04_spatial_multimodal.rds")