引言

本文介紹了如何在Seurat軟件中將查詢數(shù)據(jù)集與經(jīng)過注釋的參考數(shù)據(jù)集進(jìn)行匹配。以一個(gè)實(shí)例來說，我們把10X Genomics公司早期發(fā)布的一個(gè)包含2700個(gè)外周血單核細(xì)胞（PBMC）的單細(xì)胞RNA測(cè)序（scRNA-seq）數(shù)據(jù)集，與我們最近創(chuàng)建的一個(gè)使用228種抗體測(cè)量的、包含162,000個(gè)PBMC的CITE-seq參考數(shù)據(jù)集進(jìn)行匹配。這個(gè)例子用來說明，在參考數(shù)據(jù)集的幫助下進(jìn)行的有監(jiān)督分析，是如何幫助我們識(shí)別那些僅通過無監(jiān)督分析難以發(fā)現(xiàn)的細(xì)胞狀態(tài)。在另一個(gè)例子中，我們展示了如何將來自不同個(gè)體的人類骨髓細(xì)胞（Human BMNC）的人類細(xì)胞圖譜（Human Cell Atlas）數(shù)據(jù)集，有序地映射到一個(gè)統(tǒng)一的參考框架上。

我們之前利用參考映射的方法來標(biāo)注查詢數(shù)據(jù)集中的細(xì)胞標(biāo)簽。在Seurat v4版本中，大幅提高了執(zhí)行集成任務(wù)，包括參考映射的速度和內(nèi)存效率，并且還新增了將查詢細(xì)胞投影到之前計(jì)算好的UMAP（Uniform Manifold Approximation and Projection，均勻流形近似和投影）可視化界面的功能。

內(nèi)容

在本示例中，我們將展示如何利用一個(gè)已經(jīng)建立的參考數(shù)據(jù)集來解讀單細(xì)胞RNA測(cè)序（scRNA-seq）查詢：

1. 根據(jù)參考數(shù)據(jù)集定義的細(xì)胞狀態(tài)集，對(duì)每個(gè)查詢細(xì)胞進(jìn)行標(biāo)注。
2. 將每個(gè)查詢細(xì)胞投影到之前計(jì)算完成的UMAP可視化界面上。
3. 估算在CITE-seq參考數(shù)據(jù)集中測(cè)量到的表面蛋白的預(yù)測(cè)水平。

要運(yùn)行本示例，請(qǐng)確保安裝了Seurat v4，該軟件可在CRAN上下載。同時(shí)，您還需要安裝SeuratDisk包。

library(Seurat)
library(ggplot2)
library(patchwork)

options(SeuratData.repo.use = "http://seurat.nygenome.org")

Example 1：繪制人類外周血細(xì)胞圖譜

Data

我們從最近的論文中加載 reference，并可視化預(yù)先計(jì)算的 UMAP。

reference <- readRDS("/brahms/hartmana/vignette_data/pbmc_multimodal_2023.rds")

DimPlot(object = reference, reduction = "wnn.umap", group.by = "celltype.l2", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()

Mapping

為了演示與此多模式參考的映射，我們將使用由 10x Genomics 生成并可通過 SeuratData 獲取的 2,700 個(gè) PBMC 數(shù)據(jù)集。

library(SeuratData)
InstallData('pbmc3k')

pbmc3k <- LoadData('pbmc3k')
pbmc3k <- UpdateSeuratObject(pbmc3k)

參考集是通過 SCTransform() 函數(shù)進(jìn)行標(biāo)準(zhǔn)化處理的，因此我們同樣采用這一方法來對(duì)查詢集進(jìn)行標(biāo)準(zhǔn)化處理。

pbmc3k <- SCTransform(pbmc3k, verbose = FALSE)

接下來，我們?cè)趨⒖技c查詢集之間確立對(duì)應(yīng)關(guān)系。根據(jù)論文中的描述，本例中我們采用了預(yù)先計(jì)算的監(jiān)督主成分分析（Supervised PCA，簡(jiǎn)稱spca）變換。我們建議對(duì)CITE-seq數(shù)據(jù)集采用監(jiān)督主成分分析方法，并將在本指南的下一個(gè)部分展示如何執(zhí)行這一變換。當(dāng)然，您也可以選擇使用傳統(tǒng)的主成分分析（PCA）變換。

anchors <- FindTransferAnchors(
  reference = reference,
  query = pbmc3k,
  normalization.method = "SCT",
  reference.reduction = "spca",
  dims = 1:50
)

接下來，我們將細(xì)胞類型的標(biāo)簽和蛋白質(zhì)信息從參考集轉(zhuǎn)移到查詢集。同時(shí)，我們還把查詢集的數(shù)據(jù)映射到參考集的UMAP（均勻流形逼近和投影）結(jié)構(gòu)上。

pbmc3k <- MapQuery(
  anchorset = anchors,
  query = pbmc3k,
  reference = reference,
  refdata = list(
    celltype.l1 = "celltype.l1",
    celltype.l2 = "celltype.l2",
    predicted_ADT = "ADT"
  ),
  reference.reduction = "spca", 
  reduction.model = "wnn.umap"
)

• MapQuery 在做什么？

MapQuery() 函數(shù)封裝了三個(gè)子函數(shù)：TransferData()、IntegrateEmbeddings() 和 ProjectUMAP()。TransferData() 函數(shù)的作用是傳遞細(xì)胞類型的標(biāo)簽并估算 ADT（Ambient Dolly Technology）的數(shù)值。而 IntegrateEmbeddings() 和 ProjectUMAP() 函數(shù)則用于將查詢數(shù)據(jù)集映射到參考數(shù)據(jù)集的 UMAP（均勻流形逼近和投影）結(jié)構(gòu)上。以下是使用這些中間函數(shù)完成相同操作的代碼示例：

pbmc3k <- TransferData(
  anchorset = anchors, 
  reference = reference,
  query = pbmc3k,
  refdata = list(
    celltype.l1 = "celltype.l1",
    celltype.l2 = "celltype.l2",
    predicted_ADT = "ADT")
)
pbmc3k <- IntegrateEmbeddings(
  anchorset = anchors,
  reference = reference,
  query = pbmc3k, 
  new.reduction.name = "ref.spca"
)
pbmc3k <- ProjectUMAP(
  query = pbmc3k, 
  query.reduction = "ref.spca", 
  reference = reference, 
  reference.reduction = "spca", 
  reduction.model = "wnn.umap"
)

探索映射結(jié)果

目前，我們能夠?qū)?700個(gè)查詢細(xì)胞進(jìn)行可視化展示。這些細(xì)胞已經(jīng)被映射到由參考數(shù)據(jù)集定義的UMAP（均勻流形逼近和投影）視圖上，并且每個(gè)細(xì)胞都獲得了兩個(gè)不同詳細(xì)程度（第一級(jí)和第二級(jí)）的注釋信息。

p1 = DimPlot(pbmc3k, reduction = "ref.umap", group.by = "predicted.celltype.l1", label = TRUE, label.size = 3, repel = TRUE) + NoLegend()
p2 = DimPlot(pbmc3k, reduction = "ref.umap", group.by = "predicted.celltype.l2", label = TRUE, label.size = 3 ,repel = TRUE) + NoLegend()
p1 + p2

通過參考映射數(shù)據(jù)集，我們能夠辨識(shí)出在對(duì)查詢數(shù)據(jù)集進(jìn)行無監(jiān)督分析時(shí)難以區(qū)分的細(xì)胞類型。例如，這包括了漿細(xì)胞樣樹突狀細(xì)胞（pDC）、造血干細(xì)胞和祖細(xì)胞（HSPC）、調(diào)節(jié)性T細(xì)胞（Treg）、CD8 初始T細(xì)胞、CD56+ 自然殺傷（NK）細(xì)胞、記憶B細(xì)胞和初始B細(xì)胞以及漿細(xì)胞。對(duì)于每個(gè)預(yù)測(cè)的細(xì)胞類型，系統(tǒng)都會(huì)給出一個(gè)0到1范圍內(nèi)的評(píng)分。

FeaturePlot(pbmc3k, features = c("pDC", "CD16 Mono", "Treg"),  reduction = "ref.umap", cols = c("lightgrey", "darkred"), ncol = 3) & theme(plot.title = element_text(size = 10))

library(ggplot2)
plot <- FeaturePlot(pbmc3k, features = "CD16 Mono",  reduction = "ref.umap", cols = c("lightgrey", "darkred"))  + ggtitle("CD16 Mono") + theme(plot.title = element_text(hjust = 0.5, size = 30)) + labs(color = "Prediction Score") +  xlab("UMAP 1") + ylab("UMAP 2") + 
  theme(axis.title = element_text(size = 18), legend.text = element_text(size = 18), legend.title = element_text(size = 25)) 
ggsave(filename = "../output/images/multimodal_reference_mapping.jpg", height = 7, width = 12, plot = plot, quality = 50)

我們可以通過檢查特定標(biāo)記基因的表達(dá)模式來核實(shí)我們的預(yù)測(cè)結(jié)果。舉例來說，已有研究指出CLEC4C和LIR4是漿細(xì)胞樣樹突狀細(xì)胞（pDC）的標(biāo)志性基因，這與我們的預(yù)測(cè)相符。同樣，如果我們通過差異表達(dá)分析來篩選調(diào)節(jié)性T細(xì)胞（Treg）的標(biāo)記，我們能夠識(shí)別出一組標(biāo)準(zhǔn)標(biāo)記基因，包括RTKN2、CTLA4、FOXP3和IL2RA。

Idents(pbmc3k) <- 'predicted.celltype.l2'
VlnPlot(pbmc3k, features = c("CLEC4C", "LILRA4"), sort = TRUE) + NoLegend()

treg_markers <- FindMarkers(pbmc3k, ident.1 = "Treg", only.pos = TRUE, logfc.threshold = 0.1)
print(head(treg_markers))


##                       p_val avg_log2FC pct.1 pct.2    p_val_adj
## AC004854.4     4.795409e-26   7.800900 0.042 0.000 6.028789e-22
## RP11-297N6.4   4.795409e-26   7.800900 0.042 0.000 6.028789e-22
## RTKN2          4.795409e-26   7.800900 0.042 0.000 6.028789e-22
## CTD-2020K17.4  4.795409e-26   7.800900 0.042 0.000 6.028789e-22
## RP11-1399P15.1 2.782841e-18   4.172869 0.292 0.021 3.498588e-14
## IL2RA          1.808612e-14   4.137935 0.208 0.014 2.273787e-10

最終，我們能夠展示根據(jù) CITE-seq 參考數(shù)據(jù)推斷出的表面蛋白表達(dá)水平的可視化結(jié)果。

DefaultAssay(pbmc3k) <- 'predicted_ADT'
# see a list of proteins: rownames(pbmc3k)
FeaturePlot(pbmc3k, features = c("CD3-1", "CD45RA", "IgD"), reduction = "ref.umap", cols = c("lightgrey", "darkgreen"), ncol = 3)

計(jì)算新的 UMAP 可視化

在之前的案例中，我們將查詢細(xì)胞映射到基于參考數(shù)據(jù)集生成的 UMAP 結(jié)構(gòu)后進(jìn)行了可視化展示。維持一致的可視化方式有助于我們解讀新的數(shù)據(jù)集。然而，如果查詢數(shù)據(jù)集中包含了在參考數(shù)據(jù)集中未出現(xiàn)的細(xì)胞狀態(tài)，這些狀態(tài)將會(huì)映射到參考數(shù)據(jù)集中最為相似的細(xì)胞類型上。這是 UMAP 軟件包所預(yù)期的工作方式，但有時(shí)也可能因此忽視了查詢數(shù)據(jù)集中存在的、可能具有研究?jī)r(jià)值的新細(xì)胞類型。

在我們的論文中，我們對(duì)一個(gè)包含發(fā)育中和已分化的中性粒細(xì)胞的查詢數(shù)據(jù)集進(jìn)行了映射，這些細(xì)胞在我們的參考數(shù)據(jù)集中并未包含。我們發(fā)現(xiàn)，在將參考數(shù)據(jù)集和查詢數(shù)據(jù)集合并后，重新計(jì)算一個(gè)新的 UMAP（稱為“從頭可視化”）有助于識(shí)別這些細(xì)胞群體，這一點(diǎn)在補(bǔ)充圖 8 中有展示。在“從頭可視化”中，查詢數(shù)據(jù)集中的獨(dú)特細(xì)胞狀態(tài)仍然保持獨(dú)立。在這個(gè)例子中，2700個(gè)外周血單核細(xì)胞（PBMC）并沒有包含獨(dú)特的細(xì)胞狀態(tài)，但我們將展示如何計(jì)算這種可視化。

我們想強(qiáng)調(diào)的是，如果用戶嘗試映射的數(shù)據(jù)集樣本不是 PBMC，或者包含了參考數(shù)據(jù)集中未出現(xiàn)細(xì)胞類型，那么進(jìn)行一次“從頭可視化”是理解和分析他們數(shù)據(jù)集的一個(gè)重要步驟。

reference <- DietSeurat(reference, counts = FALSE, dimreducs = "spca")
pbmc3k <- DietSeurat(pbmc3k, counts = FALSE, dimreducs = "ref.spca")

#merge reference and query
reference$id <- 'reference'
pbmc3k$id <- 'query'
refquery <- merge(reference, pbmc3k)
refquery[["spca"]] <- merge(reference[["spca"]], pbmc3k[["ref.spca"]])
refquery <- RunUMAP(refquery, reduction = 'spca', dims = 1:50)
DimPlot(refquery, group.by = 'id', shuffle = TRUE)