當合并多個單細胞染色質(zhì)數(shù)據(jù)集時蔑滓，我們必須注意到憨闰，如果每個數(shù)據(jù)集都是獨立的進行peak calling，則它們得到的peaks可能不是完全一致的确封。Seurat在處理時會把所有不完全相同的peaks視為不同的features除呵。因此，我們在合并完多個數(shù)據(jù)集后需要創(chuàng)建一組通用的peaks爪喘。
我們可以使用GenomicRanges包中的函數(shù)創(chuàng)建統(tǒng)一的peaks集颜曾。其中，GenomicRanges中的reduce函數(shù)可以將所有相交的peaks進行合并秉剑，而disjoin函數(shù)則將創(chuàng)建不同的不重疊的peaks集泛豪。這里，我們用一個直觀的示例來說明reduce和disjoin函數(shù)之間的區(qū)別：
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("GenomicRanges")
library(GenomicRanges)
# 創(chuàng)建一個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
在本示例中侦鹏，我們將演示如何通過在每個包含一組共同peaks的對象中創(chuàng)建一個新的assay诡曙，來合并多個包含單細胞染色質(zhì)數(shù)據(jù)的Seurat對象。
這里略水，我們使用四個10x Genomics平臺產(chǎn)生的scATAC-seq的PBMC數(shù)據(jù)集進行演示：
Loading data 加載數(shù)據(jù)集
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區(qū)域可能不會完全重疊。我們可以合并所有數(shù)據(jù)集中的peaks區(qū)域以創(chuàng)建一個公共的peaks集渊涝，并且利用合并后的公共peaks集分別對每個實驗重新進行定量慎璧。
我們可以使用幾種不同的方法來創(chuàng)建一個公共的peaks集。一種是使用GenomicRanges包中的reduce或disjoin函數(shù)驶赏，另一種是直接使用Signac包中的UnifyPeaks函數(shù)從對象列表中提取peaks的坐標炸卑，并對peaks進行reduce或disjoin以創(chuàng)建單個不重疊的peaks集既鞠。
# 直接使用UnifyPeaks函數(shù)將多個不同對象的中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函數(shù)對每個數(shù)據(jù)集基于公共的peaks集重新進行計數(shù)定量，并新建一個assay存儲計數(shù)的數(shù)據(jù)嘱蛋。
# 使用FeatureMatrix函數(shù)對每個數(shù)據(jù)集重新進行計數(shù)定量
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函數(shù)新建一個assay對象存儲計數(shù)的數(shù)據(jù)
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函數(shù)將多個數(shù)據(jù)集進行合并
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"
# 對合并后的數(shù)據(jù)進行歸一化蚯姆，降維與可視化
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)
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
要創(chuàng)建合并的片段文件(fragment files)五续，我們需要先下載和解壓縮這些文件，添加相應的細胞ID并進行文件的合并龄恋，最后對合并后的文件進行壓縮并構(gòu)建索引疙驾。
這些是直接在命令行下操作的，而不是在R中執(zhí)行的郭毕，請確保您已安裝了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函數(shù)可視化peaks的信息
CoveragePlot(
object = combined,
group.by = 'dataset',
region = "chr14-99700000-99760000",
peaks = StringToGRanges(rownames(combined), sep = c(":", "-"))
)
