ATAC-seq分析：Peak Calling（8）
ATAC-seq分析：Peak Calling（8）
1. 寻找开发区域
ATACseq 的一个共同目标是识别转录因子结合和/或转录机制活跃的无核小体区域。该核小体游离信号对应于小于一个核小体的片段（如 Greenleaf 论文中定义 < 100bp）。
然而，为了识别开放的染色质，我们可以简单地使用在测序中正确配对的所有读数（< 2000bp）。
2. 无核小体区域
有许多方法可用于从 ATACseq 数据中调用无核小体区域，其中许多方法借鉴自 ChIPseq 分析。一种非常流行且标准的 ATACseq 峰值调用程序是 MAC2。
2.1. 单端数据
使用 ATACseq 的单端测序，我们不知道片段有多长。因此，与 ChIPseq 相比，要识别开放区域，MACS2 峰值调用需要一些不同的参数。采用的一种策略是将读取 5' 末端移动 -100，然后从该位置延伸 200bp。考虑到我们的无核小体片段的预期大小，这应该提供适合 MACS2 窗口大小的核小体区域的堆积。
这是我们用来执行此操作的 MACS 系统调用。
MACS2 callpeak -t singleEnd.bam --nomodel --shift -100
--extsize 200 --format BAM -g MyGenome
在 R 中，我们可以使用 Herper 运行这个系统调用，这样我们就可以访问我们安装的 MACS2。 MACS2 已安装到 ATACseq_analysis 中。所以我们可以使用 with_CondaEnv() 从 R 中使用这个环境。
with_CondaEnv("ATACseq_analysis",
system2(command="macs2",args =c("callpeak",
"-t", "singleEnd.bam",
"--nomodel",
"--shift","-100",
"--extsize", "200",
"--format", "BAM",
"-g", "hs")),
stdout = TRUE))
或者对于核小体占据的数据，我们可以调整移位和延伸以将信号集中在核小体中心（包裹在 147bp DNA 中的核小体）。
MACS2 callpeak -t singleEnd.bam --nomodel --shift 37
--extsize 73 --format BAM -g MyGenome
- R
with_CondaEnv("ATACseq_analysis", system2(command = "macs2", args = c("callpeak",
"-t", "singleEnd.bam", "--nomodel", "--shift", "37", "--extsize", "73", "--format",
"BAM", "-g", "hs")), stdout = TRUE)
2.2. 双端数据
如果我们对配对末端数据进行了测序，那么我们确实知道片段长度，并且可以向 MACS2 提供 BAM 文件，这些文件已经过预过滤以正确配对（如果您想区分核小体和无核小体区域，则提供片段大小）。
我们必须告诉 MACS2 数据是使用格式参数配对的。这很重要，因为默认情况下 MACS2 会猜测它是单端 BAM。
MACS2 callpeak -t pairedEnd.bam -f BAMPE
--outdir path/to/output/
--name pairedEndPeakName -g MyGenome
- R
with_CondaEnv("ATACseq_analysis", system2(command = "macs2", args = c("callpeak",
"-t", "pairedEnd.bam", "--format", "BAMPE", "--outdir", "/Documents/ATAC_MACS2_calls/",
"--name", "pairedEndPeakName", "-g", "hs")), stdout = TRUE)
对于此处的配对末端数据，我们从无核小体区域 BAM 文件中调用无核小体区域的峰。
MACS2 callpeak -t ~/Downloads/Sorted_ATAC_50K_2_openRegions.bam
--outdir ATAC_Data/ATAC_Peaks/ATAC_50K_2
--name Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak
-f BAMPE -g hs
- R
with_CondaEnv("ATACseq_analysis", system2(command = "macs2", args = c("callpeak",
"-t", "~/Downloads/Sorted_ATAC_50K_2_openRegions.bam", "--outdir", "ATAC_Data/ATAC_Peaks/ATAC_50K_2",
"--name", "Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak", "-f", "BAMPE", "-g",
"hs")), stdout = TRUE)
在之后，我们将获得 3 个文件。
- Name.narrowPeak – 适用于 IGV 和进一步分析的格式
- Name_peaks.xls – 适合在 excel 中查看的峰值表。（实际上不是 xls，而是 tsv）
- summits.bed – 用于查找
motifs和绘图的峰顶位置
3. QC
通常我们想做一些 QC 来检查低质量、重复和信号分布。在我们删除任何数据之前，我们可以快速评估我们的峰值读取、重复率、低质量读取和来自 ChIPQC 的伪像区域中的读取。
library(ChIPQC)
library(rtracklayer)
library(DT)
library(dplyr)
library(tidyr)
blkList <- import.bed("~/Downloads/ENCFF001TDO.bed.gz")
openRegionPeaks <- "~/Downloads/ATAC_Workshop/ATAC_Data/ATAC_Peaks/Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak"
qcRes <- ChIPQCsample("~/Downloads/ATAC_Workshop/ATAC_Data/ATAC_BAM/Sorted_ATAC_50K_2_openRegions.bam",
peaks = openRegionPeaks, annotation = "hg19", chromosomes = "chr20", blacklist = blkList,
verboseT = FALSE)
我们可以使用 ChIPQC 包来捕获我们的 ATACseq 数据的一些重要指标，例如来自 QCmetrics() 函数的峰值读取和黑名单中的读取以及来自 flagtagcounts() 函数的重复读取数。
myMetrics <- QCmetrics(qcRes)
myMetrics[c("RiBL%", "RiP%")]
flgCounts <- flagtagcounts(qcRes)
DupRate <- flgCounts["DuplicateByChIPQC"]/flgCounts["Mapped"]
DupRate * 100
4. 黑名单删除
来自测序的人工制品和不完美的基因组构建可能会混淆我们的结果。这些工件已被整理到区域的“黑名单”中。
由于列入黑名单的区域可能会混淆我们的分析，因此我们删除了在那里被调用的所有峰值。
过早删除黑名单可能会隐藏数据中的一些 QC 问题。在您的分析中应始终考虑黑名单，并建议在考虑 QC 后从这些区域中删除数据。
MacsCalls <- granges(qcRes[seqnames(qcRes) %in% "chr20"])
data.frame(Blacklisted = sum(MacsCalls %over% blkList), Not_Blacklisted = sum(!MacsCalls %over%
blkList))
MacsCalls <- MacsCalls[!MacsCalls %over% blkList]