单细胞分析(Signac)： PBMC scATAC-seq 预处理

引言

在本教学指南中，我们将探讨由10x Genomics公司提供的人类外周血单核细胞（PBMCs）的单细胞ATAC-seq数据集。本节中所涉及的所有文件均可在10x Genomics的官方网站上找到：

• 原始数据文件
• 相关的元数据文件
• 包含数据片段的文件
• 用于索引的片段文件

要下载所有必需的文件，您可以在 shell 中运行以下行：

wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz.tbi

加载包

首先加载 Signac、Seurat 和我们将用于分析人类数据的其他一些包。

if (!requireNamespace("EnsDb.Hsapiens.v75", quietly = TRUE))
    BiocManager::install("EnsDb.Hsapiens.v75")

library(Signac)
library(Seurat)
library(EnsDb.Hsapiens.v75)

library(ggplot2)
library(patchwork)

预处理工作流程

在进行染色质数据的预处理过程中，Signac 利用了两个相关输入文件的信息，这两个文件都可以通过 CellRanger 软件生成：

• 峰值/细胞矩阵。这与用于分析单细胞 RNA-seq 技术的基因表达计数矩阵相似。不同之处在于，矩阵中的每一行不是代表一个基因，而是代表基因组中的一个特定区域（称为峰值），这个区域被预测为开放染色质的区域。矩阵中的每个数值表示每个独特的条形码（也就是每个细胞）内 Tn5 整合位点的数量，这些位点都映射到特定的峰值区域内。更多详细信息可以在 10X 官网上查阅。
• 片段文件。这个文件包含了所有单细胞中所有独特片段的详尽清单。它的体积较大，处理速度较慢，并且是存储在硬盘上而非内存中。尽管如此，保留这个文件的好处在于它包含了与每个单细胞相关的所有片段信息，而不仅仅是那些映射到峰值的片段。关于片段文件的更多信息，可以在 10x Genomics 官网或 sinto 网站上找到。

我们首先利用 cellranger-atac 生成的峰值/细胞矩阵和细胞元数据来创建一个 Seurat 对象，并将硬盘上的片段文件路径信息存储在这个 Seurat 对象中：

counts <- Read10X_h5(filename = "../vignette_data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
  file = "../vignette_data/atac_v1_pbmc_10k_singlecell.csv",
  header = TRUE,
  row.names = 1
)

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

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

pbmc

## An object of class Seurat 
## 87561 features across 8728 samples within 1 assay 
## Active assay: peaks (87561 features, 0 variable features)
##  2 layers present: counts, data

• 如果我没有 H5 文件怎么办？

即便您手头没有 .h5 格式的文件，分析工作也照样可以进行。您拥有的数据可能以三个文件的形式组织：一个计数文件（.mtx 格式）、一个细胞条形码文件以及一个峰值文件。在这种情况下，您可以通过调用 Matrix::readMM() 函数来导入这些数据：

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

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

另外，如果您仅有一个片段文件，您可以利用 FeatureMatrix() 函数来生成一个计数矩阵：

fragpath <- './data/atac_v1_pbmc_10k_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-seq 数据通过一种特定的分析方法——ChromatinAssay——进行存储。这种方法支持对单细胞基因组分析（例如单细胞ATAC-seq）进行特殊功能的分析。通过展示 ChromatinAssay 的内容，我们可以获取包括基序信息、基因注释和基因组数据在内的多种补充信息。

pbmc[['peaks']]

## ChromatinAssay data with 87561 features for 8728 cells
## Variable features: 0 
## Genome: 
## Annotation present: FALSE 
## Motifs present: FALSE 
## Fragment files: 1

举例来说，我们可以对一个将 ChromatinAssay 设为当前活动分析方法的 Seurat 对象执行 granges 函数（或者直接在 ChromatinAssay 对象上执行），以此来查看该对象中每个特征所对应的基因组区域。若要进一步了解 ChromatinAssay 类的信息，请参阅关于对象交互的相关文档。

granges(pbmc)

## GRanges object with 87561 ranges and 0 metadata columns:
##           seqnames            ranges strand
##                         
##       [1]     chr1     565107-565550      *
##       [2]     chr1     569174-569639      *
##       [3]     chr1     713460-714823      *
##       [4]     chr1     752422-753038      *
##       [5]     chr1     762106-763359      *
##       ...      ...               ...    ...
##   [87557]     chrY 58993392-58993760      *
##   [87558]     chrY 58994571-58994823      *
##   [87559]     chrY 58996352-58997331      *
##   [87560]     chrY 59001782-59002175      *
##   [87561]     chrY 59017143-59017246      *
##   -------
##   seqinfo: 24 sequences from an unspecified genome; no seqlengths

此外，我们可以为人类基因组的 pbmc 对象增加基因注释信息。这样一来，后续的分析函数就能够直接从该对象获取所需的基因注释数据。

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

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

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