MIA脚本封装版

作者，Evil Genius

我们利用单细胞空间转录组数据联合分析，MIA的方法实现如下结果

就是我们需要知道空间各个区域富集了哪些细胞类型，封装版代码更新如下：

#! usr/bin/R
### zhaoyunfei
### 20230612

suppressMessages({
library(Seurat)
library(argparse)
library(dplyr)
library(ggplot2)
library(tidyverse)
})

parser = ArgumentParser()
parser$add_argument("--sc_rds", help="the sc data",required = T)
parser$add_argument("--spatial_rds", help="the sp data",required = T)
parser$add_argument("--sample", help="the sample name",required = T)
parser$add_argument("--outdir", help="the outdir",default = './')
parser$add_argument("--celltype", help="the annotation for celltype")
parser$add_argument("--normalization_method",default = 'SCT',choice = c('SCT','LogNormalize'))
parser$add_argument("-s","--sc_min_pct", help="the min pct for sc diff test",default = 0.3)
parser$add_argument("--sc_logfc", help="the threshold of logfc for sc diff test",default = 0.25)
parser$add_argument("--region", help="the spatial region annotation,eg:Barcode,Cluster")
parser$add_argument("--sp_min_pct", help="the min pct for sp diff test",default = 0.3)
parser$add_argument("--sp_logfc", help="the threshold of logfc for sp diff test",default = 0.25)
args <- parser$parse_args()
str(args)

sc_rds = args$sc_rds
spatial_rds = args$spatial_rds
outdir = args$outdir
celltype = args$celltype
normalization_method = args$normalization_method
sample = args$sample
sc_min = args$sc_min_pct
sc_logfc = args$sc_logfc
region = args$region
sp.pct = args$sp_min_pct
sp.logfc = args$sp_logfc


if (!file.exists(outdir)) {dir.create(outdir,recursive = TRUE)}


if (strsplit(sc_rds,'[.]')[[1]][length(strsplit(sc_rds,'[.]')[[1]])] == 'rds'){

cortex_sc <- readRDS(sc_rds)

if (normalization_method == 'SCT'){

cortex_sc <- SCTransform(cortex_sc, ncells = 3000, verbose = FALSE) %>% RunPCA(verbose = FALSE) %>% RunUMAP(dims = 1:20)}else{

cortex_sc <- NormalizeData(cortex_sc, normalization.method = "LogNormalize", scale.factor = 10000)

cortex_sc <- FindVariableFeatures(cortex_sc, selection.method = "vst", nfeatures = 2000)

cortex_sc <- ScaleData(cortex_sc, features = rownames(cortex_sc))

cortex_sc <- RunPCA(cortex_sc, features = VariableFeatures(object = cortex_sc))

cortex_sc <- RunUMAP(cortex_sc, dims = 1:20)

}

cortex_sc <- Seurat::FindNeighbors(cortex_sc, dims = 1:20, verbose = FALSE)

cortex_sc <- Seurat::FindClusters(cortex_sc, verbose = FALSE)

} else if (strsplit(sc_rds,'[.]')[[1]][length(strsplit(sc_rds,'[.]')[[1]])] == 'csv'){

gene_bar <- read.csv(sc_rds,header=T,row.names=1,sep=',',check.names = F)

cortex_sc <- CreateSeuratObject(counts = gene_bar, min.cells  = 3, project = 'sc')

cortex_sc <- Seurat::SCTransform(cortex_sc, verbose = FALSE , ncells = 3000)

cortex_sc <- Seurat::RunPCA(cortex_sc, verbose = FALSE)

cortex_sc <- Seurat::RunUMAP(cortex_sc, dims = 1:20, verbose = FALSE)

cortex_sc <- Seurat::FindNeighbors(cortex_sc, dims = 1:20, verbose = FALSE)

cortex_sc <- Seurat::FindClusters(cortex_sc, verbose = FALSE)
}


if ( ! is.null(celltype) ) {

anno = read.csv(celltype,header = T)

if (length(anno$Cluster) != length(rownames(cortex_sc))){print("warning ~~~ ,The length of celltype file is not match the sc data,subset will be acrry out ~~~")}

index = na.omit(match(colnames(cortex_sc),anno$Barcode))

cortex_sc = cortex_sc[,index]

Seurat::Idents(object = cortex_sc) <- anno$Cluster

cortex_sc$seurat_clusters = anno$Cluster

if (! is.null(anno$Sample)){cortex_sc$orig.ident = anno$Sample}

}else {Seurat::Idents(object = cortex_sc) <- cortex_sc$seurat_clusters

}

DefaultAssay(cortex_sc) = 'RNA'

sc.markers <- FindAllMarkers(cortex_sc, only.pos = TRUE, min.pct = sc_min, logfc.threshold = sc.logfc)

sc.markers$d = sc.markers$pct.1 - sc.markers$pct.2

sc.main.marker = subset(sc.markers , avg_log2FC > 0.25 & p_val_adj < 0.05 & d > 0.2)

sc.main.marker = sc.main.marker %>% arrange(cluster,desc(avg_log2FC))

sc.main.marker = as.data.frame(sc.main.marker)

sc.main.marker$cluster = paste('sc',sc.main.marker$cluster,sep = '_')
#####空间转录组处理

cortex_sp = readRDS(spatial_rds)

cortex_sp <- SCTransform(cortex_sp, verbose = FALSE)

cortex_sp <- RunPCA(cortex_sp, assay = "SCT", verbose = FALSE)

cortex_sp <- FindNeighbors(cortex_sp, reduction = "pca", dims = 1:20)

cortex_sp <- FindClusters(cortex_sp, verbose = FALSE)

cortex_sp <- RunUMAP(cortex_sp, reduction = "pca", dims = 1:20)

if (!is.null(region)){

sp_anno = read.csv(region,header = T)

cortex_sp = cortex_sp[,sp_anno$Barcode]

Idents(cortex_sp) = sp_anno$Cluster

cortex_sp$seurat_clusters = sp_anno$Cluster

}

region_marker=FindAllMarkers(cortex_sp,logfc.threshold = sp.logfc,only.pos = T,min.pct = sp.pct)

region_marker$d=region_marker$pct.1 - region_marker$pct.2

region_main_marker = subset(region_marker,avg_log2FC > 0.25 & p_val_adj < 0.05 & d > 0.1)

region_main_marker = region_main_marker %>% arrange(cluster,desc(avg_log2FC))

region_main_marker = as.data.frame(region_main_marker)

region_main_marker$cluster = paste('spatial',region_main_marker$cluster,sep = '_')

#####MIA

region_specific = region_main_marker[,c("cluster","gene")]

colnames(region_specific)[1]="region"

celltype_specific = sc.main.marker[,c("cluster","gene")]

colnames(celltype_specific)[1]="celltype"

source("MIA.R")

Result = zhao_MIA(region_specific,celltype_specific,sample,outdir)

source脚本如下：