作者，Evil Genius

上一篇提到了單細胞空間聯(lián)合分析的時候，如果有空間轉(zhuǎn)錄組數(shù)據(jù)，但是沒有測匹配的單細胞數(shù)據(jù)，用公共數(shù)據(jù)庫的單細胞數(shù)據(jù)，那么聯(lián)合加MIA的分析是非常合適的，今天我們來補充一下MIA的分析內(nèi)容。

MIA分析可以用來評估空間上某個region或者cluster中富集的細胞類型。需要單細胞和空間轉(zhuǎn)錄組兩種組學數(shù)據(jù)。

首先MIA的分析前提
1、單細胞空間數(shù)據(jù)，兩個組學的嚴格上講要嚴格匹配，但實際情況達不到，所以要求組織部位以及疾病類型匹配即可。
2、空間區(qū)域，MIA運用的時候空間區(qū)域可以人為劃分，也可以通過空間聚類的方法獲得。
3、空間聚類的結(jié)果其實就是匹配不同的組織區(qū)域，這個在上一篇文章中得到了印證。其中多個樣本的空間數(shù)據(jù)的整合也是需要去批次的，但是采用CCA還是harmony要看實際效果，因為有空間形態(tài)的前提下明顯可以看出來去批次的效果怎么樣。

從上圖可以得到，我們單細胞空間數(shù)據(jù)，每個聚類結(jié)果差異基因匹配的越多，越富含該細胞類型，當然，有一定的檢驗方法，MIA用到的是超幾何檢驗。

下圖為文章的示例圖

再強調(diào)一次，空間劃分區(qū)域可以是人工劃區(qū)域，也可以是空間聚類的結(jié)果，歸根結(jié)底就是不同的組織部位區(qū)域

接下來就分享一下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)

空間轉(zhuǎn)錄組流程

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數(shù)據(jù)框的方法不唯一，閾值也不唯一,但是經(jīng)過實際測試，結(jié)果差別不大。
2.第二個DEG數(shù)據(jù)框也可以是空間cluster的marker
3.MIA分析模式在單細胞和空間轉(zhuǎn)錄組場景都可以應(yīng)用，空轉(zhuǎn)場景是看細胞亞群的富集程度，單細胞場景是做細胞亞群注釋

MIA過程