hello，大家好，最近呢，10X單細(xì)胞核轉(zhuǎn)錄組 + 10X單細(xì)胞ATAC的多組學(xué)測(cè)序技術(shù)日漸成熟，因?yàn)檫@樣技術(shù)同時(shí)可以獲得一個(gè)細(xì)胞表達(dá)的轉(zhuǎn)錄組信息和ATAC的信息，所以很多學(xué)者和公司開(kāi)始關(guān)注和推廣這項(xiàng)技術(shù)，目前來(lái)看，前景非常廣闊，而且很早之前Seurat就推出了WNN（Weighted Nearest Neighbor Analysis）分析，用于分析從一個(gè)細(xì)胞得到的多組學(xué)數(shù)據(jù)，大家可以參考我之前分享的文章,Seurat4版本的WNN的運(yùn)行與原理與softmax。

注意一點(diǎn)，10X單細(xì)胞核轉(zhuǎn)錄組 + 10X單細(xì)胞ATAC技術(shù)和一個(gè)樣本分成兩份，分別測(cè)轉(zhuǎn)錄組信息和ATAC信息的聯(lián)合分析方法完全不一樣，WNN也不簡(jiǎn)單的就是一一對(duì)應(yīng)起來(lái)，而是相互輔助進(jìn)行數(shù)據(jù)分析，非常重要，一定要認(rèn)真明白其中的原理和方法，采用合適的分析策略。

示例代碼

同時(shí)測(cè)量多種模態(tài)，稱(chēng)為多模態(tài)分析，代表了單細(xì)胞基因組學(xué)的一個(gè)令人興奮的前沿領(lǐng)域，需要新的計(jì)算方法來(lái)定義基于多種數(shù)據(jù)類(lèi)型的細(xì)胞狀態(tài)。 每種模態(tài)的不同信息內(nèi)容，即使在同一數(shù)據(jù)集中的細(xì)胞之間，也代表了對(duì)多模態(tài)數(shù)據(jù)集的分析和整合的緊迫挑戰(zhàn)。 在 (Hao, Hao et al, Cell 2021) 中，我們引入了“加權(quán)最近鄰”(WNN) 分析，這是一種無(wú)監(jiān)督框架，用于學(xué)習(xí)每個(gè)單元格中每種數(shù)據(jù)類(lèi)型的相對(duì)效用，從而能夠?qū)Χ喾N模態(tài)進(jìn)行綜合分析 。（相對(duì)效用這個(gè)大家要重點(diǎn)關(guān)注一下）。

WNN analysis of 10x Multiome, RNA + ATAC

Here, we demonstrate the use of WNN analysis to a second multimodal technology, the 10x multiome RNA+ATAC kit. We use a dataset that is publicly available on the 10x website, where paired transcriptomes and ATAC-seq profiles are measured in 10,412 PBMCs.

• Create a multimodal Seurat object with paired transcriptome and ATAC-seq profiles
• Perform weighted neighbor clustering on RNA+ATAC data in single cells
• Leverage both modalities to identify putative regulators of different cell types and states

加載包

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

我們將根據(jù)基因表達(dá)數(shù)據(jù)創(chuàng)建一個(gè) Seurat 對(duì)象，然后添加 ATAC-seq 數(shù)據(jù)作為第二個(gè)檢測(cè)。

# 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

我們根據(jù)每種模式檢測(cè)到的分子數(shù)量以及線(xiàn)粒體百分比執(zhí)行基本 QC。

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

圖片.png

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

接下來(lái)，我們使用 RNA 和 ATAC-seq 數(shù)據(jù)的標(biāo)準(zhǔn)方法，獨(dú)立地對(duì)兩種檢測(cè)進(jìn)行預(yù)處理和降維。

# 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_")

We calculate a WNN graph, representing a weighted combination of RNA and ATAC-seq modalities. We use this graph for UMAP visualization and clustering

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)

We annotate the clusters below.

# 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)

我們可以可視化基于基因表達(dá)、ATAC-seq 或 WNN 分析的聚類(lèi)。 差異比之前的分析更微妙，但我們發(fā)現(xiàn) WNN 提供了最清晰的細(xì)胞狀態(tài)分離。

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))

圖片.png

例如，ATAC-seq 數(shù)據(jù)有助于分離 CD4 和 CD8 T 細(xì)胞狀態(tài)。 這是由于存在多個(gè)基因座，這些基因座在不同 T 細(xì)胞亞型之間表現(xiàn)出不同的可及性。

## 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)

圖片.png

Next, we will examine the accessible regions of each cell to determine enriched motifs. 有幾種方法可以做到這一點(diǎn)，但我們將使用 Greenleaf 實(shí)驗(yàn)室的 chromVAR 包。 這會(huì)計(jì)算已知基序的每個(gè)細(xì)胞可訪(fǎng)問(wèn)性分?jǐn)?shù)，并將這些分?jǐn)?shù)添加為 Seurat 對(duì)象中的第三個(gè)檢測(cè) (chromvar)。

加載數(shù)據(jù)

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
)

最后，我們探索多模態(tài)數(shù)據(jù)集以識(shí)別每個(gè)細(xì)胞狀態(tài)的關(guān)鍵調(diào)節(jié)器。 配對(duì)數(shù)據(jù)提供了一個(gè)獨(dú)特的機(jī)會(huì)來(lái)識(shí)別滿(mǎn)足多個(gè)標(biāo)準(zhǔn)的轉(zhuǎn)錄因子 (TF)，有助于將推定的調(diào)節(jié)器列表縮小到最有可能的候選者。 我們的目標(biāo)是識(shí)別在 RNA 測(cè)量中在多種細(xì)胞類(lèi)型中表達(dá)豐富的 TF，但在 ATAC 測(cè)量中也豐富了其基序的可及性。

作為示例和陽(yáng)性對(duì)照，CCAAT 增強(qiáng)子結(jié)合蛋白 (CEBP) 蛋白家族，包括 TF CEBPB，已反復(fù)證明在包括單核細(xì)胞和樹(shù)突細(xì)胞在內(nèi)的骨髓細(xì)胞的分化和功能中發(fā)揮重要作用。 我們可以看到 CEBPB 的表達(dá)和 MA0466.2.4 基序（編碼 CEBPB 的結(jié)合位點(diǎn)）的可及性都富含單核細(xì)胞。

#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

圖片.png

量化這種關(guān)系，并在所有細(xì)胞類(lèi)型中搜索以找到類(lèi)似的例子。 為此，我們將使用 presto 包來(lái)執(zhí)行快速差分表達(dá)。 我們進(jìn)行了兩項(xiàng)測(cè)試：一項(xiàng)使用基因表達(dá)數(shù)據(jù)，另一項(xiàng)使用 chromVAR 基序可訪(fǎng)問(wèn)性。 presto 根據(jù) Wilcox 秩和檢驗(yàn)計(jì)算 p 值，這也是 Seurat 中的默認(rèn)檢驗(yàn)，我們將搜索限制為在兩個(gè)檢驗(yàn)中都返回顯著結(jié)果的 TF。

presto 還計(jì)算“AUC”統(tǒng)計(jì)數(shù)據(jù)，它反映了每個(gè)基因（或基序）作為細(xì)胞類(lèi)型標(biāo)記的能力。 最大 AUC 值為 1 表示完美標(biāo)記。 由于基因和基序的 AUC 統(tǒng)計(jì)量在同一尺度上，我們?nèi)蓚€(gè)測(cè)試的 AUC 值的平均值，并使用它來(lái)對(duì)每種細(xì)胞類(lèi)型的 TF 進(jìn)行排序：

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)
}

We can now compute, and visualize, putative regulators for any cell type.

# 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

圖片.png

# identify top markers in pDC and visualize
##   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

圖片.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

圖片.png

不同的技術(shù)需要用到不同的方法，大家需要注意留心

生活很好，等你超越