引言
在本教程中，我们将探讨由10x Genomics公司提供的成年小鼠大脑细胞的单细胞ATAC-seq数据集。本教程中使用的所有相关文件均可在10x Genomics官方网站上获取。
本教程复现了之前在人类外周血单核细胞（PBMC）的Signac入门教程中执行的命令。我们通过在不同的系统上进行相同的分析，来展示其性能以及对不同组织类型的适用性，并提供了一个来自不同物种的示例。
实战
首先，我们需要导入Signac、Seurat等一些用于分析小鼠数据的软件包。
library(Signac)
library(Seurat)
library(EnsDb.Mmusculus.v79)
library(ggplot2)
library(patchwork)
预处理工作流程
counts <- Read10X_h5("../vignette_data/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
file = "../vignette_data/atac_v1_adult_brain_fresh_5k_singlecell.csv",
header = TRUE,
row.names = 1
)
brain_assay <- CreateChromatinAssay(
counts = counts,
sep = c(":", "-"),
genome = "mm10",
fragments = '../vignette_data/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz',
min.cells = 1
)
brain <- CreateSeuratObject(
counts = brain_assay,
assay = 'peaks',
project = 'ATAC',
meta.data = metadata
)
我们还可以向小鼠基因组的大脑对象添加基因注释。这将允许下游函数直接从对象中提取基因注释信息。
# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79)
# change to UCSC style since the data was mapped to hg19
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "mm10"
# add the gene information to the object
Annotation(brain) <- annotations
计算 QC 指标
接下来我们计算一些有用的细胞 QC 指标。
brain <- NucleosomeSignal(object = brain)
我们可以分析所有细胞的DNA片段长度的周期性变化，并根据细胞核小体信号的强弱进行分类。观察结果表明，那些在单核小体与无核小体比例上表现异常的细胞，呈现出与其他细胞不同的条带图谱。而其他细胞则显示出了一次成功的ATAC-seq实验所特有的典型模式。
brain$nucleosome_group <- ifelse(brain$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = brain, group.by = 'nucleosome_group', region = 'chr1-1-10000000')
在ATAC-seq实验中，Tn5转座酶在转录起始位点（TSS）处的整合事件的富集程度，是一个关键的质量控制指标，用于评价Tn5的定位效率。ENCODE联盟将TSS富集分数定义为TSS周围Tn5整合位点的计数与这些位点在相邻区域计数的比率。在Signac软件包中，我们可以使用TSSEnrichment()函数来为每个细胞计算这一富集分数。
brain <- TSSEnrichment(brain, fast = FALSE)
brain$high.tss <- ifelse(brain$TSS.enrichment > 2, 'High', 'Low')
TSSPlot(brain, group.by = 'high.tss') + NoLegend()
brain$pct_reads_in_peaks <- brain$peak_region_fragments / brain$passed_filters * 100
brain$blacklist_ratio <- brain$blacklist_region_fragments / brain$peak_region_fragments
VlnPlot(
object = brain,
features = c('pct_reads_in_peaks', 'peak_region_fragments',
'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
pt.size = 0.1,
ncol = 5
)
我们删除了这些 QC 指标异常值的细胞。
brain <- subset(
x = brain,
subset = peak_region_fragments > 3000 &
peak_region_fragments < 100000 &
pct_reads_in_peaks > 40 &
blacklist_ratio < 0.025 &
nucleosome_signal < 4 &
TSS.enrichment > 2
)
brain
## An object of class Seurat
## 157203 features across 3512 samples within 1 assay
## Active assay: peaks (157203 features, 0 variable features)
## 2 layers present: counts, data
归一化和线性降维
brain <- RunTFIDF(brain)
brain <- FindTopFeatures(brain, min.cutoff = 'q0')
brain <- RunSVD(object = brain)
在分析中，LSI（线性判别分析）的第一个主成分往往反映的是测序的深度（即技术层面的变异），而非生物学上的变异。如果确实如此，那么在后续的分析中应该将这一成分排除掉。为了判断是否存在这种情况，我们可以通过调用DepthCor()函数来计算每个LSI主成分与测序深度之间的相关性。
DepthCor(brain)
在这里，我们看到第一个 LSI 组件与细胞的计数总数之间存在非常强的相关性，因此我们将在没有该组件的情况下执行下游步骤。
非线性降维和聚类
细胞数据已经被嵌入到一个低维度的空间里，我们可以采用单细胞RNA测序(scRNA-seq)数据常用的分析方法，执行基于图谱的聚类分析，并通过非线性降维技术来进行数据可视化。RunUMAP()、FindNeighbors()和FindClusters()这些功能均集成在Seurat软件包中。
brain <- RunUMAP(
object = brain,
reduction = 'lsi',
dims = 2:30
)
brain <- FindNeighbors(
object = brain,
reduction = 'lsi',
dims = 2:30
)
brain <- FindClusters(
object = brain,
algorithm = 3,
resolution = 1.2,
verbose = FALSE
)
DimPlot(object = brain, label = TRUE) + NoLegend()
本文由 mdnice 多平台发布
