单细胞研究使我们能够了解细胞发育的分子过程和病理的发生机理。单细胞ATAC-seq（scATAC-seq）通过分析单细胞水平的染色质可及性来定义转录和表观遗传变化。本文将为大家介绍两种单细胞ATAC数据的读取方式。
首先，我们先来看一下进行分析的单细胞ATAC数据文件有哪些。以10x Single Cell ATAC数据为例，Cell Ranger ATAC的输出目录结构如下：
outs/
│ filtered_peak_bc_matrix.h5 过滤后hdf5格式的Peak矩阵
│ filtered_tf_bc_matrix.h5 过滤后hdf5格式的TF矩阵
│ fragments.tsv.gz fragment文件
│ fragments.tsv.gz.tbi fragment文件索引
│ peaks.bed Peak位置的bed文件
│ peak_annotation.tsv Peak对基因的注释文件
│ peak_motif_mapping.bed Peak与Motif的关联文件
│ raw_peak_bc_matrix.h5 原始hdf5格式的Peak矩阵
│ singlecell.csv 每个barcode元信息的统计文件
│ summary.csv 重要指标的统计文件
│ summary.html html格式分析报告
│ summary.json 所有数据指标汇总的json文件
├─filtered_peak_bc_matrix 过滤后mex格式的Peak矩阵
│ barcodes.tsv
│ matrix.mtx
│ peaks.bed
├─filtered_tf_bc_matrix 过滤后mex格式的TF矩阵
│ barcodes.tsv.gz
│ matrix.mtx.gz
│ motifs.tsv
└─raw_peak_bc_matrix 原始mex格式的Peak矩阵
barcodes.tsv
matrix.mtx
peaks.bed
在下游分析中用到的主要文件包括filtered_peak_bc_matrix、filtered_peak_bc_matrix.h5、fragments.tsv.gz和singlecell.csv。
下面我们就来介绍如何将上述文件导入单细胞ATAC分析软件中。
目前有两种常用的单细胞ATAC分析工具：ArchR和Signac，二者均为R包。我们首先来看如何使用ArchR读入单细胞ATAC数据。
ArchR是由斯坦福大学William J. Greenleaf实验室开发的进行单细胞ATAC数据分析的工具，但是由于实验室人员变动该软件目前已经停止更新和维护.
ArchR中分析的基本单元称为Arrow文件，每个Arrow文件存储与单个样本相关的所有数据。创建Arrow文件所需的输入文件为fragment文件。Fragment文件是tabix排序的文本文件，包含每个scATAC-seq fragment和相应的barcode，每行一个fragment。
在使用ArchR时，我们要做的第一件事是更改工作目录、设置使用的线程数和加载基因组和基因组注释。
library(ArchR)
addArchRThreads(threads = 16)
addArchRGenome("hg38") ###hg19, hg38, mm10
setwd('/path/to/outdir/')
然后，将待读入的fragment文件路径写为一个向量，可同时读入多个fragment文件。以样本名为向量的每个元素命名。
inputFiles <- c('/path/to/sample1/fragments.tsv.gz', '/path/to/sample2/fragments.tsv.gz')
names(inputFiles) <- c('sample1', 'sample2')
使用createArrowFiles函数进行数据读取。
ArrowFiles <- createArrowFiles(
inputFiles = inputFiles,
sampleNames = names(inputFiles),
filterTSS = 4,
filterFrags = 1000,
addTileMat = TRUE,
addGeneScoreMat = TRUE
)
这一步大概耗时10-15分钟，其中，对每一个样本的操作为：
- 从输入文件中读取开放fragment；
- 计算每个细胞的质控信息（即TSS富集分数和核小体信息）；
- 根据质控参数过滤细胞；
- 使用500 bp bin创建全基因组TileMatrix；
- 使用调用addArchRGenome()时定义的geneAnnotation创建GeneScoreMatrix。
运行完成后，会在输出目录下生成sample1.arrow和sample2.arrow两个文件，而ArrowFiles对象实际只是Arrow文件路径的字符向量。
ArrowFiles
[1] "sample1.arrow" "sample2.arrow"
与此同时，在输出目录下还生成一个名为QualityControl的目录，其中保存了每个样本的2个质控图片：TSS富集分数图和fragment分布图。
细心的朋友可能会发现，通过读入fragment文件得到的矩阵细胞数与html报告中的细胞数并不一致，这是因为fragments.tsv.gz相当于fragment的原始文件，包含每个cell barcode的fragment信息，而报告中的细胞数是通过骤降图截取后得出的预估细胞数。如果我们想要获取与报告中细胞数一致的矩阵，则可以按照如下方法进行处理。
ArchR中的getValidBarcodes()函数可以读取singlecell.csv文件并识别通过QC的细胞的barcode。但就如我们前面所说，ArchR已经停止更新维护，而Cell Ranger在升级至2.0后，singlecell.csv文件的格式发生了变化，导致目前的ArchR无法对v2.0以上版本的Cell Ranger产生的singlecell.csv进行提取有效barcode的处理。不过我们可以根据源码进行修改，完成这一操作。
getValidBarcodes <- function(csvFiles = NULL, sampleNames = NULL){
if(length(sampleNames) != length(csvFiles)){
stop("csvFiles and sampleNames must exist!")
}
if(!all(file.exists(csvFiles))){
stop("Not All csvFiles exists!")
}
.suppressAll <- function(expr = NULL){
suppressPackageStartupMessages(suppressMessages(suppressWarnings(expr)))
}
barcodeList <- lapply(seq_along(csvFiles), function(x){
df <- .suppressAll(data.frame(readr::read_csv(csvFiles[x])))
as.character(df[which(paste0(df$is__cell_barcode) == "1"),]$barcode)
###is__cell_barcode为singlecell.csv文件中指示细胞是否通过QC的列名
}) %>% SimpleList
names(barcodeList) <- sampleNames
barcodeList
}
singleFiles <- c('/path/to/sample1/singlecell.csv', '/path/to/sample2/singlecell.csv')
names(singleFiles) <- c('sample1', 'sample2')
validBC <- getValidBarcodes(singleFiles, sampleNames=names(singleFiles))
ArrowFiles <- createArrowFiles(
inputFiles = inputFiles,
sampleNames = names(inputFiles),
addTileMat = TRUE,
addGeneScoreMat = TRUE,
validBarcodes = validBC
)
这样，我们就使用getValidBarcodes和validBarcodes参数结合，完成了对特定barcode的fragment的提取。
通过以上方式产生的Arrow文件就可以使用ArchR继续进行后续的处理和分析啦。
参考资料：
https://www.10xgenomics.com/datasets/10k-human-pbmcs-atac-v2-chromium-x-2-standard
https://www.archrproject.com/bookdown/index.html
https://github.com/GreenleafLab/ArchR
