美文网首页单细胞学习
Seurat之整合分析(1)

Seurat之整合分析(1)

作者: jjjscuedu | 来源:发表于2022-12-15 19:20 被阅读0次

前面也说过用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

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

相关文章

网友评论

    本文标题:Seurat之整合分析(1)

    本文链接:https://www.haomeiwen.com/subject/zchlqdtx.html