python 单细胞 批次矫正 单细胞 批次效应

此文章是通过学习瑞典国家生物信息学基础设施(NBIS) 所开放的单细胞分析教程加上网上所查找的资料，自身的理解所形成的，可能会有不足之处。该部分是整合不同批次的单细胞数据并矫正批次效应。

整合方法

目前，常见的整合方法一共有四种，具体如下：

Markdown	Language	Library	Ref
CCA	R	Seurat	Cell
MNN	R/Python	Scater/Scanpy	Nat. Biotech.
Conos	R	conos	Nat. Methods
Scanorama	Python	scanorama	Nat. Biotech.

首先，我们把上一步的得到的数据读进R并且根据样本名将数据拆分成单个数据，然后对基因表达量进行标准化，并根据方差稳定性转换（variance stabilizing transformation ,vst）识别每个样本的高变基因。

suppressPackageStartupMessages({
  library(Seurat)
  library(cowplot)
  library(ggplot2)
})
alldata <- readRDS("data/results/covid_qc_dr.rds")

alldata.list <- SplitObject(alldata, split.by = "orig.ident")

for (i in 1:length(alldata.list)) {
    alldata.list[[i]] <- NormalizeData(alldata.list[[i]], verbose = FALSE)
    alldata.list[[i]] <- FindVariableFeatures(alldata.list[[i]], selection.method = "vst",
        nfeatures = 2000, verbose = FALSE)
}

hvgs_per_dataset <- lapply(alldata.list, function(x) {
    x@assays$[e-mail]
})
# venn::venn(hvgs_per_dataset,opacity = .4,zcolor = scales::hue_pal()(3),cexsn
# = 1,cexil = 1,lwd=1,col='white',frame=F,borders = NA)

temp <- unique(unlist(hvgs_per_dataset))
overlap <- sapply(hvgs_per_dataset, function(x) {
    temp %in% x
})
pheatmap::pheatmap(t(overlap * 1), cluster_rows = F, color = c("grey90", "grey20"))

Seurat4

数据整合

使用Seurat4所自带的函数对不同批次的样本进行整合,是一种基于anchor的方法。

# 首先需要找到整合数据使，作为参考的anchor
alldata.anchors <- FindIntegrationAnchors(object.list = alldata.list, dims = 1:30,
    reduction = "cca")

# 根据找到的anchor整合不同批次的单细胞数据
alldata.int <- IntegrateData(anchorset = alldata.anchors, dims = 1:30, new.assay.name = "CCA")
# 发现，现在以CCA的名字建立了一个新的基因表达矩阵，只不过现在用来计算的默认基因集由RNA变成了新的CCA
names(alldata.int@assays)
# [1] "RNA" "CCA"

[e-mail]
# [1] "CCA"

聚类差异

在进行整合之后，原始的基因表达举证依旧存在于原来的list中，所以这个时候可以比较整合前后的降维差异。

# 在新的整合数据上进行降维处理
alldata.int <- ScaleData(alldata.int, verbose = FALSE)
alldata.int <- RunPCA(alldata.int, npcs = 30, verbose = FALSE)
alldata.int <- RunUMAP(alldata.int, dims = 1:30)
alldata.int <- RunTSNE(alldata.int, dims = 1:30)

# 绘图显示差异
plot_grid(ncol = 3,
  DimPlot(alldata, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA raw_data"),
  DimPlot(alldata, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE raw_data"),
  DimPlot(alldata, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP raw_data"),
  
  DimPlot(alldata.int, reduction = "pca", group.by = "orig.ident")+NoAxes()+ggtitle("PCA integrated"),
  DimPlot(alldata.int, reduction = "tsne", group.by = "orig.ident")+NoAxes()+ggtitle("tSNE integrated"),
  DimPlot(alldata.int, reduction = "umap", group.by = "orig.ident")+NoAxes()+ggtitle("UMAP integrated")
)

PS:可以从tSNE上看到，批次效应已经得到了很好的矫正

marker基因展示

Markers	Cell Type
CD3E	T cells
CD3E CD4	CD4+ T cells
CD3E CD8A	CD8+ T cells
GNLY, NKG7	NK cells
MS4A1	B cells
CD14, LYZ, CST3, MS4A7	CD14+ Monocytes
FCGR3A, LYZ, CST3, MS4A7	FCGR3A+ Monocytes
FCER1A, CST3	DCs

myfeatures <- c("CD3E", "CD4", "CD8A", "NKG7", "GNLY", "MS4A1", "CD14", "LYZ", "MS4A7",
    "FCGR3A", "CST3", "FCER1A")
plot_list <- list()
for (i in myfeatures) {
    plot_list[[i]] <- FeaturePlot(alldata, reduction = "umap", dims = 1:2, features = i,
        ncol = 3, order = T) + NoLegend() + NoAxes() + NoGrid()
}
plot_grid(ncol = 3, plotlist = plot_list)

Harmony

Harmony是另外一个可以用来进行数据整合的R包。

# install.packages("harmony")
library(harmony)

# 不需要对数据进行拆分，可以直接整合
alldata.harmony <- RunHarmony(alldata, group.by.vars = "orig.ident", reduction = "pca",
    dims.use = 1:50, assay.use = "RNA")

# 使用Harmony计算得到的主成分进行UMAP二次降维
alldata.int[["harmony"]] <- alldata.harmony[["harmony"]]
alldata.int <- RunUMAP(alldata.int, dims = 1:50, reduction = "harmony", reduction.name = "umap_harmony")

# 展示了整合后每个样本中的特征基因的数量与细胞数量
hvgs <- unique(unlist(hvgs_per_dataset))

assaylist <- list()
genelist <- list()
for (i in 1:length(alldata.list)) {
    assaylist[[i]] <- t(as.matrix(GetAssayData(alldata.list[[i]], "data")[hvgs, ]))
    genelist[[i]] <- hvgs
}

lapply(assaylist, dim)

还有一种基于python的整合方法-scanorama，但不知到为什么，我在R中无法导入这个包，暂时也没有找到解决办法。就我目前解决单细胞批次效应，两个包已经足以，后面如果再有时间，在来继续研究这个包。

保存数据

像之前一样，保存当下的数据：

saveRDS(alldata.int, "data/results/covid_qc_dr_int.rds")