单细胞空间MIA分析思路及代码分享

作者，Evil Genius

上一篇提到了单细胞空间联合分析的时候，如果有空间转录组数据，但是没有测匹配的单细胞数据，用公共数据库的单细胞数据，那么联合加MIA的分析是非常合适的，今天我们来补充一下MIA的分析内容。

MIA分析可以用来评估空间上某个region或者cluster中富集的细胞类型。需要单细胞和空间转录组两种组学数据。

首先MIA的分析前提
1、单细胞空间数据，两个组学的严格上讲要严格匹配，但实际情况达不到，所以要求组织部位以及疾病类型匹配即可。
2、空间区域，MIA运用的时候空间区域可以人为划分，也可以通过空间聚类的方法获得。
3、空间聚类的结果其实就是匹配不同的组织区域，这个在上一篇文章中得到了印证。其中多个样本的空间数据的整合也是需要去批次的，但是采用CCA还是harmony要看实际效果，因为有空间形态的前提下明显可以看出来去批次的效果怎么样。

从上图可以得到，我们单细胞空间数据，每个聚类结果差异基因匹配的越多，越富含该细胞类型，当然，有一定的检验方法，MIA用到的是超几何检验。

下图为文章的示例图

再强调一次，空间划分区域可以是人工划区域，也可以是空间聚类的结果，归根结底就是不同的组织部位区域

接下来就分享一下MIA的代码，事先准备处理好的单细胞空间的rds

library(Seurat)
library(tidyverse)

cortex_sc <- readRDS(sc_rds)  ####单细胞的rds
cortex_sc <- SCTransform(cortex_sc, ncells = 3000, verbose = FALSE) %>% RunPCA(verbose = FALSE) %>% RunUMAP(dims = 1:20)
cortex_sc <- Seurat::FindNeighbors(cortex_sc, dims = 1:20, verbose = FALSE)
cortex_sc <- Seurat::FindClusters(cortex_sc, verbose = FALSE)

#####添加细胞类型注释信息
#####Barcode，Cluster
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
}

#######寻找每个cluster的差异基因
diff.exp <- FindAllMarkers(object = cortex_sc, only.pos = F, min.pct = 0.25)
diff.exp = diff.exp[,c(7,8,9,6,2,1,5,3,4)]  ## 7, unqiue gene name, removed
diff.sig.exp = subset(diff.exp, avg_logFC> 0.3 & p_val_adj< 0.05)
#diff.sig.exp = subset(diff.exp, abs(avg_logFC)>logfc.diff & p_val_adj<padjust.diff)
sc.marker.exp = subset(diff.exp, avg_logFC> 0.3 & p_val_adj<0.05 & pct.1>0.5 & pct.2<0.5)

空间转录组流程

cortex_sp = readRDS(spatial_rds)
### 找特异基因
sp.diff.exp = FindAllMarkers(
  object = cortex_sp, 
  only.pos = F, 
  min.pct = 0.25,
  test.use = 'wilcox',
  logfc.threshold = 0.2,
  return.thresh = 1
)

sp.diff.exp = sp.diff.exp[,c(7,8,9,6,2,1,5,3,4)]  ## 7, unqiue gene name, removed
sp.diff.sig.exp = subset(sp.diff.exp, avg_logFC> 0.3 & p_val_adj< 0.05)
#diff.sig.exp = subset(diff.exp, abs(avg_logFC)>logfc.diff & p_val_adj<padjust.diff)
sp.marker.exp = subset(sp.diff.exp, avg_logFC> 0.3 & p_val_adj<0.05 & pct.1>0.5 & pct.2<0.5)

1.上述找两个DEG数据框的方法不唯一，阈值也不唯一,但是经过实际测试，结果差别不大。
2.第二个DEG数据框也可以是空间cluster的marker
3.MIA分析模式在单细胞和空间转录组场景都可以应用，空转场景是看细胞亚群的富集程度，单细胞场景是做细胞亚群注释

MIA过程