在R语言中的 ATACseq 数据分析全流程实战（五）：peaks质控 & peaks 注释

本帖子学习资源：/

续上前面的实战：

在R语言中的 ATACseq 数据分析全流程实战（一）

在R语言中的 ATACseq 数据分析全流程实战（二）

ATACseq 数据分析全流程实战（三）：课后习题

在R语言中的 ATACseq 数据分析全流程实战（四）：Peak calling

续上面第四期~

前面在做Peak calling的时候，因为教程中主要是关注nucleosome free regions，用的bam文件为ATAC_50K_2_sorted_nucFreeRegions.bam，但是一直报错。

A common goal in ATACseq is to identify nucleosome free regions where transcription factors are binding and/or transcriptional machinery is active. This nucleosome free signal would correspond to fragments less than one nucleosome (as defined in Greenleaf paper < 100bp)

后面我改成了用全部的bam文件ATAC_50K_2_sorted.bam，先运行全部的。

代码语言：javascript代码运行次数：0运行复制

# -t:A BAM file to find enriched regions in
# –name:A Name for peak calls
# –outdir: An output folder to write peaks into
macs2 callpeak  -t GSE47753/ATAC_50K_2_sorted.bam \
                --outdir GSE47753/ \
                --name Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak \
                -f BAMPE -g hs

得到的结果如下：

• Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak_peaks.narrowPeak
• Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak_peaks.xls
• Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak_summits.bed

peaks QC

peaks鉴定出来后，需要看一下peaks结果的质量，如 low quality, duplicates and signal distribution，使用的R包为ChIPQC。

关于blacklists的概念见：.html#24，

即测序数据中经常会存在一些common artifacts，如超高信号的区域，这些区域会混淆peak识别，fragment长度估计以及QC metrics计算。因此Anshul Kundaje等定义了一个DAC blacklist 作为reference来帮助处理这些区域。

hg19的blacklist下载地址：/，

选择其中的：/@@download/ENCFF001TDO.bed.gz

我们可以使用ChIPQC R包来对ATQCseq数据来计算一些重要的指标如peaks的read数、blacklist中的reads，duplicated reads

代码语言：javascript代码运行次数：0运行复制

rm(list=ls())
# BiocManager::install("Herper")
library(Herper)
library(ChIPQC)
library(rtracklayer)
library(DT)
library(dplyr)
library(tidyr)

blkList <- import.bed("./ENCFF001TDO.bed.gz")
blkList
openRegionPeaks <- "GSE47753/Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak_peaks.narrowPeak"
qcRes <- ChIPQCsample("GSE47753/ATAC_50K_2_sorted.bam",
                      peaks = openRegionPeaks, 
                      annotation = "hg19", 
                      chromosomes = "chr20", 
                      blacklist = blkList,
                      verboseT = FALSE)


myMetrics <- QCmetrics(qcRes)
myMetrics[c("RiBL%", "RiP%")]
# RiBL%   RiP% 
# 0.573 14.700 

flgCounts <- flagtagcounts(qcRes)
DupRate <- flgCounts["DuplicateByChIPQC"]/flgCounts["Mapped"]
DupRate * 100
# 
# DuplicateByChIPQC 
# 20.49508 

去除blacklisted peaks

来自测序中的Artifacts，不完整的参考基因组会混淆peaks鉴定结果， blacklist建议在做了QC之后进行去除：

代码语言：javascript代码运行次数：0运行复制

## 去除
MacsCalls <- granges(qcRes[seqnames(qcRes) %in% "chr20"])
data.frame(Blacklisted = sum(MacsCalls %over% blkList), 
           Not_Blacklisted = sum(!MacsCalls %over% blkList))
#    Blacklisted Not_Blacklisted
# 1          11             738

MacsCalls <- MacsCalls[!MacsCalls %over% blkList]

peaks注释

peaks得到的结果一般是dna上的一个区间，我们希望知道这个区间有什么功能，处于哪个基因附近以便了解这些区域所发挥的功能。

注释的方法一般为将这些区域与他们最近的基因进行关联，或者是否在基因的转录起始位点。

使用的R包为chipseeker

代码语言：javascript代码运行次数：0运行复制

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)

## 注释
MacsCalls_Anno <- annotatePeak(MacsCalls, TxDb = TxDb.Hsapiens.UCSC.hg19.knownGene)

# 可视化注释结果
plotAnnoPie(MacsCalls_Anno)

# 提取 TSS regions (+/- 500)
MacsGR_Anno <- as.GRanges(MacsCalls_Anno)
MacsGR_TSS <- MacsGR_Anno[abs(MacsGR_Anno$distanceToTSS) < 500]
MacsGR_TSS[1, ]

peaks功能分析

这里教程用了一个没咋见过的R包，看看rGREAT

代码语言：javascript代码运行次数：0运行复制

## 功能分析
library(rGREAT)

great_Job <- submitGreatJob(MacsCalls, species = "hg19")
availableCategories(great_Job)

# GO 数据库结果
great_ResultTable = getEnrichmentTables(great_Job, category = "GO")
names(great_ResultTable)
# [1] "GO Molecular Function" "GO Biological Process" "GO Cellular Component"

temp <- great_ResultTable[["GO Biological Process"]][1:4, ]

peak差异分析

接下来，我们希望筛选到在不同分组样本中，peaks发生显著变化的区域，这里做差异分析需要用到多个样本，而从开始到现在我们只做了一个样本，所以这里我先再去找一些数据再来完成这部分好了~

下面来补充一些我在学习这个教程中的背景知识。

补充一些背景知识

了解实验背景至关重要，这可以很直观的告诉我们分析的具体数据是什么，下面是ATACseq的实验原理图，绿色的部分表示Tn5 transposase，红色和蓝色为两端的接头序列，核小体为图中的灰色部分。Tn5转座酶只插入到开放染色质的区域，所以测到的序列就是不同的核小体之间打开的那部分DNA序列。

nucleosomes

核小体(nucleosome )是组成真核细胞染色质的基本结构单位，由组蛋白和约200bp的DNA组成，形状为直径约10nm的球形小体。每个核小体包括有组蛋白八聚体及1分子组蛋白H1。

作为真核染色体DNA包装的第一级，核小体可使DNA长度压缩大约6倍。核小体之间通过连接DNA相连，在电镜下可观察到其以串珠的方式组织形成10nm纤丝，后者可进一步通过多级包装最终形成高度浓缩的染色体。相反，基因表达活跃部位以及处于复制或重组过程的区域则常伴随核小体结构的重塑或解聚，转而由多种非组蛋白结合，共同调控重要的生理过程（ATACseq测到的就是这部分序列）。

具体的组蛋白：

组装DNA过程：

ref：/

open chromatin

DNA的复制，基因的转录，是需要将DNA的高级结构解开的。这就类似于，要从压缩文件中获得信息，首先需要解压缩。

但是，这里的解压缩并不需要将整个zip包全部打开，而只需要打开一部分，即需要表达基因的区域即可。而这一个过程，主要由染色体组蛋白的修饰（尤其是乙酰化）来实现的。这部分打开的染色质，就叫开放染色质（open chromatin）。

致密的核小体结构被破坏后，这段启动子、增强子、绝缘子、沉默子等顺式调控元件和反式作用因子可以接近的区域，就是开放染色质。

而染色质一旦打开，就允许一些调控蛋白（比如转录因子和辅因子）跑过来与之相结合。而染色质的这种特性，就叫做染色质的可及性（chromatin accessibility），也理解为染色质开放性。

开放染色质与真核生物的转录调控密切相关。

那么常见的染色质开放区有哪些呢？

常见的染色质开放区主要是基因上游的启动子和远端的调控元件比如增强子和沉默子

• 启动子：是靠近转录起始点(TSS)的DNA区域，它包含转录因子的结合位点（transcription factor binding site，TFBS），所以转录因子能够结合在启动子上TFBS，招募RNA聚合酶进而转录基因。
• 增强子：一般位于启动子下游或上游1Mb的DNA区域，转录因子与增强子结合，并与启动子区域接触时，能够促进基因的转录。
• 沉默子：会减少或抑制基因的表达。

所以说，ATAC-seq可以帮助识别启动子区域、潜在的增强子或沉默子，也就是说，ATAC-seq中的peak，往往是启动子、增强子序列，以及一些反式调控因子结合的位点。

ref：

transposons

转座子是一类能够在基因组上移动其位置的DNA序列。

ref：

Tn5 transposase

ATAC-seq技术利用Tn5转座酶切割暴露的DNA，使其连接上特异性的Adapters，插入Adapters的DNA片段可被分离出来用于二代测序，从而获得大量的细胞系和组织样本基因表达、转录调控、转录修饰等的信息。

ref：

本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。原始发表：2025-03-29，如有侵权请联系 cloudcommunity@tencent 删除数据分析amp教程连接数据