Seurat软件学习7-同胞多组学结合方法-WNN

原创

小胡子刺猬的生信学习123

发布于 2022-11-04 18:55:47

9460

发布于 2022-11-04 18:55:47

文章被收录于专栏：文献分享及代码学习文献分享及代码学习

Seurat软件学习1-多个模型得数据进行整合：https://cloud.tencent.com/developer/article/2130078

Seurat软件学习2-scrna数据整合分析：https://cloud.tencent.com/developer/article/2131431

Seurat软件学习3-scrna数据整合分析注释数据集：https://cloud.tencent.com/developer/article/2133583

Seurat软件学习4-使用RPCA进行快速整合数据集：https://cloud.tencent.com/developer/article/2134684

Seurat软件学习5-scRNA-Seq和scATAC-Seq数据整合：https://cloud.tencent.com/developer/article/2136814

Seurat软件学习6-多模型参考映射的方法：https://cloud.tencent.com/developer/article/2144475

目前的单细胞组学已经不局限于在单细胞转录组学上进行相关的分析，在很多其他组学上，比如DNA层面及蛋白层面也进行了相关的发展。

因此开发seurat的作者也针对于同胞多组学的研究内容进行了算法的开发，主要是基于seurat的wnn的分析。

同时测量多个组学的数据，即所谓的多组学分析，代表了单细胞基因组学的一个令人兴奋的前沿方向，因此需要新的计算方法，可以基于多种数据类型定义细胞状态。每种模式的信息含量各不相同，甚至在同一数据集中的不同单元之间也是如此，这对多模式数据集的分析和整合提出了紧迫的挑战。在(Hao，Hao等人，Cell 2021)中，我们引入了‘加权最近邻’(WNN)分析，这是一个非监督框架，用于了解每个单元中每种数据类型的相对效用，从而实现对多个模态的综合分析。

本小节介绍了用于分析多模态单细胞数据集的WNN工作流程。该工作流程包括三个步骤：

对每种模式分别进行独立的预处理和降维处理
学习细胞特定的模式 "权重"，并构建一个整合各种模式的WNN图形
WNN图的下游分析（即可视化、聚类等）。我们展示了WNN分析在两种单细胞多模态技术中的应用。CITE-seq和10x multiome。我们在两种模式的基础上定义细胞状态，而不是单独的一种模式。

WNN analysis of CITE-seq, RNA + ADT

我们使用来自（Stuart, Butler et al, Cell 2019）的CITE-seq数据集，该数据集由30,672个scRNA-seq配置文件组成，与来自骨髓的25种抗体一起测量。该对象包含两种检测方法，RNA和抗体衍生标签（ADT）。

首先进行seurat软件的安装。

install.packages("Seurat")
library(Seurat)
library(SeuratData)
library(cowplot)
library(dplyr)
InstallData("bmcite")
bm <- LoadData(ds = "bmcite")

我们首先对两个检测项目独立进行预处理和降维。我们使用标准归一化，但你也可以使用SCTransform或任何其他方法。

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

DefaultAssay(bm) <- 'ADT'
# we will use all ADT features for dimensional reduction
# we set a dimensional reduction name to avoid overwriting the 
VariableFeatures(bm) <- rownames(bm[["ADT"]])
bm <- NormalizeData(bm, normalization.method = 'CLR', margin = 2) %>% 
  ScaleData() %>% RunPCA(reduction.name = 'apca')

对于每个细胞，我们根据RNA和蛋白质相似性的加权组合来计算其在数据集中的最接近的邻居。细胞特定的模态权重和多模态邻居是在一个函数中计算的，在这个数据集上运行需要2分钟左右。我们指定每个模态的维度（类似于指定scRNA-seq聚类中包含的PC数量），但你可以改变这些设置，看看小的变化对整体结果的影响最小。

# Identify multimodal neighbors. These will be stored in the neighbors slot, 
# and can be accessed using bm[['weighted.nn']]
# The WNN graph can be accessed at bm[["wknn"]], 
# and the SNN graph used for clustering at bm[["wsnn"]]
# Cell-specific modality weights can be accessed at bm$RNA.weight
bm <- FindMultiModalNeighbors(
  bm, reduction.list = list("pca", "apca"), 
  dims.list = list(1:30, 1:18), modality.weight.name = "RNA.weight"
)

我们现在可以使用这些结果进行下游分析，如可视化和聚类。例如，我们可以根据RNA和蛋白质数据的加权组合，创建一个UMAP的可视化数据。我们还可以在UMAP上进行基于图形的聚类，并将这些结果与一组细胞注释一起可视化。

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上直观地看到典型标记基因和蛋白质的表达，这可以帮助验证所提供的注释。

p5 <- FeaturePlot(bm, features = c("adt_CD45RA","adt_CD16","adt_CD161"),
                  reduction = 'wnn.umap', max.cutoff = 2, 
                  cols = c("lightgrey","darkgreen"), ncol = 3)
p6 <- FeaturePlot(bm, features = c("rna_TRDC","rna_MPO","rna_AVP"), 
                  reduction = 'wnn.umap', max.cutoff = 3, ncol = 3)
p5 / p6

最后，我们可以直观地看到每个细胞的模态权重。每个具有最高RNA权重的细胞群代表祖细胞，而具有最高蛋白质权重的细胞群代表T细胞。这符合我们的生物学预期，因为抗体不包含能区分不同祖细胞群的标记。

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

WNN analysis of 10x Multiome, RNA + ATAC

在此，我们展示了WNN分析在第二种多模态技术中的应用，即10x多组学的RNA+ATAC试剂盒。我们使用了10x网站上公开的数据集，其中的转录组和ATAC-seq图谱是在10,412个PBMCs中测量的。

我们使用的WNN方法与我们在前一个标签中使用的方法相同，我们将综合多模态分析应用于CITE-seq数据集。在这个例子中，我们将演示如何进行下面的分析内容：

1.用成对的转录组和ATAC-seq图谱创建一个多模式的Seurat对象

2.对单细胞中的RNA+ATAC数据进行加权邻近聚类分析

3.利用这两种方式来确定不同细胞类型和状态的调控模式

你可以从10x基因组学的网站上下载数据集，请确保下载以下文件。

1.经过过滤的特征条码矩阵（HDF5）。

2.ATAC每个片段信息文件(TSV.GZ)

3.ATAC每个片段信息索引(TSV.GZ索引)

最后，为了运行vignette，确保安装了以下软件包。

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

我们将根据基因表达数据创建一个Seurat对象，然后加入ATAC-seq数据作为第二个ChromatinAssay矩阵。你可以浏览Signac入门流程分析，了解更多关于创建和处理ChromatinAssay对象的信息。

# the 10x hdf5 file contains both data types. 
inputdata.10x <- Read10X_h5("../data/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5")

# extract RNA and ATAC data
rna_counts <- inputdata.10x$`Gene Expression`
atac_counts <- inputdata.10x$Peaks

# Create Seurat object
pbmc <- CreateSeuratObject(counts = rna_counts)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# Now add in the ATAC-seq data
# we'll only use peaks in standard chromosomes
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)
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg38"

frag.file <- "../data/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz"
chrom_assay <- CreateChromatinAssay(
   counts = atac_counts,
   sep = c(":", "-"),
   genome = 'hg38',
   fragments = frag.file,
   min.cells = 10,
   annotation = annotations
 )
pbmc[["ATAC"]] <- chrom_assay

我们根据每种模式检测到的数量以及线粒体百分比进行基本的质量控制。

VlnPlot(pbmc, features = c("nCount_ATAC", "nCount_RNA","percent.mt"), ncol = 3,
  log = TRUE, pt.size = 0) + NoLegend()

pbmc <- subset(
  x = pbmc,
  subset = nCount_ATAC < 7e4 &
    nCount_ATAC > 5e3 &
    nCount_RNA < 25000 &
    nCount_RNA > 1000 &
    percent.mt < 20
)

接下来，我们使用RNA和ATAC-seq数据的标准方法，对两个组学分别进行预处理和降维。

# RNA analysis
DefaultAssay(pbmc) <- "RNA"
pbmc <- SCTransform(pbmc, verbose = FALSE) %>% RunPCA() %>% RunUMAP(dims = 1:50, reduction.name = 'umap.rna', reduction.key = 'rnaUMAP_')

# ATAC analysis
# We exclude the first dimension as this is typically correlated with sequencing depth
DefaultAssay(pbmc) <- "ATAC"
pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
pbmc <- RunSVD(pbmc)
pbmc <- RunUMAP(pbmc, reduction = 'lsi', dims = 2:50, reduction.name = "umap.atac", reduction.key = "atacUMAP_")

我们计算出一个WNN图，代表RNA和ATAC-seq模式的加权组合。我们使用该图进行UMAP的可视化和聚类

pbmc <- FindMultiModalNeighbors(pbmc, reduction.list = list("pca", "lsi"), dims.list = list(1:50, 2:50))
pbmc <- RunUMAP(pbmc, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
pbmc <- FindClusters(pbmc, graph.name = "wsnn", algorithm = 3, verbose = FALSE)

我们在下面注释了这些簇。请注意，你也可以使用我们的监督映射管道来注释数据集，使用我们的vignette，或自动网络工具Azimuth。

# perform sub-clustering on cluster 6 to find additional structure
pbmc <- FindSubCluster(pbmc, cluster = 6, graph.name = "wsnn", algorithm = 3)
Idents(pbmc) <- "sub.cluster"
# add annotations
pbmc <- RenameIdents(pbmc, '19' = 'pDC','20' = 'HSPC','15' = 'cDC')
pbmc <- RenameIdents(pbmc, '0' = 'CD14 Mono', '9' ='CD14 Mono', '5' = 'CD16 Mono')
pbmc <- RenameIdents(pbmc, '10' = 'Naive B', '11' = 'Intermediate B', '17' = 'Memory B', '21' = 'Plasma')
pbmc <- RenameIdents(pbmc, '7' = 'NK')
pbmc <- RenameIdents(pbmc, '4' = 'CD4 TCM', '13'= "CD4 TEM", '3' = "CD4 TCM", '16' ="Treg", '1' ="CD4 Naive", '14' = "CD4 Naive")
pbmc <- RenameIdents(pbmc, '2' = 'CD8 Naive', '8'= "CD8 Naive", '12' = 'CD8 TEM_1', '6_0' = 'CD8 TEM_2', '6_1' ='CD8 TEM_2', '6_4' ='CD8 TEM_2')
pbmc <- RenameIdents(pbmc, '18' = 'MAIT')
pbmc <- RenameIdents(pbmc, '6_2' ='gdT', '6_3' = 'gdT')
pbmc$celltype <- Idents(pbmc)

我们可以根据基因表达、ATAC-seq或WNN分析来实现聚类的可视化。差异比之前的分析更微妙（你可以探索权重，它比我们的CITE-seq例子更均匀），但我们发现WNN提供了最清晰的细胞状态的分离结果。

p1 <- DimPlot(pbmc, reduction = "umap.rna", group.by = "celltype", label = TRUE, label.size = 2.5, repel = TRUE) + ggtitle("RNA")
p2 <- DimPlot(pbmc, reduction = "umap.atac", group.by = "celltype", label = TRUE, label.size = 2.5, repel = TRUE) + ggtitle("ATAC")
p3 <- DimPlot(pbmc, reduction = "wnn.umap", group.by = "celltype", label = TRUE, label.size = 2.5, repel = TRUE) + ggtitle("WNN")
p1 + p2 + p3 & NoLegend() & theme(plot.title = element_text(hjust = 0.5))

例如，ATAC-seq数据有助于分离CD4和CD8 T细胞状态。这是由于在不同的T细胞亚型之间存在着多个表现出不同可及性的基因座。例如，我们可以使用Signac可视化中的工具，将CD8A基因座的 "pseudobulk "轨迹与基因表达水平的小提琴图一起可视化。

## to make the visualization easier, subset T cell clusters
celltype.names <- levels(pbmc)
tcell.names <- grep("CD4|CD8|Treg", celltype.names,value = TRUE)
tcells <- subset(pbmc, idents = tcell.names)
CoveragePlot(tcells, region = 'CD8A', features = 'CD8A', assay = 'ATAC', expression.assay = 'SCT', peaks = FALSE)

接下来，我们将检查每个细胞的可访问区域，以确定富集的图案。正如Signac流程中所描述的，有几种方法可以做到这一点，但我们将使用Greenleaf实验室的chromVAR软件包。

要继续，请确保你已经安装了以下软件包。

1.chromVAR用于分析scATAC-seq中的主题可及性

2.presto用于快速差异表达分析。

3.TFBSTools用于TFBS分析

4.JASPAR2020用于JASPAR动机模型

5.motifmatchr用于motif匹配

6.用于chromVAR的BSgenome.Hsapiens.UCSC.hg38

library(chromVAR)
library(JASPAR2020)
library(TFBSTools)
library(motifmatchr)
library(BSgenome.Hsapiens.UCSC.hg38)

# Scan the DNA sequence of each peak for the presence of each motif, and create a Motif object
DefaultAssay(pbmc) <- "ATAC"
pwm_set <- getMatrixSet(x = JASPAR2020, opts = list(species = 9606, all_versions = FALSE))
motif.matrix <- CreateMotifMatrix(features = granges(pbmc), pwm = pwm_set, genome = 'hg38', use.counts = FALSE)
motif.object <- CreateMotifObject(data = motif.matrix, pwm = pwm_set)
pbmc <- SetAssayData(pbmc, assay = 'ATAC', slot = 'motifs', new.data = motif.object)

# Note that this step can take 30-60 minutes 
pbmc <- RunChromVAR(
  object = pbmc,
  genome = BSgenome.Hsapiens.UCSC.hg38
)

最后，我们探索多组学数据集，以确定每个细胞状态的关键调节因子。成对的数据提供了一个独特的机会来识别满足多种标准的转录因子（TFs），帮助缩小推测调节因子的范围，进行候选调节因子的确定。我们的目标是识别在RNA测量中其表达在多种细胞类型中富集的TFs，同时在ATAC测量中也有富集的可及性。

作为一个例子和阳性对照，CCAAT增强子结合蛋白（CEBP）家族的蛋白质，包括TF CEBPB，已被反复证明在髓系细胞（包括单核细胞和树突状细胞）的分化和功能中发挥重要作用。我们可以看到，CEBPB的表达和MA0466.2.4可视化（编码CEBPB的结合位点）的可及性都在富集到了单核细胞中。

#returns MA0466.2
motif.name <- ConvertMotifID(pbmc, name = 'CEBPB')
gene_plot <- FeaturePlot(pbmc, features = "sct_CEBPB", reduction = 'wnn.umap')
motif_plot <- FeaturePlot(pbmc, features = motif.name, min.cutoff = 0, cols = c("lightgrey", "darkred"), reduction = 'wnn.umap')
gene_plot | motif_plot

我们想对这种关系进行量化，并在所有的细胞类型中搜索，找到类似的结果。为此，我们将使用presto包来进行快速差异分析。对两个测试集进行分析：一个使用基因表达数据，另一个使用chromVAR可及性数据。presto根据Wilcox秩和测试计算p值，这也是Seurat中的默认测试，我们将搜索在两个测试中结果显著的TFs。

presto还计算了一个 "AUC "统计量，它反映了每个基因（或图案）作为细胞类型标记的能力。最大AUC值为1表示一个完美的标记。由于AUC统计量对基因和图案都是在同一尺度上，我们从两个测试中取AUC值的平均值，并以此对每个细胞类型的TFs进行排序。

markers_rna <- presto:::wilcoxauc.Seurat(X = pbmc, group_by = 'celltype', assay = 'data', seurat_assay = 'SCT')
markers_motifs <- presto:::wilcoxauc.Seurat(X = pbmc, group_by = 'celltype', assay = 'data', seurat_assay = 'chromvar')
motif.names <- markers_motifs$feature
colnames(markers_rna) <- paste0("RNA.", colnames(markers_rna))
colnames(markers_motifs) <- paste0("motif.", colnames(markers_motifs))
markers_rna$gene <- markers_rna$RNA.feature
markers_motifs$gene <- ConvertMotifID(pbmc, id = motif.names)
# a simple function to implement the procedure above
topTFs <- function(celltype, padj.cutoff = 1e-2) {
  ctmarkers_rna <- dplyr::filter(
    markers_rna, RNA.group == celltype, RNA.padj < padj.cutoff, RNA.logFC > 0) %>% 
    arrange(-RNA.auc)
  ctmarkers_motif <- dplyr::filter(
    markers_motifs, motif.group == celltype, motif.padj < padj.cutoff, motif.logFC > 0) %>% 
    arrange(-motif.auc)
  top_tfs <- inner_join(
    x = ctmarkers_rna[, c(2, 11, 6, 7)], 
    y = ctmarkers_motif[, c(2, 1, 11, 6, 7)], by = "gene"
  )
  top_tfs$avg_auc <- (top_tfs$RNA.auc + top_tfs$motif.auc) / 2
  top_tfs <- arrange(top_tfs, -avg_auc)
  return(top_tfs)
}

我们现在可以计算出任何细胞类型的假设TF，并使之可视化。我们复现了成熟的regulators，包括NK细胞的TBX21，浆细胞的IRF4，造血祖细胞的SOX4，B细胞的EBF1和PAX5，pDC的IRF8和TCF4。我们相信，类似的方法可以用来获得不同条件下的一系列的putative regulators。

# identify top markers in NK and visualize
head(topTFs("NK"), 3)
##   RNA.group  gene   RNA.auc      RNA.pval motif.group motif.feature motif.auc
## 1        NK TBX21 0.7264660  0.000000e+00          NK      MA0690.1 0.9223858
## 2        NK EOMES 0.5895889 1.552097e-100          NK      MA0800.1 0.9263286
## 3        NK RUNX3 0.7722418 7.131401e-122          NK      MA0684.2 0.7083570
##      motif.pval   avg_auc
## 1 2.343664e-211 0.8244259
## 2 2.786462e-215 0.7579587
## 3  7.069675e-53 0.7402994
motif.name <- ConvertMotifID(pbmc, name = 'TBX21')
gene_plot <- FeaturePlot(pbmc, features = "sct_TBX21", reduction = 'wnn.umap')
motif_plot <- FeaturePlot(pbmc, features = motif.name, min.cutoff = 0, cols = c("lightgrey", "darkred"), reduction = 'wnn.umap')
gene_plot | motif_plot

# identify top markers in pDC and visualize
head(topTFs("pDC"), 3)
##   RNA.group gene   RNA.auc      RNA.pval motif.group motif.feature motif.auc
## 1       pDC TCF4 0.9998833 3.347777e-163         pDC      MA0830.2 0.9959622
## 2       pDC IRF8 0.9905372 2.063258e-124         pDC      MA0652.1 0.8814713
## 3       pDC SPIB 0.9114648  0.000000e+00         pDC      MA0081.2 0.8997955
##     motif.pval   avg_auc
## 1 2.578226e-69 0.9979228
## 2 9.702602e-42 0.9360043
## 3 1.130040e-45 0.9056302
motif.name <- ConvertMotifID(pbmc, name = 'TCF4')
gene_plot <- FeaturePlot(pbmc, features = "sct_TCF4", reduction = 'wnn.umap')
motif_plot <- FeaturePlot(pbmc, features = motif.name, min.cutoff = 0, cols = c("lightgrey", "darkred"), reduction = 'wnn.umap')
gene_plot | motif_plot

# identify top markers in HSPC and visualize
head(topTFs("CD16 Mono"),3)
##   RNA.group  gene   RNA.auc      RNA.pval motif.group motif.feature motif.auc
## 1 CD16 Mono  SPI1 0.8764099 4.116679e-291   CD16 Mono      MA0080.5 0.8831213
## 2 CD16 Mono CEBPB 0.8675114 8.321489e-292   CD16 Mono      MA0466.2 0.7859496
## 3 CD16 Mono MEF2C 0.7132221  4.210640e-79   CD16 Mono      MA0497.1 0.7986104
##      motif.pval   avg_auc
## 1 3.965856e-188 0.8797656
## 2 1.092808e-105 0.8267305
## 3 4.462855e-115 0.7559162
motif.name <- ConvertMotifID(pbmc, name = 'SPI1')
gene_plot <- FeaturePlot(pbmc, features = "sct_SPI1", reduction = 'wnn.umap')
motif_plot <- FeaturePlot(pbmc, features = motif.name, min.cutoff = 0, cols = c("lightgrey", "darkred"), reduction = 'wnn.umap')
gene_plot | motif_plot

# identify top markers in other cell types
head(topTFs("Naive B"), 3)
##   RNA.group   gene   RNA.auc      RNA.pval motif.group motif.feature motif.auc
## 1   Naive B  MEF2C 0.8814368 1.002848e-183     Naive B      MA0497.1 0.6511956
## 2   Naive B   PAX5 0.9437146  0.000000e+00     Naive B      MA0014.3 0.5634996
## 3   Naive B POU2F2 0.6022295  4.067907e-13     Naive B      MA0507.1 0.8773954
##      motif.pval   avg_auc
## 1  3.903712e-23 0.7663162
## 2  3.175138e-05 0.7536071
## 3 5.457378e-135 0.7398125
head(topTFs("HSPC"), 3)
##   RNA.group  gene   RNA.auc     RNA.pval motif.group motif.feature motif.auc
## 1      HSPC  SOX4 0.9869221 7.766427e-71        HSPC      MA0867.2 0.7104016
## 2      HSPC GATA2 0.7884615 0.000000e+00        HSPC      MA0036.3 0.8308485
## 3      HSPC MEIS1 0.8830582 0.000000e+00        HSPC      MA0498.2 0.6781207
##     motif.pval   avg_auc
## 1 2.059705e-04 0.8486618
## 2 5.336146e-09 0.8096550
## 3 1.677267e-03 0.7805894

head(topTFs("Plasma"), 3)

##   RNA.group  gene   RNA.auc     RNA.pval motif.group motif.feature motif.auc
## 1    Plasma  IRF4 0.8189901 5.206829e-35      Plasma      MA1419.1 0.9782567
## 2    Plasma MEF2C 0.9070644 4.738760e-12      Plasma      MA0497.1 0.7687288
## 3    Plasma  TCF4 0.8052937 5.956053e-12      Plasma      MA0830.2 0.8046950
##     motif.pval   avg_auc
## 1 2.180043e-12 0.8986234
## 2 7.951876e-05 0.8378966
## 3 7.678507e-06 0.8049943