1. 尋找開發(fā)區(qū)域
ATACseq 的一個共同目標是識別轉錄因子結合和/或轉錄機制活躍的無核小體區(qū)域陨收。該核小體游離信號對應于小于一個核小體的片段（如 Greenleaf 論文中定義 < 100bp）篓足。
然而踱稍，為了識別開放的染色質品山，我們可以簡單地使用在測序中正確配對的所有讀數(shù)（< 2000bp）即彪。
2. 無核小體區(qū)域
有許多方法可用于從 ATACseq 數(shù)據(jù)中調(diào)用無核小體區(qū)域，其中許多方法借鑒自 ChIPseq 分析一疯。一種非常流行且標準的 ATACseq 峰值調(diào)用程序是 MAC2。
2.1. 單端數(shù)據(jù)
使用 ATACseq 的單端測序夺姑，我們不知道片段有多長墩邀。因此，與 ChIPseq 相比盏浙，要識別開放區(qū)域眉睹，MACS2 峰值調(diào)用需要一些不同的參數(shù)。采用的一種策略是將讀取 5' 末端移動 -100废膘，然后從該位置延伸 200bp竹海。考慮到我們的無核小體片段的預期大小丐黄，這應該提供適合 MACS2 窗口大小的核小體區(qū)域的堆積斋配。
這是我們用來執(zhí)行此操作的 MACS 系統(tǒng)調(diào)用。
MACS2 callpeak -t singleEnd.bam --nomodel --shift -100
--extsize 200 --format BAM -g MyGenome
在 R 中灌闺，我們可以使用 Herper 運行這個系統(tǒng)調(diào)用艰争，這樣我們就可以訪問我們安裝的 MACS2。 MACS2 已安裝到 ATACseq_analysis 中菩鲜。所以我們可以使用 with_CondaEnv() 從 R 中使用這個環(huán)境园细。
with_CondaEnv("ATACseq_analysis",
system2(command="macs2",args =c("callpeak",
"-t", "singleEnd.bam",
"--nomodel",
"--shift","-100",
"--extsize", "200",
"--format", "BAM",
"-g", "hs")),
stdout = TRUE))
或者對于核小體占據(jù)的數(shù)據(jù)，我們可以調(diào)整移位和延伸以將信號集中在核小體中心（包裹在 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. 雙端數(shù)據(jù)
如果我們對配對末端數(shù)據(jù)進行了測序猛频，那么我們確實知道片段長度，并且可以向 MACS2 提供 BAM 文件蛛勉，這些文件已經(jīng)過預過濾以正確配對（如果您想?yún)^(qū)分核小體和無核小體區(qū)域鹿寻，則提供片段大小）诽凌。
我們必須告訴 MACS2 數(shù)據(jù)是使用格式參數(shù)配對的毡熏。這很重要，因為默認情況下 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)
對于此處的配對末端數(shù)據(jù)痢法，我們從無核小體區(qū)域 BAM 文件中調(diào)用無核小體區(qū)域的峰。
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 來檢查低質量尖奔、重復和信號分布。在我們刪除任何數(shù)據(jù)之前，我們可以快速評估我們的峰值讀取提茁、重復率淹禾、低質量讀取和來自 ChIPQC 的偽像區(qū)域中的讀取。
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 數(shù)據(jù)的一些重要指標茴扁，例如來自 QCmetrics() 函數(shù)的峰值讀取和黑名單中的讀取以及來自 flagtagcounts() 函數(shù)的重復讀取數(shù)铃岔。
myMetrics <- QCmetrics(qcRes)
myMetrics[c("RiBL%", "RiP%")]
flgCounts <- flagtagcounts(qcRes)
DupRate <- flgCounts["DuplicateByChIPQC"]/flgCounts["Mapped"]
DupRate * 100
4. 黑名單刪除
來自測序的人工制品和不完美的基因組構建可能會混淆我們的結果。這些工件已被整理到區(qū)域的“黑名單”中丹弱。
由于列入黑名單的區(qū)域可能會混淆我們的分析德撬，因此我們刪除了在那里被調(diào)用的所有峰值铲咨。
過早刪除黑名單可能會隱藏數(shù)據(jù)中的一些 QC 問題躲胳。在您的分析中應始終考慮黑名單，并建議在考慮 QC 后從這些區(qū)域中刪除數(shù)據(jù)纤勒。
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]
本文由mdnice多平臺發(fā)布
