Seurat4.0系列教程14:整合scRNA-seq and scATAC-seq数据

生信技能树jimmy

发布于 2022-01-10 08:57:46

3K1

发布于 2022-01-10 08:57:46

文章被收录于专栏：单细胞天地

单细胞转录学改变了我们描述细胞状态的能力，但深入的生物学解释需要的不仅仅是分群。随着测量不同细胞模式的新方法的出现，一个关键的分析挑战是整合这些数据集，以更好地了解细胞身份和功能。例如，用户可以在同一生物系统上执行 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"))