方法更新
作者,Evil Genius
今天我们分享一个分析方法思路。
在HD 16um(8um通常不采用)、Stereo-seq bin50以下(25um以下),在做单细胞空间联合分析的时候,采用的是RCTD的double模式,那么在这种模式下,细胞类型之间的空间关系该如何分析呢?
这内容本来放在2025培训班上讲的,不过既然粉丝们要求了,今天就分享出来。
其实在很早之前(2022年)进行空间培训的时候谈到过这个问题,大家可以看一下视频空间转录组高级个性化数据分析第一期
相信大家也看了很多的文献,关于double模式下的分析其实参考很少,Stereo-seq发表的文章大多也在bin60以上,采用更多的类似full模式的分析。
今天我提一篇文献,很有参考价值,很早之前分享过,不过我发现认真好好学习的人,真的很少。
其实最关键的地方,大家看一下
采用的是bin20进行分析(10-15um),接近HD的16um精度。
那么注释的方式是什么呢?简而言之,pattern模式。
就是在double的模式下,将注释分为两层,first pattern和second pattern。
那么在这种思路下,空间细胞类型的关系问题就演化成了两层pattern的细胞关系问题。
那么细胞邻域富集,细胞邻域的差异分析等均可以采用squidpy进行分析了,不过个人推荐cellcharter。
其中最为关键的是,此时的空间配受体分析直接与细胞类型相关联。
分享到这里,大家应该知道怎么分析double模式了吧。
最后写上封装好的RCTD的代码,其实脚本封装大家一定要学会。
代码语言:javascript代码运行次数:0运行复制# 创建解析器
parser <- ArgumentParser(description = "RCTD analysis")
# 添加参数
# 创建解析器
parser <- ArgumentParser(description = "RCTD analysis")
# 添加参数
parser$add_argument("-s", "--scrna", required = TRUE, help = "scrna file")
parser$add_argument("-t", "--type", required = TRUE, help = "scrna file type",choices = c('csv','rds'))
parser$add_argument("-a", "--anno",help = "scrna file type;csv")
parser$add_argument("-p", "--spatial", required = TRUE, help = "spatial rds file")
parser$add_argument("-o", "--output", required = TRUE, help = "Output directory for saving results")
parser$add_argument("-m", "--mode", choices = c("doublet", "full","multi"), required = TRUE, help = "RCTD mode")
parser$add_argument("-n", "--name", required = TRUE, help = "sample name")
scrna = args$scrna
type = args$type
anno = args$anno
spatial = args$spatial
outdir = args$output
mode = args$mode
sample = args$name
完整脚本如下:
代码语言:javascript代码运行次数:0运行复制#! usr/bin/R
####zhaoyunfei
####20241202
library(argparse)
library(Seurat)
library(ggplot2)
library(patchwork)
library(dplyr)
library(spacexr)
library(ggplot2)
# 创建解析器
parser <- ArgumentParser(description = "RCTD and CC analysis")
# 添加参数
parser$add_argument("-s", "--scrna", required = TRUE, help = "scrna file")
parser$add_argument("-t", "--type", required = TRUE, help = "scrna file type",choices = c('csv','rds'))
parser$add_argument("-a", "--anno",help = "scrna file type;csv")
parser$add_argument("-p", "--spatial", required = TRUE, help = "spatial rds file")
parser$add_argument("-o", "--output", required = TRUE, help = "Output directory for saving results")
parser$add_argument("-m", "--mode", choices = c("doublet", "full","multi"), required = TRUE, help = "RCTD mode")
parser$add_argument("-n", "--name", required = TRUE, help = "sample name")
# 解析命令行参数
args <- parser$parse_args()
scrna = args$scrna
type = args$type
anno = args$anno
spatial = args$spatial
outdir = args$output
mode = args$mode
sample = args$name
if(! file.exists(outdir)){dir.create(outdir,recursive = TRUE)}
if (type == 'csv'){
ref = read.table(scrna,header = T,row.names = 1,check.names = F,sep = '\t')
counts <- ref
if (!is.null(anno)){
anno = read.table(anno,header = T,sep = '\t',row.names = 1)
anno = anno[colnames(counts),]
cluster <- as.factor(anno$CellType)
names(cluster) <- colnames(ref)
}else{stop("错误:必须提供 --anno 参数", call. = FALSE)}
nUMI <- colSums(ref)
names(nUMI) <- colnames(ref)
reference <- Reference(counts, cluster, nUMI)
}else{ref = readRDS(scrna)
counts <- ref@assays$RNA@counts
cluster = ref$celltype
names(cluster) <- colnames(ref)
nUMI <- colSums(counts)
names(nUMI) <- colnames(ref)
reference <- Reference(counts, cluster, nUMI)
}
spatial = Load10X_Spatial(data.dir = spatial, bin.size = c(16))
spatial <- subset(spatial, subset = nFeature_Spatial.016um > 200)
counts <- spatial@assays$Spatial.016um@layers$counts
colnames(counts) = colnames(spatial)
rownames(counts) = rownames(spatial)
coords <- GetTissueCoordinates(spatial)
coords = coords[,c(1,2)]
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 = mode)
spatial@meta.data$first_type = RCTD@results$results_df$first_type
spatial@meta.data$second_type = RCTD@results$results_df$second_type
pdf(paste(outdir,paste(sample,'RCTD.first_type.pdf',sep = '.'),sep = '/'))
SpatialDimPlot(spatial, group.by = "first_type", label = T, repel = T, label.size = 4)
dev.off()
pdf(paste(outdir,paste(sample,'RCTD.second_type.pdf',sep = '.'),sep = '/'))
SpatialDimPlot(spatial, group.by = "second_type", label = T, repel = T, label.size = 4)
dev.off()
saveRDS(spatial,file = paste(outdir,paste(sample,'HD.RCTD.rds',sep = '.'),sep = '/'))
生活很好,有你更好。
发布者:admin,转转请注明出处:http://www.yc00.com/web/1747699877a4682868.html
评论列表(0条)