Seurat之整合分析（1）

前面也说过用CCA，SCTransform和Harmony等可以消除多个数据的批次效应，来实现多组数据的整合。但有时候我们会遇到两个datasets只有部分重叠，这和之前介绍的方法就有一点不同了。特别是，在标准工作流下，识别存在于多个数据集中的基因可能存在问题。Seurat v4 包括一组方法，以匹配（或"对齐"）跨数据集共同的基因。这些方法首先识别处于匹配生物状态的交叉数据集对（"锚点"），既可用于纠正数据集之间的技术差异（即批量效应校正），又可用于对整个实验条件进行比较scRNA-seq分析，当然也可以实现跨物种的整合，例如我们利用跨物种的同源基因。

原理

Seurat整合数据集的过程首先是对不同的数据集进行降维，在低维空间中寻找成对的anchors，然后再以anchor作为锚点，对其他的点的坐标进行调整，从而将整个数据集整合在一起。

https://www.jianshu.com/p/ebc328f9fb73

Integration

官方教程中，在进行数据整合之前，先将所有的数据整合在一起，再split成一个object list，在object list中又进行了标准化等操作。真正进行整合的核心代码其实是FindIntegationAnchors以及IntegrateData两个函数来实现的。

下面，我们以在ctrl或干扰素刺激状态下对人体免疫细胞 （PBMC） 进行比较分析。

library(Seurat)

library(SeuratData)

library(patchwork)

install.packages("ifnb.SeuratData")  #我是直接下载下来自己安装的测试数据，也可以通过命令安装：InstallData("ifnb")

library(ifnb.SeuratData)

LoadData("ifnb")

> ifnb

An object of class Seurat

14053 features across 13999 samples within 1 assay

Active assay: RNA (14053 features, 0 variable features)

#把数据集分成2个 (stim和CTRL)

ifnb.list <- SplitObject(ifnb, split.by = "stim")

#对每个数据集进行标准化，提取细胞间变异系数较大的top 2000个基因

ifnb.list <- lapply(X = ifnb.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 = ifnb.list)

执行整合

然后，我们使用该功能识别锚点，以 Seurat 对象列表作为输入，并使用这些锚点将两个数据集集成在一起。其实，这个工作真正进行整合的核心代码是FindIntegationAnchors以及IntegrateData两个函数来实现的。

immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features)

immune.combined <- IntegrateData(anchorset = immune.anchors)

FindIntegrationAnchors：

这个函数主要完成以下几个步骤：

1. 降维

对list中的数据集分别进行降维，以便下一步在低维空间中寻找细胞的nearest neighbors。这一步可以选择不同的方法，包括cca(canonical correlation analysis) rpca(reciprocal pca) rlsi三种，可以用reduction进行设置，默认方法为cca。

cca是针对两个不同数据集，寻找两个数据集之间相关性最大的成分的分析方法。

2. 寻找anchors

在低维空间中，对于一个数据集中的每个细胞，寻找它在另一个数据集中距离最近的k个细胞（k nearest neighbors）。如果X数据集中的细胞a在Y数据集中找到了b1 b2 … bk这几个kNN，而Y数据集中的b1细胞在X数据集中找到的kNN又恰巧包括a，那么a和b1就是一对nearest neighbor，又称为mutual nearest neighbor(MNN)。在这一步中，可以设置的参数为k.anchor，就是在相对数据集中寻找到的neighbor的个数，默认值为5.

3. 过滤anchors

这一步将在未降维的原始数据集中，提取与CCV(canonical correlation vectors)相关性最高的n个features。对其进行normalization，然后在这些feaure构成的空间中，对上一步找到的nearest neighbor pairs进行筛选。比如，若上一步找到的a和b1这个pair，在这个高位空间中，a的前k个nearest neighbor中不包含b1，那么这个nearest neighbor pair将会被过滤掉，不会再被使用。在这一步中，可以自行设置的参数包括max.features，即选取多少个最相关基因构成高维空间，默认为200；k.filter，即在相对数据集中寻找nearest neighbor的个数，默认为200。需要注意的是，这一步寻找NN的范围设置较宽，因为不是为了寻找新的pair，而是为了过滤掉在上一步中找到的不可靠的Pair。

4. 对anchors进行打分

如果X数据集的细胞a与Y数据集中的细胞b相似，那么细胞a应该也与Y中与b相似的那些细胞相似，细胞b应该也与X中与a相似的那些细胞相似。基于这种原理对anchor进行打分。对于一个X Y数据及中的一对anchor pair a和b，在X数据及中找到与a最相似的k个细胞，同时在Y数据集中也找到与其最相似的k个细胞，同时也对b做相同的事情。然后计算他们又多少neighbor是共享的，即shared nearest neighbors(SNN)。根据SNN的比例对这个anchor pair进行打分。SNN越多，打分越高，意味着这个paire越可靠，在后续的整合中这个pair就会有更高的权重。涉及这一步的参数是k.weight，即在不同数据集中寻找的neighbor的个数，默认为30。

经过这几步，FindIntegrationAnchors就会找到两个数据集之间的较为可靠的anchors，并且对可靠性有一个评估。

IntegrateData：

这个函数主要完成对于每个细胞校准。

对于每个细胞，找到与其距离最近的k个anchors，然后根据细胞与anchor的距离，以及上一步对anchor可靠性的打分，决定每个anchors对这个细胞的权重。k值由k.weight设置，默认值为100。然后利用anchors的坐标以及不同anchors所占的权重，对每个细胞的表达数据进行校准。

IntegrateData返回的是与SCTransform或NormalizeData相同格式的数据slot。如果normalization.method设置为LogNormalize，则返回数据为normalized slot，储存在@data中，如果设置为SCT，则返回数据为完成中心化的数据，储存在@scale.data中。

执行整合分析：

下面，我们可以对所有细胞进行单次整合分析，像以前的seurat一样。

# unmodified data still resides in the 'RNA' assay

DefaultAssay(immune.combined) <- "integrated"

immune.combined <- ScaleData(immune.combined, verbose = FALSE)

immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)

immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)

immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)

immune.combined <- FindClusters(immune.combined, resolution = 0.5)

p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")

p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE, repel = TRUE)

p1 + p2

DimPlot(immune.combined, reduction = "umap", split.by = "stim")

鉴定在不同的组或数据集里面保守的细胞标记基因：

在以前的seurat分析中，我们是利用FindMarkers等来鉴定不同cluster之间的细胞标记基因，现在我们可以利用FindConservedMarkers鉴定特定cluster在不同组之间都存在的，也就是保守的marker基因。例如：我们可以计算出无论刺激条件如何，第6组（NK细胞）中保守标记的基因。

DefaultAssay(immune.combined) <- "RNA"

nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 6, grouping.var = "stim", verbose = FALSE)

head(nk.markers)

从结果可以看出，前半部分和普通的seurat的分析类似，只不过是分了2组而已，在2个组中都是marker基因。

画基因的表达图，也和前面讲seurat类似：

FeaturePlot(immune.combined, features = c("CD3D", "SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A", "CCL2", "PPBP"), min.cutoff = "q9")

immune.combined <- RenameIdents(immune.combined, '0' = "CD14 Mono", '1' = "CD4 Naive T", '2' = "CD4 Memory T",'3' = "CD16 Mono", '4' = "B", '5' = "CD8 T", '6' = "NK", '7' = "T activated", '8' = "DC", '9' = "B Activated",'10' = "Mk", '11' = "pDC", '12' = "Eryth", '13' = "Mono/Mk Doublets", '14' = "HSPC")

DimPlot(immune.combined, label = TRUE)

Idents(immune.combined) <- factor(Idents(immune.combined), levels = c("HSPC", "Mono/Mk Doublets", "pDC", "Eryth", "Mk", "DC", "CD14 Mono", "CD16 Mono", "B Activated", "B", "CD8 T", "NK", "T activated", "CD4 Naive T", "CD4 Memory T"))

markers.to.plot <- c("CD3D", "CREM", "HSPH1", "SELL", "GIMAP5", "CACYBP", "GNLY", "NKG7", "CCL5",

    "CD8A", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A", "VMO1", "CCL2", "S100A9", "HLA-DQA1",

    "GPR183", "PPBP", "GNG11", "HBA2", "HBB", "TSPAN13", "IL3RA", "IGJ", "PRSS57")

DotPlot(immune.combined, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "stim") +

    RotatedAxis()

两组数据对比，识别差异基因：

因为我们有2组数据：STIM和CTRL，也就是刺激和对照细胞，我们可以开始做比较分析，看看刺激引起的差异。观察这些变化的一种方法是绘制受刺激细胞和对照细胞的平均表达，并寻找在散点图上离群的基因。在这里，我们采取刺激和对照组的幼稚T细胞和CD14+单核细胞群的平均表达，并产生散点图，突出显示干扰素刺激剧烈反应的基因。

library(ggplot2)

library(cowplot)

theme_set(theme_cowplot())

t.cells <- subset(immune.combined, idents = "CD4 Naive T")

Idents(t.cells) <- "stim"

avg.t.cells <- as.data.frame(log1p(AverageExpression(t.cells, verbose = FALSE)$RNA))

avg.t.cells$gene <- rownames(avg.t.cells)

cd14.mono <- subset(immune.combined, idents = "CD14 Mono")

Idents(cd14.mono) <- "stim"

avg.cd14.mono <- as.data.frame(log1p(AverageExpression(cd14.mono, verbose = FALSE)$RNA))

avg.cd14.mono$gene <- rownames(avg.cd14.mono)

genes.to.label = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1", "IFIT2", "IFIT1", "CXCL10", "CCL8")

p1 <- ggplot(avg.t.cells, aes(CTRL, STIM)) + geom_point() + ggtitle("CD4 Naive T Cells")

p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)

p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point() + ggtitle("CD14 Monocytes")

p2 <- LabelPoints(plot = p2, points = genes.to.label, repel = TRUE)

p1 + p2

这样，我们可以清晰的发现那些基因在细胞水平有明显的差异。

我们当然，也可以挑选相同细胞类型在不同组中显著差异的细胞。还是利用Seurat中常用的FindMarkers函数。只不过要把细胞类型重新命名，加上组别的信息。

immune.combined$celltype.stim <- paste(Idents(immune.combined), immune.combined$stim, sep = "_")

immune.combined$celltype <- Idents(immune.combined)

Idents(immune.combined) <- "celltype.stim"

b.interferon.response <- FindMarkers(immune.combined, ident.1 = "B_STIM", ident.2 = "B_CTRL", verbose = FALSE)

head(b.interferon.response, n = 15)

我们也可以对于给定基因列表的特征图，按分组变量（此处为刺激对照）进行拆分。

FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3,cols = c("grey", "red"))

plots <- VlnPlot(immune.combined, features = c("LYZ", "ISG15", "CXCL10"), split.by = "stim", group.by = "celltype",pt.size = 0, combine = FALSE)

wrap_plots(plots = plots, ncol = 1)

当然，网上还建议了另外一个改进版的整合功能，是使用sctransform来整合数据

有几个关键差异：

1. 通过SCTransform()单独（而不是在整合之前）实现数据集标准化；

2. 在识别锚点之前运行：[PrepSCTIntegration()]

3. 运行[FindIntegrationAnchors()]和[IntegrateData()]时，设置参数到normalization.method为SCT

4. 当运行基于 sctransform 的工作流（包括整合）时，不运行该功能[ScaleData()]

大致流程如下：（但是具体那种好，我现在也不是很清楚）

LoadData("ifnb")

ifnb.list <- SplitObject(ifnb, split.by = "stim")

ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)

features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)

ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)

immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",

    anchor.features = features)

immune.combined.sct <- IntegrateData(anchorset = immune.anchors, normalization.method = "SCT")

immune.combined.sct <- RunPCA(immune.combined.sct, verbose = FALSE)

immune.combined.sct <- RunUMAP(immune.combined.sct, reduction = "pca", dims = 1:30)

p1 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "stim")

p2 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "seurat_annotations", label = TRUE,repel = TRUE)

p1 + p2

我们后面会用自己的数据集试一下看看。