library(Seurat) 
library(tidyverse)
library(sloop)


## 10x Visium
seu_vis <- Load10X_Spatial(data.dir = "data/human_lymph_node/",
                           filename = "V1_Human_Lymph_Node_filtered_feature_bc_matrix.h5")

image <- Read10X_Image(image.dir = "data/human_lymph_node/spatial/", image.name = "tissue_hires_image.png")
seu_vis_hires <- Load10X_Spatial(data.dir = "data/human_lymph_node/", image = image,
                                 filename = "V1_Human_Lymph_Node_filtered_feature_bc_matrix.h5")

class(seu_vis[["slice1"]])
head(seu_vis[["slice1"]]@coordinates)
class(seu_vis[["slice1"]]@image)
dim(seu_vis[["slice1"]]@image)

head(GetTissueCoordinates(seu_vis))
head(GetTissueCoordinates(seu_vis[["slice1"]]))
s3_dispatch(GetTissueCoordinates(seu_vis[["slice1"]]))

ScaleFactors(seu_vis[["slice1"]])

head(slot(seu_vis[["slice1"]], name = "coordinates"))

View(slot(seu_vis[["slice1"]], name = "coordinates"))

data.use <- seu_vis[["slice1"]]@coordinates
p1 <- ggplot(data.use, aes(row, col)) + geom_point(size = 2)
p2 <- ggplot(data.use, aes(imagerow, imagecol)) + geom_point(size = 2)
p1 + p2

data.use <- GetTissueCoordinates(seu_vis)
p3 <- ggplot(data.use, aes(imagerow, imagecol)) + geom_point(size = 2)

p1 + p2 + p3



library(grid)
library(gridExtra)
slice1 <- seu_vis[["slice1"]]@image

grid.raster(slice1)

slice1.R <- slice1[,,1]
slice1.G <- slice1[,,2]
slice1.B <- slice1[,,3]

img1 = rasterGrob(slice1.R)
img2 = rasterGrob(slice1.G)
img3 = rasterGrob(slice1.B)

grid.arrange(img1, img2, img3, nrow=1)

df = data.frame(
  red = matrix(slice1[,,1], ncol=1),
  green = matrix(slice1[,,2], ncol=1),
  blue = matrix(slice1[,,3], ncol=1)
)

K = kmeans(df,4)
df$label = K$cluster

C1 = matrix(df$label==1, nrow=dim(slice1)[1])
C2 = matrix(df$label==2, nrow=dim(slice1)[1])
C3 = matrix(df$label==3, nrow=dim(slice1)[1])
C4 = matrix(df$label==4, nrow=dim(slice1)[1])

img0 <- rasterGrob(slice1)
img1 <- rasterGrob(C1)
img2 <- rasterGrob(C2)
img3 <- rasterGrob(C3)
img4 <- rasterGrob(C4)
grid.arrange(img0, img1, img2, img3, img4, nrow=1)

colors = data.frame(
  label = 1:nrow(K$centers),
  R = K$centers[,"red"],
  G = K$centers[,"green"],
  B = K$centers[,"blue"]
)

df$order = 1:nrow(df)
df = merge(df, colors)
df = df[order(df$order),]
df$order = NULL

# get mean color channel values for each row of the df.
R = matrix(df$R, nrow=dim(slice1)[1])
G = matrix(df$G, nrow=dim(slice1)[1])
B = matrix(df$B, nrow=dim(slice1)[1])

slice1_segmented = array(dim=dim(slice1))
slice1_segmented[,,1] = R
slice1_segmented[,,2] = G
slice1_segmented[,,3] = B


img1 = rasterGrob(slice1)
img2 = rasterGrob(slice1_segmented)
grid.arrange(img1, img2, nrow=1)


img0 <- rasterGrob(slice1_segmented)
img1 <- rasterGrob(slice1)
grid.arrange(img0, img1, nrow=1)

library(Seurat) ## V4
library(tidyverse)

seu_vis <- Load10X_Spatial(data.dir = "data/human_lymph_node/",
                           filename = "V1_Human_Lymph_Node_filtered_feature_bc_matrix.h5")

p1 <- VlnPlot(seu_vis, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
p2 <- SpatialFeaturePlot(seu_vis, features = "nCount_Spatial") + theme(legend.position = "right")
p1 | p2


seu_vis <- subset(seu_vis, nCount_Spatial > 5000)
seu_vis <- subset(seu_vis, nCount_Spatial < 35000)

seu_vis <- NormalizeData(seu_vis)
seu_vis <- FindVariableFeatures(seu_vis)
seu_vis <- ScaleData(seu_vis)

seu_vis <- RunPCA(seu_vis)
seu_vis <- RunUMAP(seu_vis, reduction = "pca", dims = 1:50)

seu_vis <- FindNeighbors(seu_vis, reduction = "pca", dims=1:50)
seu_vis <- FindClusters(seu_vis, resolution = 0.8)

DimPlot(seu_vis, reduction = "umap", group.by = "Spatial_snn_res.0.8")

SpatialDimPlot(seu_vis, group.by = "Spatial_snn_res.0.8") + ggsci::scale_fill_d3()

Idents(seu_vis) <- seu_vis$Spatial_snn_res.0.8
all.markers <- FindAllMarkers(seu_vis, only.pos = T)
head(all.markers)

cells <- rownames(subset(seu_vis@meta.data, Spatial_snn_res.0.8 == 0))
p1 <- DimPlot(seu_vis, reduction = "umap", cells.highlight = cells)
p2 <- FeaturePlot(seu_vis, reduction = "umap", features = "CCL14") + theme(legend.position = "right")
p1 | p2


cells <- rownames(subset(seu_vis@meta.data, Spatial_snn_res.0.8 == 0))
p1 <- SpatialDimPlot(seu_vis, cells.highlight = cells)
p2 <- SpatialFeaturePlot(seu_vis, features = "CST3") + theme(legend.position = "right")
p1 | p2


SpatialDimPlot(seu_vis, group.by = "Spatial_snn_res.0.8") + ggsci::scale_fill_d3()

seu_vis2 <- subset(seu_vis, Spatial_snn_res.0.8 %in% c(2,4))
SpatialDimPlot(seu_vis2, group.by = "Spatial_snn_res.0.8") + ggsci::scale_fill_d3()

library(spacexr)
library(Seurat)
library(tidyverse)

## 导入参考数据集
ref <- sceasy::convertFormat("data/immune_reference/sc.h5ad", from = "anndata", to = "seurat")


unique(ref$Subset)

Idents(ref) <- ref$Subset
plot <- DimPlot(ref, reduction = "umap")
LabelClusters(plot = plot, id = "ident")

table(ref$Subset) %>% sort()

ref.ds <- subset(ref, downsample = 350)
table(ref.ds$Subset) %>% sort()

ref.ds <- subset(ref.ds, Subset != "Mast")

# extract information to pass to the RCTD Reference function
counts <- ref.ds[["RNA"]]@counts
cluster <- ref.ds$Subset # should be named factors
cluster <- droplevels(cluster)
nUMI <- colSums(counts)
reference <- Reference(counts, cluster, nUMI) # only support <= 10,000 cells
class(reference)

# set up query with the RCTD function SpatialRNA
seu_vis <- Load10X_Spatial(data.dir = "data/human_lymph_node/",
                           filename = "V1_Human_Lymph_Node_filtered_feature_bc_matrix.h5")

counts <- seu_vis[["Spatial"]]@counts
coords <- GetTissueCoordinates(seu_vis)
colnames(coords) <- c("x", "y")
query <- SpatialRNA(coords, counts, colSums(counts))
class(query)

# deconvolution
RCTD <- create.RCTD(query, reference, max_cores = 8)
RCTD <- run.RCTD(RCTD, doublet_mode = "doublet")

dir.create("tmp")
qs::qsave(RCTD, "tmp/RCTD.qs")

seu_vis <- AddMetaData(seu_vis, metadata = RCTD@results$results_df)
Idents(seu_vis) <- seu_vis$first_type

table(is.na(Idents(seu_vis)))

qs::qsave(seu_vis, "tmp/seu_vis.RCTD.qs")

SpatialDimPlot(seu_vis, cells.highlight = CellsByIdentities(object = seu_vis, idents = levels(Idents(seu_vis))),
               facet.highlight = T, ncol = 10)

weights.df <- as.data.frame(as.matrix(RCTD@results$weights))

seu_vis <- AddMetaData(seu_vis, metadata = weights.df)
colnames(seu_vis@meta.data)[13:45]
SpatialFeaturePlot(seu_vis, features = colnames(seu_vis@meta.data)[13:23], ncol = 5)
SpatialFeaturePlot(seu_vis, features = colnames(seu_vis@meta.data)[24:33], ncol = 5)
SpatialFeaturePlot(seu_vis, features = colnames(seu_vis@meta.data)[34:45], ncol = 5)
SpatialFeaturePlot(seu_vis, features = "T_CD4._TfH_GC")

## correlation
pheatmap::pheatmap(cor(weights.df))

SpatialFeaturePlot(seu_vis, features = c("FDC", "B_GC_LZ", "B_GC_prePB", "T_CD4._TfH_GC"))
seu_vis$Endo[is.na(seu_vis$Endo)] <- 0
seu_vis$VSMC[is.na(seu_vis$VSMC)] <- 0
SpatialFeaturePlot(seu_vis, features = c("Endo", "VSMC"), max.cutoff = "q99")

## NMF
K = 12
model <- RcppML::nmf(t(weights.df), k = K)
w <- model$w
h <- model$h
dim(w)
dim(h)

colnames(w) <- paste0("factor.", 1:K)
rownames(w) <- colnames(weights.df)

pheatmap::pheatmap(w)

w <- t(w) %>% as.data.frame()
w <- w %>% mutate(celltype = rownames(.), .before=1)
w <- w %>% pivot_longer(cols = 2:ncol(.), names_to = "factors", values_to = "weight")
ggplot(w, aes(factors, celltype)) +
  geom_point(aes(color = weight, size = weight)) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_color_viridis_c()