单细胞ATAC数据的两种读取方式（下）

本期我们将来介绍如何使用Signac读入单细胞ATAC数据。
Signac是由Seurat同一团队开发的用于分析单细胞染色质数据的R包，能够分析包括scATAC-seq、scCUT&Tag和 scNTT-seq以及多模态数据集在内的染色质数据。
在上一期中，我们提到在单细胞ATAC下游分析中用到的主要文件包括filtered_peak_bc_matrix、filtered_peak_bc_matrix.h5、fragments.tsv.gz和singlecell.csv。下面，我们就来看一下如何使用Signac读取标准的Cell Ranger的输出文件。
Signac使用Cell Ranger ATAC生成的峰/细胞矩阵和细胞元数据创建Seurat对象，并将磁盘上fragment文件的路径存储在Seurat对象中。

library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v86)
library(ggplot2)
library(patchwork)

counts <- Read10X_h5(filename = "/path/to/filtered_peak_bc_matrix.h5")
metadata <- read.csv(
  file = "/path/to/singlecell.csv",
  header = TRUE,
  row.names = 1
)

chrom_assay <- CreateChromatinAssay(
  counts = counts,
  sep = c(":", "-"),
  fragments = '/path/to/fragments.tsv.gz',
  min.cells = 10,
  min.features = 200
)

pbmc <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = metadata
)

如果我们想要复现一个公开的数据，而这个数据中又没有matrix.h5文件，该如何读入ATAC数据呢？此时我们就要使用filtered_peak_bc_matrix文件夹和Matrix::readMM()函数加载数据。

counts <- Matrix::readMM("/path/to/filtered_peak_bc_matrix/matrix.mtx")
barcodes <- readLines("/path/to/filtered_peak_bc_matrix/barcodes.tsv")
peaks <- read.table("/path/to/filtered_peak_bc_matrix/peaks.bed", sep="\t")
peaknames <- paste(peaks$V1, peaks$V2, peaks$V3, sep="-")

colnames(counts) <- barcodes
rownames(counts) <- peaknames

如果我们只有fragments.tsv.gz文件，又该如何读入呢？在这种情况下可以使用FeatureMatrix()函数创建计数矩阵。

fragpath <- '/path/to/fragments.tsv.gz'

# Define cells
# If you already have a list of cell barcodes to use you can skip this step
total_counts <- CountFragments(fragpath)
cutoff <- 1000 # Change this number depending on your dataset!
barcodes <- total_counts[total_counts$frequency_count > cutoff, ]$CB

# Create a fragment object
frags <- CreateFragmentObject(path = fragpath, cells = barcodes)

# First call peaks on the dataset
# If you already have a set of peaks you can skip this step
peaks <- CallPeaks(frags)

# Quantify fragments in each peak
counts <- FeatureMatrix(fragments = frags, features = peaks, cells = barcodes)

经由以上三种方式，我们就可以得到一个包含ATAC peak数据和fragment数据的Seurat对象。ATAC-seq数据存储在特定的ChromatinAssay中，其中还包含了motif、基因注释和基因组等额外信息。
我们还可以将基因注释信息添加到Seurat对象中，以便下游函数直接从Seurat对象中提取基因注释信息。

# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)

# change to UCSC style since the data was mapped to hg19
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg38"

# add the gene information to the object
Annotation(pbmc) <- annotations

这样，我们就使用Signac完成了对单细胞ATAC数据的读入。

参考资料：
https://www.10xgenomics.com/datasets/10k-human-pbmcs-atac-v2-chromium-x-2-standard
https://stuartlab.org/signac/articles/pbmc_vignette

相关阅读：
单细胞ATAC数据的两种读取方式（上）