在本教程中，我們將學(xué)習(xí)使用Seurat3對(duì)scATAC-seq和scRNA-seq的數(shù)據(jù)進(jìn)行整合分析格嗅，使用一種新的數(shù)據(jù)轉(zhuǎn)移映射方法，將scATAC-seq的數(shù)據(jù)基于scRNA-seq數(shù)據(jù)聚類的結(jié)果進(jìn)行細(xì)胞分群汉匙，并進(jìn)行整合分析涤姊。
對(duì)于scATAC-seq，我們還開(kāi)發(fā)了一個(gè)獨(dú)立的Signac包霜定，可以用于單細(xì)胞ATAC-seq的數(shù)據(jù)分析及其整合分析档悠，詳細(xì)的使用方法可以參考Signac包官網(wǎng)的說(shuō)明文檔。
本示例中望浩，我們將使用10X genomics測(cè)序平臺(tái)產(chǎn)生的PBMC數(shù)據(jù)集（包含大約10K個(gè)細(xì)胞的scRNA-seq和scATAC-seq的測(cè)序數(shù)據(jù)）進(jìn)行分析演示辖所，整體的分析流程主要包括以下幾個(gè)步驟：
- Estimate RNA-seq levels from ATAC-seq (quantify gene expression ‘a(chǎn)ctivity’ from ATAC-seq reads)
- Learn the internal structure of the ATAC-seq data on its own (accomplished using LSI)
- Identify ‘a(chǎn)nchors’ between the ATAC-seq and RNA-seq datasets
- Transfer data between datasets (either transfer labels for classification, or impute RNA levels in the ATAC-seq data to enable co-embedding)
基因表達(dá)活性的定量
首先，我們加載scATAC-seq數(shù)據(jù)產(chǎn)生的一個(gè)peak matrix磨德，并將其整合合并成一個(gè)gene activity matrix缘回。我們基于這樣一個(gè)簡(jiǎn)單的假設(shè)：基因的表達(dá)活性可以簡(jiǎn)單的通過(guò)基因上下游2kb范圍內(nèi)覆蓋的reads數(shù)的加和進(jìn)行定量，最后獲得一個(gè)gene x cell的表達(dá)矩陣典挑。
# 加載所需的R包
library(Seurat)
library(ggplot2)
library(patchwork)
# 加載peak矩陣文件
peaks <- Read10X_h5("../data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
# create a gene activity matrix from the peak matrix and GTF, using chromosomes 1:22, X, and Y.
# Peaks that fall within gene bodies, or 2kb upstream of a gene, are considered
# 使用CreateGeneActivityMatrix函數(shù)構(gòu)建基因表達(dá)活性矩陣
activity.matrix <- CreateGeneActivityMatrix(peak.matrix = peaks, annotation.file = "../data/Homo_sapiens.GRCh37.82.gtf", seq.levels = c(1:22, "X", "Y"), upstream = 2000, verbose = TRUE)
構(gòu)建Seurat對(duì)象
接下來(lái)酥宴，我們創(chuàng)建一個(gè)seurat對(duì)象存儲(chǔ)scATAC-seq的數(shù)據(jù)，其中原始的peak counts存儲(chǔ)在"ATAC"這個(gè)Assay中您觉，基因表達(dá)活性的矩陣存儲(chǔ)在"RNA"這個(gè)Assay中拙寡。對(duì)于原始的count，我們將進(jìn)行QC過(guò)濾掉那些總counts數(shù)低于5k的細(xì)胞琳水。
# 創(chuàng)建seurat對(duì)象
pbmc.atac <- CreateSeuratObject(counts = peaks, assay = "ATAC", project = "10x_ATAC")
pbmc.atac[["ACTIVITY"]] <- CreateAssayObject(counts = activity.matrix)
# meta信息
meta <- read.table("../data/atac_v1_pbmc_10k_singlecell.csv", sep = ",", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
meta <- meta[colnames(pbmc.atac), ]
pbmc.atac <- AddMetaData(pbmc.atac, metadata = meta)
# 數(shù)據(jù)過(guò)濾
pbmc.atac <- subset(pbmc.atac, subset = nCount_ATAC > 5000)
pbmc.atac$tech <- "atac"
數(shù)據(jù)預(yù)處理
# 對(duì)基因表達(dá)活性矩陣進(jìn)行標(biāo)準(zhǔn)化
DefaultAssay(pbmc.atac) <- "ACTIVITY"
pbmc.atac <- FindVariableFeatures(pbmc.atac)
pbmc.atac <- NormalizeData(pbmc.atac)
pbmc.atac <- ScaleData(pbmc.atac)
Here we perform latent semantic indexing to reduce the dimensionality of the scATAC-seq data. This procedure learns an ‘internal’ structure for the scRNA-seq data, and is important when determining the appropriate weights for the anchors when transferring information. We utilize Latent Semantic Indexing (LSI) to learn the structure of ATAC-seq data, as proposed in Cusanovich et al, Science 2015 and implemented in the RunLSI function. LSI is implemented here by performing computing the term frequency-inverse document frequency (TF-IDF) followed by SVD.
# 對(duì)原始的peak矩陣進(jìn)行降維處理
DefaultAssay(pbmc.atac) <- "ATAC"
VariableFeatures(pbmc.atac) <- names(which(Matrix::rowSums(pbmc.atac) > 100))
# 使用LSI方法進(jìn)行降維
pbmc.atac <- RunLSI(pbmc.atac, n = 50, scale.max = NULL)
pbmc.atac <- RunUMAP(pbmc.atac, reduction = "lsi", dims = 1:50)
# 加載scRNA-seq數(shù)據(jù)
pbmc.rna <- readRDS("../data/pbmc_10k_v3.rds")
pbmc.rna$tech <- "rna"
# 數(shù)據(jù)可視化
p1 <- DimPlot(pbmc.atac, reduction = "umap") + NoLegend() + ggtitle("scATAC-seq")
p2 <- DimPlot(pbmc.rna, group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend() + ggtitle("scRNA-seq")
p1 + p2
接下來(lái)肆糕，我們將對(duì)scATAC-seq和scRNA-seq的數(shù)據(jù)進(jìn)行整合分析般堆。
# 識(shí)別整合的anchors，使用scRNA-seq的數(shù)據(jù)作為參照
transfer.anchors <- FindTransferAnchors(reference = pbmc.rna, query = pbmc.atac, features = VariableFeatures(object = pbmc.rna), reference.assay = "RNA", query.assay = "ACTIVITY", reduction = "cca")
# 使用TransferData函數(shù)進(jìn)行數(shù)據(jù)轉(zhuǎn)移映射
celltype.predictions <- TransferData(anchorset = transfer.anchors, refdata = pbmc.rna$celltype, weight.reduction = pbmc.atac[["lsi"]])
# 添加預(yù)測(cè)信息
pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)
然后擎宝，我們可以檢查預(yù)測(cè)分?jǐn)?shù)的分布郁妈，并可以選擇過(guò)濾掉那些得分較低的細(xì)胞。在這里绍申，我們發(fā)現(xiàn)超過(guò)95％的細(xì)胞得分為0.5或更高噩咪。
hist(pbmc.atac$prediction.score.max)
abline(v = 0.5, col = "red")
table(pbmc.atac$prediction.score.max > 0.5)
##
## FALSE TRUE
## 290 7576
同樣的，我們可以通過(guò)UMAP降維可視化scATAC-seq數(shù)據(jù)預(yù)測(cè)的細(xì)胞類型极阅，并發(fā)現(xiàn)所轉(zhuǎn)移的標(biāo)簽與scRNA-seq數(shù)據(jù)UMAP結(jié)果的結(jié)構(gòu)高度一致胃碾。
# 數(shù)據(jù)過(guò)濾
pbmc.atac.filtered <- subset(pbmc.atac, subset = prediction.score.max > 0.5)
pbmc.atac.filtered$predicted.id <- factor(pbmc.atac.filtered$predicted.id, levels = levels(pbmc.rna)) # to make the colors match
# 數(shù)據(jù)可視化
p1 <- DimPlot(pbmc.atac.filtered, group.by = "predicted.id", label = TRUE, repel = TRUE) + ggtitle("scATAC-seq cells") + NoLegend() + scale_colour_hue(drop = FALSE)
p2 <- DimPlot(pbmc.rna, group.by = "celltype", label = TRUE, repel = TRUE) + ggtitle("scRNA-seq cells") + NoLegend()
p1 + p2
共嵌入聯(lián)合可視化分析
最后，為了將所有的細(xì)胞一起可視化筋搏，我們可以將scRNA-seq和scATAC-seq細(xì)胞共同嵌入到同一低維空間中仆百。在這里，我們首先使用先前識(shí)別出的相同anchors來(lái)轉(zhuǎn)移映射細(xì)胞類型標(biāo)簽奔脐，以估算填充scATAC-seq細(xì)胞的RNA-seq值俄周。然后，我們合并測(cè)得的和估算填充的scRNA-seq數(shù)據(jù)髓迎，并運(yùn)行標(biāo)準(zhǔn)UMAP進(jìn)行降維可視化所有的細(xì)胞峦朗。請(qǐng)注意，此步驟僅用于可視化目的排龄，不是數(shù)據(jù)轉(zhuǎn)移映射分析的必要部分波势。
# 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"]])
# this line adds the imputed data matrix to the pbmc.atac object
pbmc.atac[["RNA"]] <- imputation
# 進(jìn)行merge合并
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)
coembed$celltype <- ifelse(!is.na(coembed$celltype), coembed$celltype, coembed$predicted.id)
Here we plot all cells colored by either their assigned cell type (from the 10K dataset) or their predicted cell type from the data transfer procedure.
p1 <- DimPlot(coembed, group.by = "tech")
p2 <- DimPlot(coembed, group.by = "celltype", label = TRUE, repel = TRUE)
p1 + p2
How can I further investigate unmatched groups of cells that only appear in one assay?
在查看UMAP嵌和數(shù)據(jù)可視化的結(jié)果中，我們發(fā)現(xiàn)有幾種細(xì)胞似乎僅在一種測(cè)序類型數(shù)據(jù)中存在橄维。例如尺铣，血小板細(xì)胞僅出現(xiàn)在scRNA-seq數(shù)據(jù)中。這些細(xì)胞被認(rèn)為正在經(jīng)歷從巨核細(xì)胞到血小板細(xì)胞分化的過(guò)程中争舞，因此要么完全缺乏核物質(zhì)凛忿，要么它們的染色質(zhì)狀態(tài)與其轉(zhuǎn)錄的過(guò)程不一致。因此竞川，我們不能期望它們?cè)诖朔治鲋袝?huì)保持一致店溢。
DimPlot(coembed, split.by = "tech", group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend()
此外，我們還發(fā)現(xiàn)在B細(xì)胞祖細(xì)胞旁邊似乎也有一個(gè)細(xì)胞類群流译，該類群只存在于scATAC-seq數(shù)據(jù)中逞怨，并且與scRNA-seq數(shù)據(jù)的細(xì)胞整合得不好者疤。對(duì)這些細(xì)胞的元數(shù)據(jù)信息進(jìn)一步檢查發(fā)現(xiàn)福澡，大量的reads映射到black list區(qū)域（由10x Genomics QC指標(biāo)提供）。這表明這些細(xì)胞可能代表了死亡或垂死的細(xì)胞驹马，環(huán)境中的DNA或其他人工的污染革砸，導(dǎo)致其不存在于scRNA-seq數(shù)據(jù)集中除秀。
coembed$blacklist_region_fragments[is.na(coembed$blacklist_region_fragments)] <- 0
FeaturePlot(coembed, features = "blacklist_region_fragments", max.cutoff = 500)
How can I evaluate the success of cross-modality integration?
- Overall, the prediction scores are high which suggest a high degree of confidence in our cell type assignments.
- We observe good agreement between the scATAC-seq only dimensional reduction (UMAP plot above) and the transferred labels.
- The co-embedding based on the same set of anchors gives good mixing between the two modalities.
- When we collapse the ATAC-seq reads on a per-cluster basis into “pseudo bulk” profiles, we observe good concordance with the chromatin patterns observed in bulk data
參考來(lái)源：https://satijalab.org/seurat/v3.1/atacseq_integration_vignette.html
