10月20日，Satija Lab發(fā)布了最新教程：Weighted Nearest Neighbor Analysis
在Seurat V4中使用加權(quán)最近鄰法（weighted nearest neighbor, WNN）分析多模態(tài)單細(xì)胞數(shù)據(jù)：
聯(lián)合分析CITE-seq（RNA +蛋白質(zhì)）或10x multiome（RNA + ATAC）數(shù)據(jù)
執(zhí)行多模態(tài)聚類(lèi)
構(gòu)建聯(lián)合可視化
多模態(tài)數(shù)據(jù)整合和WNN算法原理可以查詢預(yù)印本文章：
Integrated analysis of multimodal single-cell data
https://www.biorxiv.org/content/10.1101/2020.10.12.335331v1
該Vignette有BMNC - RNA & ADT和PBMC - RNA & ATAC兩個(gè)教程，本文介紹了第一個(gè)BMNC - RNA & ADT，整合分析單細(xì)胞轉(zhuǎn)錄組和蛋白組數(shù)據(jù)。
BMNC - RNA & ADT
本段介紹了多模態(tài)單細(xì)胞數(shù)據(jù)集分析的WNN工作流
工作流由三個(gè)步驟組成：
每個(gè)模態(tài)獨(dú)立的預(yù)處理和降維
學(xué)習(xí)細(xì)胞特定的模態(tài)“權(quán)重”，并構(gòu)建一個(gè)集成了這些模態(tài)的WNN圖
WNN圖的下游分析(如可視化、聚類(lèi)等)
我們使用（Stuart, Butler et al, Cell 2019）的CITE-seq數(shù)據(jù)集，該數(shù)據(jù)集包含30,672個(gè)scRNA-seq圖譜和25個(gè)抗體，包含兩個(gè)assay：RNA和抗體衍生標(biāo)簽（ADT）
1. R包安裝和數(shù)據(jù)下載
在運(yùn)行前先安裝Seurat V4，在作者github頁(yè)面上有beta版本。
remotes::install_github("satijalab/seurat", ref = "release/4.0.0")
結(jié)果報(bào)錯(cuò)：
1.jpg
報(bào)錯(cuò)namespace 'sctransform' 0.2.1 is being loaded, but >= 0.3.1 is required 顯示sctransform包有問(wèn)題，需要更新最新版0.3.1，然后本地安裝：
install.packages("D:/Program Files/R/R-3.6.1/library/sctransform_0.3.1.tar.gz",repos = NULL)
sctransform安裝成功：
2.png
然后，再次安裝Seurat v4：
remotes::install_github("satijalab/seurat", ref = "release/4.0.0")
Seurat v4安裝成功！
3.png
接下來(lái)開(kāi)始一系列操作，載入必要的包、載入數(shù)據(jù)：
library(Seurat)
library(SeuratData)
library(cowplot)
library(dplyr)
InstallData("bmcite") # 99.7 MB
bm <- LoadData(ds = "bmcite")
2. 預(yù)處理和降維
首先分別對(duì)RNA和ADT這兩種assays進(jìn)行預(yù)處理和降維。這里使用標(biāo)準(zhǔn)歸一化流程，但也可以使用SCTransform或任何替代方法。
DefaultAssay(bm) <- 'RNA'
bm <- NormalizeData(bm) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()
DefaultAssay(bm) <- 'ADT's
# 使用所有ADT 特征進(jìn)行降維
# 這里設(shè)置了一個(gè)降維名，防止改寫(xiě)
VariableFeatures(bm) <- rownames(bm[["ADT"]])
bm <- NormalizeData(bm, normalization.method = 'CLR', margin = 2) %>% # CLR：中心對(duì)數(shù)比 margin：當(dāng)執(zhí)行CLR時(shí)，對(duì)特征（1）或細(xì)胞（2）進(jìn)行標(biāo)準(zhǔn)化
ScaleData() %>%
RunPCA(reduction.name = 'apca')
3. 計(jì)算多模態(tài)近鄰
我們根據(jù)RNA和蛋白質(zhì)相似度的加權(quán)組合來(lái)計(jì)算每個(gè)細(xì)胞在數(shù)據(jù)集中的最近鄰。利用FindMultiModalNeighbors函數(shù)中計(jì)算細(xì)胞特異性模態(tài)權(quán)重和多模態(tài)近鄰，在此數(shù)據(jù)集上運(yùn)行約 2 min。我們需要指定每個(gè)模態(tài)的維數(shù)（類(lèi)似于指定scRNA-seq集群中的PC數(shù)量)，并且可以改變這些設(shè)置，會(huì)發(fā)現(xiàn)小改動(dòng)對(duì)總體結(jié)果的影響極小。
bm <- FindMultiModalNeighbors( bm
, reduction.list = list("pca", "apca")
, dims.list = list(1:30, 1:18)
, modality.weight.name = "RNA.weight"
)
識(shí)別的多模態(tài)近鄰將存儲(chǔ)于neighbors slot中，可通過(guò)
bm[['weighted.nn']]訪問(wèn)
通過(guò)bm[["wknn"]]訪問(wèn)WNN圖；通過(guò)bm[[“wsnn”]]訪問(wèn)用于聚類(lèi)的SNN圖
通過(guò)bm$RNA.weight訪問(wèn)細(xì)胞特異性模態(tài)權(quán)重
4. 下游分析：可視化和聚類(lèi)
4.1 多個(gè)模態(tài)數(shù)據(jù)聚類(lèi)和可視化
我們現(xiàn)在可以使用這些結(jié)果進(jìn)行下游分析，比如可視化和聚類(lèi)。
例如，可以基于RNA和蛋白質(zhì)數(shù)據(jù)的加權(quán)組合創(chuàng)建數(shù)據(jù)的UMAP可視化。
還可以執(zhí)行基于圖形的聚類(lèi)，并在UMAP上可視化這些結(jié)果，以及細(xì)胞注釋。
（類(lèi)似于Seurat V3提出的共嵌入分析，將兩個(gè)模態(tài)數(shù)據(jù)一同展示和聚類(lèi)）
bm <- RunUMAP(bm, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
bm <- FindClusters(bm, graph.name = "wsnn", algorithm = 3, resolution = 2, verbose = FALSE)
p1 <- DimPlot(bm, reduction = 'wnn.umap', label = TRUE, repel = TRUE, label.size = 2.5) + NoLegend()
p2 <- DimPlot(bm, reduction = 'wnn.umap', group.by = 'celltype.l2', label = TRUE, repel = TRUE, label.size = 2.5) + NoLegend()
p1 + p2
多模態(tài)數(shù)據(jù)UMAP圖P1按聚類(lèi)結(jié)果顯示，P2按細(xì)胞類(lèi)型顯示
4.2 單模態(tài)數(shù)據(jù)的可視化和比較
我們也可以對(duì)RNA和蛋白質(zhì)數(shù)據(jù)分別UMAP可視化，并進(jìn)行比較。
我們發(fā)現(xiàn)，RNA分析相比于ADT分析在識(shí)別祖細(xì)胞狀態(tài)方面能夠提供更多的信息，表現(xiàn)在
ADT面板包含分化細(xì)胞的markers，而T細(xì)胞狀態(tài)則相反（ADT分析優(yōu)于RNA分析)。
bm <- RunUMAP(bm, reduction = 'pca', dims = 1:30, assay = 'RNA',
reduction.name = 'rna.umap', reduction.key = 'rnaUMAP_')
bm <- RunUMAP(bm, reduction = 'apca', dims = 1:18, assay = 'ADT',
reduction.name = 'adt.umap', redruction.key = 'adtUMAP_')
p3 <- DimPlot(bm, reduction = 'rna.umap', group.by = 'celltype.l2', label = TRUE,
repel = TRUE, label.size = 2.5) + NoLegend()
p4 <- DimPlot(bm, reduction = 'adt.umap', group.by = 'celltype.l2', label = TRUE,
repel = TRUE, label.size = 2.5) + NoLegend()
p3 + p4
RNA（左）和ADT（右）的UMAP圖
4.3 marker基因和蛋白表達(dá)的可視化
我們還可以在多模態(tài)UMAP圖上可視化一些典型marker gene和蛋白的表達(dá)，有助于驗(yàn)證所提供的注釋?zhuān)?/p>
p5 <- FeaturePlot(bm, features = c("adt_CD45RA","adt_CD16","adt_CD161"),
reduction = 'wnn.umap', max.cutoff = 2,
colrs = c("lightgrey","red"), ncol = 3)
p6 <- FeaturePlot(bm, features = c("rna_TRDC","rna_MPO","rna_AVP"),
reduction = 'wnn.umap', max.cutoff = 3, ncol = 3)
p5 / p6
eg：Naive CD4+ T cell, B cell—CD45RA；NK cell—CD16；NKT cell—TRDC；GMP cell（Granulocyte-monocyte progenitor）—MPOThe expression of canonical marker genes and proteins on the multimodal UMAP
以上我是通過(guò)細(xì)胞標(biāo)志物數(shù)據(jù)庫(kù)——CellMarker查詢的。
4.4 模態(tài)權(quán)重的可視化
最后，我們可視化每個(gè)細(xì)胞的模態(tài)權(quán)重。RNA權(quán)重最高的群體代表祖細(xì)胞，而蛋白質(zhì)權(quán)重最高的群體代表T細(xì)胞。這符合我們的生物學(xué)預(yù)期，因?yàn)榭贵w中不包含可以區(qū)分不同祖細(xì)胞群體的marker。
VlnPlot(bm, features = "RNA.weight", group.by = 'celltype.l2', sort = TRUE, pt.size = 0.1) + NoLegend()
模態(tài)權(quán)重的可視化
從最初的CCA到CCA+Anchor到4.0版本的WNN，從多源的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)整合到多個(gè)模態(tài)的單細(xì)胞數(shù)據(jù)整合，真的是十分敬佩Satija實(shí)驗(yàn)室，年年高產(chǎn)，驚喜不斷。
