20迫靖、第二十計(jì) 混水摸魚(yú)
比喻趁混亂時(shí)機(jī)攫取不正當(dāng)?shù)睦嫜咽病Ｒ沧鳌皽喫~(yú)”测秸。
此計(jì)是說(shuō)打仗時(shí)要得于抓住敵方的可乘之隙，借機(jī)行事灾常，使亂順我之意霎冯，我便亂中取利。
單細(xì)胞轉(zhuǎn)錄組學(xué)已經(jīng)改變了我們表征細(xì)胞狀態(tài)的能力钞瀑，但是對(duì)生物學(xué)的深入了解不僅需要分類學(xué)的分類列表沈撞。隨著測(cè)量不同細(xì)胞形態(tài)的新方法的出現(xiàn)，一個(gè)關(guān)鍵的分析挑戰(zhàn)是整合這些數(shù)據(jù)集以更好地理解細(xì)胞的特性和功能雕什。例如缠俺，用戶可以在同一生物系統(tǒng)上執(zhí)行scRNA-seq和scATAC-seq實(shí)驗(yàn)，并用相同的細(xì)胞類型標(biāo)簽集一致地注釋兩個(gè)數(shù)據(jù)集贷岸。由于scATAC-seq數(shù)據(jù)集難以注釋壹士，因?yàn)橐詥渭?xì)胞分辨率收集的基因組數(shù)據(jù)稀疏，并且在scRNA-seq數(shù)據(jù)中缺乏可解釋的基因標(biāo)記偿警，因此該分析尤其具有挑戰(zhàn)性躏救。
在Butler *等人的Stuart *中，2019年户敬，我們介紹了整合從同一生物系統(tǒng)收集的scRNA-seq和scATAC-seq數(shù)據(jù)集的方法落剪，并在此小插圖中演示了這些方法。特別是尿庐，我們演示了以下分析：
- 如何使用帶注釋的scRNA-seq數(shù)據(jù)集標(biāo)記來(lái)自scATAC-seq實(shí)驗(yàn)的細(xì)胞
- 如何從scRNA-seq和scATAC-seq共同可視化（共嵌入）細(xì)胞
- 如何將scATAC-seq細(xì)胞投影到源自scRNA-seq實(shí)驗(yàn)的UMAP
該插圖廣泛使用了Signac軟件包忠怖，該軟件包是最近開(kāi)發(fā)的，用于分析以單細(xì)胞分辨率（包括scATAC-seq）收集的染色質(zhì)數(shù)據(jù)集抄瑟。請(qǐng)參閱西涅克網(wǎng)站凡泣，了解更多的護(hù)身符和文檔分析scATAC-seq的數(shù)據(jù)。
我們使用來(lái)自10x Genomics的?12,000個(gè)人PBMC'multiome'數(shù)據(jù)集公開(kāi)展示了這些方法皮假。在此數(shù)據(jù)集中鞋拟，scRNA-seq和scATAC-seq圖譜同時(shí)收集在同一細(xì)胞中。出于本插圖的目的惹资，我們將數(shù)據(jù)集視為源自兩個(gè)不同的實(shí)驗(yàn)贺纲，并將其集成在一起。由于它們最初是在同一單元中測(cè)量的褪测，因此提供了一個(gè)基本事實(shí)猴誊，我們可以用來(lái)評(píng)估積分的準(zhǔn)確性潦刃。我們強(qiáng)調(diào)，此處我們將多基因組數(shù)據(jù)集用于演示和評(píng)估目的懈叹，并且用戶應(yīng)將這些方法應(yīng)用于分別收集的scRNA-seq和scATAC-seq數(shù)據(jù)集乖杠。我們提供單獨(dú)的加權(quán)最近鄰小插圖（WNN） 描述了多組單細(xì)胞數(shù)據(jù)的分析策略。
加載數(shù)據(jù)并分別處理每個(gè)模式
PBMC多基因組數(shù)據(jù)集可從10x基因組學(xué)獲得澄成。為了方便加載和瀏覽胧洒，它也可以作為我們的SeuratData軟件包的一部分提供。我們分別加載RNA和ATAC數(shù)據(jù)墨状，并假設(shè)這些配置文件是在單獨(dú)的實(shí)驗(yàn)中測(cè)量的卫漫。我們?cè)?a target="_blank">WNN小插圖中對(duì)這些單元進(jìn)行了注釋，并且注釋也包含在SeuratData中歉胶。
library(SeuratData)
# install the dataset and load requirements
InstallData("pbmcMultiome")
library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)
library(ggplot2)
library(cowplot)
# load both modalities
pbmc.rna <- LoadData("pbmcMultiome", "pbmc.rna")
pbmc.atac <- LoadData("pbmcMultiome", "pbmc.atac")
# repeat QC steps performed in the WNN vignette
pbmc.rna <- subset(pbmc.rna, seurat_annotations != "filtered")
pbmc.atac <- subset(pbmc.atac, seurat_annotations != "filtered")
# Perform standard analysis of each modality independently RNA analysis
pbmc.rna <- NormalizeData(pbmc.rna)
pbmc.rna <- FindVariableFeatures(pbmc.rna)
pbmc.rna <- ScaleData(pbmc.rna)
pbmc.rna <- RunPCA(pbmc.rna)
pbmc.rna <- RunUMAP(pbmc.rna, dims = 1:30)
# ATAC analysis add gene annotation information
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- "UCSC"
genome(annotations) <- "hg38"
Annotation(pbmc.atac) <- annotations
# We exclude the first dimension as this is typically correlated with sequencing depth
pbmc.atac <- RunTFIDF(pbmc.atac)
pbmc.atac <- FindTopFeatures(pbmc.atac, min.cutoff = "q0")
pbmc.atac <- RunSVD(pbmc.atac)
pbmc.atac <- RunUMAP(pbmc.atac, reduction = "lsi", dims = 2:30, reduction.name = "umap.atac", reduction.key = "atacUMAP_")
現(xiàn)在汛兜，我們繪制兩種模態(tài)的結(jié)果。先前已根據(jù)轉(zhuǎn)錄組狀態(tài)對(duì)細(xì)胞進(jìn)行了注釋通今。我們將預(yù)測(cè)scATAC-seq單元的注釋。
p1 <- DimPlot(pbmc.rna, group.by = "seurat_annotations", label = TRUE) + NoLegend() + ggtitle("RNA")
p2 <- DimPlot(pbmc.atac, group.by = "orig.ident", label = FALSE) + NoLegend() + ggtitle("ATAC")
p1 + p2
識(shí)別scRNA-seq和scATAC-seq數(shù)據(jù)集之間的錨點(diǎn)
為了識(shí)別scRNA-seq和scATAC-seq實(shí)驗(yàn)之間的“錨點(diǎn)”肛根，我們首先使用GeneActivity()功能通過(guò)對(duì)2kb上游區(qū)域和基因體中的ATAC-seq計(jì)數(shù)進(jìn)行量化辫塌，來(lái)生成每個(gè)基因的轉(zhuǎn)錄活性的粗略估計(jì)在Signac包中。然后將來(lái)自scATAC-seq數(shù)據(jù)的隨后的基因活性得分與來(lái)自scRNA-seq的基因表達(dá)定量一起用作規(guī)范相關(guān)分析的輸入派哲。我們對(duì)從scRNA-seq數(shù)據(jù)集確定為高度可變的所有基因進(jìn)行此量化臼氨。
# quantify gene activity
gene.activities <- GeneActivity(pbmc.atac, features = VariableFeatures(pbmc.rna))
# add gene activities as a new assay
pbmc.atac[["ACTIVITY"]] <- CreateAssayObject(counts = gene.activities)
# normalize gene activities
DefaultAssay(pbmc.atac) <- "ACTIVITY"
pbmc.atac <- NormalizeData(pbmc.atac)
pbmc.atac <- ScaleData(pbmc.atac, features = rownames(pbmc.atac))
# Identify anchors
transfer.anchors <- FindTransferAnchors(reference = pbmc.rna, query = pbmc.atac, features = VariableFeatures(object = pbmc.rna),
reference.assay = "RNA", query.assay = "ACTIVITY", reduction = "cca")
通過(guò)標(biāo)簽轉(zhuǎn)移注釋scATAC-seq細(xì)胞
識(shí)別錨點(diǎn)后，我們可以將注釋從scRNA-seq數(shù)據(jù)集轉(zhuǎn)移到scATAC-seq細(xì)胞上芭届。注釋存儲(chǔ)在seurat_annotations字段中储矩，并作為refdata參數(shù)的輸入提供。輸出將包含一個(gè)矩陣褂乍，其中包含每個(gè)ATAC序列信元的預(yù)測(cè)和置信度分?jǐn)?shù)持隧。
celltype.predictions <- TransferData(anchorset = transfer.anchors, refdata = pbmc.rna$seurat_annotations,
weight.reduction = pbmc.atac[["lsi"]], dims = 2:30)
pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)
執(zhí)行轉(zhuǎn)移后，ATAC-seq單元具有預(yù)測(cè)的注釋（從scRNA-seq數(shù)據(jù)集轉(zhuǎn)移）存儲(chǔ)在predicted.id字段中逃片。由于這些細(xì)胞是用多基因組試劑盒測(cè)量的屡拨，因此我們還有一個(gè)可以用于評(píng)估的真實(shí)注釋。您可以看到預(yù)測(cè)的注釋和實(shí)際的注釋非常相似褥实。
pbmc.atac$annotation_correct <- pbmc.atac$predicted.id == pbmc.atac$seurat_annotations
p1 <- DimPlot(pbmc.atac, group.by = "predicted.id", label = TRUE) + NoLegend() + ggtitle("Predicted annotation")
p2 <- DimPlot(pbmc.atac, group.by = "seurat_annotations", label = TRUE) + NoLegend() + ggtitle("Ground-truth annotation")
p1 | p2
在此示例中呀狼，約90％的時(shí)間通過(guò)scRNA-seq整合正確預(yù)測(cè)了scATAC-seq譜圖的注釋。另外损离，該prediction.score.max領(lǐng)域量化了與我們預(yù)測(cè)的注釋相關(guān)的不確定性哥艇。我們可以看到正確注釋的單元格通常與較高的預(yù)測(cè)分?jǐn)?shù)（> 90％）相關(guān)聯(lián)，而錯(cuò)誤注釋的單元格與較低的預(yù)測(cè)分?jǐn)?shù)（<50％）相關(guān)聯(lián)僻澎。錯(cuò)誤的分配也傾向于反映密切相關(guān)的細(xì)胞類型（即中級(jí)與原始B細(xì)胞）貌踏。
predictions <- table(pbmc.atac$seurat_annotations, pbmc.atac$predicted.id)
predictions <- predictions/rowSums(predictions) # normalize for number of cells in each cell type
predictions <- as.data.frame(predictions)
p1 <- ggplot(predictions, aes(Var1, Var2, fill = Freq)) + geom_tile() + scale_fill_gradient(name = "Fraction of cells",
low = "#ffffc8", high = "#7d0025") + xlab("Cell type annotation (RNA)") + ylab("Predicted cell type label (ATAC)") +
theme_cowplot() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
correct <- length(which(pbmc.atac$seurat_annotations == pbmc.atac$predicted.id))
incorrect <- length(which(pbmc.atac$seurat_annotations != pbmc.atac$predicted.id))
data <- FetchData(pbmc.atac, vars = c("prediction.score.max", "annotation_correct"))
p2 <- ggplot(data, aes(prediction.score.max, fill = annotation_correct, colour = annotation_correct)) +
geom_density(alpha = 0.5) + theme_cowplot() + scale_fill_discrete(name = "Annotation Correct",
labels = c(paste0("FALSE (n = ", incorrect, ")"), paste0("TRUE (n = ", correct, ")"))) + scale_color_discrete(name = "Annotation Correct",
labels = c(paste0("FALSE (n = ", incorrect, ")"), paste0("TRUE (n = ", correct, ")"))) + xlab("Prediction Score")
p1 + p2
共嵌入scRNA-seq和scATAC-seq數(shù)據(jù)集
除了跨模式轉(zhuǎn)移標(biāo)記外瓮增，還可以在同一圖中可視化scRNA-seq和scATAC-seq細(xì)胞。我們強(qiáng)調(diào)此步驟主要用于可視化哩俭，并且是可選步驟绷跑。通常，當(dāng)我們?cè)趕cRNA-seq和scATAC-seq數(shù)據(jù)集之間進(jìn)行整合分析時(shí)凡资，我們主要集中在如上所述的標(biāo)簽轉(zhuǎn)移上砸捏。我們?cè)谙旅嬲故玖擞糜诠睬度氲墓ぷ髁鞒蹋⒃俅螐?qiáng)調(diào)了這是出于演示目的隙赁，尤其是在這種特殊情況下垦藏，實(shí)際上在同一單元中測(cè)量了scRNA-seq圖譜和scATAC-seq圖譜。
為了執(zhí)行共嵌入伞访，我們首先根據(jù)先前計(jì)算的錨點(diǎn)將RNA表達(dá)“插入”到scATAC-seq細(xì)胞中掂骏，然后合并數(shù)據(jù)集。
# note that we restrict the imputation to variable genes from scRNA-seq, but could impute the
# full transcriptome if we wanted to
genes.use <- VariableFeatures(pbmc.rna)
refdata <- GetAssayData(pbmc.rna, assay = "RNA", slot = "data")[genes.use, ]
# refdata (input) contains a scRNA-seq expression matrix for the scRNA-seq cells. imputation
# (output) will contain an imputed scRNA-seq matrix for each of the ATAC cells
imputation <- TransferData(anchorset = transfer.anchors, refdata = refdata, weight.reduction = pbmc.atac[["lsi"]],
dims = 2:30)
pbmc.atac[["RNA"]] <- imputation
coembed <- merge(x = pbmc.rna, y = pbmc.atac)
# Finally, we run PCA and UMAP on this combined object, to visualize the co-embedding of both
# datasets
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)
coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)
DimPlot(coembed, group.by = c("orig.ident", "seurat_annotations"))
