空间转录组数据注释分析：SPOTlight反卷积

# 设置一下镜像
## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
options("repos"=c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"))

#install.packages("devtools")
library(devtools)
install_github("https://github.com/MarcElosua/SPOTlight")

rm(list=ls())
library(ggplot2)
library(SPOTlight)
library(SingleCellExperiment)
library(SpatialExperiment)
library(scater)
library(scran)

# 数据准备
# 空转数据
library(TENxVisiumData)
spe <- MouseKidneyCoronal()
spe
# 将spe的行名转为基因名字
head(rownames(spe))
rownames(spe) <- rowData(spe)$symbol
# 先保存一下，下次可以直接导入
library(qs)
qsave(spe, file="spe.qs")


# 单细胞数据
library(TabulaMurisSenisData)
sce <- TabulaMurisSenisDroplet(tissues = "Kidney")$Kidney
sce
# 先保存一下，下次可以直接导入
qsave(sce, file="sce.qs")
table(sce$free_annotation, sce$age)

# 选择同一个时间段的细胞：18m mice
sce <- sce[, sce$age == "18m"]
# Keep cells with clear cell type annotations
sce <- sce[, !sce$free_annotation %in% c("nan", "CD45")]
sce

# Feature selection
sce <- logNormCounts(sce)

# Variance modelling
# 首先排除 ribosomal or mitochondrial 基因
genes <- !grepl(pattern = "^Rp[l|s]|Mt", x = rownames(sce))
table(genes)
dec <- modelGeneVar(sce, subset.row = genes)
plot(dec$mean, dec$total, xlab = "Mean log-expression", ylab = "Variance")
curve(metadata(dec)$trend(x), col = "blue", add = TRUE)
# Get the top 3000 genes.
hvg <- getTopHVGs(dec, n = 3000)
hvg
grep("CD3", hvg, ignore.case = T,value = T)

# 2.细胞类型的特征基因
colLabels(sce) <- colData(sce)$free_annotation
# Compute marker genes
mgs <- scoreMarkers(sce, subset.row = genes)

mgs_fil <- lapply(names(mgs), function(i) {
  x <- mgs[[i]]
# Filter and keep relevant marker genes, those with AUC > 0.8
  x <- x[x$mean.AUC > 0.8, ]
# Sort the genes from highest to lowest weight
  x <- x[order(x$mean.AUC, decreasing = TRUE), ]
# Add gene and cluster id to the dataframe
  x$gene <- rownames(x)
  x$cluster <- i
  data.frame(x)
})
mgs_df <- do.call(rbind, mgs_fil)

# 3.Cell Downsampling
# 生成一个每种细胞类型的list对象
idx <- split(seq(ncol(sce)), sce$free_annotation)
# 降采样到每种细胞20个细胞，真实数据中需要使用100个以上
n_cells <- 20
cs_keep <- lapply(idx, function(i) {
  n <- length(i)
  if (n < n_cells)
    n_cells <- n
  sample(i, n_cells)
})
sce <- sce[, unlist(cs_keep)]
sce
table(sce$free_annotation)

## 反卷积
res <- SPOTlight(
  x = sce,
  y = spe,
  groups = as.character(sce$free_annotation),
  mgs = mgs_df,
  hvg = hvg,
  weight_id = "mean.AUC",
  group_id = "cluster",
  gene_id = "gene")


# 提取反卷积矩阵：每一行为一个spot，每一列为一种细胞类型，值为细胞类型在spot中的相对百分比
# 行和一定为1
head(mat <- res$mat)[, 1:3]
dim(mat)
gead(rowSums(mat))
# 提取NMF模型
mod <- res$NMF

## 可视化
plotTopicProfiles(
  x = mod,
  y = sce$free_annotation,
  facet = FALSE,
  min_prop = 0.01,
  ncol = 1) +
  theme(aspect.ratio = 1)

plotTopicProfiles(
  x = mod,
  y = sce$free_annotation,
  facet = TRUE,
  min_prop = 0.01,
  ncol = 6)

library(NMF)
sign <- basis(mod)
colnames(sign) <- paste0("Topic", seq_len(ncol(sign)))
head(sign)
# This can be dynamically visualized with DT as shown below
DT::datatable(sign, fillContainer = TRUE, filter = "top")

plotCorrelationMatrix(mat)

plotInteractions(mat, which = "heatmap", metric = "prop")
plotInteractions(mat, which = "heatmap", metric = "jaccard")
plotInteractions(mat, which = "network")

## 饼图
ct <- colnames(mat)
# 占比小于0.1的不展示
mat[mat < 0.1] <- 0

# 颜色设置
paletteMartin <- c(
"#000000", "#004949", "#009292", "#ff6db6", "#ffb6db", 
"#490092", "#006ddb", "#b66dff", "#6db6ff", "#b6dbff", 
"#920000", "#924900", "#db6d00", "#24ff24", "#ffff6d")
pal <- colorRampPalette(paletteMartin)(length(ct))
names(pal) <- ct
pal

plotSpatialScatterpie(
  x = spe,
  y = mat,
  cell_types = colnames(mat),
  img = FALSE,
  scatterpie_alpha = 1,
  pie_scale = 0.4) +
  scale_fill_manual(
    values = pal,
    breaks = names(pal))