在之前的文章中，我們解讀過一篇引用達(dá)2000多次的ATAC經(jīng)典文獻(xiàn)

引用2115次的ATAC經(jīng)典論文解讀

這篇文章中，使用了Genrich這個(gè)軟件來進(jìn)行peak calling。該軟件適用于chip_seq, DNase_seq, ATAC_seq等多種文庫的peak calling，源代碼保存在github上，鏈接如下

https://github.com/jsh58/Genrich

該軟件用法簡單，只需要輸入排序之后的bam文件即可，其peak calling的基本算法如下所示

1. 對原始bam文件進(jìn)行過濾，該軟件只針對雙端測序的bam文件，默認(rèn)情況下只保留雙端同時(shí)比對到基因組上的序列。和macs2不同，該軟件支持multimapping reads,對于比對到基因組多個(gè)位置的reads, 每個(gè)位置會計(jì)算一個(gè)權(quán)重。同時(shí)還可以過濾PCR重復(fù)序列，去除指定染色體上的序列

2. 如果提供了control樣本，首先利用control樣本構(gòu)建一個(gè)背景

3. 對于treat樣本，利用log-normal分布來計(jì)算每個(gè)位點(diǎn)的p值

4. 多重假設(shè)檢驗(yàn)的校正，將p值轉(zhuǎn)換為q值

5. 計(jì)算統(tǒng)計(jì)學(xué)顯著，即p值或者q值小于0.05的區(qū)域，對應(yīng)的曲線下的面積，即AUC, 設(shè)定AUC的閾值，如果一個(gè)區(qū)域的AUC大于閾值，則定義為peak

該軟件有兩種運(yùn)行模式，第一種適用于chip_seq, 第二種適用于ATAC。對于ATAC的reads，需要對reads比對位置進(jìn)行shift, 所以peak calling的模式也會有所不同，圖示如下

默認(rèn)情況下為第一種模式，通過-j參數(shù)可以調(diào)整為ATAC模式。默認(rèn)模式下，將fragment兩個(gè)末端對應(yīng)的區(qū)域直接作為peak, 以chip_seq為例，抗體抓下來的reads就是蛋白結(jié)合區(qū)域的reads。而ATAC文庫中Tn5切割的序列是在結(jié)合區(qū)域的兩側(cè)，所以在ATAC模式中，以fragment的中心點(diǎn)為基準(zhǔn)，來判斷真實(shí)的peak區(qū)域。

對于ATAC數(shù)據(jù)的peak calling而言，該軟件的用法如下

-t參數(shù)指定輸入的bam文件，-o參數(shù)指定輸出的peak文件，-f參數(shù)指定log文件的名稱。-r參數(shù)表示去除PCR重復(fù)，-x參數(shù)表示保留只有一端比對上基因組的reads，-e參數(shù)剔除指定染色體上的序列，-E參數(shù)剔除bed文件中指定的基因組區(qū)域。-q參數(shù)指定qvalue的閾值，-a參數(shù)指定AUC面積的閾值，

對于ATAC的peak calling而言，除了MACS2, 該軟件也是一個(gè)不錯(cuò)的選擇。