使用Signac包进行单细胞ATAC-seq数据分析（三）：sc
2020-04-30 本文已影响0人
Davey1220
在本教程中，我们将学习使用Signac包对多样本的scATAC-seq数据进行整合分析。这里，我们对来自10x Genomics和sci-ATAC-seq技术测序的成年小鼠大脑的多个单细胞ATAC-seq数据集进行了整合分析。
其中，10x Genomics平台产生的原始数据可从官网下载：
- The Raw data
- The Metadata
- The fragments file
- The fragments file index
sci-ATAC-seq技术产生的数据集由Cusanovich和Hill等人生成。原始数据可从作者的网站下载：
我们将演示使用Seurat v3中的数据整合方法（dataset integration and label transfer）对多个scATAC-seq数据集进行整合分析，以及使用harmony包进行数据整合。
安装并加载所需的R包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("rtracklayer")
library(devtools)
install_github("immunogenomics/harmony")
library(Signac)
library(Seurat)
library(patchwork)
set.seed(1234)
加载数据集并构建Seurat对象
# this object was created following the mouse brain vignette
# 加载10x Genomics的小鼠大脑的scATAC-seq数据
tenx <- readRDS(file = "/home/dongwei/scATAC-seq/brain/10x/adult_mouse_brain.rds")
tenx$tech <- '10x'
tenx$celltype <- Idents(tenx)
# 加载sci-ATAC-seq的小鼠大脑数据
sci.metadata <- read.table(
file = "/home/dongwei/scATAC-seq/brain/sci/cell_metadata.txt",
header = TRUE,
row.names = 1,
sep = "\t"
)
# subset to include only the brain data
sci.metadata <- sci.metadata[sci.metadata$tissue == 'PreFrontalCortex', ]
sci.counts <- readRDS(file = "/home/dongwei/scATAC-seq/brain/sci/atac_matrix.binary.qc_filtered.rds")
sci.counts <- sci.counts[, rownames(x = sci.metadata)]
数据预处理
在上述两个数据集中，sci-ATAC-seq数据是比对到小鼠mm9参考基因组的，而10x的数据是比对到小鼠mm10参考基因组的，因此这两个数据集中peaks的基因组坐标信息是不同的。我们可以使用rtracklayer包将mm9参考基因组的坐标信息转换到mm10中，并使用mm10的坐标更换sci-ATAC-seq数据中peaks的坐标，其中liftover转换的chain文件可以从UCSC官网进行下载。
# 将peaks坐标信息转换成GRanges格式
sci_peaks_mm9 <- StringToGRanges(regions = rownames(sci.counts), sep = c("_", "_"))
# 导入mm9ToMm10.over.chain文件
mm9_mm10 <- rtracklayer::import.chain("/home/dongwei/data/liftover/mm9ToMm10.over.chain")
# 使用rtracklayer包中的liftOver函数转换坐标信息
sci_peaks_mm10 <- rtracklayer::liftOver(x = sci_peaks_mm9, chain = mm9_mm10)
names(sci_peaks_mm10) <- rownames(sci.counts)
# discard any peaks that were mapped to >1 region in mm10
correspondence <- S4Vectors::elementNROWS(sci_peaks_mm10)
sci_peaks_mm10 <- sci_peaks_mm10[correspondence == 1]
sci_peaks_mm10 <- unlist(sci_peaks_mm10)
sci.counts <- sci.counts[names(sci_peaks_mm10), ]
# rename peaks with mm10 coordinates
rownames(sci.counts) <- GRangesToString(grange = sci_peaks_mm10)
# create Seurat object and perform some basic QC filtering
# 构建Seurat对象
sci <- CreateSeuratObject(
counts = sci.counts,
meta.data = sci.metadata,
assay = 'peaks',
project = 'sci'
)
# 数据过滤
sci.ds <- sci[, sci$nFeature_peaks > 2000 & sci$nCount_peaks > 5000 & !(sci$cell_label %in% c("Collisions", "Unknown"))]
sci$tech <- 'sciATAC'
# 使用RunTFIDF函数进行数据归一化
sci <- RunTFIDF(sci)
# 使用FindTopFeatures函数提取高变异的peaks
sci <- FindTopFeatures(sci, min.cutoff = 50)
# 使用RunSVD函数进行线性降维
sci <- RunSVD(sci, n = 30, reduction.name = 'lsi', reduction.key = 'LSI_')
# 使用RunUMAP函数进行非线性降维
sci <- RunUMAP(sci, reduction = 'lsi', dims = 2:30)
现在，我们构建好了两个scATAC-seq对象，并且它们都含有基于相同的mm10参考基因组坐标系统得到的peaks信息。但是，由于这两个实验都单独进行了peak calling，因此这两个数据集中得到的peaks坐标不太可能完全重叠。为了在我们要整合的数据集中具有共同的特征，我们可以基于10x Genomics数据集对sci-ATAC-seq中peaks的reads进行计数，并使用这些计数创建一个新的assay。
# find peaks that intersect in both datasets
# 使用GetIntersectingFeatures函数提取两个数据集中重叠的peak区域
intersecting.regions <- GetIntersectingFeatures(
object.1 = sci,
object.2 = tenx,
sep.1 = c("-", "-"),
sep.2 = c(":", "-")
)
# choose a subset of intersecting peaks
peaks.use <- sample(intersecting.regions[[1]], size = 10000, replace = FALSE)
# count fragments per cell overlapping the set of peaks in the 10x data
# 使用FeatureMatrix函数对peaks中的reads进行计数
sci_peaks_tenx <- FeatureMatrix(
fragments = GetFragments(object = tenx, assay = 'peaks'),
features = StringToGRanges(peaks.use),
cells = colnames(tenx)
)
# create a new assay and add it to the 10x dataset
# 使用CreateAssayObject函数新建一个assay对象
tenx[['sciPeaks']] <- CreateAssayObject(counts = sci_peaks_tenx, min.cells = 1)
# 数据归一化
tenx <- RunTFIDF(object = tenx, assay = 'sciPeaks')
多样本scATAC-seq数据集的整合
在进行数据整合之前，我们最好先检查下是否存在数据集特异的差异，并将其删除。如果没有，我们可以简单地将多个对象进行合并而不执行整合。在本示例中，由于使用不同的测序技术，两个数据集之间存在很大的差异。我们可以使用Seurat v3中的数据整合方法来消除这种影响。
# 先简单的将两个数据集进行合并看一下聚类的效果
# Look at the data without integration first
# 使用MergeWithRegions函数将两个数据对象进行合并
unintegrated <- MergeWithRegions(
object.1 = sci,
object.2 = tenx,
assay.1 = 'peaks',
assay.2 = 'sciPeaks',
sep.1 = c("-", "-"),
sep.2 = c("-", "-")
)
# 对合并后数据进行归一化，特征选择和降维可视化
unintegrated <- RunTFIDF(unintegrated)
unintegrated <- FindTopFeatures(unintegrated, min.cutoff = 50)
unintegrated <- RunSVD(unintegrated, n = 30, reduction.name = 'lsi', reduction.key = 'LSI_')
unintegrated <- RunUMAP(unintegrated, reduction = 'lsi', dims = 2:30)
p1 <- DimPlot(unintegrated, group.by = 'tech', pt.size = 0.1) + ggplot2::ggtitle("Unintegrated")
# 使用Seurat v3的数据整合方法进行数据集的整合
# find integration anchors between 10x and sci-ATAC
# 使用FindIntegrationAnchors函数识别共享的整合anchors
anchors <- FindIntegrationAnchors(
object.list = list(tenx, sci),
anchor.features = peaks.use,
assay = c('sciPeaks', 'peaks'),
k.filter = NA
)
# integrate data and create a new merged object
# 使用IntegrateData函数根据识别的anchors进行数据整合
integrated <- IntegrateData(
anchorset = anchors,
weight.reduction = sci[['lsi']],
dims = 2:30,
preserve.order = TRUE
)
# we now have a "corrected" TF-IDF matrix, and can run LSI again on this corrected matrix
# 对整合后的数据进行降维可视化
integrated <- RunSVD(
object = integrated,
n = 30,
reduction.name = 'integratedLSI'
)
integrated <- RunUMAP(
object = integrated,
dims = 2:30,
reduction = 'integratedLSI'
)
p2 <- DimPlot(integrated, group.by = 'tech', pt.size = 0.1) + ggplot2::ggtitle("Integrated")
p1 + p2
image
Label transfer
我们还可以使用Seurat v3中的Label transfer方法进行数据集的整合，它将数据从一个query数据集映射到另一个reference数据集中。在这里，我们通过将细胞类型标签从10x Genomics scATAC-seq数据映射到到sci-ATAC-seq数据中。
# 使用FindTransferAnchors函数识别整合的anchors
transfer.anchors <- FindTransferAnchors(
reference = tenx,
query = sci,
reference.assay = 'sciPeaks',
query.assay = 'peaks',
reduction = 'cca',
features = peaks.use,
k.filter = NA
)
# 使用TransferData函数基于识别好的anchors进行数据映射整合
predicted.id <- TransferData(
anchorset = transfer.anchors,
refdata = tenx$celltype,
weight.reduction = sci[['lsi']],
dims = 2:30
)
sci <- AddMetaData(
object = sci,
metadata = predicted.id
)
sci$predicted.id <- factor(sci$predicted.id, levels = levels(tenx$celltype)) # to make the colors match
# 数据可视化
p3 <- DimPlot(tenx, group.by = 'celltype', label = TRUE) + NoLegend() + ggplot2::ggtitle("Celltype labels 10x scATAC-seq")
p4 <- DimPlot(sci, group.by = 'predicted.id', label = TRUE) + NoLegend() + ggplot2::ggtitle("Predicted labels sci-ATAC-seq")
p3 + p4
image
Integration with Harmony使用Harmony包进行数据整合
Harmony需要一个对象作为输入，因此这里我们使用MergeWithRegions函数以coordinate-aware的方式将sci-ATAC-seq数据集和10x的scATAC-seq数据集进行合并，然后对合并后的对象进行LSI线性降维。数据降维后，我们可以使用RunHarmony函数调用Harmony的方法进行数据的整合，并提供用作分组变量的技术来消除sci-ATAC-seq和10x Genomics scATAC-seq数据集之间的批次差异。这会产生一组“校正”的LSI嵌入，可以进一步用于UMAP或tSNE降维，并进行细胞聚类分群。
# 加载harmony包
library(harmony)
# 使用RunHarmony函数进行数据整合
hm.integrated <- RunHarmony(
object = unintegrated,
group.by.vars = 'tech',
reduction = 'lsi',
assay.use = 'peaks',
project.dim = FALSE
)
# re-compute the UMAP using corrected LSI embeddings
# 数据降维可视化
hm.integrated <- RunUMAP(hm.integrated, dims = 2:30, reduction = 'harmony')
p5 <- DimPlot(hm.integrated, group.by = 'tech', pt.size = 0.1) + ggplot2::ggtitle("Harmony integration")
p1 + p5
image
