單細(xì)胞分析：多模態(tài) reference mapping (2)

引言

本文介紹了如何在Seurat軟件中將查詢數(shù)據(jù)集與經(jīng)過注釋的參考數(shù)據(jù)集進(jìn)行匹配。我們展示了如何將來自不同個(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 2：繪制人類骨髓細(xì)胞圖譜

Data

例如，我們將由人類細(xì)胞圖譜項(xiàng)目生成的，來自八位不同捐獻(xiàn)者的人類骨髓單核細(xì)胞（BMNC）數(shù)據(jù)集進(jìn)行了映射。我們以之前使用加權(quán)最近鄰分析（WNN）方法分析過的人類BMNC的CITE-seq參考集作為比對(duì)標(biāo)準(zhǔn)。

本文除了展示與之前PBMC案例相同的參考映射功能外，還進(jìn)一步介紹了：

  • 如何構(gòu)建一個(gè)監(jiān)督的主成分分析（sPCA）轉(zhuǎn)換。
  • 如何將多個(gè)不同的數(shù)據(jù)集依次映射到同一個(gè)參考集上。
  • 采取哪些優(yōu)化措施來提高映射過程的速度。
# Both datasets are available through SeuratData
library(SeuratData)
#load reference data
InstallData("bmcite")
bm <- LoadData(ds = "bmcite")
#load query data
InstallData('hcabm40k')
hcabm40k <- LoadData(ds = "hcabm40k")

參考數(shù)據(jù)集構(gòu)建了一個(gè)加權(quán)最近鄰（WNN）圖，該圖體現(xiàn)了在本次CITE-seq實(shí)驗(yàn)中RNA和蛋白質(zhì)數(shù)據(jù)的加權(quán)整合情況。

基于這個(gè)WNN圖，我們可以生成一個(gè)UMAP（Uniform Manifold Approximation and Projection）的可視化表示。在計(jì)算過程中，我們?cè)O(shè)置參數(shù)return.model為TRUE，這樣就可以將待查詢的數(shù)據(jù)集映射到這個(gè)UMAP可視化空間中。

bm <- RunUMAP(bm, nn.name = "weighted.nn", reduction.name = "wnn.umap", 
              reduction.key = "wnnUMAP_", return.model = TRUE)
DimPlot(bm, group.by = "celltype.l2", reduction = "wnn.umap")

計(jì)算 sPCA 變換

如我們?cè)谡撐闹兴觯覀兪紫葓?zhí)行一個(gè)“監(jiān)督式”的主成分分析（PCA）。該分析旨在找出轉(zhuǎn)錄組數(shù)據(jù)的最佳轉(zhuǎn)換方式，以最準(zhǔn)確地反映加權(quán)最近鄰（WNN）圖中的結(jié)構(gòu)特征。通過這種方法，我們可以將蛋白質(zhì)和RNA的測(cè)量值進(jìn)行加權(quán)組合，以“指導(dǎo)”PCA的計(jì)算過程，從而凸顯出數(shù)據(jù)中最為重要的變異因素。一旦計(jì)算出這種轉(zhuǎn)換，就可以將其應(yīng)用到任何查詢數(shù)據(jù)集上。盡管我們也可以計(jì)算并應(yīng)用傳統(tǒng)的PCA投影，但在處理通過WNN分析構(gòu)建的多模態(tài)參考數(shù)據(jù)時(shí)，我們更推薦使用監(jiān)督式PCA（sPCA）。

sPCA的計(jì)算過程只需進(jìn)行一次，之后就可以快速地將其應(yīng)用到每一個(gè)查詢數(shù)據(jù)集上。

bm <- ScaleData(bm, assay = 'RNA')
bm <- RunSPCA(bm, assay = 'RNA', graph = 'wsnn')

計(jì)算緩存的鄰居索引

鑒于我們需要將多個(gè)查詢樣本與同一個(gè)參考集進(jìn)行比對(duì)，我們可以對(duì)那些僅與參考集相關(guān)的特定步驟進(jìn)行緩存處理。這個(gè)步驟雖然是可選的，但在處理多個(gè)樣本的映射時(shí)，它可以有效提升運(yùn)算速度。

我們首先在參考集的監(jiān)督式PCA（sPCA）空間內(nèi)計(jì)算出前50個(gè)最近鄰。然后，我們將這些信息保存在Seurat對(duì)象的spca.annoy.neighbors屬性中，并通過設(shè)置cache.index = TRUE來緩存annoy索引數(shù)據(jù)結(jié)構(gòu)。

bm <- FindNeighbors(
  object = bm,
  reduction = "spca",
  dims = 1:50,
  graph.name = "spca.annoy.neighbors", 
  k.param = 50,
  cache.index = TRUE,
  return.neighbor = TRUE,
  l2.norm = TRUE
)
  • 如何保存和加載緩存的煩惱索引？

如果您需要保存或加載一個(gè)利用 "annoy" 方法和啟用了緩存索引（通過設(shè)置 cache.index = TRUE）創(chuàng)建的 Neighbor 對(duì)象的緩存索引，可以使用 SaveAnnoyIndex() 和 LoadAnnoyIndex() 這兩個(gè)函數(shù)來完成。需要注意的是，這個(gè)索引不能通過常規(guī)方式保存到 RDS 或 RDA 文件，這意味著它不會(huì)在 R 會(huì)話重新啟動(dòng)或使用 saveRDS/readRDS 函數(shù)保存和讀取包含該索引的 Seurat 對(duì)象時(shí)被正確保留。因此，每次當(dāng) R 重新啟動(dòng)或者您從 RDS 文件加載參考 Seurat 對(duì)象時(shí)，都需要使用 LoadAnnoyIndex() 函數(shù)來重新將 Annoy 索引加載到 Neighbor 對(duì)象中。SaveAnnoyIndex() 函數(shù)生成的文件可以與參考 Seurat 對(duì)象一起分發(fā)，以便在需要時(shí)將其添加到參考對(duì)象中的 Neighbor 對(duì)象里。

bm[["spca.annoy.neighbors"]]

## A Neighbor object containing the 50 nearest neighbors for 30672 cells

SaveAnnoyIndex(object = bm[["spca.annoy.neighbors"]], file = "/brahms/shared/vignette-data/reftmp.idx")
bm[["spca.annoy.neighbors"]] <- LoadAnnoyIndex(object = bm[["spca.annoy.neighbors"]], file = "/brahms/shared/vignette-data/reftmp.idx")

查詢數(shù)據(jù)集預(yù)處理

本節(jié)我們將展示如何將來自多位捐獻(xiàn)者的骨髓樣本與一個(gè)多模態(tài)骨髓參考集進(jìn)行比對(duì)。這些待查詢的數(shù)據(jù)集來源于人類細(xì)胞圖譜（Human Cell Atlas，HCA）的免疫細(xì)胞圖譜中的骨髓數(shù)據(jù)集，可以通過SeuratData包訪問。提供的數(shù)據(jù)集是一個(gè)合并后的對(duì)象，涵蓋了8位捐獻(xiàn)者的數(shù)據(jù)。我們首先需要將這些數(shù)據(jù)拆分成8個(gè)獨(dú)立的Seurat對(duì)象，對(duì)應(yīng)每位捐獻(xiàn)者，然后分別進(jìn)行映射分析。

library(dplyr)
library(SeuratData)
InstallData('hcabm40k')
hcabm40k.batches <- SplitObject(hcabm40k, split.by = "orig.ident")

接下來，我們按照參考數(shù)據(jù)集的處理方式對(duì)查詢數(shù)據(jù)集進(jìn)行標(biāo)準(zhǔn)化處理。具體來說，參考數(shù)據(jù)集是通過NormalizeData()函數(shù)采用對(duì)數(shù)標(biāo)準(zhǔn)化的方法進(jìn)行處理的。如果參考數(shù)據(jù)集是利用SCTransform()函數(shù)進(jìn)行標(biāo)準(zhǔn)化的，那么查詢數(shù)據(jù)集同樣需要應(yīng)用SCTransform()函數(shù)來進(jìn)行標(biāo)準(zhǔn)化處理。

hcabm40k.batches <- lapply(X = hcabm40k.batches, FUN = NormalizeData, verbose = FALSE)

Mapping

接下來，我們?cè)诿课痪璜I(xiàn)者的數(shù)據(jù)集與多模態(tài)參考集之間確定錨點(diǎn)。為了縮短映射時(shí)間，我們采用了一種優(yōu)化的命令，該命令通過輸入預(yù)先計(jì)算好的參考鄰居集合，并關(guān)閉錨點(diǎn)篩選功能來實(shí)現(xiàn)效率提升。

anchors <- list()
for (i in 1:length(hcabm40k.batches)) {
  anchors[[i]] <- FindTransferAnchors(
    reference = bm,
    query = hcabm40k.batches[[i]],
    k.filter = NA,
    reference.reduction = "spca", 
    reference.neighbors = "spca.annoy.neighbors", 
    dims = 1:50
  )
}

然后我們單獨(dú)映射每個(gè)數(shù)據(jù)集。

for (i in 1:length(hcabm40k.batches)) {
  hcabm40k.batches[[i]] <- MapQuery(
    anchorset = anchors[[i]], 
    query = hcabm40k.batches[[i]],
    reference = bm, 
    refdata = list(
      celltype = "celltype.l2", 
      predicted_ADT = "ADT"),
    reference.reduction = "spca",
    reduction.model = "wnn.umap"
  )
}

探索映射結(jié)果

現(xiàn)在映射已完成，我們可以可視化各個(gè)對(duì)象的結(jié)果

p1 <- DimPlot(hcabm40k.batches[[1]], reduction = 'ref.umap', group.by = 'predicted.celltype', label.size = 3)
p2 <- DimPlot(hcabm40k.batches[[2]], reduction = 'ref.umap', group.by = 'predicted.celltype', label.size = 3)
p1 + p2 + plot_layout(guides = "collect")

我們還可以把所有的數(shù)據(jù)對(duì)象合并成一個(gè)統(tǒng)一的數(shù)據(jù)集。需要注意的是，這些數(shù)據(jù)對(duì)象都已經(jīng)通過參考集被整合到了一個(gè)共同的分析空間中。之后，我們就能夠?qū)⑦@些數(shù)據(jù)的分析結(jié)果一并展現(xiàn)出來。

# Merge the batches 
hcabm40k <- merge(hcabm40k.batches[[1]], hcabm40k.batches[2:length(hcabm40k.batches)], merge.dr = "ref.umap")
DimPlot(hcabm40k, reduction = "ref.umap", group.by =  "predicted.celltype", label = TRUE, repel = TRUE, label.size = 3) + NoLegend()

我們可以對(duì)查詢細(xì)胞中的基因表達(dá)模式、聚類預(yù)測(cè)得分以及（估算得到的）表面蛋白水平進(jìn)行可視化展示：

p3 <- FeaturePlot(hcabm40k, features = c("rna_TRDC", "rna_MPO", "rna_AVP"), reduction = 'ref.umap', 
                  max.cutoff = 3, ncol = 3)

# cell type prediction scores
DefaultAssay(hcabm40k) <- 'prediction.score.celltype'
p4 <- FeaturePlot(hcabm40k, features = c("CD16 Mono", "HSC", "Prog-RBC"), ncol = 3, 
                  cols = c("lightgrey", "darkred"))

# imputed protein levels
DefaultAssay(hcabm40k) <- 'predicted_ADT'
p5 <- FeaturePlot(hcabm40k, features = c("CD45RA", "CD16", "CD161"), reduction = 'ref.umap',
                  min.cutoff = 'q10', max.cutoff = 'q99', cols = c("lightgrey", "darkgreen") ,
                  ncol = 3)
p3 / p4 / p5
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容