数据整合概述
对两个或多个单细胞数据集的联合分析带来了独特的挑战。特别是，在标准工作流程下，识别存在于多个数据集中的细胞群可能会有问题。Seurat有包括一组跨数据集match（或“align”）共有细胞群的方法。这些方法首先识别处于匹配生物学状态（“anchors”）的跨数据集细胞对，既可用于校正数据集之间的技术差异（即批量效应校正），也可用于跨实验条件的比较scRNA-seq分析。
整合目的
以下教程旨在向您概述使用 Seurat 积分程序可以对复杂细胞类型进行的比较分析。在这里，我们解决了几个关键目标：
- 为下游分析创建“整合”的数据分析
- 识别两个数据集中存在的细胞类型
- 获取在对照细胞和刺激细胞中都保守的细胞类型标记物
- 比较数据集以找到细胞类型对刺激的特定反应
构建seurat 对象
seurat官网给出的数据是 SeuratData 内置数据，但是使用 SeuratData 下载会提示相应的问题，也有通过数据连接直接下载。这里我使用个人数据，特点是：已经通过 cellranger 处理了 FASTQ 文件，得到了 feature， barcode， matrix 文件， 如下
上一级文件如下，
层级显示关系如下
然后需要将这些文件统一读入到一个列表中，操作如下
# 创建一个空列表来存储数据块
seu.list <- list()
# 指定文件夹路径
folder_path <- "./HRA000212_raw_feature_bc_matrix"
# 获取文件夹中的文件列表
file_list <- list.files(path = folder_path, full.names = TRUE)
# 遍历文件列表并读入数据
for (file_path in file_list) {
# 提取文件名（不包括路径和扩展名）
file_name <- tools::file_path_sans_ext(basename(file_path))
file_name <- substr(basename(file_name), 1, 9)
# 读入文件数据（这里假设文件是文本文件）
data <- Read10X(file_path)
scobj <- CreateSeuratObject(data, project = file_name, min.cells = 3, min.features = 200)
# 将数据添加到列表中，以文件名作为标识
seu.list[[file_name]] <- scobj
# return(seu.list)
}
# 打印列表中的数据
print(seu.list)
这样就得到了一个list，seu.list，当中把seutat处理过的几个样本都存放到了一起。
这里是引用
这一步就是和seurat官网中的得到列表的操作
# split the dataset into a list of two seurat objects (stim and CTRL)
ifnb.list <- SplitObject(ifnb, split.by = "stim")
之所以得到列表，是因为后续的 FindIntegrationAnchors() 函数使用的是列表读入
读入之后，可以对数据进行标准化、高变基因操作。类似于单个数据处理。这一步可以在这里处理，也可以整合之后再处理，都可。
seu.list <- lapply(X = seu.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})
识别在多个数据集之间具有一致性变化的特征
# select features that are repeatedly variable across datasets for integration
features <- SelectIntegrationFeatures(object.list = seu.list)
SelectIntegrationFeatures 是 Seurat 包中的一个函数，用于选择在单细胞RNA测序数据整合分析中用于集成的特征（基因）。在整合多个单细胞RNA测序数据集时，通常不是所有的基因都对整合有用。因此，这个函数的主要目的是识别在多个数据集之间具有一致性变化的特征，以便将它们用于整合。
具体来说，SelectIntegrationFeatures 函数执行以下任务：
对于每个单细胞RNA测序数据集，它计算每个基因的方差。方差衡量了基因在该数据集中的变化程度。通常，对于一致性变化的基因，它们在不同数据集之间的方差较小。
然后，函数根据计算的方差值对基因进行排序，选择具有最高方差的一组基因。这些高方差的基因被认为是用于整合的候选特征。
最后，用户可以根据需要设置一个特定的阈值，从排序后的候选特征中选择前几个用于整合。
通过选择具有一致性变化的高方差基因，SelectIntegrationFeatures
函数有助于确保在整合分析中使用的特征具有生物学上的相关性，从而提高了整合的准确性和有效性。在调用这个函数时，你可以通过设置不同的参数来控制如何选择整合特征，例如选择特定数量的特征或根据阈值选择特征。这些参数可以根据你的数据集和研究问题进行调整。
整合数据
然后使用 FindIntegrationAnchors() 函数识别anchors，这个函数将 Seurat 对象列表作为输入，并使用这些anchors通过 IntegrateData() 将两个数据集集成在一起
anchors <- FindIntegrationAnchors(object.list = seu.list, anchor.features = features)
FindIntegrationAnchors 是 Seurat 包中的一个函数，用于在单细胞RNA测序数据的集成分析中寻找集成锚点（integration anchors）。这个函数的主要目的是将多个单细胞RNA测序数据集整合在一起，以便在一个统一的数据空间中进行分析，例如聚类、降维、可视化等。
在单细胞RNA测序研究中，可能存在多个不同的实验批次或样本来源，每个批次或样本都有其独特的特点和噪声。为了将这些数据集成在一起，可以使用集成分析的方法。FindIntegrationAnchors
函数的作用是寻找可以用于整合不同数据集的共同特征，这些共同特征被称为 “锚点”。具体来说，FindIntegrationAnchors 函数会对输入的单细胞RNA测序数据集执行以下步骤：
首先，对每个数据集进行降维操作，通常使用主成分分析（PCA）或其他降维技术，以减少数据的维度并捕捉关键的变化。
然后，函数会计算不同数据集之间的共同变化，以识别集成锚点。这些共同变化通常通过比较降维后的数据空间中的数据点来确定。
最后，找到的锚点可以用于整合不同数据集，将它们映射到一个共享的数据空间中，从而允许在整合后的数据上执行一致的分析。
整合分析有助于减轻不同数据集之间的批次效应（batch effects）并提高数据的可比性。一旦找到锚点，你可以使用
IntegrateData 函数来执行实际的数据整合操作。这些功能在单细胞RNA测序研究中特别有用，因为它们允许将来自不同来源的单细胞数据整合在一起，以便更全面地理解细胞类型和状态。
这一步会特别耗时，所以又学习了一下 seurat 联合 feature，多线程处理的方法，进行如下修改
# Perform integration -------------------------------------------------------------------------
library(future)
# check the current active plan
plan()
# change the current plan to access parallelization
plan("multiprocess", workers = 15)
anchors <- FindIntegrationAnchors(object.list = seu.list, anchor.features = features)
这样调动了多线程处理，但是如果数据比较大的话，（大概率比较大）有可能出现这样的报错
anchors <- FindIntegrationAnchors(object.list = seu.list, anchor.features = features)
Scaling features for provided objects
Finding all pairwise anchors
Error in getGlobalsAndPackages(expr, envir = envir, globals = globals) :
The total size of the 25 globals exported for future expression (‘FUN()’) is 8.60 GiB.. This exceeds the maximum allowed size of 500.00 MiB (option 'future.globals.maxSize'). The three largest globals are ‘object.list’ (8.60 GiB of class ‘list’), ‘FUN’ (323.28 KiB of class ‘function’) and ‘ReciprocalProject’ (258.20 KiB of class ‘function’)
对于某些函数，每个工作线程都需要访问某些全局变量。如果这些限制大于默认限制，您将看到此错误。要解决此问题，可以设置 options（future.globals.maxSize = X），其中 X 是允许的最大大小（以字节为单位）。因此，要将其设置为 1GB，您需要运行 options（future.globals.maxSize = 1000 * 1024^2）。这个操作将增加RAM 使用率。
根据个人的内存情况去设置就可以
library(future)
# check the current active plan
plan()
# change the current plan to access parallelization
plan("multiprocess", workers = 15)
options(future.globals.maxSize = 10000 * 1024^2)
anchors <- FindIntegrationAnchors(object.list = seu.list, anchor.features = features)
然后就可以跑起来了
> anchors
An AnchorSet object containing 579144 anchors between 11 Seurat objects
This can be used as input to IntegrateData.
还是要忍不住吐槽一下seurat这个整合过程，特别耗费内存，特别慢
然后整合数据
# this command creates an 'integrated' data assay
combined <- IntegrateData(anchorset = anchors)
得到的数据
> combined
An object of class Seurat
32922 features across 152190 samples within 2 assays
Active assay: integrated (2000 features, 2000 variable features)
1 other assay present: RNA
到这里，就算是把数据去批次，合并完成，可以看一下都有哪些变化
> slotNames(combined)
[1] "assays" "meta.data" "active.assay" "active.ident" "graphs" "neighbors" "reductions" "images" "project.name"
[10] "misc" "version" "commands" "tools"
> slotNames(anchors)
[1] "object.list" "reference.cells" "reference.objects" "query.cells" "anchors" "offsets" "anchor.features"
[8] "neighbors" "command"
> slotNames(scobj)
[1] "assays" "meta.data" "active.assay" "active.ident" "graphs" "neighbors" "reductions" "images" "project.name"
[10] "misc" "version" "commands" "tools"
可以看出，合并后的数据combined 和 scobj 这两个是一致的，也就是都整合成了seurat对象，中间部分，寻找anchors的步骤，有 “reference.cells” “reference.objects” “query.cells” “anchors” “offsets” “anchor.features” " neighbors" ，可以看出，大概就是以某一个位参考，然后依次对比添加的过程。
上面这几个步骤都十分耗时，建议使用system.time()函数评估一下时间
既然这么难得到数据，在这个步骤还是应该保存一下，方便下次读取
qs::qsave(combined, file = 'outputs_data/seurat_combined.qs')
保存方式这里，对比了一下qssave函数和saveRDS这两种方法
saveRDS(combined, file = 'outputs_data/processing_data/seurat_combined.rds')
combined是10G大小，压缩后，qssave大小为3.4G， saveRDS压缩后的大小为3.2G，压缩效率上saveRDS更好一些，没有使用system.time比较时间，但是直观感受还是qssave要快很多。
整合后的分析
这里的分析和单个样本的是类似的
# specify that we will perform downstream analysis on the corrected data note that the
# original unmodified data still resides in the 'RNA' assay
DefaultAssay(combined) <- "integrated"
在 Seurat 包中，DefaultAssay() 函数用于设置默认的数据存储"assay"，以便在后续分析中使用。DefaultAssay() 函数允许你指定在 Seurat 对象中的哪个 “assay”
存储了你希望在分析中使用的主要数据。通常，这对于存储和管理不同处理步骤或数据类型的数据非常有用。在你提供的示例代码 DefaultAssay(combined) <- “integrated” 中，它的作用是将 combined中的默认 “assay” 设置为 “integrated”。这意味着在后续的分析中，Seurat 将默认使用名为 "integrated"的数据存储进行操作，除非另有指定。这通常用于整合多个数据集或不同处理步骤的数据，以便在一个统一的数据空间中进行进一步的分析和可视化。
设置默认 “assay” 有助于简化代码，因为你可以在后续的分析中省略指定要使用哪个 “assay”。如果你没有明确指定默认"assay"，则在某些情况下需要在函数中明确指定要使用的 “assay”，这可能会使代码更复杂。
# Run the standard workflow for visualization and clustering
combined <- ScaleData(combined, verbose = FALSE)
combined <- RunPCA(combined, npcs = 30, verbose = FALSE)
combined <- RunUMAP(combined, reduction = "pca", dims = 1:30)
combined <- RunTSNE(combined, reduction = "pca", dims = 1:30)
combined <- FindNeighbors(combined, reduction = "pca", dims = 1:30)
combined <- FindClusters(combined, resolution = 0.5)
展示去批次
# Visualization
p1 <- DimPlot(combined, reduction = "umap", group.by = "orig.ident")
p2 <- DimPlot(combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2
换个颜色也还是可以的
p1 <- DimPlot(combined, reduction = "umap", group.by = "orig.ident")+
ggsci::scale_color_d3("category20")
p2 <- DimPlot(combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2
可以看出合并去批次之后的效果还是有的，样本都聚合都一起。
seurat package可以进行数据合并和整合，进行批次矫正，但是缺点是速度有些慢。为了知识的完整性，还是要把这个整合方法进行学习和记录。
参考文章
Introduction to scRNA-seq integration
Batch Effect Correction
Single-cell integration benchmarking
扫一扫
693
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
