簡介

在單細(xì)胞基因組學(xué)領(lǐng)域，將新數(shù)據(jù)集映射到已建立的參考數(shù)據(jù)集上的能力，與讀取映射工具變革基因組序列分析的方式如出一轍。

一個關(guān)鍵的挑戰(zhàn)是將這種參考映射框架擴(kuò)展到不測量基因表達(dá)的技術(shù)，即使底層參考是基于單細(xì)胞RNA測序（scRNA-seq）的。在Hao et al, Nat Biotechnol 2023中，介紹了“橋接整合”（bridge integration），它使得將補(bǔ)充技術(shù)（如單細(xì)胞ATAC-seq（scATAC-seq）、單細(xì)胞DNA甲基化（scDNAme）、細(xì)胞因子分析（CyTOF））映射到scRNA-seq參考數(shù)據(jù)上成為可能，使用一個“多組學(xué)”數(shù)據(jù)集作為分子橋接。

在本文中，我們將展示如何將人類PBMC的scATAC-seq數(shù)據(jù)集映射到我們之前構(gòu)建的PBMC參考數(shù)據(jù)。我們使用一個公開可用的10x多組學(xué)數(shù)據(jù)集作為橋接數(shù)據(jù)集，該數(shù)據(jù)集在同一細(xì)胞中同時測量基因表達(dá)和染色質(zhì)可及性。

我們展示了以下內(nèi)容：

• 加載和預(yù)處理scATAC-seq、多組學(xué)和scRNA-seq參考數(shù)據(jù)集
• 通過橋接整合映射scATAC-seq數(shù)據(jù)集
• 探索和評估所得注釋

Azimuth ATAC用于橋接整合

用戶現(xiàn)在可以在Azimuth網(wǎng)站或R中使用新發(fā)布的Azimuth ATAC工作流程自動運(yùn)行PBMC和骨髓scATAC-seq查詢的橋接整合。有關(guān)在R中本地運(yùn)行的更多詳細(xì)信息，請參見此vignette中的ATAC數(shù)據(jù)部分。

library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)
library(dplyr)
library(ggplot2)

加載橋接、查詢和參考數(shù)據(jù)集

我們首先加載一個包含約12,000個健康捐獻(xiàn)者的PBMC的10x多組學(xué)數(shù)據(jù)集。該數(shù)據(jù)集在同一細(xì)胞中測量RNA-seq和ATAC-seq，并可從10x Genomics 此處下載。我們遵循Signac包vignette中的加載說明。請注意，使用Signac時，請確保您使用的是最新版本的Bioconductor，因?yàn)?a target="_blank">用戶在使用較舊的BioC版本時報告了錯誤。

加載和設(shè)置10x多組學(xué)對象

# 10x hdf5文件包含兩種數(shù)據(jù)類型。
inputdata.10x <- Read10X_h5("/brahms/hartmana/vignette_data/pbmc_cellranger_arc_2/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5")

# 提取RNA和ATAC數(shù)據(jù)
rna_counts <- inputdata.10x$`Gene Expression`
atac_counts <- inputdata.10x$Peaks

# 創(chuàng)建Seurat對象
obj.multi <- CreateSeuratObject(counts = rna_counts)

# 獲取線粒體基因的百分比
obj.multi[["percent.mt"]] <- PercentageFeatureSet(obj.multi, pattern = "^MT-")

# 添加ATAC-seq檢測
grange.counts <- StringToGRanges(rownames(atac_counts), sep = c(":", "-"))
grange.use <- seqnames(grange.counts) %in% standardChromosomes(grange.counts)
atac_counts <- atac_counts[as.vector(grange.use), ]

# 獲取基因注釋
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)

# 切換到UCSC風(fēng)格
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg38"

# 包含ATAC每個片段信息的文件
frag.file <- "/brahms/hartmana/vignette_data/pbmc_cellranger_arc_2/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz"

# 將ATAC-seq數(shù)據(jù)作為ChromatinAssay對象添加
chrom_assay <- CreateChromatinAssay(
  counts = atac_counts,
  sep = c(":", "-"),
  genome = 'hg38',
  fragments = frag.file,
  min.cells = 10,
  annotation = annotations)

# 將ATAC檢測添加到多組學(xué)對象
obj.multi[["ATAC"]] <- chrom_assay

# 根據(jù)QC指標(biāo)過濾ATAC數(shù)據(jù)
obj.multi <- subset(
  x = obj.multi,
  subset = nCount_ATAC < 7e4 &
    nCount_ATAC > 5e3 &
    nCount_RNA < 25000 &
    nCount_RNA > 1000 &
    percent.mt < 20)

scATAC-seq 查詢數(shù)據(jù)集代表來自健康捐贈者的約 10,000 個 PBMC。我們加載峰/細(xì)胞矩陣，存儲片段文件的路徑，并向?qū)ο筇砑踊蜃⑨專裱嘟M實(shí)驗(yàn)中 ATAC 數(shù)據(jù)的步驟。

加載并設(shè)置 10x scATAC-seq 查詢

我們從最近的論文中加載 reference:

obj.rna <- readRDS("/brahms/haoy/seurat4_pbmc/pbmc_multimodal_2023.rds")

映射scATAC-seq數(shù)據(jù)集使用橋接整合

現(xiàn)在我們已經(jīng)設(shè)置了參考、查詢和橋接數(shù)據(jù)集，我們可以開始整合。橋接數(shù)據(jù)集使得scRNA-seq參考和scATAC-seq查詢之間的轉(zhuǎn)換成為可能，有效地擴(kuò)展了參考數(shù)據(jù)集，使其能夠映射新的數(shù)據(jù)類型。

我們稱這個擴(kuò)展的參考為“擴(kuò)展參考”，首先設(shè)置它。請注意，您可以保存此函數(shù)的結(jié)果，而無需重新運(yùn)行即可映射多個scATAC-seq數(shù)據(jù)集。

首先，我們丟棄ATAC降維的第一個維度。

dims.atac <- 2:50
dims.rna <- 1:50
DefaultAssay(obj.multi) <- "RNA"
DefaultAssay(obj.rna) <- "SCT"
obj.rna.ext <- PrepareBridgeReference(
  reference = obj.rna,
  bridge = obj.multi,
  reference.reduction = "spca",
  reference.dims = dims.rna,
  normalization.method = "SCT")

現(xiàn)在，我們可以直接在擴(kuò)展參考和查詢對象之間找到錨點(diǎn)。我們使用FindBridgeTransferAnchors函數(shù)，它使用與轉(zhuǎn)換參考相同的字典來轉(zhuǎn)換查詢數(shù)據(jù)集，然后在該空間中識別錨點(diǎn)。該函數(shù)旨在模仿我們的FindTransferAnchors函數(shù)，但是要識別跨模態(tài)的對應(yīng)關(guān)系。

bridge.anchor <- FindBridgeTransferAnchors(
  extended.reference = obj.rna.ext,
  query = obj.atac,
  reduction = "lsiproject",
  dims = dims.atac)

一旦我們識別了錨點(diǎn)，我們就可以將查詢數(shù)據(jù)集映射到參考上。MapQuery函數(shù)與我們之前介紹的用于參考映射的函數(shù)相同。它從參考數(shù)據(jù)集中轉(zhuǎn)移細(xì)胞注釋，并且還在先前計算的UMAP嵌入上可視化查詢數(shù)據(jù)集。由于我們的參考數(shù)據(jù)集包含三個分辨率級別的細(xì)胞類型注釋（l1 - l3），我們可以將每個級別轉(zhuǎn)移到查詢數(shù)據(jù)集中。

obj.atac <- MapQuery(
  anchorset = bridge.anchor,
  reference = obj.rna.ext,
  query = obj.atac,
  refdata = list(
    l1 = "celltype.l1",
    l2 = "celltype.l2",
    l3 = "celltype.l3"
  ),
  reduction.model = "wnn.umap")

現(xiàn)在我們可以可視化結(jié)果，在參考UMAP嵌入上繪制基于其預(yù)測注釋的scATAC-seq細(xì)胞。您可以看到，每個scATAC-seq細(xì)胞都根據(jù)scRNA-seq定義的細(xì)胞本體學(xué)被分配了一個細(xì)胞名稱。

DimPlot(
  obj.atac, group.by = "predicted.l2",
  reduction = "ref.umap", label = TRUE
) + ggtitle("ATAC") + NoLegend()

評估映射

為了評估映射和細(xì)胞類型預(yù)測，我們首先看看預(yù)測的細(xì)胞類型標(biāo)簽是否與scATAC-seq數(shù)據(jù)集的無監(jiān)督分析一致。我們遵循scATAC-seq數(shù)據(jù)的標(biāo)準(zhǔn)無監(jiān)督處理工作流程：

obj.atac <- FindTopFeatures(obj.atac, min.cutoff = "q0")
obj.atac <- RunSVD(obj.atac)
obj.atac <- RunUMAP(obj.atac, reduction = "lsi", dims = 2:50)

現(xiàn)在，我們在無監(jiān)督UMAP嵌入上可視化預(yù)測的簇標(biāo)簽。我們可以看到，預(yù)測的簇標(biāo)簽（來自scRNA-seq參考）與scATAC-seq數(shù)據(jù)的結(jié)構(gòu)一致。然而，有些細(xì)胞類型（例如Treg）在無監(jiān)督分析中似乎沒有分開。這些可能是預(yù)測錯誤，或者是參考映射提供額外分辨率的情況。

DimPlot(obj.atac, group.by = "predicted.l2", reduction = "umap", label = FALSE)

最后，我們通過檢查scATAC-seq數(shù)據(jù)在典型位點(diǎn)的染色質(zhì)可及性剖面來驗(yàn)證預(yù)測的細(xì)胞類型。我們使用CoveragePlot函數(shù)在CD8A、FOXP3和RORC處可視化可及性模式，按照其預(yù)測標(biāo)簽對細(xì)胞進(jìn)行分組。我們在每種情況下都看到了預(yù)期的模式。例如，PAX5位點(diǎn)顯示出僅在B細(xì)胞中可訪問的峰，CD8A位點(diǎn)在CD8 T細(xì)胞亞群中也是如此。同樣，預(yù)測的Tregs中FOXP3的可及性為我們的預(yù)測準(zhǔn)確性提供了強(qiáng)有力的支持。

CoveragePlot(
  obj.atac, region = "PAX5", group.by = "predicted.l1",
  idents = c("B", "CD4 T", "Mono", "NK"), window = 200,
  extend.upstream = -150000
)

CoveragePlot(
  obj.atac, region = "CD8A", group.by = "predicted.l2",
  idents = c("CD8 Naive", "CD4 Naive", "CD4 TCM", "CD8 TCM"),
  extend.downstream = 5000, extend.upstream = 5000
)

CoveragePlot(
  obj.atac, region = "FOXP3", group.by = "predicted.l2",
  idents = c("CD4 Naive", "CD4 TCM", "CD4 TEM", "Treg"),
  extend.downstream = 0, extend.upstream = 0
)

CoveragePlot(
  obj.atac, region = "RORC", group.by = "predicted.l2",
  idents = c("CD8 Naive", "CD8 TEM", "CD8 TCM", "MAIT"),
  extend.downstream = 5000, extend.upstream = 5000
)

未完待續(xù)，歡迎關(guān)注！