作者,Evil Genius
关于空间转录组的整合分析已经是非常常见了,昨天分享的是python版本,随着数据量不断地增大,R版本的计算效率明显下降,但是作为一种更加具有接受度的语言,R语言的整合分析也是必要的,也是去年空转系列课程的主要整合分析方法。
从算法上来看,R和python的内部逻辑是一样的,分析拿到的是一模一样的结果,如果不一样,说明在处理上产生了问题。
![](https://img-blog.csdnimg.cn/img_convert/442dc27420fd453d987d2638654b6323.webp?x-oss-process=image/format,png)
R脚本封装版如下:
首先准备文件
![](https://img-blog.csdnimg.cn/img_convert/4f0ddda5244804a65e395a861d25c1dc.webp?x-oss-process=image/format,png)
然后进行R版本的整合分析,注意以10X为例
library(argparse)
parser = ArgumentParser(description="Seurat analysis for spatial integration")
parser$add_argument("--inputs", help="input files or directories. Required.eg:sample,file", required=T)
parser$add_argument("--outdir", required=T)
parser$add_argument("--resolution", default = 0.5)
parser$add_argument("--harmony", help="Perform harmony batch effect analysis.", action='store_true')
qc_group$add_argument("--normalize-method", help="normalize method, default: LogNormalize", choices=c('LogNormalize','CLR','RC','SCT'), default='LogNormalize')
args = parser$parse_args()
str(args)
input_files = args$inputs
outdir = args$outdir
resolution = args$resolution
is_harmony = args$harmony
normalize-method = args$normalize-method
dir.create(outdir, recursive = TRUE)
library(grid)
library(dplyr)
library(RColorBrewer)
library(Seurat)
library(foreach)
library(reshape2)
library(ggplot2)
library(psych)
library(pheatmap)
library(reticulate)
library(cowplot)
library(plotly)
library(RColorBrewer)
sample_cols = unique(c(brewer.pal(8,"Dark2"),brewer.pal(9,"Set1"),brewer.pal(8,"Set2"),brewer.pal(12,"Set3")))
defined_cols = c('#e6194b', '#3cb44b', '#ffe119', '#4363d8', '#f58231', '#911eb4', '#46f0f0', '#f032e6', '#bcf60c', '#fabebe', '#008080', '#e6beff', '#9a6324', '#fffac8', '#800000', '#aaffc3', '#808000', '#ffd8b1', '#000075', '#808080', sample_cols)
names(defined_cols) = 1:length(defined_cols) - 1
use_colors = args$defined_colors
samples = read.csv(input_files, header = T)
SP.list = list()
for (i in dim(samples)[1]){
sptmp = Load10X_Spatial(samples[i,2])
sptmp$sample = samples[i,1]
cat("normalize ... \n")
if (normalize_method != "SCT"){
scale_factor = median(sctmp@meta.data[,feature.names.tmp[1]])
print(paste("scale_factor = ", scale_factor))
sptmp = NormalizeData(sctmp, assay='Spatial', normalization.method=normalize-method, scale.factor=scale_factor)
sptmp = FindVariableFeatures(object = sctmp,selection.method = 'vst',assay = 'Spatial',mean.function = ExpMean,dispersion.function = LogVMR,mean.cutoff = c(x_low_cutoff, x_high_cutoff),dispersion.cutoff = c(dispersion_cutoff, Inf),nfeatures = 2000)
}else{
sptmp = SCTransform(sptmp, assay = 'Spatial', verbose = FALSE, variable.features.n = 2000)
}
HVGs = c(HVGs, VariableFeatures(object = sptmp))
SP.list[[samples[i,1]]] = sptmp
}
if(normalize_method == "SCT"){active.assay = "SCT"}
SP.anchors <- FindIntegrationAnchors(object.list = SP.list, dims = 1:20)
Spatial <- IntegrateData(anchorset = SP.anchors, dims = 1:20)
DefaultAssay(Spatial) = active.assay
Spatial$integrated = NULL
VariableFeatures(Spatial) = unique(HVGs)
Spatial = ScaleData(object = Spatial, assay = active.assay)
Spatial = RunPCA(object = Spatial, assay = active.assay, do.print = F, npcs = pc_dim)
if(is_harmony){
library(harmony)
Spatial = RunHarmony(Spatial, group.by.vars = 'sample')
reduction = "harmony"
}
# cluster
cat("cluster ...\n")
Spatial = FindNeighbors(Spatial, reduction = reduction dims = 1:pc.num, annoy.metric = annoy_metric)
Spatial = FindClusters(Spatial, resolution = resolution)
cluster.ids = levels(Spatial@meta.data$seurat_clusters)
# tSNE / UMAP
cat("tsne / umap ...\n")
Spatial = RunTSNE(object = Spatial, reduction = reduction, dims.use = 1:pc.num, do.fast = TRUE)
Spatial = RunUMAP(object = Spatial, reduction = reduction, dims = 1:pc.num, umap.method = "umap-learn",metric = "correlation")
sample = 'combined'
saveRDS(Spatial, file=file.path(outdir,paste(sample,".seurat.rds",sep="")))
umapplot = DimPlot(SCrna, reduction = "umap",label = TRUE, cols=defined_cols)
![](https://img-blog.csdnimg.cn/img_convert/1d77cfa996e6732052e36fbe9c296fd0.webp?x-oss-process=image/format,png)
关于HD数据,截止到这一篇文章发表的时候,没有接到要整合分析的任何需求。
生活很好,有你更好