admin管理员组文章数量:1032889
一个成年小鼠大脑的scATAC
前面在写ATACseq学习笔记的时候，后台有很多粉丝询问scATAC-seq有没有相关资料，现在就来学习，开个新专辑《scATAC-seq分析》~
scATAC-seq用的最多的包就是Signac包了，这个包来自Seurat包的扩展，分析步骤等跟Seurat非常像，下面来看看。
与Seurat软件一样，官方给了很多示例，前面我们学习了PBMC组织的scATAC-seq分析：一个PBMC的scATAC-seq基础分析：Signac
今天再来学习一个对小鼠大脑组织的scATAC-seq分析，Signac可以分析各种组织类型的数据：
安装
代码语言：javascript代码运行次数：0运行复制## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror=";)
options("repos"=c(CRAN="/"))
setRepositories(ind=1:3)
install.packages("Signac")
rm(list=ls())
## 加载包
library(Signac)
library(Seurat)
library(EnsDb.Mmusculus.v79)
library(ggplot2)
library(patchwork)
0.示例数据
这个数据来自10x Genomics，为组织类型为： mouse brain。
下载地址为：
代码语言：javascript代码运行次数：0运行复制# bash命令，或者自己选择用浏览器，迅雷都可以
wget .1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5
wget .1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_singlecell.csv
wget .1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz
wget .1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz.tbi
1.数据预处理
1.1 读取数据并创建seurat对象：
这里的代码跟前面差不多：
代码语言：javascript代码运行次数：0运行复制# 1.读取数据：
counts <- Read10X_h5(filename = "atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
file = "atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_singlecell.csv",
header = TRUE,
row.names = 1
)
# 简单的探索一下数据
head(metadata)
dim(counts)
# 每一行为一个peak区间，每一列为一个细胞
counts[1:5,1:5]
# 创建对象
brain_assay <- CreateChromatinAssay(
counts = counts,
sep = c(":", "-"),
genome = "mm10",
fragments = 'atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz',
min.cells = 1
)
# 生成seurat对象
brain <- CreateSeuratObject(
counts = brain_assay,
assay = 'peaks',
project = 'ATAC',
meta.data = metadata
)
brain
brain[['peaks']]
granges(brain)
总共5337个细胞：
代码语言：javascript代码运行次数：0运行复制An object of class Seurat
157203 features across 5337 samples within 1 assay
Active assay: peaks (157203 features, 0 variable features)
2 layers present: counts, data
ATACseq数据存储在ChromatinAssay里面：
brain[['peaks']]
# ChromatinAssay data with 157203 features for 5337 cells
# Variable features: 0
# Genome: mm10
# Annotation present: FALSE
# Motifs present: FALSE
# Fragment files: 1
# 查看基因组范围
granges(brain)
# GRanges object with 157203 ranges and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] chr1 3094708-3095552 *
# [2] chr1 3119414-3121782 *
# [3] chr1 3204809-3205178 *
# [4] chr1 3217330-3217359 *
# [5] chr1 3228123-3228463 *
# ... ... ... ...
# [157199] chrY 90761049-90762169 *
# [157200] chrY 90800417-90800831 *
# [157201] chrY 90804515-90805440 *
# [157202] chrY 90808545-90809111 *
# [157203] chrY 90810741-90810960 *
# -------
# seqinfo: 21 sequences from an unspecified genome; no seqlengths
1.2 添加基因信息
这里需要注意与上游分析使用的参考基因需要是同一个版本，这个数据的10x Genomics用的版本为：GRCh38-2020-A，对应Ensembl v98 patch release，他们的对应关系可以从这里获取：.html
#### 添加参考基因组信息
# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79)
head(annotations)
# 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
参考基因组信息：
2.计算QC指标
ATAC-seq测序有一些数据质量评估指标，下面来看看~
2.1 核小体信号模式
Nucleosome banding pattern：DNA片段大小的直方图应该显示出一个强烈的核小体条带模式，这对应于缠绕在单个核小体周围的DNA长度。我们针对每个单细胞计算这一模式，并定量估算单核小体片段与无核小体片段的近似比例（存储为
nucleosome_signal）；
我们可以观察所有细胞的片段长度周期性，并根据核小体信号强度的高低对细胞进行分组。可以看到，单核小体/无核小体比例异常的细胞具有不同的条带模式。其余细胞表现出典型的成功ATAC-seq实验的模式。
代码语言：javascript代码运行次数：0运行复制## 2.计算QC指标
brain <- NucleosomeSignal(object = brain)
# 分组
brain$nucleosome_group <- ifelse(brain$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
# 可视化
FragmentHistogram(object = brain, group.by = 'nucleosome_group', region = 'chr1-1-10000000')
2.2 TSS富集打分
Tn5整合事件在转录起始位点（TSS）的富集也可以作为评估ATAC-seq实验中Tn5靶向性的一个重要质量控制指标。ENCODE联盟定义了一个TSS富集评分，即TSS周围的Tn5整合位点数量与侧翼区域的Tn5整合位点数量的比值。关于TSS富集评分的更多信息，请参阅ENCODE文档：/。我们可以使用Signac包中的TSSEnrichment()函数为每个细胞计算TSS富集评分。
# TSS富集打分
brain <- TSSEnrichment(brain)
brain$pct_reads_in_peaks <- brain$peak_region_fragments / brain$passed_filters * 100
brain$blacklist_ratio <- brain$blacklist_region_fragments / brain$peak_region_fragments
head(brain@meta.data)
VlnPlot(
object = brain,
features = c('pct_reads_in_peaks', 'peak_region_fragments',
'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
pt.size = 0.1,
ncol = 5
)
根据上面的数据分布，过滤条件如下，如果是自己的数据需要根据项目背景进行过滤cutoff调整：
代码语言：javascript代码运行次数：0运行复制# 过滤
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
还剩下 3512 多个细胞。
3.标准化和线性降维
分为三个小步骤：
- 标准化：Signac进行
TF-IDF标准化 - 特征选择：scATAC-seq数据的低动态范围使得我们难以像处理scRNA-seq数据那样进行可变特征选择。但可以选择仅使用排名top n%的特征（peaks）来进行降维，或者使用
FindTopFeatures()函数移除出现在少于n个细胞中的特征 - 降维：在上述选择的特征（峰值）所构成的TF-IDF矩阵上运行奇异值分解（SVD）
## 3.标准化，线性降维
brain <- RunTFIDF(brain)
brain <- FindTopFeatures(brain, min.cutoff = 'q0')
brain <- RunSVD(object = brain)
DepthCor(brain)
第一个 LSI 成分通常捕获测序深度（技术变异）而非生物学变异。如果情况如此，那么该成分应从下游分析中移除。可以使用DepthCor()函数来评估每个LSI成分与测序深度之间的相关性：
在这里，我们看到第一个 LSI 成分与细胞的总计数数量之间存在非常强的相关性（图中看到是就是前10个LSI，除了第一个基本上都在0.3以内）。我们将在后续步骤中不使用这个成分，dim为 2:30。
4.非线性标准化和降维
现在细胞已经被嵌入到低维空间中，我们可以使用通常用于scRNA-seq数据分析的方法来进行基于图的聚类和非线性降维以用于可视化。RunUMAP()、FindNeighbors()和FindClusters()这些函数都来自Seurat包。
## 4.非线性标准化聚类
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()
聚类结果如下：
5.构建基因活性矩阵
通过评估与基因相关的染色质可及性来量化基因组中每个基因的活性，并从scATAC-seq数据中创建一个新的基因活性检测方法。这些步骤可以通过GeneActivity()函数自动完成：
## 5.构建基因活性矩阵
# compute gene activities
gene.activities <- GeneActivity(brain)
dim(gene.activities)
gene.activities[1:5,1:5]
# add the gene activity matrix to the Seurat object as a new assay
brain[['RNA']] <- CreateAssayObject(counts = gene.activities)
brain <- NormalizeData(
object = brain,
assay = 'RNA',
normalization.method = 'LogNormalize',
scale.factor = median(brain$nCount_RNA)
)
现在可以可视化典型标记基因的活性，以帮助解释ATAC-seq聚类结果。以下是这些基因在小鼠脑部组织中的细胞类型及详细解释的表格：
基因名 | 细胞类型（中文/英文） | 详细解释 |
|---|---|---|
Sst | 胆囊收缩素神经元/Somatostatin neuron | Sst基因标记的神经元是一种抑制性神经元，在小鼠大脑皮层中存在多种亚型，包括建立长程投射的长投射神经元和建立局部投射的中间神经元。 |
Pvalb | 副钙蛋白神经元/Parvalbumin neuron | Pvalb基因标记的神经元是一种快速放电的抑制性中间神经元，参与调节神经网络的同步活动。 |
Gad2 | 谷氨酸脱羧酶神经元/Glutamate Decarboxylase neuron | Gad2基因编码的谷氨酸脱羧酶是合成GABA的关键酶，Gad2标记的神经元是抑制性神经元。 |
Neurod6 | 神经分化因子6神经元/Neurod6 neuron | Neurod6基因在神经分化过程中起重要作用，其标记的神经元可能与神经系统的发育和功能有关。 |
Rorb | 视网膜孤儿受体β神经元/Retinoid-related orphan receptor beta neuron | Rorb基因在视网膜神经节细胞和某些中间神经元中表达，参与调节神经元的发育和功能。 |
Syt6 | 突触融合蛋白6神经元/Synaptotagmin 6 neuron | Syt6基因编码的突触融合蛋白6参与神经递质的释放过程，其标记的神经元可能与神经信号的传递有关。 |
需要注意的是，部分基因的细胞类型解释是基于其在神经系统中的已知功能和表达模式，具体的细胞类型可能因研究方法和实验条件的不同而有所差异。
代码语言：javascript代码运行次数：0运行复制# 可视化
DefaultAssay(brain) <- 'RNA'
FeaturePlot(
object = brain,
features = c('Sst','Pvalb',"Gad2","Neurod6","Rorb","Syt6"),
pt.size = 0.1,
max.cutoff = 'q95',
ncol = 3
)
6.与单细胞数据进行整合
为了帮助解释scATAC-seq数据，我们可以基于同一生物系统（成年小鼠大脑）的scRNA-seq实验来对细胞进行分类。利用跨模态整合和标签转移的方法，旨在识别基因活性矩阵和单细胞RNA-seq数据集中共享的相关性模式，以识别两种模态之间的匹配生物学状态。这一过程将为每个细胞返回每个基于单细胞RNA-seq定义的聚类标签的分类得分。
可以从艾伦研究所网站下载该实验的原始数据：，并在GitHub上查看用于构建该对象的代码。或者，也可以在此处下载预处理过的Seurat对象：.rds。
代码语言：javascript代码运行次数：0运行复制## 6.与单细胞数据整合
# Load the pre-processed scRNA-seq data
allen_rna <- readRDS("atac_v1_adult_brain_fresh_5k/allen_brain.rds")
allen_rna <- UpdateSeuratObject(allen_rna)
allen_rna <- FindVariableFeatures( object = allen_rna, nfeatures = 5000 )
Idents(allen_rna) <- "subclass"
DimPlot(object = allen_rna, label = TRUE) + NoLegend()
已经注释好的单细胞：
找到 anchors：
代码语言：javascript代码运行次数：0运行复制transfer.anchors <- FindTransferAnchors(
reference = pbmc_rna,
query = pbmc,
reduction = 'cca'
)
transfer.anchors
predicted.labels <- TransferData(
anchorset = transfer.anchors,
refdata = pbmc_rna$celltype,
weight.reduction = pbmc[['lsi']],
dims = 2:30
)
可视化预测结果：
代码语言：javascript代码运行次数：0运行复制## 添加到对象中
brain <- AddMetaData(object = brain, metadata = predicted.labels)
plot1 <- DimPlot(allen_rna, group.by = 'subclass', label = TRUE, repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')
plot2 <- DimPlot(brain, group.by = 'predicted.id', label = TRUE, repel = TRUE) + NoLegend() + ggtitle('scATAC-seq')
plot1 + plot2
现在这个图就是很多文献中scRNA-seq+scATAC-seq联合分析的结果啦：
7.细胞类型差异peaks分析
在这里，我们寻找皮层不同层之间兴奋性神经元的差异可及性区域：
代码语言：javascript代码运行次数：0运行复制## 7.差异peaks分析
# switch back to working with peaks instead of gene activities
DefaultAssay(brain) <- 'peaks'
Idents(brain) <- "predicted.id"
da_peaks <- FindMarkers(
object = brain,
ident.1 = c("L2/3 IT"),
ident.2 = c("L4", "L5 IT", "L6 IT"),
test.use = 'wilcox',
min.pct = 0.1
)
head(da_peaks)
选一个peaks可视化：
代码语言：javascript代码运行次数：0运行复制# 可视化差异peaks
rownames(da_peaks)[1]
plot1 <- VlnPlot(
object = brain,
features = rownames(da_peaks)[1],
pt.size = 0.1,
idents = c("L4","L5 IT","L2/3 IT")
)
plot2 <- FeaturePlot(
object = brain,
features = rownames(da_peaks)[1],
pt.size = 0.1,
max.cutoff = 'q95'
)
plot1 | plot2
对peaks进行注释，与基因关联：
代码语言：javascript代码运行次数：0运行复制# 对peaks进行注释，与基因关联：
open_l23 <- rownames(da_peaks[da_peaks$avg_log2FC > 3, ])
open_l456 <- rownames(da_peaks[da_peaks$avg_log2FC < 3, ])
closest_l23 <- ClosestFeature(brain, open_l23)
closest_l456 <- ClosestFeature(brain, open_l456)
head(closest_l23)
8.绘制 genomic regions
我们还可以使用CoveragePlot()函数，针对对象中存储的任何基因组区域，按簇、细胞类型或其他元数据对细胞进行分组，绘制覆盖度图。这些图表示伪批量可及性轨迹，其中来自一个组内所有细胞的信号被平均在一起，以可视化一个区域内的DNA可及性。
# show cell types with at least 50 cells
idents.plot <- names(which(table(Idents(brain)) > 50))
brain <- SortIdents(brain)
CoveragePlot(
object = brain,
region = c("Neurod6", "Gad2"),
idents = idents.plot,
extend.upstream = 1000,
extend.downstream = 1000,
ncol = 1
)
到此为止，我们已经学习了两个组织的scATAC-Seq基础分析，对这个流程也有了一定的认知，下期我们来寻找一些文献看看常规的个性化~
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-04-02，如有侵权请联系 cloudcommunity@tencent 删除数据object对象函数可视化一个成年小鼠大脑的scATAC
前面在写ATACseq学习笔记的时候，后台有很多粉丝询问scATAC-seq有没有相关资料，现在就来学习，开个新专辑《scATAC-seq分析》~
scATAC-seq用的最多的包就是Signac包了，这个包来自Seurat包的扩展，分析步骤等跟Seurat非常像，下面来看看。
与Seurat软件一样，官方给了很多示例，前面我们学习了PBMC组织的scATAC-seq分析：一个PBMC的scATAC-seq基础分析：Signac
今天再来学习一个对小鼠大脑组织的scATAC-seq分析，Signac可以分析各种组织类型的数据：
安装
代码语言：javascript代码运行次数：0运行复制## 使用西湖大学的 Bioconductor镜像
options(BioC_mirror=";)
options("repos"=c(CRAN="/"))
setRepositories(ind=1:3)
install.packages("Signac")
rm(list=ls())
## 加载包
library(Signac)
library(Seurat)
library(EnsDb.Mmusculus.v79)
library(ggplot2)
library(patchwork)
0.示例数据
这个数据来自10x Genomics，为组织类型为： mouse brain。
下载地址为：
代码语言：javascript代码运行次数：0运行复制# bash命令，或者自己选择用浏览器，迅雷都可以
wget .1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5
wget .1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_singlecell.csv
wget .1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz
wget .1.0/atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz.tbi
1.数据预处理
1.1 读取数据并创建seurat对象：
这里的代码跟前面差不多：
代码语言：javascript代码运行次数：0运行复制# 1.读取数据：
counts <- Read10X_h5(filename = "atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
file = "atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_singlecell.csv",
header = TRUE,
row.names = 1
)
# 简单的探索一下数据
head(metadata)
dim(counts)
# 每一行为一个peak区间，每一列为一个细胞
counts[1:5,1:5]
# 创建对象
brain_assay <- CreateChromatinAssay(
counts = counts,
sep = c(":", "-"),
genome = "mm10",
fragments = 'atac_v1_adult_brain_fresh_5k/atac_v1_adult_brain_fresh_5k_fragments.tsv.gz',
min.cells = 1
)
# 生成seurat对象
brain <- CreateSeuratObject(
counts = brain_assay,
assay = 'peaks',
project = 'ATAC',
meta.data = metadata
)
brain
brain[['peaks']]
granges(brain)
总共5337个细胞：
代码语言：javascript代码运行次数：0运行复制An object of class Seurat
157203 features across 5337 samples within 1 assay
Active assay: peaks (157203 features, 0 variable features)
2 layers present: counts, data
ATACseq数据存储在ChromatinAssay里面：
brain[['peaks']]
# ChromatinAssay data with 157203 features for 5337 cells
# Variable features: 0
# Genome: mm10
# Annotation present: FALSE
# Motifs present: FALSE
# Fragment files: 1
# 查看基因组范围
granges(brain)
# GRanges object with 157203 ranges and 0 metadata columns:
# seqnames ranges strand
# <Rle> <IRanges> <Rle>
# [1] chr1 3094708-3095552 *
# [2] chr1 3119414-3121782 *
# [3] chr1 3204809-3205178 *
# [4] chr1 3217330-3217359 *
# [5] chr1 3228123-3228463 *
# ... ... ... ...
# [157199] chrY 90761049-90762169 *
# [157200] chrY 90800417-90800831 *
# [157201] chrY 90804515-90805440 *
# [157202] chrY 90808545-90809111 *
# [157203] chrY 90810741-90810960 *
# -------
# seqinfo: 21 sequences from an unspecified genome; no seqlengths
1.2 添加基因信息
这里需要注意与上游分析使用的参考基因需要是同一个版本，这个数据的10x Genomics用的版本为：GRCh38-2020-A，对应Ensembl v98 patch release，他们的对应关系可以从这里获取：.html
#### 添加参考基因组信息
# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Mmusculus.v79)
head(annotations)
# 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
参考基因组信息：
2.计算QC指标
ATAC-seq测序有一些数据质量评估指标，下面来看看~
2.1 核小体信号模式
Nucleosome banding pattern：DNA片段大小的直方图应该显示出一个强烈的核小体条带模式，这对应于缠绕在单个核小体周围的DNA长度。我们针对每个单细胞计算这一模式，并定量估算单核小体片段与无核小体片段的近似比例（存储为
nucleosome_signal）；
我们可以观察所有细胞的片段长度周期性，并根据核小体信号强度的高低对细胞进行分组。可以看到，单核小体/无核小体比例异常的细胞具有不同的条带模式。其余细胞表现出典型的成功ATAC-seq实验的模式。
代码语言：javascript代码运行次数：0运行复制## 2.计算QC指标
brain <- NucleosomeSignal(object = brain)
# 分组
brain$nucleosome_group <- ifelse(brain$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
# 可视化
FragmentHistogram(object = brain, group.by = 'nucleosome_group', region = 'chr1-1-10000000')
2.2 TSS富集打分
Tn5整合事件在转录起始位点（TSS）的富集也可以作为评估ATAC-seq实验中Tn5靶向性的一个重要质量控制指标。ENCODE联盟定义了一个TSS富集评分，即TSS周围的Tn5整合位点数量与侧翼区域的Tn5整合位点数量的比值。关于TSS富集评分的更多信息，请参阅ENCODE文档：/。我们可以使用Signac包中的TSSEnrichment()函数为每个细胞计算TSS富集评分。
# TSS富集打分
brain <- TSSEnrichment(brain)
brain$pct_reads_in_peaks <- brain$peak_region_fragments / brain$passed_filters * 100
brain$blacklist_ratio <- brain$blacklist_region_fragments / brain$peak_region_fragments
head(brain@meta.data)
VlnPlot(
object = brain,
features = c('pct_reads_in_peaks', 'peak_region_fragments',
'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
pt.size = 0.1,
ncol = 5
)
根据上面的数据分布，过滤条件如下，如果是自己的数据需要根据项目背景进行过滤cutoff调整：
代码语言：javascript代码运行次数：0运行复制# 过滤
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
还剩下 3512 多个细胞。
3.标准化和线性降维
分为三个小步骤：
- 标准化：Signac进行
TF-IDF标准化 - 特征选择：scATAC-seq数据的低动态范围使得我们难以像处理scRNA-seq数据那样进行可变特征选择。但可以选择仅使用排名top n%的特征（peaks）来进行降维，或者使用
FindTopFeatures()函数移除出现在少于n个细胞中的特征 - 降维：在上述选择的特征（峰值）所构成的TF-IDF矩阵上运行奇异值分解（SVD）
## 3.标准化，线性降维
brain <- RunTFIDF(brain)
brain <- FindTopFeatures(brain, min.cutoff = 'q0')
brain <- RunSVD(object = brain)
DepthCor(brain)
第一个 LSI 成分通常捕获测序深度（技术变异）而非生物学变异。如果情况如此，那么该成分应从下游分析中移除。可以使用DepthCor()函数来评估每个LSI成分与测序深度之间的相关性：
在这里，我们看到第一个 LSI 成分与细胞的总计数数量之间存在非常强的相关性（图中看到是就是前10个LSI，除了第一个基本上都在0.3以内）。我们将在后续步骤中不使用这个成分，dim为 2:30。
4.非线性标准化和降维
现在细胞已经被嵌入到低维空间中，我们可以使用通常用于scRNA-seq数据分析的方法来进行基于图的聚类和非线性降维以用于可视化。RunUMAP()、FindNeighbors()和FindClusters()这些函数都来自Seurat包。
## 4.非线性标准化聚类
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()
聚类结果如下：
5.构建基因活性矩阵
通过评估与基因相关的染色质可及性来量化基因组中每个基因的活性，并从scATAC-seq数据中创建一个新的基因活性检测方法。这些步骤可以通过GeneActivity()函数自动完成：
## 5.构建基因活性矩阵
# compute gene activities
gene.activities <- GeneActivity(brain)
dim(gene.activities)
gene.activities[1:5,1:5]
# add the gene activity matrix to the Seurat object as a new assay
brain[['RNA']] <- CreateAssayObject(counts = gene.activities)
brain <- NormalizeData(
object = brain,
assay = 'RNA',
normalization.method = 'LogNormalize',
scale.factor = median(brain$nCount_RNA)
)
现在可以可视化典型标记基因的活性，以帮助解释ATAC-seq聚类结果。以下是这些基因在小鼠脑部组织中的细胞类型及详细解释的表格：
基因名 | 细胞类型（中文/英文） | 详细解释 |
|---|---|---|
Sst | 胆囊收缩素神经元/Somatostatin neuron | Sst基因标记的神经元是一种抑制性神经元，在小鼠大脑皮层中存在多种亚型，包括建立长程投射的长投射神经元和建立局部投射的中间神经元。 |
Pvalb | 副钙蛋白神经元/Parvalbumin neuron | Pvalb基因标记的神经元是一种快速放电的抑制性中间神经元，参与调节神经网络的同步活动。 |
Gad2 | 谷氨酸脱羧酶神经元/Glutamate Decarboxylase neuron | Gad2基因编码的谷氨酸脱羧酶是合成GABA的关键酶，Gad2标记的神经元是抑制性神经元。 |
Neurod6 | 神经分化因子6神经元/Neurod6 neuron | Neurod6基因在神经分化过程中起重要作用，其标记的神经元可能与神经系统的发育和功能有关。 |
Rorb | 视网膜孤儿受体β神经元/Retinoid-related orphan receptor beta neuron | Rorb基因在视网膜神经节细胞和某些中间神经元中表达，参与调节神经元的发育和功能。 |
Syt6 | 突触融合蛋白6神经元/Synaptotagmin 6 neuron | Syt6基因编码的突触融合蛋白6参与神经递质的释放过程，其标记的神经元可能与神经信号的传递有关。 |
需要注意的是，部分基因的细胞类型解释是基于其在神经系统中的已知功能和表达模式，具体的细胞类型可能因研究方法和实验条件的不同而有所差异。
代码语言：javascript代码运行次数：0运行复制# 可视化
DefaultAssay(brain) <- 'RNA'
FeaturePlot(
object = brain,
features = c('Sst','Pvalb',"Gad2","Neurod6","Rorb","Syt6"),
pt.size = 0.1,
max.cutoff = 'q95',
ncol = 3
)
6.与单细胞数据进行整合
为了帮助解释scATAC-seq数据，我们可以基于同一生物系统（成年小鼠大脑）的scRNA-seq实验来对细胞进行分类。利用跨模态整合和标签转移的方法，旨在识别基因活性矩阵和单细胞RNA-seq数据集中共享的相关性模式，以识别两种模态之间的匹配生物学状态。这一过程将为每个细胞返回每个基于单细胞RNA-seq定义的聚类标签的分类得分。
可以从艾伦研究所网站下载该实验的原始数据：，并在GitHub上查看用于构建该对象的代码。或者，也可以在此处下载预处理过的Seurat对象：.rds。
代码语言：javascript代码运行次数：0运行复制## 6.与单细胞数据整合
# Load the pre-processed scRNA-seq data
allen_rna <- readRDS("atac_v1_adult_brain_fresh_5k/allen_brain.rds")
allen_rna <- UpdateSeuratObject(allen_rna)
allen_rna <- FindVariableFeatures( object = allen_rna, nfeatures = 5000 )
Idents(allen_rna) <- "subclass"
DimPlot(object = allen_rna, label = TRUE) + NoLegend()
已经注释好的单细胞：
找到 anchors：
代码语言：javascript代码运行次数：0运行复制transfer.anchors <- FindTransferAnchors(
reference = pbmc_rna,
query = pbmc,
reduction = 'cca'
)
transfer.anchors
predicted.labels <- TransferData(
anchorset = transfer.anchors,
refdata = pbmc_rna$celltype,
weight.reduction = pbmc[['lsi']],
dims = 2:30
)
可视化预测结果：
代码语言：javascript代码运行次数：0运行复制## 添加到对象中
brain <- AddMetaData(object = brain, metadata = predicted.labels)
plot1 <- DimPlot(allen_rna, group.by = 'subclass', label = TRUE, repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')
plot2 <- DimPlot(brain, group.by = 'predicted.id', label = TRUE, repel = TRUE) + NoLegend() + ggtitle('scATAC-seq')
plot1 + plot2
现在这个图就是很多文献中scRNA-seq+scATAC-seq联合分析的结果啦：
7.细胞类型差异peaks分析
在这里，我们寻找皮层不同层之间兴奋性神经元的差异可及性区域：
代码语言：javascript代码运行次数：0运行复制## 7.差异peaks分析
# switch back to working with peaks instead of gene activities
DefaultAssay(brain) <- 'peaks'
Idents(brain) <- "predicted.id"
da_peaks <- FindMarkers(
object = brain,
ident.1 = c("L2/3 IT"),
ident.2 = c("L4", "L5 IT", "L6 IT"),
test.use = 'wilcox',
min.pct = 0.1
)
head(da_peaks)
选一个peaks可视化：
代码语言：javascript代码运行次数：0运行复制# 可视化差异peaks
rownames(da_peaks)[1]
plot1 <- VlnPlot(
object = brain,
features = rownames(da_peaks)[1],
pt.size = 0.1,
idents = c("L4","L5 IT","L2/3 IT")
)
plot2 <- FeaturePlot(
object = brain,
features = rownames(da_peaks)[1],
pt.size = 0.1,
max.cutoff = 'q95'
)
plot1 | plot2
对peaks进行注释，与基因关联：
代码语言：javascript代码运行次数：0运行复制# 对peaks进行注释，与基因关联：
open_l23 <- rownames(da_peaks[da_peaks$avg_log2FC > 3, ])
open_l456 <- rownames(da_peaks[da_peaks$avg_log2FC < 3, ])
closest_l23 <- ClosestFeature(brain, open_l23)
closest_l456 <- ClosestFeature(brain, open_l456)
head(closest_l23)
8.绘制 genomic regions
我们还可以使用CoveragePlot()函数，针对对象中存储的任何基因组区域，按簇、细胞类型或其他元数据对细胞进行分组，绘制覆盖度图。这些图表示伪批量可及性轨迹，其中来自一个组内所有细胞的信号被平均在一起，以可视化一个区域内的DNA可及性。
# show cell types with at least 50 cells
idents.plot <- names(which(table(Idents(brain)) > 50))
brain <- SortIdents(brain)
CoveragePlot(
object = brain,
region = c("Neurod6", "Gad2"),
idents = idents.plot,
extend.upstream = 1000,
extend.downstream = 1000,
ncol = 1
)
到此为止，我们已经学习了两个组织的scATAC-Seq基础分析，对这个流程也有了一定的认知，下期我们来寻找一些文献看看常规的个性化~
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-04-02，如有侵权请联系 cloudcommunity@tencent 删除数据object对象函数可视化本文标签： 一个成年小鼠大脑的scATAC
版权声明：本文标题：一个成年小鼠大脑的scATAC 内容由热心网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://it.en369.cn/jiaocheng/1747984333a2237480.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论