整合受刺激与对照PBMC数据集，以了解细胞类型特异性反应
教程背景
本教程介绍了康等人在2017年提出的两组pbmc。本实验将PBMCs分为受刺激组和对照组，受刺激组给予干扰素治疗。对干扰素的反应引起了细胞类型特异性基因表达的变化，这使得所有数据的联合分析变得困难，因为细胞会因刺激条件和细胞类型而聚集。这里，我们展示了斯图亚特和巴特勒等人(2018年)所述的整合策略，用于执行整合分析，以促进常见细胞类型的识别并进行比较分析。虽然这个例子演示了两个数据集(条件)的整合，但这些方法已经扩展到多个数据集。这个工作流提供了一个整合四个胰岛数据集的示例。
背景知识
interferon beta
干扰素是一类糖蛋白，它具有高度的种属特异性，故动物的干扰素对人无效，干扰素具有抗病毒、抑制细胞增殖、调节免疫及抗肿瘤作用。干扰素 (IFN)是一种广谱抗病毒剂，并不直接杀伤或抑制病毒，而主要是通过细胞表面受体作用使细胞产生抗病毒蛋白，从而抑制乙肝病毒的复制，其类型分为三类，α-(白细胞)型、 β-(成纤维细胞)型，γ-(淋巴细胞)型;同时还可增强自然杀伤细胞(NK细胞)、巨噬细胞和T淋巴细胞的活力，从而起到免疫调节作用，并增强抗病毒能力。干扰素是一组具有多种功能的活性蛋白质(主要是糖蛋白)，是一种由单核细胞和淋巴细胞产生的细胞因子。它们在同种细胞上具有广谱的抗病毒、影响细胞生长，以及分化、调节免疫功能等多种生物活性。
整合的目标
下面的教程旨在向您概述使用Seurat整合过程可以对复杂细胞类型进行的各种比较分析。在此，我们提出三个主要目标:
- 标识出现在两个数据集中的细胞类型
- 获得在对照细胞和受刺激细胞中保存的细胞类型标记
- 比较数据集来发现细胞类型对刺激的特定反应
1.设置Seurat对象
为了方便，我们通过SeuratData包分发这个数据集。
rm(list=ls())
options(stringsAsFactors = F)
# 加载包
library(Seurat)
library(SeuratData)
library(cowplot)
library(patchwork)
## 1.设置Seurat对象
InstallData("ifnb")
data("ifnb")
class(ifnb)
ifnb.list <- SplitObject(ifnb, split.by = "stim")
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
x <- NormalizeData(x)
x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
print(x) # 显示进程以缓解焦虑
})
2.执行整合
然后，我们使用FindIntegrationAnchors函数来标识锚，该函数接受一列Seurat对象作为输入，并使用这些锚将两个数据集与IntegrateData集成在一起。
## 2.执行整合
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, dims = 1:20)
immune.combined <- IntegrateData(anchorset = immune.anchors, dims = 1:20)
3.执行综合分析
现在我们可以在所有细胞上运行一次整合分析!
## 3.进行整合分析
DefaultAssay(immune.combined) <- "integrated"
# 运行可视化和聚类的标准工作流
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
# t-SNE和聚类
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:20)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:20)
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)
p1
p2
plot_grid(p1, p2)
绘制图结果如下：
Snipaste_2020-09-22_19-53-45.png
为了将这两个条件并行地形象化，我们可以使用分割。参数来显示按集群着色的每个条件
DimPlot(immune.combined, reduction = "umap", split.by = "stim")
Snipaste_2020-09-22_19-56-41.png
4.识别保守细胞类型标记
为了识别不同条件下保守的典型细胞类型标记基因，我们提供了findconservedmarker功能。该函数对每个数据集/组执行差异基因表达测试，并使用MetaDE R包中的元分析方法组合p值。例如，我们可以计算出在cluster 6 (NK细胞)中保守标记的基因，而不考虑刺激条件。
## 4.识别保守细胞类型标记
DefaultAssay(immune.combined) <- "RNA"
nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 6, grouping.var = "stim", verbose = T)
head(nk.markers)
CTRL_p_val CTRL_avg_logFC CTRL_pct.1 CTRL_pct.2 CTRL_p_val_adj STIM_p_val STIM_avg_logFC
GNLY 0 4.179377 0.944 0.045 0 0.000000e+00 4.077669
NKG7 0 3.163193 0.954 0.084 0 0.000000e+00 2.901784
GZMB 0 2.926436 0.841 0.043 0 0.000000e+00 3.124419
CLIC3 0 2.400545 0.599 0.024 0 0.000000e+00 2.446643
FGFBP2 0 2.242469 0.503 0.020 0 9.638295e-158 1.476625
CTSW 0 2.077683 0.533 0.029 0 0.000000e+00 2.180716
STIM_pct.1 STIM_pct.2 STIM_p_val_adj max_pval minimump_p_val
GNLY 0.951 0.059 0.00000e+00 0.000000e+00 0
NKG7 0.957 0.080 0.00000e+00 0.000000e+00 0
GZMB 0.895 0.060 0.00000e+00 0.000000e+00 0
CLIC3 0.620 0.031 0.00000e+00 0.000000e+00 0
FGFBP2 0.256 0.016 1.35447e-153 9.638295e-158 0
CTSW 0.596 0.035 0.00000e+00 0.000000e+00 0
我们可以为每个细胞簇探索这些标记基因，并使用它们注释我们的细胞簇作为特定的细胞类型。
FeaturePlot(immune.combined, features = c("CD3D", "SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A", "CCL2", "PPBP"), min.cutoff = "q9")
Snipaste_2020-09-23_00-04-09.png
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")
DimPlot(immune.combined, label = TRUE)
Snipaste_2020-09-23_00-05-40.png
带有拆分的DotPlot函数。参数可用于观察不同条件下的保守细胞类型标记，显示表达任何给定基因的细胞群的表达水平和百分比。在这里，我们为13个聚类分别标记2-3个强标记基因。
Idents(immune.combined) <- factor(Idents(immune.combined), levels = c("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")
DotPlot(immune.combined, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "stim") + RotatedAxis()
Snipaste_2020-09-23_00-07-49.png
5.识别不同条件下的差异表达基因
现在我们已经排列好了受刺激细胞和对照细胞，我们可以开始做比较分析，看看刺激所引起的差异。一种广泛观察这些变化的方法是绘制受刺激细胞和对照细胞的平均表达，并在散点图中寻找视觉异常的基因。在这里，我们取受刺激和对照幼稚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 <- 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 <- 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)
plot_grid(p1, p2)
Snipaste_2020-09-23_00-12-15.png
正如你所看到的，许多相同的基因在这两种细胞类型中上调，可能代表了保守的干扰素反应途径。
因为我们有信心在不同条件下识别出共同的细胞类型，我们可以问在不同条件下相同类型的细胞有什么基因改变。首先，我们在元数据中创建一个列。数据槽，以保存单元类型和刺激信息，并将当前ident切换到该列。然后我们使用findmarker来发现受刺激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_logFC pct.1 pct.2 p_val_adj
ISG15 9.008631e-168 3.2061225 0.998 0.235 1.265983e-163
IFIT3 1.566158e-161 3.1240285 0.961 0.049 2.200921e-157
ISG20 1.206176e-158 2.0549983 1.000 0.662 1.695038e-154
IFI6 1.544804e-158 2.9139826 0.959 0.077 2.170914e-154
IFIT1 1.531805e-144 2.8529052 0.899 0.031 2.152645e-140
MX1 3.490304e-129 2.2712525 0.902 0.113 4.904924e-125
LY6E 3.320706e-127 2.1767293 0.897 0.146 4.666588e-123
TNFSF10 2.108526e-114 2.5996223 0.773 0.021 2.963111e-110
IFIT2 7.373988e-113 2.5114519 0.781 0.035 1.036267e-108
B2M 1.570724e-101 0.4194467 1.000 1.000 2.207338e-97
PLSCR1 8.204551e-101 1.9635626 0.792 0.115 1.152986e-96
IRF7 1.517587e-100 1.8203577 0.837 0.181 2.132665e-96
CXCL10 7.783861e-92 3.6869595 0.660 0.012 1.093866e-87
UBE2L6 1.659777e-88 1.4951634 0.855 0.296 2.332484e-84
EPSTI1 1.324448e-82 1.7540298 0.717 0.103 1.861247e-78
另一种可视化这些基因表达变化的有效方法是分裂。可以选择FeaturePlot或VlnPlot函数。这将显示给定基因列表的特征图，由分组变量(这里是刺激条件)分割。CD3D和GNLY等基因是典型的细胞类型标记(对于T细胞和NK/CD8 T细胞)，它们几乎不受干扰素刺激的影响，在对照组和受刺激组显示相似的基因表达模式。另一方面，IFI6和ISG15是核心干扰素应答基因，在所有细胞类型中都相应上调。最后，CD14和CXCL10是显示细胞类型特异性干扰素反应的基因。CD14单核细胞刺激后CD14表达下降，这可能导致在监督分析框架下的错误分类，强调了综合分析的价值。在干扰素刺激后，CXCL10在单核细胞和B细胞中显示明显的上调，但在其他细胞类型中没有。
FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3,
cols = c("grey", "red"))
Snipaste_2020-09-23_00-15-47.png
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)
Snipaste_2020-09-23_00-16-43.png
网友评论