Seurat软件学习5-scRNA-Seq和scATAC-Seq数据整合
原创Seurat软件学习5-scRNA-Seq和scATAC-Seq数据整合
原创
Seurat软件学习1-多个模型得数据进行整合：https://cloud.tencent.com/developer/article/2130078
Seurat软件学习2-scrna数据整合分析：https://cloud.tencent.com/developer/article/2131431
Seurat软件学习3-scrna数据整合分析注释数据集：https://cloud.tencent.com/developer/article/2133583
Seurat软件学习4-使用RPCA进行快速整合数据集：https://cloud.tencent.com/developer/article/2134684
单细胞转录学改变了我们描述细胞状态的能力，但深入的生物学理解需要的不仅仅是簇的分类列表。随着测量不同细胞形态的新方法的出现，一个关键的分析挑战是整合这些数据集，以更好地了解细胞的身份和功能。例如，用户可以在相同的生物系统上执行scRNA-seq和scatac-seq实验，并用相同的细胞类型标签集一致地注释这两个数据集。这一分析特别具有挑战性，因为scatac-seq数据集很难注释，因为以单细胞分辨率收集的基因组数据稀少，而且scRNA-seq数据中缺乏可解释的基因标记。
在Stuart，Butler等人，2019年，我们介绍了集成从同一生物系统收集的scRNA-seq和scATAC-seq数据集的方法，并在本章中演示了这些方法。特别是，我们演示了以下分析：
1.如何使用带注释scRNA-seq数据集来标记来自scatac-seq实验的细胞？
2.如何从scRNA-seq和scatac-seq共同显示(共同嵌入)细胞？
3.如何将scatac-seq细胞投影到从scRNA-seq实验派生的UMAP上？
这个Vignette广泛使用了Signac软件包，该软件包是最近开发的，用于分析以单细胞分辨率收集的染色质数据集，包括scATAC-seq。有关分析scatac-seq数据的其他文档，请参阅Signac网站。
我们使用从10x基因组公司公开获得的约12,000个人外周血单核细胞‘多组’数据集来演示这些方法。在该数据集中，在同一细胞中同时收集了scRNA-seq和scATAC-seq图谱。出于本分析的目的，我们将数据集视为来自两个不同的实验，并将它们集成在一起。由于它们最初是在相同的单元格中测量的，这提供了一个基本事实，我们可以用它来评估积分的准确性。我们强调，我们在这里使用多组数据集是为了演示和评估目的，用户应该将这些方法应用于分别收集的scRNA-seq和scATAC-seq数据集。我们提供了一个单独的加权最近邻域(WNN)来描述多组单细胞数据的分析策略。
数据加载并单独处理每一种模式
PBMC多组数据集可从10x基因组学获得。为了便于加载和探索，它还作为我们的SeuratData包的一部分提供。我们分别加载RNA和ATAC数据，并假装这些配置文件是在单独的实验中测量的。我们在我们的WNN小节中注释了这些单元格，这些注释也包含在SeuratData中。
library(SeuratData)
# install the dataset and load requirements
InstallData("pbmcMultiome")
library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)
library(ggplot2)
library(cowplot)
# load both modalities
pbmc.rna <- LoadData("pbmcMultiome", "pbmc.rna")
pbmc.atac <- LoadData("pbmcMultiome", "pbmc.atac")
# repeat QC steps performed in the WNN vignette
pbmc.rna <- subset(pbmc.rna, seurat_annotations != "filtered")
pbmc.atac <- subset(pbmc.atac, seurat_annotations != "filtered")
# Perform standard analysis of each modality independently RNA analysis
pbmc.rna <- NormalizeData(pbmc.rna)
pbmc.rna <- FindVariableFeatures(pbmc.rna)
pbmc.rna <- ScaleData(pbmc.rna)
pbmc.rna <- RunPCA(pbmc.rna)
pbmc.rna <- RunUMAP(pbmc.rna, dims = 1:30)
# ATAC analysis add gene annotation information
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- "UCSC"
genome(annotations) <- "hg38"
Annotation(pbmc.atac) <- annotations
# We exclude the first dimension as this is typically correlated with sequencing depth
pbmc.atac <- RunTFIDF(pbmc.atac)
pbmc.atac <- FindTopFeatures(pbmc.atac, min.cutoff = "q0")
pbmc.atac <- RunSVD(pbmc.atac)
pbmc.atac <- RunUMAP(pbmc.atac, reduction = "lsi", dims = 2:30, reduction.name = "umap.atac", reduction.key = "atacUMAP_")现在我们把两种模式的结果都画出来。细胞之前已经根据转录组状态进行了注释。我们将预测scATAC-seq细胞的注释。
p1 <- DimPlot(pbmc.rna, group.by = "seurat_annotations", label = TRUE) + NoLegend() + ggtitle("RNA")
p2 <- DimPlot(pbmc.atac, group.by = "orig.ident", label = FALSE) + NoLegend() + ggtitle("ATAC")
p1 + p2
识别scRNA-seq和scATAC-seq数据集之间的锚点
为了确定scRNA-seq和scATAC-seq实验之间的 "锚"，我们首先使用Signac软件包中的GeneActivity()函数对2kb上游区域和基因体的ATAC-seq计数进行量化，从而产生对每个基因转录活性的粗略估计。随后，来自scATAC-seq数据的基因活性得分与来自scRNA-seq的基因表达量化一起被用作典型相关分析的输入。我们对所有从scRNA-seq数据集中确定为高变量的基因进行这种量化。
# quantify gene activity
gene.activities <- GeneActivity(pbmc.atac, features = VariableFeatures(pbmc.rna))
# add gene activities as a new assay
pbmc.atac[["ACTIVITY"]] <- CreateAssayObject(counts = gene.activities)
# normalize gene activities
DefaultAssay(pbmc.atac) <- "ACTIVITY"
pbmc.atac <- NormalizeData(pbmc.atac)
pbmc.atac <- ScaleData(pbmc.atac, features = rownames(pbmc.atac))通过标签转移对scATAC-seq细胞进行注释
在确定了锚之后，我们可以将注释从scRNA-seq数据集转移到scATAC-seq细胞上。这些注释被储存在seurat_annotations字段中，并被作为输入提供给refdata参数。输出将包含一个矩阵，其中有每个ATAC-seq细胞的预测和置信度分数。
celltype.predictions <- TransferData(anchorset = transfer.anchors, refdata = pbmc.rna$seurat_annotations,
weight.reduction = pbmc.atac[["lsi"]], dims = 2:30)
pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)Why do you choose different (non-default) values for reduction and weight.reduction?
在FindTransferAnchors()中，当在scRNA-seq数据集之间转移时，我们通常将参考文献中的PCA结构投影到查询中。然而，当跨模式转移时，我们发现CCA能更好地捕捉到共享的特征相关结构，因此在这里设置减少='cca'。此外，在TransferData()中，我们默认使用相同的投影PCA结构来计算影响每个细胞预测的锚点的局部邻域的权重。在scRNA-seq向scATAC-seq转移的情况下，我们使用通过计算ATAC-seq数据的LSI所学到的低维空间来计算这些权重，因为这更好地捕捉了ATAC-seq数据的内部结构。
执行转移后，ATAC-seq细胞的预测注释（从scRNA-seq数据集转移而来）存储在predicted.id字段中。由于这些细胞是用multiome试剂盒测量的，我们也有一个可用于评估的基础真实注释。你可以看到预测的和实际的注释是极其相似的。
pbmc.atac$annotation_correct <- pbmc.atac$predicted.id == pbmc.atac$seurat_annotations
p1 <- DimPlot(pbmc.atac, group.by = "predicted.id", label = TRUE) + NoLegend() + ggtitle("Predicted annotation")
p2 <- DimPlot(pbmc.atac, group.by = "seurat_annotations", label = TRUE) + NoLegend() + ggtitle("Ground-truth annotation")
p1 | p2
在这个例子中，通过scATAC-seq整合，scATAC-seq图谱的注释有90%的时间是正确预测的。此外，prediction.score.max字段量化了与我们预测的注释相关的不确定性。我们可以看到，被正确注释的细胞通常与高预测分数（>90%）有关，而被错误注释的细胞则与急剧降低的预测分数（<50%）有关。不正确的分配也倾向于反映密切相关的细胞类型（Intermediate vs. Naive B cells）。
predictions <- table(pbmc.atac$seurat_annotations, pbmc.atac$predicted.id)
predictions <- predictions/rowSums(predictions) # normalize for number of cells in each cell type
predictions <- as.data.frame(predictions)
p1 <- ggplot(predictions, aes(Var1, Var2, fill = Freq)) + geom_tile() + scale_fill_gradient(name = "Fraction of cells",
low = "#ffffc8", high = "#7d0025") + xlab("Cell type annotation (RNA)") + ylab("Predicted cell type label (ATAC)") +
theme_cowplot() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
correct <- length(which(pbmc.atac$seurat_annotations == pbmc.atac$predicted.id))
incorrect <- length(which(pbmc.atac$seurat_annotations != pbmc.atac$predicted.id))
data <- FetchData(pbmc.atac, vars = c("prediction.score.max", "annotation_correct"))
p2 <- ggplot(data, aes(prediction.score.max, fill = annotation_correct, colour = annotation_correct)) +
geom_density(alpha = 0.5) + theme_cowplot() + scale_fill_discrete(name = "Annotation Correct",
labels = c(paste0("FALSE (n = ", incorrect, ")"), paste0("TRUE (n = ", correct, ")"))) + scale_color_discrete(name = "Annotation Correct",
labels = c(paste0("FALSE (n = ", incorrect, ")"), paste0("TRUE (n = ", correct, ")"))) + xlab("Prediction Score")
p1 + p2predictions <- table(pbmc.atac$seurat_annotations, pbmc.atac$predicted.id)
predictions <- predictions/rowSums(predictions) # normalize for number of cells in each cell type
predictions <- as.data.frame(predictions)
p1 <- ggplot(predictions, aes(Var1, Var2, fill = Freq)) + geom_tile() + scale_fill_gradient(name = "Fraction of cells",
low = "#ffffc8", high = "#7d0025") + xlab("Cell type annotation (RNA)") + ylab("Predicted cell type label (ATAC)") +
theme_cowplot() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
correct <- length(which(pbmc.atac$seurat_annotations == pbmc.atac$predicted.id))
incorrect <- length(which(pbmc.atac$seurat_annotations != pbmc.atac$predicted.id))
data <- FetchData(pbmc.atac, vars = c("prediction.score.max", "annotation_correct"))
p2 <- ggplot(data, aes(prediction.score.max, fill = annotation_correct, colour = annotation_correct)) +
geom_density(alpha = 0.5) + theme_cowplot() + scale_fill_discrete(name = "Annotation Correct",
labels = c(paste0("FALSE (n = ", incorrect, ")"), paste0("TRUE (n = ", correct, ")"))) + scale_color_discrete(name = "Annotation Correct",
labels = c(paste0("FALSE (n = ", incorrect, ")"), paste0("TRUE (n = ", correct, ")"))) + xlab("Prediction Score")
p1 + p2
联合嵌入scRNA-seq和scATAC-seq数据集
除了跨模态转移标签外，还可以在同一图上可视化scrna-seq和scatac-seq细胞。 我们强调，此步骤主要用于可视化，并且是一个可选步骤。 通常，当我们在SCRNA-SEQ和SCATAC-SEQ数据集之间执行综合分析时，我们主要关注上述标签传输。 我们证明了下面共同插入的工作流程，并再次强调这是出于演示目的，尤其是在这种特殊情况下，SCRNA-SEQ轮廓和SCATAC-SEQ剖面实际上实际上在同一细胞中测量了。
为了进行共同嵌入，我们首先根据先前计算的锚点将RNA表达 "归入 "scATAC-seq细胞，然后合并数据集。
# note that we restrict the imputation to variable genes from scRNA-seq, but could impute the
# full transcriptome if we wanted to
genes.use <- VariableFeatures(pbmc.rna)
refdata <- GetAssayData(pbmc.rna, assay = "RNA", slot = "data")[genes.use, ]
# refdata (input) contains a scRNA-seq expression matrix for the scRNA-seq cells. imputation
# (output) will contain an imputed scRNA-seq matrix for each of the ATAC cells
imputation <- TransferData(anchorset = transfer.anchors, refdata = refdata, weight.reduction = pbmc.atac[["lsi"]],
dims = 2:30)
pbmc.atac[["RNA"]] <- imputation
coembed <- merge(x = pbmc.rna, y = pbmc.atac)
# Finally, we run PCA and UMAP on this combined object, to visualize the co-embedding of both
# datasets
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)
coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)
DimPlot(coembed, group.by = c("orig.ident", "seurat_annotations"))
总结
作者也提到了在这节的处理中是将来源于同一barcode的数据进行了不同的整合处理，可以参考wnn的处理方法，目前也会在后面对这个内容进行解析。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。
原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。
如有侵权，请联系 cloudcommunity@tencent.com 删除。