第五计 趁火打劫
本指趁人家失火的时候去抢东西。现比喻乘人之危，捞一把。是以“刚”喻己，以“柔”喻敌，言乘敌之危，就势而取胜。
scRNA-seq整合简介
对两个或多个单细胞数据集的联合分析提出了独特的挑战。特别是，在标准工作流程下，识别多个数据集中存在的细胞群体可能会成问题。Seurat v4包括一组用于匹配（或“对齐”）跨数据集的共享细胞群体的方法。这些方法首先确定处于匹配生物学状态（“锚”）的细胞的跨数据集对，既可以用于校正数据集之间的技术差异（即批效应校正），也可以用于对基因组进行比较性scRNA-seq分析跨实验条件。
下面，我们展示了Stuart *，Butler *等人，2019中所述的scRNA-seq整合方法，以对处于静止或干扰素刺激状态的人免疫细胞（PBMC）进行比较分析。
整合目标
以下教程旨在概述使用Seurat集成过程可能进行的复杂细胞类型的比较分析。在这里，我们解决了一些关键目标：
- 创建“集成”数据分析以进行下游分析
- 识别两个数据集中都存在的单元格类型
- 获得在对照和刺激细胞中均保守的细胞类型标记
- 比较数据集以找到对刺激的细胞类型特异性反应
设置Seurat对象
为了方便起见，我们通过SeuratData软件包分发此数据集。
library(Seurat)
library(SeuratData)
library(patchwork)
# install dataset
InstallData("ifnb")
# load dataset
LoadData("ifnb")
# split the dataset into a list of two seurat objects (stim and CTRL)
ifnb.list <- SplitObject(ifnb, split.by = "stim")
# normalize and identify variable features for each dataset independently
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)
执行整合
然后，我们使用FindIntegrationAnchors()函数来识别锚点，该函数将Seurat对象的列表作为输入，并使用这些锚点将两个数据集与集成在一起IntegrateData()`。
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features)
# this command creates an 'integrated' data assay
immune.combined <- IntegrateData(anchorset = immune.anchors)
进行综合分析
现在，我们可以在所有单元上运行单个集成分析！
# specify that we will perform downstream analysis on the corrected data note that the original
# unmodified data still resides in the 'RNA' assay
DefaultAssay(immune.combined) <- "integrated"
# Run the standard workflow for visualization and clustering
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)
# Visualization
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2
为了并排可视化这两个条件，我们可以使用split.by参数来显示每个以聚类着色的条件。
DimPlot(immune.combined, reduction = "umap", split.by = "stim")
识别保守的细胞类型标记
为了鉴定在各种条件下保守的规范细胞类型标记基因，我们提供了该FindConservedMarkers()`功能。此功能对每个数据集/组执行差异基因表达测试，并使用MetaDE R软件包中的荟萃分析方法组合p值。例如，无论簇6中的刺激条件如何，我们都可以计算出保守标记的基因（NK细胞）。
# For performing differential expression after integration, we switch back to the original data
DefaultAssay(immune.combined) <- "RNA"
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 6, grouping.var = "stim", verbose = FALSE)
head(nk.markers)
## CTRL_p_val CTRL_avg_log2FC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj
## GNLY 0 6.006422 0.944 0.045 0
## FGFBP2 0 3.223246 0.503 0.020 0
## CLIC3 0 3.466418 0.599 0.024 0
## PRF1 0 2.654683 0.424 0.017 0
## CTSW 0 2.991829 0.533 0.029 0
## KLRD1 0 2.781453 0.497 0.019 0
## STIM_p_val STIM_avg_log2FC STIM_pct.1 STIM_pct.2 STIM_p_val_adj
## GNLY 0.000000e+00 5.853573 0.956 0.060 0.000000e+00
## FGFBP2 7.275492e-161 2.200379 0.260 0.016 1.022425e-156
## CLIC3 0.000000e+00 3.549919 0.627 0.031 0.000000e+00
## PRF1 0.000000e+00 4.102686 0.862 0.057 0.000000e+00
## CTSW 0.000000e+00 3.139620 0.596 0.035 0.000000e+00
## KLRD1 0.000000e+00 2.880055 0.558 0.027 0.000000e+00
## max_pval minimump_p_val
## GNLY 0.000000e+00 0
## FGFBP2 7.275492e-161 0
## CLIC3 0.000000e+00 0
## PRF1 0.000000e+00 0
## CTSW 0.000000e+00 0
## KLRD1 0.000000e+00 0
我们可以为每个簇探索这些标记基因，并使用它们将我们的簇注释为特定的细胞类型。
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)
DotPlot()带有split.by`参数的函数可用于查看各种条件下的保守细胞类型标记，显示表达水平和表达任何给定基因的簇中细胞的百分比。在这里，我们为14个簇中的每个簇绘制2-3个强标记基因。
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()
确定跨条件的差异表达基因
现在，我们已经排列了刺激细胞和对照细胞，我们可以开始进行比较分析，并观察刺激引起的差异。广泛观察这些变化的一种方法是绘制受刺激细胞和对照细胞的平均表达，并在散点图上寻找视觉异常值的基因。在这里，我们采用受刺激的和对照的原始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
如您所见，许多相同的基因在这两种细胞类型中均被上调，可能代表保守的干扰素应答途径。
因为我们有信心确定出跨条件的常见细胞类型，所以我们可以询问相同条件下不同条件下哪些基因会发生变化。首先，我们在meta.data插槽中创建一列，以保存细胞类型和刺激信息，并将当前标识切换到该列。然后，我们用于FindMarkers()`查找受激B细胞和对照B细胞之间不同的基因。请注意，此处显示的许多顶级基因与我们之前绘制的核心干扰素应答基因相同。此外，我们看到的像CXCL10的基因对单核细胞和B细胞干扰素的反应也具有特异性，在该列表中也显示出很高的意义。
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)
## p_val avg_log2FC pct.1 pct.2 p_val_adj
## ISG15 5.398167e-155 4.5889194 0.998 0.240 7.586044e-151
## IFIT3 2.209577e-150 4.5032297 0.964 0.052 3.105118e-146
## IFI6 7.060888e-150 4.2375542 0.969 0.080 9.922666e-146
## ISG20 7.147214e-146 2.9387415 1.000 0.672 1.004398e-141
## IFIT1 7.650201e-138 4.1295888 0.914 0.032 1.075083e-133
## MX1 1.124186e-120 3.2883709 0.905 0.115 1.579819e-116
## LY6E 2.504364e-118 3.1297866 0.900 0.152 3.519383e-114
## TNFSF10 9.454398e-110 3.7783774 0.791 0.025 1.328627e-105
## IFIT2 1.672384e-105 3.6569980 0.783 0.035 2.350201e-101
## B2M 5.564362e-96 0.6100242 1.000 1.000 7.819599e-92
## PLSCR1 1.128239e-93 2.8205802 0.796 0.117 1.585514e-89
## IRF7 6.602529e-92 2.5832239 0.838 0.190 9.278534e-88
## CXCL10 4.402118e-82 5.2406913 0.639 0.010 6.186297e-78
## UBE2L6 2.995453e-81 2.1487435 0.852 0.300 4.209510e-77
## PSMB9 1.755809e-76 1.6379482 0.940 0.573 2.467438e-72
可视化基因表达中这些变化的另一种有用方法是split.by选择FeaturePlot()或VlnPlot()功能。这将显示给定基因列表的FeaturePlots，并按分组变量（此处为刺激条件）进行划分。诸如CD3D和GNLY之类的基因是典型的细胞类型标记（对于T细胞和NK / CD8 T细胞），实际上不受干扰素刺激的影响，并且在对照组和受刺激组中显示出相似的基因表达模式。另一方面，IFI6和ISG15是核心干扰素反应基因，因此在所有细胞类型中均被上调。最后，CD14和CXCL10是显示细胞类型特异性干扰素应答的基因。CD14单核细胞受刺激后，CD14表达下降，这可能导致在有监督的分析框架中进行错误分类，从而强调了整合分析的价值。
FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3,
cols = c("grey", "red"))
[图片上传失败...(image-275df6-1615650992246)]
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)
