美文网首页
单细胞ATAC数据的两种读取方式(上)

单细胞ATAC数据的两种读取方式(上)

作者: Cherrieg | 来源:发表于2024-08-11 17:16 被阅读0次

单细胞研究使我们能够了解细胞发育的分子过程和病理的发生机理。单细胞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_matrixfiltered_peak_bc_matrix.h5fragments.tsv.gzsinglecell.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分钟,其中,对每一个样本的操作为:

  1. 从输入文件中读取开放fragment;
  2. 计算每个细胞的质控信息(即TSS富集分数和核小体信息);
  3. 根据质控参数过滤细胞;
  4. 使用500 bp bin创建全基因组TileMatrix;
  5. 使用调用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
)

这样,我们就使用getValidBarcodesvalidBarcodes参数结合,完成了对特定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

相关文章

网友评论

      本文标题:单细胞ATAC数据的两种读取方式(上)

      本文链接:https://www.haomeiwen.com/subject/mwhrkjtx.html