批次校正-单细胞转录组分析批次效应的校正scRNA-02

背景：2013年到2023年是单细胞转录组分析的第一个十年。批次校正在单细胞转录组分析中至关重要，根据不同样本的基因表达模式可以分析样本与样本之间的生物学差异。然而基因表达模式可能会因为实验方法、技术方法等原因引入偶然误差，例如明明是同一个细胞类型的细胞，却在不同的样本间无法聚类在一起。

基因表达模式简介
什么是单细胞转录组分析的基因表达模式呢？举例来说，如果我们要描述不同的人，可以通过身高，体重，肤色等特征来说明：

同学A [身高：180，体重：50，肤色：白色...]
同学B [身高：170，体重：55，肤色：黄色...]
同学C [身高：175，体重:60，肤色：白色...]
同学D [身高：165，体重：50，肤色：黑色...]

那么我们对这些同学进行分类（单细胞转录组细胞分群），就是根据他们的这些特征情况，计算每个同学之间的相似性，比如肤色相同的同学更相近，如果在肤色相同的情况下，他们的体重也相近，那么他们的关系就更加接近...对于单细胞转录组分析来说，每个细胞就是每个同学，每个特征就是每个基因。

于是，构建同学的信息矩阵如下：

#'同学的信息矩阵
        同学A 同学B 同学C 同学D
身高   180    170   175  165
体重    50    55    60    50
肤色   白色   黄色  白色  黑色
...
#'单细胞转录组的counts矩阵
        细胞A 细胞B 细胞C 细胞D
基因1   0    1    1    0
基因2   1    0    0    1
基因3   1   1     0    0
...

这其实是和单细胞转录组的原始counts矩阵类似，那么单细胞转录组后续的分析就是算法在这种矩阵信息的基础上，根据每个细胞的所有基因表达情况对细胞进行基因表达模式的分析，将相似的基因表达模式的细胞聚类在一起。

而批次效应，就是因为每个批次的实验方法等不同造成同一个类型的细胞的某些（非细胞类型的关键基因）基因表达量存在较大的差异，使得在聚类可视化的时候，明明是同一个类型的细胞，却因为批次没有聚在一起。例如下图红框中的细胞，根据细胞类型注释它们都是CD14+细胞，然而在最左边的可视化图中，却因为batch不同形成了分群。

batch.png

样本的批次无法消除，要数据量多就需要尽可能使用不同实验室、不同机构测出来的数据。去除批次效应，只能通过各种算法，重新计算细胞之间的相似度（根据细胞类型的关键基因），将同一个类型的细胞的基因表达量校正回来。例如在肿瘤分析中，尽量减少不同病人之间同一细胞类型的误差。

类型一：使用R包Seurat方法

library(Seurat)
library(dplyr)
sc<-readRDS("sc.merge.rds")
data<-NormalizeData(sc)%>% FindVariableFeatures(nfeatures=2000) %>% 
      ScaleData() %>% RunPCA() %>% FindNeighbors(dims = 1:30) %>%
      RunUMAP(dims = 1:30) %>% FindClusters(resolution=0.05)
saveRDS(data,"sc.rds")
#'未去除批次时候降维聚类

seurat_clusters.png

samples.png

方法一 Seurat-integration

library(Seurat)
sc<-"" #'Seurat对象
split.by<-"" #'需要去除批次的批次标签名字
method<-"" #'标准化的方法 log或者sct或者harmony
nfeatures<-3500 #'高变基因保留的数目，默认3500
dims<-30 #'降维的维度
resolution<-0.5 #'聚类分辨率
sc.list<-SplitObject(sc, split.by = split.by)
if (method == "log") {
  sc.list = lapply(sc.list, function(x) {
    x = NormalizeData(x) %>% FindVariableFeatures(nfeatures = nfeatures) %>% 
      ScaleData()
    return(x)
  })
  sc.anchors <- FindIntegrationAnchors(object.list = sc.list, 
                                        dims = 1:30)
  sc.combined <- IntegrateData(anchorset = sc.anchors, 
                                dims = 1:30)
  DefaultAssay(sc.combined) <- "integrated"
  sc.combined <- sc.combined %>% ScaleData() %>% 
    RunPCA() %>% FindNeighbors(dims = 1:dims) %>% RunUMAP(dims = 1:dims) %>% 
    FindClusters(resolution = resolution)
  data <- sc.combined
}
pdf("umap_seurat_clusters.pdf",height=7,width=10)
DimPlot(data,reduction = "umap",group.by="seurat_clusters")
dev.off()
pdf("umap_groups.pdf",height=7,width=10)
DimPlot(data,reduction = "umap",group.by="groups")
dev.off()
pdf("umap_samples.pdf",height=7,width=10)
DimPlot(data,reduction = "umap",group.by="orig.ident")
dev.off()
#'data就是去除批次后的数据

integrated_seurat_clusters.png

integrated_samples.png

方法二 Seurat-SCTransform

library(Seurat)
sc<-"" #'Seurat对象
split.by<-"" #'需要去除批次的批次标签名字
method<-"" #'标准化的方法 log或者sct或者harmony
nfeatures<-3500 #'高变基因保留的数目，默认3500
dims<-30 #'降维的维度
resolution<-0.5 #'聚类分辨率
sc.list<-SplitObject(sc, split.by = split.by)
if (method == "sct") {
sc.list <- lapply(X = sc.list, FUN = function(x){
    x <- SCTransform(x,assay = Assays(x)[1])
    return(x)
  })
  features <- SelectIntegrationFeatures(object.list = sc.list,nfeatures=nfeatures)
  sc.list <- PrepSCTIntegration(object.list = sc.list, anchor.features = features) 
  sc.list <- FindIntegrationAnchors(object.list = sc.list, normalization.method = "SCT",anchor.features = features)
  sc.combined <- IntegrateData(anchorset = sc.list, normalization.method = "SCT")
  sc.combined <- RunPCA(sc.combined, verbose = FALSE)
  sc.combined <- RunUMAP(sc.combined, reduction = "pca", dims = 1:dims)
  sc.combined <- FindNeighbors(sc.combined, reduction = "pca", dims = 1:dims)
  sc.combined <- FindClusters(sc.combined, resolution = resolution)
  data<-sc.combined
}

SCT_seurat_clusters.png

SCT_samples.png

方法三 Seurat-harmony

library(Seurat)
library(harmony)
sc<-"" #'Seurat对象
split.by<-"" #'需要去除批次的批次标签名字
method<-"" #'标准化的方法 log或者sct或者harmony
nfeatures<-3500 #'高变基因保留的数目，默认3500
dims<-30 #'降维的维度
resolution<-0.5 #'聚类分辨率
sc.list<-SplitObject(sc, split.by = split.by)
if(method=="harmony"){
  sc <- NormalizeData(sc, normalization.method = "LogNormalize", scale.factor = 10000)
  sc <- FindVariableFeatures(sc, selection.method = "vst", nfeatures = nfeatures)
  all.genes <- rownames(sc) #'ScaleData使用的基因，否则后面可能报错
  sc <- ScaleData(sc , features = all.genes, vars.to.regress = "nCount_RNA")
  sc <- RunPCA(sc, features = VariableFeatures(object = sc))
  sc <- RunHarmony(sc, split.by , plot_convergence = F,dims.use = 1:dims)
  sc.combine <- sc
  sc.combine = RunTSNE(sc.combine, reduction = "harmony", dims = 1:dims)
  sc.combine = RunUMAP(sc.combine, reduction = "harmony", dims = 1:dims)
  sc.combine = FindNeighbors(sc.combine, reduction = "harmony",dims = 1:dims)
  sc.combine = FindClusters(sc.combine, resolution = resolution)
  data<-sc.combine
}

harmony_seurat_clusters.png

samples.png

类型二：使用python中的包（待整理）

安装需要的python包：

conda install -c conda-forge scanpy
pip install scib
conda install -c conda-forge omicverse
conda install -c bioconda scvi

首先需要将Seurat对象转换为.h5ad文件：

library(SeuratDisk)
library(Seurat)
sc<-readRDS("sc.rds") #'读入RDS文件
SaveH5Seurat(sc,filename="sc.h5seurat", overwrite = TRUE)
Convert("sc.h5seurat", dest = "h5ad", overwrite = TRUE) #'转换为当前路径下的.h5ad文件

未去批次时的降维聚类可视化：

import omicverse as ov
import scanpy as sc
adata=sc.read('sc.h5ad')
adata=ov.pp.qc(adata,
              tresh={'mito_perc': 0.2, 'nUMIs': 500, 'detected_genes': 250},
              batch_key='groups') #'质控
ov.utils.store_layers(adata,layers='counts')
adata=ov.pp.preprocess(adata,mode='shiftlog|pearson',
                       n_HVGs=2000,batch_key='groups') #'高变基因
#adata
adata.raw = adata
adata = adata[:, adata.var.highly_variable_features]
ov.pp.scale(adata)
ov.pp.pca(adata,layer='scaled',n_pcs=30) #PC选择30个维度
adata.obsm["X_mde_pca"] = ov.utils.mde(adata.obsm["scaled|original|X_pca"])
ov.utils.embedding(adata,
                basis='X_mde_pca',frameon='small',
                color=['groups','seurat_clusters'],show=False)
#'未去批次时的降维聚类

方法一 combat

adata_combat=ov.single.batch_correction(adata,batch_key='groups',
                                        methods='combat',n_pcs=30)
adata.obsm["X_mde_combat"] = ov.utils.mde(adata.obsm["X_combat"])
ov.utils.embedding(adata,
                basis='X_mde_combat',frameon='small',
                color=['groups','seurat_clusters'],show=False)

方法二 scanorama

adata_scanorama=ov.single.batch_correction(adata,batch_key='groups',
                                        methods='scanorama',n_pcs=30)
adata.obsm["X_mde_scanorama"] = ov.utils.mde(adata.obsm["X_scanorama"])
ov.utils.embedding(adata,
                basis='X_mde_scanorama',frameon='small',
                color=['groups','seurat_clusters'],show=False)

总结：批次校正是一个相当复杂的环节，对不同的样本数据有着不同的方法应用效果，具体使用哪个方法需要根据校正结果决定。此外，批次校正方法众多，还有很多新方法可以尝试。