在R语言中的 ATACseq 数据分析全流程实战（四）：Peak calling
在R语言中的 ATACseq 数据分析全流程实战（四）：Peak calling
ATAC-seq（转座酶可及染色质测序）利用转座酶，提供了一种从单个样本中同时提取转录因子结合位点和核小体位置信号的方法。
这意味着我们的数据可能包含多种信号类型的混合：
- 会检测到无核小体区域以及转录因子周围的信号（即较短的片段）；
- 还有一部分信号来自开放染色质中的核小体周围（较长的片段）。
我们所有的数据都来自开放染色质区域，因为这些区域是转座酶能够进入并发挥作用的地方：
评估TSS信号
如果较短的片段代表转录因子和转录机制周围的开放区域，那么预期会在转录起始位点看到信号。
较长的片段将代表核小体周围的信号，因此信号应该出现在转录起始位点之外，并且更多地集中在 +1 和 -1 核小体位置。
TBP相关因子(TAFs) ：
TATA box结合蛋白(TBP)，TFIID共有13个TBP相关因子，即TAFs。TAFs的功能主要表现在两个方面：与TATA box以外的启动子元件(Inr和DCE等)结合，以及与激活因子结合。有些没有TATA box的基因，必须以来TAFs来启动转录。
Phased nucleosomes：
“Phased nucleosomes”可以翻译为“相位核小体”或“规则排列的核小体”。在生物学中，核小体的相位排列（phasing）是指核小体在基因组上的有序排列，通常与转录起始位点（TSS）或其他特定的基因组特征相关。这种排列方式在基因表达调控、染色质结构维持等方面具有重要作用。
我们可以创建一个覆盖所有转录起始位点（TSS）区域的meta-plot，以展示核小体缺失区域和核小体占据区域的信号分布最为集中的位置。
meta-plot通过对一组区域的信号进行平均或求和，来识别数据中的趋势。
绘制上面的图，我们可以使用soGGi包：输入数据为bam 以及参考基因组区域。
bam为：在R语言中的 ATACseq 数据分析全流程实战（二）这一步骤中生成的：ATAC_50K_2_sorted.bam
或者教程中有直接提供的处理好的：
# bam文件
https://s3.amazonaws.com/rubioinformatics/ATAC_Workshop/ATAC_Data/ATAC_BAM/Sorted_ATAC_50K_2.bam
https://s3.amazonaws.com/rubioinformatics/ATAC_Workshop/ATAC_Data/ATAC_BAM/Sorted_ATAC_50K_2.bam.bai
在转录起始位点（TSS）区域进行绘图，因此首先需要为 hg19 基因组生成一个包含 TSS 位置的 GRanges 对象：
format：提供输入文件格式的信息paired：说明数据是否为配对数据style：指定绘图类型minFragmentLength和maxFragmentLength：指定用于绘图的配对读取片段的最小和最大长度。这样能够仅选择核小体缺失区域的信号（<100碱基对），以在转录起始位点（TSS）区域生成metaplot。
rm(list=ls())
# BiocManager::install("soGGi")
library(soGGi)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
TxDb.Hsapiens.UCSC.hg19.knownGene
# extract gene locations (TSS to TTS) using the genes() function and our TxDb object
genesLocations <- genes(TxDb.Hsapiens.UCSC.hg19.knownGene)
genesLocations
# 使用 `resize()` 函数以考虑链方向的方式提取每个基因的起始位置（即转录起始位点TSS）。
tssLocations <- resize(genesLocations, fix = "start", width = 1)
tssLocations
# 当我们创建索引时，我们将基因组限定在主要染色体上。
# 我们也可以对我们的TSS（转录起始位点）GRanges对象进行同样的操作，并更新其水平（分类）层级。
# 这意味着BAM文件和GRanges对象能够更好地兼容和协同工作。
mainChromosomes <- paste0("chr", c(1:21, "X", "Y", "M"))
mainChromosomes
seqnames(tssLocations)
myindex <- (match(seqnames(tssLocations), mainChromosomes))
myindex <- myindex[!is.na(myindex)]
as.numeric(myindex)
tssLocations <- tssLocations[as.numeric(myindex)]
seqlevels(tssLocations) <- mainChromosomes
tssLocations
library(soGGi)
sortedBAM <- "GSE47753/ATAC_50K_2_sorted.bam"
library(Rsamtools)
# Nucleosome free
allSignal <- regionPlot(bamFile = sortedBAM, testRanges = tssLocations)
allSignal
# 设置其它参数
nucFree <- regionPlot(bamFile = sortedBAM, testRanges = tssLocations, style = "point",
format = "bam", paired = TRUE, minFragmentLength = 0, maxFragmentLength = 100,
forceFragment = 50)
class(nucFree)
plotRegion(nucFree)
在这里，我们看到了在转录起始位点（TSS）区域上方的核小体缺失区域的预期信号峰值：
我们可以通过调整minFragmentLength和maxFragmentLength参数为核小体长度片段的预期范围（此处为180到240），来创建单核小体信号的图：
# mono-nucleosome signal
monoNuc <- regionPlot(bamFile = sortedBAM, testRanges = tssLocations, style = "point",
format = "bam", paired = TRUE, minFragmentLength = 180, maxFragmentLength = 240,
forceFragment = 80)
plotRegion(monoNuc)
在这个图表中，我们可以清晰地看到预期的+1核小体信号峰值，以及其他几个核小体信号峰值：
Peak calling
找到染色质开放区域
在 ATAC-seq 中，一个常见的目标是识别转录因子结合和/或转录机制活跃的无核小体区域。这种无核小体信号对应于小于一个核小体的片段（如 Greenleaf 文章中定义的 <100bp）。然而，为了识别开放染色质，我们其实可以简单地使用所有在测序中正确配对的reads（<2000bp）。
后面聚焦于分析 ATAC-seq 数据中无核小体的部分，bam为在R语言中的 ATACseq 数据分析全流程实战（二）中的ATAC_50K_2_sorted.bam文件。
先得到 无核小体区域的bam：
# 得到 ATAC_50K_2_sorted_nucFreeRegions.bam
##
if(F) {
################### 得到 ATAC_50K_2_sorted_nucFreeRegions.bam
library(GenomicAlignments)
sortedBAM <- "GSE47753/ATAC_50K_2_sorted.bam"
flags = scanBamFlag(isProperPair = TRUE)
myParam = ScanBamParam(flag = flags, what = c("qname", "mapq", "isize"))
myParam
atacReads <- readGAlignmentPairs(sortedBAM, param = myParam)
class(atacReads)
atacReads[1:2, ]
# 使用first()和second()来访问GAlignments对象，以分别获取第一条或第二条读取的信息
read1 <- GenomicAlignments::first(atacReads)
read2 <- GenomicAlignments::second(atacReads)
read1[1, ]
############### 提取插入片段大小
atacReads_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)
################ 创建代表无核小体、单核小体和双核小体的读取的BAM文件
atacReads_NucFree <- atacReads[insertSizes < 100, ]
atacReads_MonoNuc <- atacReads[insertSizes > 180 & insertSizes < 240, ]
atacReads_diNuc <- atacReads[insertSizes > 315 & insertSizes < 437, ]
nucFreeRegionBam <- gsub("\\.bam", "_nucFreeRegions\\.bam", sortedBAM);nucFreeRegionBam
monoNucBam <- gsub("\\.bam", "_monoNuc\\.bam", sortedBAM);monoNucBam
diNucBam <- gsub("\\.bam", "_diNuc\\.bam", sortedBAM);diNucBam
library(rtracklayer)
export(atacReads_NucFree, nucFreeRegionBam, format = "bam")
export(atacReads_MonoNuc, monoNucBam, format = "bam")
export(atacReads_diNuc, diNucBam, format = "bam")
}
有许多方法可用于从 ATAC-seq 数据中识别无核小体区域，其中许多方法是从 ChIP-seq 分析中借鉴而来的。MACS2 是一种非常流行且标准的 ATAC-seq peaks 识别工具。
MACS2 还发布了一个对应的R包：https://bioconductor.org/packages/release/bioc/vignettes/MACSr/inst/doc/MACSr.html（下次学习）。
这次先用bash版本的。
R语言中使用conda来运行bash命令，来看看~ 需要借助R包Herper。
bash命令：
macs2 callpeak -t GSE47753/ATAC_50K_2_sorted_nucFreeRegions.bam \
--outdir GSE47753/ \
--name GSE47753/Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak \
-f BAMPE -g hs
R版本命令：
# 使用已有conda环境 chip, 在chip中安装macs2
salmon_paths <- install_CondaTools(tools = "masc2", env = "chip",
updateEnv = T,
pathToMiniConda = "/nas2/zhangj/biosoft/miniconda3/")
salmon_paths
# 运行
with_CondaEnv(new = "chip", # conda 的环境名称
pathToMiniConda = "/nas2/zhangj/biosoft/miniconda3/", # conda的路径
code = system2(command = "macs2", args = c("callpeak",
"-t", "GSE47753/ATAC_50K_2_sorted_nucFreeRegions.bam",
"--outdir", "GSE47753/",
"--name", "Sorted_ATAC_50K_2_Small_Paired_peaks.narrowPeak",
"-f", "BAMPE",
"-g","hs"))
)
上面的结果运行还要一点时间，就让它运行着吧~
peak calling会得到3个结果文件：
- Name.narrowPeak -- Narrow peak格式，用于IGV可视化和其他分析
- Name_peaks.xls -- Peak 结果表格，tsv格式，可以使用表格打开
- summits.bed -- Summit positions for peaks useful for finding motifs and plotting
腾讯云开发者
