Seurat软件学习2-scrna数据整合分析
Seurat软件学习2-scrna数据整合分析
Seurat软件学习1-多个模型得数据进行整合：https://cloud.tencent.com/developer/article/2130078
scRNA-seq整合介绍
两个或多个单细胞数据集的联合分析带来了独特的挑战。特别是，在标准工作流程下，识别存在于多个数据集的细胞群可能是个问题。Seurat v4包括一套方法来匹配（或 "对齐"）跨数据集的共享细胞群。这些方法首先识别处于匹配生物状态的跨数据集的细胞对（"锚"），可用于校正数据集之间的技术差异（即批处理效应校正），并对不同的实验条件进行scRNA-seq比较分析。
下面，我们展示了Stuart, Butler等人在2019年描述的scRNA-seq整合方法，对静止状态或干扰素刺激状态下的人类免疫细胞（PBMC）进行比较分析。
整合的目标
下面的教程旨在让你了解使用Seurat整合程序可以对复杂细胞类型进行的各种比较分析。在这里，我们解决几个关键问题。
1.为下游分析创建一个 "整合 "的数据检测方法
2.识别两个数据集中都有的细胞类型
3.获得在对照组和刺激组细胞中都保守的细胞类型标志物
4.比较数据集以找到细胞类型对刺激的具体反应
设置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)
##SplitObject：基于单个属性将对象拆分为子集合对象列表，每个子集合对象对应于属性的每个级别。例如，用于获取包含来自许多患者的细胞的对象，并将其细分为特定于患者的对象
ifnb.list <- SplitObject(ifnb, split.by = "stim")
# normalize and identify variable features for each dataset independently
##lapply：可以循环处理列表中的每一个元素
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对象的列表作为输入，并使用这些锚点将两个数据集整合在一起。
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.006173 0.944 0.045 0
## FGFBP2 0 3.243588 0.505 0.020 0
## CLIC3 0 3.461957 0.597 0.024 0
## PRF1 0 2.650548 0.422 0.017 0
## CTSW 0 2.987507 0.531 0.029 0
## KLRD1 0 2.777231 0.495 0.019 0
## STIM_p_val STIM_avg_log2FC STIM_pct.1 STIM_pct.2 STIM_p_val_adj
## GNLY 0.000000e+00 5.858634 0.954 0.059 0.000000e+00
## FGFBP2 3.408448e-165 2.191113 0.261 0.015 4.789892e-161
## CLIC3 0.000000e+00 3.536367 0.623 0.030 0.000000e+00
## PRF1 0.000000e+00 4.094579 0.862 0.057 0.000000e+00
## CTSW 0.000000e+00 3.128054 0.592 0.035 0.000000e+00
## KLRD1 0.000000e+00 2.863797 0.552 0.027 0.000000e+00
## max_pval minimump_p_val
## GNLY 0.000000e+00 0
## FGFBP2 3.408448e-165 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)
带有split.by参数的DotPlot()函数对于查看不同条件下的保守细胞类型标记是很有用的，它既能显示表达水平，又能显示一个簇中表达任何特定基因的细胞的百分比。这里我们为我们的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细胞之间不同的基因。注意到这里出现的许多顶级基因与我们之前绘制的核心干扰素反应基因是一样的。此外，像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 1.212995e-155 4.5997247 0.998 0.239 1.704622e-151
## IFIT3 4.743486e-151 4.5017769 0.964 0.052 6.666020e-147
## IFI6 1.680324e-150 4.2361116 0.969 0.080 2.361359e-146
## ISG20 1.595574e-146 2.9452675 1.000 0.671 2.242260e-142
## IFIT1 3.499460e-137 4.1278656 0.910 0.032 4.917791e-133
## MX1 8.571983e-121 3.2876616 0.904 0.115 1.204621e-116
## LY6E 1.359842e-117 3.1251242 0.895 0.152 1.910986e-113
## TNFSF10 4.454596e-110 3.7816677 0.790 0.025 6.260044e-106
## IFIT2 1.290640e-106 3.6584511 0.787 0.035 1.813736e-102
## B2M 2.019314e-95 0.6073495 1.000 1.000 2.837741e-91
## PLSCR1 1.464429e-93 2.8195675 0.794 0.117 2.057961e-89
## IRF7 3.893097e-92 2.5867694 0.837 0.190 5.470969e-88
## CXCL10 1.624151e-82 5.2608266 0.640 0.010 2.282419e-78
## UBE2L6 2.482113e-81 2.1450306 0.852 0.299 3.488114e-77
## PSMB9 5.977328e-77 1.6457686 0.940 0.571 8.399938e-73另一个可视化这些基因表达变化的有用方法是使用FeaturePlot()或VlnPlot()函数的split.by选项。这将显示给定基因列表的FeaturePlot，通过分组变量（这里是刺激条件）进行分割。CD3D和GNLY等基因是典型的细胞类型标记（用于T细胞和NK/CD8 T细胞），几乎不受干扰素刺激的影响，在对照组和刺激组中显示类似的基因表达模式。另一方面，IFI6和ISG15是干扰素反应的核心基因，在所有细胞类型中都相应地被上调。最后，CD14和CXCL10是显示细胞类型特定的干扰素反应的基因。CD14单核细胞受刺激后表达量下降，这可能导致监督分析框架下的错误分类，强调了综合分析的价值。CXCL10在干扰素刺激后在单核细胞和B细胞中显示出明显的上调，但在其他细胞类型中没有。
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规范化的数据集进行整合
在Hafemeister和Satija, 2019年，我们介绍了一种改进的scRNA-seq归一化方法，基于正则化的负二项回归。该方法被命名为 "sctransform"，并避免了标准归一化工作流程的一些缺陷，包括添加伪计数和对数转换。你可以在手稿或我们的SCTransform说明中阅读更多关于sctransform的信息。
下面，我们将演示如何修改Seurat整合工作流程，以处理已经用sctransform工作流程规范化了的数据集。这些命令大体上是相似的，但有几个关键的区别。
1.在整合前通过SCTransform()单独规范数据集，而不是NormalizeData()
2.正如我们在SCTransform中进一步讨论的那样，我们通常使用3000个或更多的特征来分析SCTransform的下游。
3.在识别锚点之前运行PrepSCTIntegration()函数
4.当运行FindIntegrationAnchors()和IntegrateData()时，将normalization.method参数设定为SCT值。
5.当运行基于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
现在，数据集已被整合，你可以按照中前述步骤确定细胞类型和细胞类型的特定反应。
总结
针对一个样本的不同处理是十分重要的，有可能处理之后值后差别很大，然后导致数据出现一定的问题，因此引入别人的数据集进行整合分析，除了可以看自己的处理的效果怎么样，还可以更精的一步找到不一样的东西。作为单细胞的初学者，确实是很有必要把官网作者推荐的东西都看完。
本文系外文翻译，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文系外文翻译，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。