原理：使用NMFreg的方法。

install.packages("devtools") devtools::install_github("https://github.com/MarcElosua/SPOTlight") library(SPOTlight) library(Seurat) library(dplyr)library(ggplot2) library(SPOTlight) library(SingleCellExperiment) install.packages("SpatialExperiment") library(SpatialExperiment) library(scater) library(scran)BiocManager::install("TENxVisiumData") BiocManager::install("TabulaMurisSenisData") BiocManager::install("SpatialExperiment") library(TENxVisiumData) spe <- MouseKidneyCoronal() # Use symbols instead of Ensembl IDs as feature names rownames(spe) <- rowData(spe)$symbollibrary(TabulaMurisSenisData) sce <- TabulaMurisSenisDroplet(tissues = "Kidney")$Kidney scetable(sce$free_annotation, sce$age) # Keep cells from 18m mice sce <- sce[, sce$age == "18m"] # Keep cells with clear cell type annotations sce <- sce[, !sce$free_annotation %in% c("nan", "CD45")]sce <- logNormCounts(sce)dec <- modelGeneVar(sce) plot(dec$mean, dec$total, xlab = "Mean log-expression", ylab = "Variance") curve(metadata(dec)$trend(x), col = "blue", add = TRUE)hvg <- getTopHVGs(dec, n = 3000)colLabels(sce) <- colData(sce)$free_annotation# Get vector indicating which genes are neither ribosomal or mitochondrial genes <- !grepl(pattern = "^Rp[l|s]|Mt", x = rownames(sce))# 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.8x <- x[x$mean.AUC > 0.8, ]# Sort the genes from highest to lowest weightx <- x[order(x$mean.AUC, decreasing = TRUE), ]# Add gene and cluster id to the dataframex$gene <- rownames(x)x$cluster <- idata.frame(x) }) mgs_df <- do.call(rbind, mgs_fil)#downsample # split cell indices by identity idx <- split(seq(ncol(sce)), sce$free_annotation) # downsample to at most 20 per identity & subset # We are using 5 here to speed up the process but set to 75-100 for your real # life analysis n_cells <- 5 cs_keep <- lapply(idx, function(i) {n <- length(i)if (n < n_cells)n_cells <- nsample(i, n_cells) }) sce <- sce[, unlist(cs_keep)] scetable(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")head(mat <- res$mat)[, seq_len(3)]# Extract NMF model fit mod <- res$NMF sum(mat[2,])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)plotCorrelationMatrix(mat)plotInteractions(mat, which = "heatmap", metric = "prop")ct <- colnames(mat) mat[mat < 0.1] <- 0# Define color palette # (here we use 'paletteMartin' from the 'colorBlindness' package) paletteMartin <- c("#000000", "#004949", "#009292", "#ff6db6", "#ffb6db", "#490092", "#006ddb", "#b66dff", "#6db6ff", "#b6dbff", "#920000", "#924900", "#db6d00", "#24ff24", "#ffff6d")pal <- colorRampPalette(paletteMartin)(length(ct)) names(pal) <- ctplotSpatialScatterpie(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)) saveRDS(mat,"mat.rds") saveRDS(mod,"mod.rds")

矩陣在mat里面 每行加和為1 不同spot的細胞類型比例

總結

以上是生活随笔為你收集整理的空间转录组整合方法SPOTlight的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。