合并包含单细胞染色质数据的多个 Seurat 对象。如果对每个数据集独立进行 peak calling，则峰值不太可能完全相同。因此，需要在要合并的数据集上创建一组共同的峰值。
创建统一的峰的方法：
- 使用 GenomicRanges 包中的 reduce 函数整合所有相交的峰值；
-
使用 disjoin 函数创建不同的非重叠峰集；
不同方法示意图.png
gr <- GRanges(seqnames = "chr1", ranges = IRanges(start = c(20, 70, 300), end = c(120, 200, 400)))
gr
## GRanges object with 3 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 20-120 *
## [2] chr1 70-200 *
## [3] chr1 300-400 *
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggbio package.
## Please report the issue at <https://github.com/lawremi/ggbio/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
创建共同的峰集
在每个实验中单独识别峰值，那么它们可能不会完全重叠。可以合并来自所有数据集的峰值以创建一个公共峰值集，并在合并对象之前在每个实验中量化该峰值集。
首先，加载每个实验的峰值坐标并将其转换为基因组范围，然后使用 GenomicRanges::reduce函数创建一组公共峰值以在每个数据集中进行量化。
library(Signac)
library(Seurat)
library(GenomicRanges)
library(future)
plan("multicore", workers = 4)
options(future.globals.maxSize = 50000 * 1024^2) # for 50 Gb RAM
# read in peak sets
peaks.500 <- read.table(
file = "pbmc500/atac_pbmc_500_nextgem_peaks.bed",
col.names = c("chr", "start", "end")
)
peaks.1k <- read.table(
file = "pbmc1k/atac_pbmc_1k_nextgem_peaks.bed",
col.names = c("chr", "start", "end")
)
peaks.5k <- read.table(
file = "pbmc5k/atac_pbmc_5k_nextgem_peaks.bed",
col.names = c("chr", "start", "end")
)
peaks.10k <- read.table(
file = "pbmc10k/atac_pbmc_10k_nextgem_peaks.bed",
col.names = c("chr", "start", "end")
)
# convert to genomic ranges
gr.500 <- makeGRangesFromDataFrame(peaks.500)
gr.1k <- makeGRangesFromDataFrame(peaks.1k)
gr.5k <- makeGRangesFromDataFrame(peaks.5k)
gr.10k <- makeGRangesFromDataFrame(peaks.10k)
# Create a unified set of peaks to quantify in each dataset
combined.peaks <- reduce(x = c(gr.500, gr.1k, gr.5k, gr.10k))
# Filter out bad peaks based on length
peakwidths <- width(combined.peaks)
combined.peaks <- combined.peaks[peakwidths < 10000 & peakwidths > 20]
combined.peaks
## GRanges object with 89951 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 565153-565499 *
## [2] chr1 569185-569620 *
## [3] chr1 713551-714783 *
## [4] chr1 752418-753020 *
## [5] chr1 762249-763345 *
## ... ... ... ...
## [89947] chrY 23422151-23422632 *
## [89948] chrY 23583994-23584463 *
## [89949] chrY 23602466-23602779 *
## [89950] chrY 28816593-28817710 *
## [89951] chrY 58855911-58856251 *
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengths
创建 Fragment 对象
对每个数据创建一个 Fragment 对象，去量化峰值组合。Fragment 对象是在 Signac 中定义的专门用于保存与单个片段文件相关的所有信息。
首先，加载每个数据的细胞元数据，了解细胞条形码；然后，使用 CreateFragmentObject 函数创建 Fragment 对象。
CreateFragmentObject 函数检查以确保文件存在于磁盘上，并被压缩和索引，计算文件和tabix索引的MD5和，以便我们可以判断文件是否在任何时候被修改，并检查文件中是否存在预期的单元格。
# load metadata
md.500 <- read.table(
file = "pbmc500/atac_pbmc_500_nextgem_singlecell.csv",
stringsAsFactors = FALSE,
sep = ",",
header = TRUE,
row.names = 1
)[-1, ] # remove the first row
md.1k <- read.table(
file = "pbmc1k/atac_pbmc_1k_nextgem_singlecell.csv",
stringsAsFactors = FALSE,
sep = ",",
header = TRUE,
row.names = 1
)[-1, ]
md.5k <- read.table(
file = "pbmc5k/atac_pbmc_5k_nextgem_singlecell.csv",
stringsAsFactors = FALSE,
sep = ",",
header = TRUE,
row.names = 1
)[-1, ]
md.10k <- read.table(
file = "pbmc10k/atac_pbmc_10k_nextgem_singlecell.csv",
stringsAsFactors = FALSE,
sep = ",",
header = TRUE,
row.names = 1
)[-1, ]
# perform an initial filtering of low count cells
md.500 <- md.500[md.500$passed_filters > 500, ]
md.1k <- md.1k[md.1k$passed_filters > 500, ]
md.5k <- md.5k[md.5k$passed_filters > 500, ]
md.10k <- md.10k[md.10k$passed_filters > 1000, ] # sequenced deeper so set higher cutoff
# create fragment objects
frags.500 <- CreateFragmentObject(
path = "pbmc500/atac_pbmc_500_nextgem_fragments.tsv.gz",
cells = rownames(md.500)
)
frags.1k <- CreateFragmentObject(
path = "pbmc1k/atac_pbmc_1k_nextgem_fragments.tsv.gz",
cells = rownames(md.1k)
)
frags.5k <- CreateFragmentObject(
path = "pbmc5k/atac_pbmc_5k_nextgem_fragments.tsv.gz",
cells = rownames(md.5k)
)
frags.5k <- CreateFragmentObject(
path = "pbmc5k/atac_pbmc_5k_nextgem_fragments.tsv.gz",
cells = rownames(md.5k)
)
frags.10k <- CreateFragmentObject(
path = "pbmc10k/atac_pbmc_10k_nextgem_fragments.tsv.gz",
cells = rownames(md.10k)
)
量化每个数据集中的峰值
使用 FeatureMatrix 函数为每个样本创建一个 peaks x cell 矩阵。该函数使用 future 包并行加速。
pbmc500.counts <- FeatureMatrix(
fragments = frags.500,
features = combined.peaks,
cells = rownames(md.500)
)
pbmc1k.counts <- FeatureMatrix(
fragments = frags.1k,
features = combined.peaks,
cells = rownames(md.1k)
)
pbmc5k.counts <- FeatureMatrix(
fragments = frags.5k,
features = combined.peaks,
cells = rownames(md.5k)
)
pbmc10k.counts <- FeatureMatrix(
fragments = frags.10k,
features = combined.peaks,
cells = rownames(md.10k)
)
创建对象
使用量化矩阵为每个数据集创建一个 Seurat 对象，将每个数据集的 Fragment 对象存储在 assay 中。
pbmc500_assay <- CreateChromatinAssay(pbmc500.counts, fragments = frags.500)
pbmc500 <- CreateSeuratObject(pbmc500_assay, assay = "ATAC", meta.data=md.500)
pbmc1k_assay <- CreateChromatinAssay(pbmc1k.counts, fragments = frags.1k)
pbmc1k <- CreateSeuratObject(pbmc1k_assay, assay = "ATAC", meta.data=md.1k)
pbmc5k_assay <- CreateChromatinAssay(pbmc5k.counts, fragments = frags.5k)
pbmc5k <- CreateSeuratObject(pbmc5k_assay, assay = "ATAC", meta.data=md.5k)
pbmc10k_assay <- CreateChromatinAssay(pbmc10k.counts, fragments = frags.10k)
pbmc10k <- CreateSeuratObject(pbmc10k_assay, assay = "ATAC", meta.data=md.10k)
合并对象
现在每个对象都包含一个相同特征集的 assay，使用标准 merge 函数合并对象。这也将合并所有片段对象，以便我们保留最终合并对象中的每个细胞的片段信息。
# add information to identify dataset of origin
pbmc500$dataset <- 'pbmc500'
pbmc1k$dataset <- 'pbmc1k'
pbmc5k$dataset <- 'pbmc5k'
pbmc10k$dataset <- 'pbmc10k'
# merge all datasets, adding a cell ID to make sure cell names are unique
combined <- merge(
x = pbmc500,
y = list(pbmc1k, pbmc5k, pbmc10k),
add.cell.ids = c("500", "1k", "5k", "10k")
)
combined[["ATAC"]]
## ChromatinAssay data with 89951 features for 21688 cells
## Variable features: 0
## Genome:
## Annotation present: FALSE
## Motifs present: FALSE
## Fragment files: 4
标准分析流程：
combined <- RunTFIDF(combined)
combined <- FindTopFeatures(combined, min.cutoff = 20)
combined <- RunSVD(combined)
combined <- RunUMAP(combined, dims = 2:50, reduction = 'lsi')
DimPlot(combined, group.by = 'dataset', pt.size = 0.1)
image.png
合并的对象包含四个 fragment 对象，也包含每个片段文件中的细胞的内部映射，以便不修改细胞名而检索片段信息。通过绘制基因组的一个区域来检查从片段文件中提取出数据是否正常。
CoveragePlot(
object = combined,
group.by = 'dataset',
region = "chr14-99700000-99760000"
)
不同数据UMAP.png
没有公共特征集的合并
上述方法要求能够访问每个数据集的片段文件。若没有这些数据，仍可以创建一个合并对象，但合并后的计数矩阵可能不那么准确。
Signac 中用于合并 chroinassay对象的 merge 函数，将考虑重叠的峰，并调整峰所跨越的基因组范围，使每个被合并对象的特征变得等效。注意，这可能导致计数矩阵不准确，因为一些峰值将扩展到覆盖最初未量化的区域。
# load the count matrix for each object that was generated by cellranger
counts.500 <- Read10X_h5("pbmc500/atac_pbmc_500_nextgem_filtered_peak_bc_matrix.h5")
counts.1k <- Read10X_h5("pbmc1k/atac_pbmc_1k_nextgem_filtered_peak_bc_matrix.h5")
counts.5k <- Read10X_h5("pbmc5k/atac_pbmc_5k_nextgem_filtered_peak_bc_matrix.h5")
counts.10k <- Read10X_h5("pbmc10k/atac_pbmc_10k_nextgem_filtered_peak_bc_matrix.h5")
# create objects
pbmc500_assay <- CreateChromatinAssay(counts = counts.500, sep = c(":", "-"), min.features = 500)
pbmc500 <- CreateSeuratObject(pbmc500_assay, assay = "peaks")
pbmc1k_assay <- CreateChromatinAssay(counts = counts.1k, sep = c(":", "-"), min.features = 500)
pbmc1k <- CreateSeuratObject(pbmc1k_assay, assay = "peaks")
pbmc5k_assay <- CreateChromatinAssay(counts = counts.5k, sep = c(":", "-"), min.features = 500)
pbmc5k <- CreateSeuratObject(pbmc5k_assay, assay = "peaks")
pbmc10k_assay <- CreateChromatinAssay(counts = counts.10k, sep = c(":", "-"), min.features = 1000)
pbmc10k <- CreateSeuratObject(pbmc10k_assay, assay = "peaks")
# add information to identify dataset of origin
pbmc500$dataset <- 'pbmc500'
pbmc1k$dataset <- 'pbmc1k'
pbmc5k$dataset <- 'pbmc5k'
pbmc10k$dataset <- 'pbmc10k'
# merge
combined <- merge(
x = pbmc500,
y = list(pbmc1k, pbmc5k, pbmc10k),
add.cell.ids = c("500", "1k", "5k", "10k")
)
# process
combined <- RunTFIDF(combined)
combined <- FindTopFeatures(combined, min.cutoff = 20)
combined <- RunSVD(combined)
combined <- RunUMAP(combined, dims = 2:50, reduction = 'lsi')
DimPlot(combined, group.by = 'dataset', pt.size = 0.1)
image.png
