课前准备----空间转录组R版本的harmony整合（封装版）-CSDN博客

本文链接：https://blog.csdn.net/weixin_53637133/article/details/140210918

作者，Evil Genius

关于空间转录组的整合分析已经是非常常见了，昨天分享的是python版本，随着数据量不断地增大，R版本的计算效率明显下降，但是作为一种更加具有接受度的语言，R语言的整合分析也是必要的，也是去年空转系列课程的主要整合分析方法。

从算法上来看，R和python的内部逻辑是一样的，分析拿到的是一模一样的结果，如果不一样，说明在处理上产生了问题。

R脚本封装版如下：

首先准备文件

然后进行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)