Seurat v4教程：Weighted Nearest Neighbor Analysis（一）

10月20日，Satija Lab发布了最新教程：Weighted Nearest Neighbor Analysis
在Seurat V4中使用加权最近邻法（weighted nearest neighbor, WNN）分析多模态单细胞数据：

联合分析CITE-seq（RNA +蛋白质）或10x multiome（RNA + ATAC）数据

执行多模态聚类

构建联合可视化

多模态数据整合和WNN算法原理可以查询预印本文章：
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两个教程，本文介绍了第一个BMNC - RNA & ADT，整合分析单细胞转录组和蛋白组数据。

BMNC - RNA & ADT

本段介绍了多模态单细胞数据集分析的WNN工作流
工作流由三个步骤组成：

每个模态独立的预处理和降维

学习细胞特定的模态“权重”，并构建一个集成了这些模态的WNN图

WNN图的下游分析(如可视化、聚类等)

我们使用（Stuart, Butler et al, Cell 2019）的CITE-seq数据集，该数据集包含30,672个scRNA-seq图谱和25个抗体，包含两个assay：RNA和抗体衍生标签（ADT）

1. R包安装和数据下载

在运行前先安装Seurat V4，在作者github页面上有beta版本。

remotes::install_github("satijalab/seurat", ref = "release/4.0.0")

结果报错：

1.jpg

报错namespace 'sctransform' 0.2.1 is being loaded, but >= 0.3.1 is required 显示sctransform包有问题，需要更新最新版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

接下来开始一系列操作，载入必要的包、载入数据：

library(Seurat)
library(SeuratData)
library(cowplot)
library(dplyr)
InstallData("bmcite") # 99.7 MB
bm <- LoadData(ds = "bmcite")

2. 预处理和降维

首先分别对RNA和ADT这两种assays进行预处理和降维。这里使用标准归一化流程，但也可以使用SCTransform或任何替代方法。

DefaultAssay(bm) <- 'RNA'
bm <- NormalizeData(bm) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()

DefaultAssay(bm) <- 'ADT's
# 使用所有ADT 特征进行降维
# 这里设置了一个降维名，防止改写 
VariableFeatures(bm) <- rownames(bm[["ADT"]])
bm <- NormalizeData(bm, normalization.method = 'CLR', margin = 2) %>%   # CLR：中心对数比  margin：当执行CLR时，对特征（1）或细胞（2）进行标准化
      ScaleData() %>% 
      RunPCA(reduction.name = 'apca')

3. 计算多模态近邻

我们根据RNA和蛋白质相似度的加权组合来计算每个细胞在数据集中的最近邻。利用FindMultiModalNeighbors函数中计算细胞特异性模态权重和多模态近邻，在此数据集上运行约 2 min。我们需要指定每个模态的维数（类似于指定scRNA-seq集群中的PC数量)，并且可以改变这些设置，会发现小改动对总体结果的影响极小。

bm <- FindMultiModalNeighbors( bm
      , reduction.list = list("pca", "apca")
      , dims.list = list(1:30, 1:18)
      , modality.weight.name = "RNA.weight"
      )

识别的多模态近邻将存储于neighbors slot中，可通过 bm[['weighted.nn']]访问
通过bm[["wknn"]]访问WNN图；通过bm[[“wsnn”]]访问用于聚类的SNN图
通过bm$RNA.weight访问细胞特异性模态权重

4. 下游分析：可视化和聚类

4.1 多个模态数据聚类和可视化

我们现在可以使用这些结果进行下游分析，比如可视化和聚类。
例如，可以基于RNA和蛋白质数据的加权组合创建数据的UMAP可视化。
还可以执行基于图形的聚类，并在UMAP上可视化这些结果，以及细胞注释。
（类似于Seurat V3提出的共嵌入分析，将两个模态数据一同展示和聚类）

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

多模态数据UMAP图

P1按聚类结果显示，P2按细胞类型显示

4.2 单模态数据的可视化和比较

我们也可以对RNA和蛋白质数据分别UMAP可视化，并进行比较。
我们发现，RNA分析相比于ADT分析在识别祖细胞状态方面能够提供更多的信息，表现在
ADT面板包含分化细胞的markers，而T细胞状态则相反（ADT分析优于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基因和蛋白表达的可视化

我们还可以在多模态UMAP图上可视化一些典型marker gene和蛋白的表达，有助于验证所提供的注释：

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

The expression of canonical marker genes and proteins on the multimodal UMAP

eg：Naive CD4+ T cell, B cell—CD45RA；NK cell—CD16；NKT cell—TRDC；GMP cell（Granulocyte-monocyte progenitor）—MPO
以上我是通过细胞标志物数据库——CellMarker查询的。

4.4 模态权重的可视化

最后，我们可视化每个细胞的模态权重。RNA权重最高的群体代表祖细胞，而蛋白质权重最高的群体代表T细胞。这符合我们的生物学预期，因为抗体中不包含可以区分不同祖细胞群体的marker。

 VlnPlot(bm, features = "RNA.weight", group.by = 'celltype.l2', sort = TRUE, pt.size = 0.1) + NoLegend()

模态权重的可视化

从最初的CCA到CCA+Anchor到4.0版本的WNN，从多源的单细胞转录组数据整合到多个模态的单细胞数据整合，真的是十分敬佩Satija实验室，年年高产，惊喜不断。