Seurat4.0系列教程14:整合scRNA-seq and scATAC-seq数据
单细胞转录学改变了我们描述细胞状态的能力，但深入的生物学解释需要的不仅仅是分群。随着测量不同细胞模式的新方法的出现，一个关键的分析挑战是整合这些数据集，以更好地了解细胞身份和功能。例如，用户可以在同一生物系统上执行 scRNA-seq 和 scATAC-seq 实验，并一致地用同一组细胞类型标签对两个数据集进行注释。这种分析尤其具有挑战性，因为 scATAC-seq 数据集难以注释，单细胞分辨率收集的基因组数据很少，而且 scRNA-seq 数据中缺乏可解释的基因标记。
我们介绍了一种方法，以整合从同一生物系统收集的 scRNA-seq 和 scATAC-seq 数据集，特别演示了以下分析：
如何使用带注释的 scRNA-seq 数据集来标记 scATAC-seq 实验中的细胞 如何共同可视化从 scRNA-seq 和 scATAC-seq 中得到的细胞 如何将 scATAC-seq 细胞投射到从 scRNA-seq 实验中提取的 UMAP 上
我们使用公开提供的 12，000 个人类 PBMC"多组学"数据集来演示这些方法。在此数据集中，scRNA-seq 和 scATAC-seq 配置文件同时收集在同一个细胞中。我们把数据集视为来自两个不同的实验，并把它们整合在一起。由于它们最初是在同一细胞中测量的，我们可以用它来评估整合的准确性。
加载数据并单独处理每个模式
PBMC多组数据集可从10x Genomics[1].装载。我们单独加载RNA和ATAC数据，并假设这些配置文件是在单独的实验中测量的。
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 实验之间的"锚点"，我们首先利用 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))# Identify anchors
transfer.anchors <- FindTransferAnchors(reference = pbmc.rna, query = pbmc.atac, features = VariableFeatures(object = pbmc.rna),
reference.assay = "RNA", query.assay = "ACTIVITY", reduction = "cca")
通过标签转移注释 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)
执行转移后，ATAC-seq 细胞有了预测的注释，存储在该predicted.id中。
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
在此示例中，通过 scRNA-seq 整合正确预测了大约 90%scATAC-seq 配置文件的注释 。此外，prediction.score.max量化了与我们预测的注释相关的不确定性。我们可以看到正确注释的细胞通常与高预测分数（>90%）相关，而错误注释的预测分数较低。不正确的分配也往往反映在密切相关的细胞类型（如中间态细胞与Naive B细胞）。
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 + 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"))
参考资料
10x Genomics: https://support.10xgenomics.com/single-cell-multiome-atac-gex/datasets/1.0.0/pbmc_granulocyte_sorted_10k