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
image.png目前的单细胞组学已经不局限于在单细胞转录组学上进行相关的分析，在很多其他组学上，比如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软件的安装。
代码语言：javascript复制install.packages("Seurat")
library(Seurat)
library(SeuratData)
library(cowplot)
library(dplyr)
InstallData("bmcite")
bm <- LoadData(ds = "bmcite")我们首先对两个检测项目独立进行预处理和降维。我们使用标准归一化，但你也可以使用SCTransform或任何其他方法。
代码语言：javascript复制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数量），但你可以改变这些设置，看看小的变化对整体结果的影响最小。
代码语言：javascript复制# 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上进行基于图形的聚类，并将这些结果与一组细胞注释一起可视化。
代码语言：javascript复制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
wumap.plot-1.png我们可以在多模态UMAP上直观地看到典型标记基因和蛋白质的表达，这可以帮助验证所提供的注释。
代码语言：javascript复制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
ftplot-1.png最后，我们可以直观地看到每个细胞的模态权重。每个具有最高RNA权重的细胞群代表祖细胞，而具有最高蛋白质权重的细胞群代表T细胞。这符合我们的生物学预期，因为抗体不包含能区分不同祖细胞群的标记。
代码语言：javascript复制 VlnPlot(bm, features = "RNA.weight", group.by = 'celltype.l2', sort = TRUE, pt.size = 0.1)
NoLegend()
plotwts-1.pngWNN 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，确保安装了以下软件包。
代码语言：javascript复制library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)
library(dplyr)
library(ggplot2)我们将根据基因表达数据创建一个Seurat对象，然后加入ATAC-seq数据作为第二个ChromatinAssay矩阵。你可以浏览Signac入门流程分析，了解更多关于创建和处理ChromatinAssay对象的信息。
代码语言：javascript复制# 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我们根据每种模式检测到的数量以及线粒体百分比进行基本的质量控制。
代码语言：javascript复制VlnPlot(pbmc, features = c("nCount_ATAC", "nCount_RNA","percent.mt"), ncol = 3,
log = TRUE, pt.size = 0) NoLegend()
QCObject-1.pngpbmc <- subset(
x = pbmc,
subset = nCount_ATAC < 7e4 &
nCount_ATAC > 5e3 &
nCount_RNA < 25000 &
nCount_RNA > 1000 &
percent.mt < 20
)接下来，我们使用RNA和ATAC-seq数据的标准方法，对两个组学分别进行预处理和降维。
代码语言：javascript复制# 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的可视化和聚类
代码语言：javascript复制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。
代码语言：javascript复制# 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提供了最清晰的细胞状态的分离结果。
代码语言：javascript复制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))
UMAPs-1.png例如，ATAC-seq数据有助于分离CD4和CD8 T细胞状态。这是由于在不同的T细胞亚型之间存在着多个表现出不同可及性的基因座。例如，我们可以使用Signac可视化中的工具，将CD8A基因座的 "pseudobulk "轨迹与基因表达水平的小提琴图一起可视化。
代码语言：javascript复制## 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)
coverageplotcd8-1.png接下来，我们将检查每个细胞的可访问区域，以确定富集的图案。正如Signac流程中所描述的，有几种方法可以做到这一点，但我们将使用Greenleaf实验室的chromVAR软件包。
要继续，请确保你已经安装了以下软件包。
1.chromVAR用于分析scATAC-seq中的主题可及性
2.presto用于快速差异表达分析。
3.TFBSTools用于TFBS分析
4.JASPAR2020用于JASPAR动机模型
5.motifmatchr用于motif匹配
6.用于chromVAR的BSgenome.Hsapiens.UCSC.hg38
代码语言：javascript复制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的结合位点）的可及性都在富集到了单核细胞中。
代码语言：javascript复制#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
CEBPB-1.png我们想对这种关系进行量化，并在所有的细胞类型中搜索，找到类似的结果。为此，我们将使用presto包来进行快速差异分析。对两个测试集进行分析：一个使用基因表达数据，另一个使用chromVAR可及性数据。presto根据Wilcox秩和测试计算p值，这也是Seurat中的默认测试，我们将搜索在两个测试中结果显著的TFs。
presto还计算了一个 "AUC "统计量，它反映了每个基因（或图案）作为细胞类型标记的能力。最大AUC值为1表示一个完美的标记。由于AUC统计量对基因和图案都是在同一尺度上，我们从两个测试中取AUC值的平均值，并以此对每个细胞类型的TFs进行排序。
代码语言：javascript复制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。
代码语言：javascript复制# 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
NK-1.png# 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
pDC-1.png# 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
CD16Mono-1-1.png# 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总结
目前针对于同胞多组学的数据还是相对比较少的，在一个细胞的层面上对两种组学进行分析，减少了后续分析的假阳性的出现，也是研究发育生物学很好的一个研究手段之一。
