单细胞转录组｜scATAC-seq 数据整合

引言

本文在此展示了如何将多个源自人类外周血单核细胞的单细胞染色质数据集进行整合。其中一个数据集是通过10x Genomics的多组学技术获得的，它涵盖了每个细胞的DNA可及性和基因表达数据。另一个数据集则是通过10x Genomics的单细胞ATAC测序(scATAC-seq)技术得到的，仅包含DNA可及性信息。

将利用Seurat软件包中的工具，基于两者共有的DNA可及性检测方法，将这两个数据集合并。此外，本文还将展示如何将连续变量（如基因表达）和分类变量（如细胞标签）的信息，从一个参考数据集转移到另一个查询数据集中。

预处理

在这里，将加载预处理的 PBMC 多组数据以及 PBMC scATAC-seq 数据：

library(Signac)
library(Seurat)
library(ggplot2)

# load the pre-processed multiome data
pbmc.multi <- readRDS("pbmc_multiome/pbmc_multiomic.rds")

# load the pre-processed atac data
pbmc.atac <- readRDS("pbmc_vignette/pbmc.rds")

进行单细胞染色质数据综合分析时，关键的第一步是确认每个数据集都包含了相同的特征。为此，对ATAC数据集中的多组学峰值进行了量化，以确保两个数据集之间存在共通之处。更多关于如何合并染色质检测的信息，请参阅合并指南。

# quantify multiome peaks in the scATAC-seq dataset
counts <- FeatureMatrix(
  fragments = Fragments(pbmc.atac),
  features = granges(pbmc.multi),
  cells = colnames(pbmc.atac)
)

# add new assay with multiome peaks
pbmc.atac[['ATAC']] <- CreateChromatinAssay(
  counts = counts,
  fragments = Fragments(pbmc.atac)
)

# compute LSI
DefaultAssay(pbmc.atac) <- "ATAC"
pbmc.atac <- FindTopFeatures(pbmc.atac, min.cutoff = 10)
pbmc.atac <- RunTFIDF(pbmc.atac)
pbmc.atac <- RunSVD(pbmc.atac)

然后，可以将多组学数据和单细胞ATAC数据集合并，注意到它们之间存在差异，这种差异似乎源于批次效应（即由实验和特定技术引起的变异）。

# first add dataset-identifying metadata
pbmc.atac$dataset <- "ATAC"
pbmc.multi$dataset <- "Multiome"

# merge
pbmc.combined <- merge(pbmc.atac, pbmc.multi)

# process the combined dataset
pbmc.combined <- FindTopFeatures(pbmc.combined, min.cutoff = 10)
pbmc.combined <- RunTFIDF(pbmc.combined)
pbmc.combined <- RunSVD(pbmc.combined)
pbmc.combined <- RunUMAP(pbmc.combined, reduction = "lsi", dims = 2:30)
p1 <- DimPlot(pbmc.combined, group.by = "dataset")

整合

为了在两个数据集之间建立整合的基准点，需要将它们映射到一个共同的低维空间。将通过设置reduction参数为"rlsi"，使用互反最小二乘投影（即让每个数据集映射到对方的最小二乘空间）来实现这一点。

与通常对单细胞RNA测序(scRNA-seq)数据进行的标准化数据矩阵整合不同，将使用IntegrateEmbeddings()函数，对跨数据集的低维细胞嵌入（即最小二乘坐标）进行整合。

# find integration anchors
integration.anchors <- FindIntegrationAnchors(
  object.list = list(pbmc.multi, pbmc.atac),
  anchor.features = rownames(pbmc.multi),
  reduction = "rlsi",
  dims = 2:30
)

# integrate LSI embeddings
integrated <- IntegrateEmbeddings(
  anchorset = integration.anchors,
  reductions = pbmc.combined[["lsi"]],
  new.reduction.name = "integrated_lsi",
  dims.to.integrate = 1:30
)

# create a new UMAP using the integrated embeddings
integrated <- RunUMAP(integrated, reduction = "integrated_lsi", dims = 2:30)
p2 <- DimPlot(integrated, group.by = "dataset")

最终，对比了合并与整合后的数据集结果，发现整合过程成功地去除了数据集中由技术差异引起的变异，同时保留了与细胞类型相关的（生物学上的）变异。

(p1 + ggtitle("Merged")) | (p2 + ggtitle("Integrated"))

在这里，演示了使用两个数据集的整合方法，但可以应用相同的工作流程来整合任意数量的数据集。

Reference mapping

当手头有一个数据量庞大且质量上乘的数据集，或者是一个包含了其他数据集中没有的独特信息（比如细胞类型的标注或额外的数据形式）的数据集时，通常希望将其作为基准数据集，并把查询数据集映射到它上面，这样就能在已有基准的背景下对查询数据集进行解读。

为了展示如何利用单细胞染色质的基准数据集和查询数据集来实现这一点，将把PBMC多组学数据集作为基准，并利用Seurat软件包中的FindTransferAnchors()和MapQuery()函数，将scATAC-seq数据集映射到这个基准数据集上。

# compute UMAP and store the UMAP model
pbmc.multi <- RunUMAP(pbmc.multi, reduction = "lsi", dims = 2:30, return.model = TRUE)

# find transfer anchors
transfer.anchors <- FindTransferAnchors(
  reference = pbmc.multi,
  query = pbmc.atac,
  reference.reduction = "lsi",
  reduction = "lsiproject",
  dims = 2:30
)

# map query onto the reference dataset
pbmc.atac <- MapQuery(
  anchorset = transfer.anchors,
  reference = pbmc.multi,
  query = pbmc.atac,
  refdata = pbmc.multi$predicted.id,
  reference.reduction = "lsi",
  new.reduction.name = "ref.lsi",
  reduction.model = 'umap'
)

• MapQuery() 函数的作用是什么？

MapQuery() 是一个辅助函数，它自动执行 TransferData()、IntegrateEmbeddings() 和 ProjectUMAP() 这三个步骤，用于处理查询数据集，并根据生成锚点对象的方式自动设定合适的默认参数。如果你需要更细致地调整这些函数的参数，你可以通过 MapQuery() 传递参数到每个具体的函数，使用 MapQuery() 提供的 transferdata.args、integrateembeddings.args 和 projectumap.args 参数，或者你也可以直接单独运行这些函数。例如：

pbmc.atac <- TransferData(
  anchorset = transfer.anchors, 
  reference = pbmc.multi,
  weight.reduction = "lsiproject",
  query = pbmc.atac,
  refdata = list(
    celltype = "predicted.id",
    predicted_RNA = "RNA")
)
pbmc.atac <- IntegrateEmbeddings(
  anchorset = transfer.anchors,
  reference = pbmc.multi,
  query = pbmc.atac, 
  reductions = "lsiproject",
  new.reduction.name = "ref.lsi"
)
pbmc.atac <- ProjectUMAP(
  query = pbmc.atac, 
  query.reduction = "ref.lsi",
  reference = pbmc.multi, 
  reference.reduction = "lsi",
  reduction.model = "umap"
)

执行了 MapQuery() 函数后，成功地将单细胞ATAC测序(scATAC-seq)数据集与多模态基准数据集进行了映射，并且实现了从基准数据集到查询数据集的细胞类型标签传递。现在，可以将这些映射结果进行可视化，并查看与 scATAC-seq 数据集现在关联的细胞类型标签：

p1 <- DimPlot(pbmc.multi, reduction = "umap", group.by = "predicted.id", label = TRUE, repel = TRUE) + NoLegend() + ggtitle("Reference")
p2 <- DimPlot(pbmc.atac, reduction = "ref.umap", group.by = "predicted.id", label = TRUE, repel = TRUE) + NoLegend() + ggtitle("Query")

p1 | p2

RNA 数据填充

之前已经将分类信息（即细胞标签）转移，并把查询数据映射到了已有的参考 UMAP 上。同样，也可以将连续数据从参考数据集转移到查询数据集中。这里，展示了如何将基因表达数据从 PBMC 多组学数据集（该数据集在同一细胞中同时测量了 DNA 可及性和基因表达）转移到只测量了 DNA 可及性的 PBMC scATAC-seq 数据集中。需要注意的是，也可以通过在上述 MapQuery() 函数调用中设置 refdata 参数为一个值列表，来实现这些数据的转移。

# predict gene expression values
rna <- TransferData(
  anchorset = transfer.anchors,
  refdata = LayerData(pbmc.multi, assay = "SCT", layer = "data"),
  weight.reduction = pbmc.atac[["lsi"]],
  dims = 2:30
)

# add predicted values as a new assay
pbmc.atac[["predicted"]] <- rna

可以查看一些免疫标记基因，发现预测的表达模式与基于已知表达模式的预期相符。

DefaultAssay(pbmc.atac) <- "predicted"

FeaturePlot(
  object = pbmc.atac,
  features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),
  pt.size = 0.1,
  max.cutoff = 'q95',
  reduction = "ref.umap",
  ncol = 3
)