当合并多个单细胞染色质数据集时，我们必须注意到，如果每个数据集都是独立的进行peak calling，则它们得到的peaks可能不是完全一致的。Seurat在处理时会把所有不完全相同的peaks视为不同的features。因此，我们在合并完多个数据集后需要创建一组通用的peaks。
我们可以使用GenomicRanges包中的函数创建统一的peaks集。其中，GenomicRanges中的reduce函数可以将所有相交的peaks进行合并，而disjoin函数则将创建不同的不重叠的peaks集。这里，我们用一个直观的示例来说明reduce和disjoin函数之间的区别：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("GenomicRanges")
library(GenomicRanges)
# 创建一个GRanges格式的peaks文件
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
reduce(gr)
GRanges object with 2 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 20-200 *
[2] chr1 300-400 *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
disjoin(gr)
GRanges object with 4 ranges and 0 metadata columns:
seqnames ranges strand
<Rle> <IRanges> <Rle>
[1] chr1 20-69 *
[2] chr1 70-120 *
[3] chr1 121-200 *
[4] chr1 300-400 *
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
image
在本示例中，我们将演示如何通过在每个包含一组共同peaks的对象中创建一个新的assay，来合并多个包含单细胞染色质数据的Seurat对象。
这里，我们使用四个10x Genomics平台产生的scATAC-seq的PBMC数据集进行演示：
Loading data 加载数据集
library(Signac)
library(Seurat)
# define a convenient function to load all the data and create a Seurat object
create_obj <- function(dir) {
count.path <- list.files(path = dir, pattern = "*_filtered_peak_bc_matrix.h5", full.names = TRUE)
fragment.path <- list.files(path = dir, pattern = "*_fragments.tsv.gz", full.names = TRUE)[1]
counts <- Read10X_h5(count.path)
md.path <- list.files(path = dir, pattern = "*_singlecell.csv", full.names = TRUE)
md <- read.table(file = md.path, stringsAsFactors = FALSE, sep = ",", header = TRUE, row.names = 1)
obj <- CreateSeuratObject(counts = counts, assay = "ATAC", meta.data = md)
obj <- SetFragments(obj, file = fragment.path)
return(obj)
}
pbmc500 <- create_obj("/home/dongwei/scATAC-seq/data/pbmc500")
pbmc1k <- create_obj("/home/dongwei/scATAC-seq/data/pbmc1k")
pbmc5k <- create_obj("/home/dongwei/scATAC-seq/data/pbmc5k")
pbmc10k <- create_obj("/home/dongwei/scATAC-seq/data/pbmc10k")
Creating a common peak set
如果每个实验中都单独进行了peak calling，那么它们最终得到的peaks区域可能不会完全重叠。我们可以合并所有数据集中的peaks区域以创建一个公共的peaks集，并且利用合并后的公共peaks集分别对每个实验重新进行定量。
我们可以使用几种不同的方法来创建一个公共的peaks集。一种是使用GenomicRanges包中的reduce或disjoin函数，另一种是直接使用Signac包中的UnifyPeaks函数从对象列表中提取peaks的坐标，并对peaks进行reduce或disjoin以创建单个不重叠的peaks集。
# 直接使用UnifyPeaks函数将多个不同对象的中peaks进行合并
combined.peaks <- UnifyPeaks(object.list = list(pbmc500, pbmc1k, pbmc5k, pbmc10k), mode = "reduce")
# 查看合并后的peaks集
combined.peaks
## GRanges object with 90694 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 *
## ... ... ... ...
## [90690] chrY 23583994-23584463 *
## [90691] chrY 23602466-23602779 *
## [90692] chrY 23899155-23899164 *
## [90693] chrY 28816593-28817710 *
## [90694] chrY 58855911-58856251 *
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengths
Quantify peaks in each dataset
得到合并好的公共peaks集后，我们可以使用FeatureMatrix函数对每个数据集基于公共的peaks集重新进行计数定量，并新建一个assay存储计数的数据。
# 使用FeatureMatrix函数对每个数据集重新进行计数定量
pbmc500.counts <- FeatureMatrix(
fragments = GetFragments(pbmc500),
features = combined.peaks,
sep = c(":", "-"),
cells = colnames(pbmc500)
)
pbmc1k.counts <- FeatureMatrix(
fragments = GetFragments(pbmc1k),
features = combined.peaks,
sep = c(":", "-"),
cells = colnames(pbmc1k)
)
pbmc5k.counts <- FeatureMatrix(
fragments = GetFragments(pbmc5k),
features = combined.peaks,
sep = c(":", "-"),
cells = colnames(pbmc5k)
)
pbmc10k.counts <- FeatureMatrix(
fragments = GetFragments(pbmc10k),
features = combined.peaks,
sep = c(":", "-"),
cells = colnames(pbmc10k)
)
# 使用CreateAssayObject函数新建一个assay对象存储计数的数据
pbmc500[['peaks']] <- CreateAssayObject(counts = pbmc500.counts)
pbmc1k[['peaks']] <- CreateAssayObject(counts = pbmc1k.counts)
pbmc5k[['peaks']] <- CreateAssayObject(counts = pbmc5k.counts)
pbmc10k[['peaks']] <- CreateAssayObject(counts = pbmc10k.counts)
Merge objects
Now that the objects each contain an assay with the same set of features, we can use the standard merge function from Seurat to merge the objects.
# 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
# 直接使用merge函数将多个数据集进行合并
combined <- merge(x = pbmc500, y = list(pbmc1k, pbmc5k, pbmc10k), add.cell.ids = c("500", "1k", "5k", "10k"))
# make sure to change to the assay containing common peaks
DefaultAssay(combined) <- "peaks"
# 对合并后的数据进行归一化，降维与可视化
combined <- RunTFIDF(combined)
combined <- FindTopFeatures(combined, min.cutoff = 20)
combined <- RunSVD(
combined,
reduction.key = 'LSI_',
reduction.name = 'lsi',
irlba.work = 400
)
combined <- RunUMAP(combined, dims = 2:30, reduction = 'lsi')
DimPlot(combined, group.by = 'dataset', pt.size = 0.1)
image
At this stage it is possible to proceed with all downstream analysis without creating a merged fragment file if you have computed quality control metrics and gene activities for each object individually prior to merging the datasets, and don’t need to plot coverage tracks with the merged data.
Merge fragment files
要创建合并的片段文件(fragment files)，我们需要先下载和解压缩这些文件，添加相应的细胞ID并进行文件的合并，最后对合并后的文件进行压缩并构建索引。
这些是直接在命令行下操作的，而不是在R中执行的，请确保您已安装了tabix和bgzip程序。
# decompress files and add the same cell prefix as was added to the Seurat object
gzip -dc atac_pbmc_500_nextgem_fragments.tsv.gz | awk 'BEGIN {FS=OFS="\t"} {print $1,$2,$3,"500_"$4,$5}' - > pbmc500_fragments.tsv
gzip -dc atac_pbmc_1k_nextgem_fragments.tsv.gz | awk 'BEGIN {FS=OFS="\t"} {print $1,$2,$3,"1k_"$4,$5}' - > pbmc1k_fragments.tsv
gzip -dc atac_pbmc_5k_nextgem_fragments.tsv.gz | awk 'BEGIN {FS=OFS="\t"} {print $1,$2,$3,"5k_"$4,$5}' - > pbmc5k_fragments.tsv
gzip -dc atac_pbmc_10k_nextgem_fragments.tsv.gz | awk 'BEGIN {FS=OFS="\t"} {print $1,$2,$3,"10k_"$4,$5}' - > pbmc10k_fragments.tsv
# merge files (avoids having to re-sort)
sort -m -k 1,1V -k2,2n pbmc500_fragments.tsv pbmc1k_fragments.tsv pbmc5k_fragments.tsv pbmc10k_fragments.tsv > fragments.tsv
# block gzip compress the merged file
bgzip -@ 4 fragments.tsv # -@ 4 uses 4 threads
# index the bgzipped file
tabix -p bed fragments.tsv.gz
# remove intermediate files
rm pbmc500_fragments.tsv pbmc1k_fragments.tsv pbmc5k_fragments.tsv pbmc10k_fragments.tsv
接下来，我们可以将合并后的片段文件的路径添加到合并后的Seurat对象中。
combined <- SetFragments(combined, "/home/dongwei/scATAC-seq/data/pbmc_combined/fragments.tsv.gz")
# 使用CoveragePlot函数可视化peaks的信息
CoveragePlot(
object = combined,
group.by = 'dataset',
region = "chr14-99700000-99760000",
peaks = StringToGRanges(rownames(combined), sep = c(":", "-"))
)
image
