8.4 真实数据集的整合实例
有很多关于整合标准的文章发表，最详细的文章（Tran，2020）使用不同大小和复杂程度的多个模拟和真实数据集比较了14种scRNA-seq数据集整合方法。根据这些文章的测试，Harmony、LIGER和Seurat（v3）的表现最佳。我们将在两个任务中说明这三种方法的性能：1）整合高度相关的3’和5’ PBMC数据集；2）仅整合部分重叠的数据集，即全血（包括红细胞和中性粒细胞）和3’ PBMC数据集。
让我们加载所有必要的库：
> library(Seurat)
> library(SeuratDisk)
> library(SeuratWrappers)
> library(patchwork)
> library(harmony)
> library(rliger)
> library(reshape2)
> library(RColorBrewer)
> library(dplyr)
另外，让我们导入一个自定义函数（https://github.com/cellgeni/scRNA.seq.course/blob/master/course_files/utils/custom_seurat_functions.R），以可视化每个cluster中不同数据集的细胞分布以及cluster大小：
> source("utils/custom_seurat_functions.R")
8.5 Seurat v3, 3’ vs 5’ 10k PBMC
加载从10x Genomics数据库下载的经过过滤的Cell Ranger h5矩阵，也可以使用以下命令下载：
download.file("https://cf.10xgenomics.com/samples/cell-exp/4.0.0/Parent_NGSC3_DI_PBMC/Parent_NGSC3_DI_PBMC_filtered_feature_bc_matrix.h5",
destfile = "3p_pbmc10k_filt.h5")
download.file("https://cf.10xgenomics.com/samples/cell-vdj/5.0.0/sc5p_v2_hs_PBMC_10k/sc5p_v2_hs_PBMC_10k_filtered_feature_bc_matrix.h5",
destfile = "5p_pbmc10k_filt.h5")
读取数据并创建Seurat对象。请注意，5’数据集另有VDJ数据。
> matrix_3p <- Read10X_h5("data/update/3p_pbmc10k_filt.h5",use.names = T)
> matrix_5p <- Read10X_h5("data/update/5p_pbmc10k_filt.h5",use.names = T)$`Gene Expression`
> srat_3p <- CreateSeuratObject(matrix_3p,project = "pbmc10k_3p")
> srat_5p <- CreateSeuratObject(matrix_5p,project = "pbmc10k_5p")
删除矩阵以节省内存：
> rm(matrix_3p)
> rm(matrix_5p)
计算线粒体基因和核糖体蛋白的比例，并对数据集进行快速过滤：
> srat_3p[["percent.mt"]] <- PercentageFeatureSet(srat_3p, pattern = "^MT-")
> srat_3p[["percent.rbp"]] <- PercentageFeatureSet(srat_3p, pattern = "^RP[SL]")
> srat_5p[["percent.mt"]] <- PercentageFeatureSet(srat_5p, pattern = "^MT-")
> srat_5p[["percent.rbp"]] <- PercentageFeatureSet(srat_5p, pattern = "^RP[SL]")
> VlnPlot(srat_3p, features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rbp"), ncol = 4)
> VlnPlot(srat_5p, features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rbp"), ncol = 4)
现在，让我们比较基因名称，结果发现它们是相同的：都使用了最新的Cell Ranger注释GRCh38-2020A。
> table(rownames(srat_3p) %in% rownames(srat_5p))
TRUE
36601
快速过滤数据集去除死细胞和双胞：
> srat_3p <- subset(srat_3p, subset = nFeature_RNA > 500 & nFeature_RNA < 5000 & percent.mt < 15)
> srat_5p <- subset(srat_5p, subset = nFeature_RNA > 500 & nFeature_RNA < 5000 & percent.mt < 10)
按照Seurat示例进行整合。为此，我们需要为这两个对象创建一个简单的R列表，并为每个对象标准化/查找HVG：
> pbmc_list <- list()
> pbmc_list[["pbmc10k_3p"]] <- srat_3p
> pbmc_list[["pbmc10k_5p"]] <- srat_5p
> for (i in 1:length(pbmc_list)) {
pbmc_list[[i]] <- NormalizeData(pbmc_list[[i]], verbose = F)
pbmc_list[[i]] <- FindVariableFeatures(pbmc_list[[i]],
selection.method = "vst",
nfeatures = 2000,
verbose = F)
}
使用以下两个Seurat命令来查找整合锚点并执行整合。这大约需要10分钟：
> pbmc_anchors <- FindIntegrationAnchors(object.list = pbmc_list, dims = 1:30)
> pbmc_seurat <- IntegrateData(anchorset = pbmc_anchors, dims = 1:30)
删除所有未使用的数据结构来节省RAM：
> rm(pbmc_list)
> rm(pbmc_anchors)
Seurat整合创建了一个统一的对象，其中包含原始数据（“RNA” assay）和整合数据（“integrated” assay）。让我们将assay设置为RNA，可视化整合前的数据集。
> DefaultAssay(pbmc_seurat) <- "RNA"
对未整合（RNA）的数据进行标准化、HVG查找、Scale、PCA和UMAP：
> pbmc_seurat <- NormalizeData(pbmc_seurat, verbose = F)
> pbmc_seurat <- FindVariableFeatures(pbmc_seurat, selection.method = "vst", nfeatures = 2000, verbose = F)
> pbmc_seurat <- ScaleData(pbmc_seurat, verbose = F)
> pbmc_seurat <- RunPCA(pbmc_seurat, npcs = 30, verbose = F)
> pbmc_seurat <- RunUMAP(pbmc_seurat, reduction = "pca", dims = 1:30, verbose = F)
整合之前数据集的UMAP图显示出明显的分离。可以将patchwork语法与Seurat绘图函数一起使用：
> DimPlot(pbmc_seurat,reduction = "umap") +
plot_annotation(title = "10k 3' PBMC and 10k 5' PBMC cells, before integration")
现在让我们将assay改为“integrated”，并进行同样的处理（它已经标准化并且筛选了HVG）：
> DefaultAssay(pbmc_seurat) <- "integrated"
> pbmc_seurat <- ScaleData(pbmc_seurat, verbose = F)
> pbmc_seurat <- RunPCA(pbmc_seurat, npcs = 30, verbose = F)
> pbmc_seurat <- RunUMAP(pbmc_seurat, reduction = "pca", dims = 1:30, verbose = F)
最后，绘制整合后的UMAP：
> DimPlot(pbmc_seurat, reduction = "umap") +
plot_annotation(title = "10k 3' PBMC and 10k 5' PBMC cells, after integration (Seurat 3)")
数据明显整合得非常好。让我们将图像分离以便于对比：
> DimPlot(pbmc_seurat, reduction = "umap", split.by = "orig.ident") + NoLegend()
对整合后的矩阵进行聚类，并查看聚类在两个数据集之间的分布情况：
> pbmc_seurat <- FindNeighbors(pbmc_seurat, dims = 1:30, k.param = 10, verbose = F)
> pbmc_seurat <- FindClusters(pbmc_seurat, verbose = F)
> DimPlot(pbmc_seurat,label = T) + NoLegend()
计算3’或5’数据集中每个cluster的细胞数：
> count_table <- table(pbmc_seurat@meta.data$seurat_clusters, pbmc_seurat@meta.data$orig.ident)
> count_table
pbmc10k_3p pbmc10k_5p
0 1313 2140
1 1427 945
2 1214 994
3 898 1110
4 612 769
5 576 437
6 401 603
7 338 646
8 311 403
9 359 223
10 251 292
11 347 131
12 202 240
13 136 288
14 251 61
15 109 190
16 118 143
17 76 92
18 93 61
19 55 87
20 35 35
21 24 43
22 13 40
23 15 10
24 14 5
25 2 12
使用自定义函数绘制聚类之间的分布：
> plot_integrated_clusters(pbmc_seurat)
rm(pbmc_seurat)
