单细胞空间转录组RCTD去卷积分析学习和整理

rm(list = ls())
library(spacexr)
library(Matrix)
library(Seurat)
library(qs)
library(BiocParallel)
register(MulticoreParam(workers = 8, progressbar = TRUE)) 

# 10x Visium
sce_s <- Load10X_Spatial(data.dir = "./RawData/GSE6716963/",
                           filename = "GSM6716963_19G081_filtered_feature_bc_matrix.h5")

# reference
sc_dataset <- qread("sc_dataset.qs")

# 过滤
sce_s <- subset(sce_s, nCount_Spatial > 5000)
sce_s <- subset(sce_s, nCount_Spatial < 35000)
# 归一化/找高变/ScaleData
sce_s <- NormalizeData(sce_s)
sce_s <- FindVariableFeatures(sce_s)
sce_s <- ScaleData(sce_s)
# PCA/UMAP
sce_s <- RunPCA(sce_s)
sce_s <- RunUMAP(sce_s, reduction = "pca", dims = 1:50)
# FindNeihbors/Cluster
sce_s <- FindNeighbors(sce_s, reduction = "pca", dims=1:50)
sce_s <- FindClusters(sce_s, resolution = 0.5, verbose = FALSE)

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

# reference处理##################
# reference细胞数量少一点，多了就好慢
ref <- subset(sc_dataset,downsample = 2000)
ref <- SCTransform(sc_dataset,ncells = 3000, 
                         verbose = FALSE) %>%
    RunPCA(verbose = FALSE) %>%
    RunUMAP(dims = 1:30)

# 建议大家都不要用特殊符号
# R语言不喜欢这些特殊符号
a <- ref@meta.data
a$celltype <- sub("/","_",a$celltype)
ref@meta.data <- a

DimPlot(ref, group.by = "celltype", label = TRUE)

# 如果抽样之后的细胞数太少(25个？)需要去除
ref <- subset(ref,celltype !="C3"&celltype !="C5")

# extract information to pass to the RCTD Reference function
counts <- ref[["RNA"]]$counts # 官网使用counts
cluster <- as.factor(ref$celltype)
cluster <- droplevels(cluster) # 如果去除了某些celltype
names(cluster) <- colnames(ref)
nUMI <- ref$nCount_RNA 
names(nUMI) <- colnames(ref)
reference <- Reference(counts, cluster, nUMI)

# set up query with the RCTD function SpatialRNA
counts <- sce_s[["Spatial"]]@counts
coords <- GetTissueCoordinates(sce_s)
colnames(coords) <- c("x", "y")
coords[is.na(colnames(coords))] <- NULL
query <- SpatialRNA(coords, counts, colSums(counts))

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

sce_s <- AddMetaData(sce_s, metadata = RCTD@results$results_df)
# 通常关注这个SPOT丰度最高
p1 <- SpatialDimPlot(sce_s, group.by = "first_type") 
# SPOT丰度第二高
p2 <- SpatialDimPlot(sce_s, group.by = "second_type")
p1+p2

barcodes <- colnames(RCTD@spatialRNA@counts)
weights <- RCTD@results$weights
norm_weights <- normalize_weights(weights)

# observe weight values
celltypes <- c('pericyte_SMC', 'endothelial cells','mast cells')
print(head(norm_weights[,celltypes])) 
head(norm_weights)

# plot
p <- plot_puck_continuous(RCTD@spatialRNA,
                        barcodes,
                        norm_weights[,'pericyte_SMC'],
                        ylimit=c(0,0.5),
                        title='plotofDentateweights',
                        size=4.5,alpha=0.8);p
ggsave("Spaital_weights.pdf",width=8,height=6,plot=p,bg="white")

library(STdeconvolve)
library(ggplot2)
library(ggsci)
packageVersion("STdeconvolve")

m <- as.matrix(norm_weights)
p <- coords

plt <- vizAllTopics(theta = m,
             pos = p,
             topicOrder=seq(ncol(m)),
             topicCols=rainbow(ncol(m)),
             groups = NA,
             group_cols = NA,
             r = 3, # size of scatterpies; adjust depending on the coordinates of the pixels
             lwd = 0.3,
             showLegend = TRUE,
             plotTitle = "scatterpies")

## function returns a `ggplot2` object, so other aesthetics can be added on:
plt <- plt + ggplot2::guides(fill=ggplot2::guide_legend(ncol=2));plt
ggsave("Spaital_scatterpies.pdf", width=9, height=6, plot=plt, bg="white")