關(guān)于ATAC-seq分析，在網(wǎng)上看到兩篇關(guān)于同一片綜述的翻譯回懦，寫的很好：
ATAC-seq數(shù)據(jù)分析工具的比較和推薦（Genome Biology綜述）
ATAC-seq生信分析綜述翻譯(Genome Biology)
這一篇文獻(xiàn)是對現(xiàn)有的ATAC-sea分析給出的一個(gè)總結(jié)性的介紹堂鲤。里面提到了一個(gè)分析步驟，稱為footprinting。
什么是轉(zhuǎn)錄因子的footprinting演怎？
真核生物基因表達(dá)是一個(gè)復(fù)雜而有序的過程，它是眾多反式作用因子和順式作用元件之間相互作用的結(jié)果避乏。反式作用因子是指能直接或間接識別和結(jié)合在順式作用元件上爷耀，調(diào)控靶基因表達(dá)的蛋白質(zhì)因子，一般也稱為轉(zhuǎn)錄因子（transcriptional factor,TF）拍皮，轉(zhuǎn)錄因子結(jié)合位點(diǎn)((Transcription factor binding site, TFBS)是與轉(zhuǎn)錄因子結(jié)合的DNA 序列歹叮。確定TFBS是理解轉(zhuǎn)錄調(diào)控機(jī)制, 建立轉(zhuǎn)錄調(diào)控網(wǎng)絡(luò)的關(guān)鍵問題。
（參考：http://www.ebiotrade.com/custom/Genomics/180613/index.html）
ATAC-seq footprints可以幫助我們查看轉(zhuǎn)錄因子在全基因組上結(jié)合的狀態(tài)铆帽。ATAC-seq中的足跡是指活性TF與DNA結(jié)合咆耿，從而阻止Tn5在結(jié)合位點(diǎn)切割。所以客觀上形成一個(gè)被保護(hù)的區(qū)域爹橱，一個(gè)低coverage區(qū)域萨螺。比如上面這張圖。在轉(zhuǎn)錄因子結(jié)合的位置是一個(gè)“低谷”的形狀。（ATAC-seq數(shù)據(jù)分析）
目前有哪些軟件可以進(jìn)行footprinting分析屑迂？
在本文開頭提到的綜述里有一個(gè)表浸策，列出了目前所有可以進(jìn)行footprinting分析的軟件：
在這個(gè)表里的倒數(shù)第二列，標(biāo)出了這些軟件是否可以用來分析ATAC-seq數(shù)據(jù)惹盼。而在這些軟件里庸汗，HINT-ATAC考慮了ATAC-seq特有的偏向性(對鏈特異性的Tn5切割偏好進(jìn)行校正)，而且在低質(zhì)量和發(fā)現(xiàn)新的motif方面有優(yōu)勢手报。
所以綜合看來HINT-ATAC是個(gè)不錯(cuò)的選擇蚯舱。那么就來學(xué)習(xí)一下如何利用HINT-ATAC做footprinting分析。
HINT-ATAC是附屬在RGT軟件里的一個(gè)功能掩蛤。RGT全稱Regulatory Genomics Toolbox(網(wǎng)站：https://www.regulatory-genomics.org/)枉昏，顧名思義它是一個(gè)工具箱，里面很多工具揍鸟，而HINT-ATAC就是其中之一兄裂。RGT是一個(gè)開源的python庫，所以你需要確保你的電腦里安裝了python阳藻。需要注意的是晰奖，這個(gè)軟件目前只支持python2，不支持python3!
HINT-ATAC的官方教程：
http://www.regulatory-genomics.org/hint/introduction/
http://www.regulatory-genomics.org/hint/tutorial/
在本教程中腥泥，將展示如何使用HINT-ATAC來比較活化的轉(zhuǎn)錄因子的足跡變化匾南。使用在HINT-ATAC論文中提供的數(shù)據(jù)，對比兩個(gè)樹突細(xì)胞的footprint蛔外。從小鼠骨髓中提取并培養(yǎng)經(jīng)典的樹突細(xì)胞1型(cDC1)和漿細(xì)胞樣樹突細(xì)胞 (pDC)蛆楞，進(jìn)行了Omni ATAC-seq實(shí)驗(yàn)（原始fastq文件：here）。
這個(gè)教程里可供下載的文件是已經(jīng)進(jìn)行了比對和peak calling之后的文件夹厌。
NOTE:這個(gè)軟件的最后一步分析豹爹，一定要用服務(wù)器運(yùn)行，普通電腦是根本帶不起來的尊流。
（一）安裝RGT
首先你要確保你安裝了pip帅戒，然后用pip一鍵安裝：
$ pip install --user RGT
$ rgt-hint ＃查看幫助信息
usage: rgt-hint [-h] [--version]
{footprinting,differential,plotting,training,estimation,evaluation,evidence,tracks}
...
positional arguments:
{footprinting,differential,plotting,training,estimation,evaluation,evidence,tracks}
Commands:
footprinting detect footprints based on reads.bam and regions.bed
differential perform differential analysis based on footprints of
two conditions
plotting generate plots based on input
training train Hidden Markov models
estimation estimate sequence specific bias
evaluation evaluate the predicted footprints based on ChIP-seq
data
evidence create evidence file based on motif matching results
and ChIP-seq data
tracks create wig track file for visualization
optional arguments:
-h, --help show this help message and exit
--version show program's version number and exit
（二）下載示例數(shù)據(jù)
比較大崖技，一共十幾個(gè)G：
$ wget --no-check-certificate https://costalab.ukaachen.de/open_data/hint/tutorial/DendriticCells.tar.gz
$ tar xvfz DendriticCells.tar.gz
$ cd DendriticCells
（三）處理基因組
參考官網(wǎng)步驟：http://www.regulatory-genomics.org/rgt/rgt-data-folder/逻住，這個(gè)軟件支持的基因組文件有：hg19, hg38, mm9, mm10, zv9, and zv10。
這里使用的基因組不是我們平時(shí)下載的基因組的fa或者gtf文件迎献，是你在上面安裝RGT的時(shí)候瞎访，里面自帶的基因組信息，然后需要用python進(jìn)行配置（configuration）：
$ cd ~/rgtdata
$ python setupGenomicData.py --mm9 #這一步會(huì)運(yùn)行幾分鐘
（四）對兩種細(xì)胞進(jìn)行call footprint
NOTE：從這一步你需要注意的是吁恍，如果你用的是linux系統(tǒng)扒秸，你先需要保證上面提到的rgtdata這個(gè)文件夾放在home目錄里播演。因?yàn)檫@個(gè)軟件是高度依賴python的，我的conda是安裝在了home里伴奥。如果你的rgtdata這個(gè)文件夾沒有拷貝到home里写烤，這里會(huì)報(bào)錯(cuò)，例如：
FileNotFoundError: [Errno 2] No such file or directory: '/home/yf/rgtdata/data.config'
這一步的時(shí)候會(huì)比較費(fèi)時(shí)（如果你的電腦配置不是很高）：
$ rgt-hint footprinting --atac-seq --paired-end --organism=mm9 --output-location=(PATH) --output-prefix=cDC1 cDC1.bam cDC1_peaks.narrowPeak
[total time: 2h 15m 8s]＃運(yùn)行后會(huì)顯示一共花費(fèi)了多少時(shí)間拾徙，我的電腦是i7,8G內(nèi)存
$ rgt-hint footprinting --atac-seq --paired-end --organism=mm9 --output-location=(PATH) --output-prefix=pDC pDC.bam pDC_peaks.narrowPeak
[total time: 1h 22m 23s]
Hint需要一個(gè)bam文件(indexed and sorted)和narrowpeak文件作為輸入文件洲炊。這個(gè)軟件只考慮peak區(qū)域進(jìn)行footprinting，以加快分析速度尼啡。每個(gè)命令運(yùn)行后暂衡，將在當(dāng)前文件夾中生成一個(gè)輸出文件，其中包含footprinting的基因組位置崖瞭。它還生成一個(gè)以“.info”結(jié)尾的文件狂巢，其中包含來自庫的統(tǒng)計(jì)信息，即總讀取次數(shù)等等书聚。當(dāng)指定paired-end時(shí)唧领，會(huì)顯著提高footprinting的預(yù)測準(zhǔn)確度。另外你還需要指定帶有標(biāo)記-有機(jī)體的基因組/有機(jī)體寺惫。
上面兩行代碼得到４個(gè)文件：
打開一個(gè)info文件疹吃，里面是這樣的信息：
Number of reads: <map object at 0x7fb5ed18c160>
Number of peaks: 94222 #顯示這個(gè)bam文件里一共多少個(gè)peak
Number of tag counts within peaks: 79092202.0 #在這些peak里有多少個(gè)count數(shù)
Number of footprints: 513503 #有多少個(gè)footprints，也就是說在這些peak里有多少個(gè)TF結(jié)合位點(diǎn)
（五）生成bw文件
HINT還可以輸出用于可視化peak的信號（比如IGV）西雀。所以你還可以用下面的命令生成基因組配置文件(BigWig文件)：
$ rgt-hint tracks --bc --bigWig --organism=mm9 cDC1.bam cDC1_peaks.narrowPeak --output-prefix=cDC1_BC
[total time: 0h 14m 21s]
$ rgt-hint tracks --bc --bigWig --organism=mm9 pDC.bam pDC_peaks.narrowPeak --output-prefix=pDC_BC
[total time: 0h 9m 23s]
生成的bigwig文件包含ATAC-seq的read數(shù)量，這個(gè)數(shù)量是在信號標(biāo)準(zhǔn)化和偏差校正后歉摧，由HINT-ATAC估計(jì)的在每個(gè)基因組位置的read數(shù)艇肴。要注意的是，上面的代碼目前只支持偏差校正和ATAC-seq的移位叁温。這比簡單地查看bam文件的覆蓋率概要更精確再悼。
用IGV打開我們上面生成的bed文件(footprint)和bw文件，需要注意的是你要選擇mm9作為參考基因組膝但，因?yàn)閷?shí)驗(yàn)是基于小鼠基因組mm9的冲九。比如查看基因Zbtb46周圍的peak情況，這個(gè)基因在cDC里表達(dá)跟束，在pDC不表達(dá)莺奸。但是在IGV里兩個(gè)細(xì)胞的bed文件在這個(gè)基因附近都有peak，不同的是cDC1細(xì)胞的bed文件里有一個(gè)特定的區(qū)域有footprint冀宴，說明這個(gè)footprint有TF的結(jié)合灭贷，有可能促進(jìn)這個(gè)基因的表達(dá):
（六）預(yù)測細(xì)胞特異性的轉(zhuǎn)錄因子
HINT-ATAC的另一個(gè)重要的應(yīng)用是根據(jù)細(xì)胞特異性尋找轉(zhuǎn)錄因子。RGT還提供了一個(gè)工具略贮，利用JASPAR, UNIPROBE或者HOCOMOCO數(shù)據(jù)庫里的motifs尋找motif結(jié)合位點(diǎn)（MPBS）甚疟。
比如下面的代碼仗岖，使用JASPAR數(shù)據(jù)庫中的motifs作為默認(rèn)值進(jìn)行motif匹配:
$ rgt-motifanalysis matching --organism=mm9 --input-files pDC.bed cDC1.bed
>> output location: /media/yf/ATAC/Footprinting/hint/DendriticCells/match
>> genome: mm9
>> pseudocounts: 1.0
>> fpr threshold: 0.0001
>>> input file pDC loaded: 312544 regions
>>> input file cDC1 loaded: 513503 regions
>> all files loaded
>> used database(s): jaspar_vertebrates
>> motifs loaded: 746
>> matching [pDC], 312544 regions... [297.790 seconds]
>> matching [cDC1], 513503 regions... [485.401 seconds]
[total time: 0h 13m 32s]
在上面的第四步里，生成了兩個(gè)細(xì)胞的bed文件览妖，這一步的input文件就是bed文件轧拄。根據(jù)你的bed文件里footprint數(shù)量的多少和motif數(shù)據(jù)庫的大小，這一步需要的時(shí)間不同讽膏。運(yùn)行后會(huì)生成一個(gè)新文件夾檩电，名字與bed文件相匹配，并且產(chǎn)生一個(gè)新的bed文件桅打，包含每一個(gè)footprint區(qū)域匹配上的motif是嗜。
其中第四列是motif的名字，第五列是motif匹配的bit-score挺尾。你可以閱讀這個(gè)網(wǎng)頁了解更多的信息：Basic Introduction
$ less cDC1_mpbs.bed
chr1 4759959 4759968 MA0077.1.SOX9 10.85478888971908 +
chr1 4759957 4759968 MA0143.4.SOX2 13.36651132464563 -
chr1 4759957 4759968 MA1120.1.SOX13 12.576069892112077 -
chr1 4759959 4759969 MA1152.1.SOX15 11.928506391469964 +
chr1 4760005 4760015 MA0100.3.MYB 12.303021888774307 +
（七）Generate average ATAC-seq profiles around binding sites of particular TF
接下來鹅搪，我們使用HINT生成在特定TF結(jié)合位點(diǎn)周圍的平均ATAC-seq profile。這個(gè)分析使得我們檢查每個(gè)特定TF的染色質(zhì)可接近性遭铺。此外丽柿，通過比較兩個(gè)ATAC-seq文件(cDC1細(xì)胞和pDC細(xì)胞)，我們可以了解兩個(gè)細(xì)胞中TF結(jié)合的變化魂挂。
下面這行代碼將讀取上面第四步里生成的motif 匹配文件和原始的bam文件甫题。其中代碼里的-bc的意思是使用bias矯正信號（目前只支持ATAC-seq數(shù)據(jù)）。-nc是cpu數(shù)量：
#這一行代碼是教程里的原代碼
$ rgt-hint differential --organism=mm9 --bc --nc 30 --mpbs-files=./match/cDC1_mpbs.bed,./match/pDC_mpbs.bed --reads-files=cDC1.bam,pDC.bam --conditions=cDC1,pDC --output-location=cDC1_pDC
我在服務(wù)器里的使用的代碼：
#你可以看到這里我調(diào)用了128G的內(nèi)存涂召，之前也試過8G和16G坠非，完全跑不起來
$ srun --mem=128G --cpus-per-task=20 --time=6:00:00 --pty rgt-hint differential --organism=mm9 --bc --nc 20 --mpbs-files=./match/cDC1_mpbs.bed,./match/pDC_mpbs.bed --reads-files=cDC1.bam,pDC.bam --conditions=cDC1,pDC --output-location=cDC1_pDC
srun: job 7612760 queued and waiting for resources
srun: job 7612760 has been allocated resources
generating signal for each motif and condition...
signal generation is done!
generating line plot for each motif...
[total time: 1h 3m 15s]
這里我用了128G的內(nèi)存來運(yùn)行，還要1個(gè)小時(shí)果正。炎码。。
運(yùn)行后秋泳，產(chǎn)生一個(gè)新的文件夾名為“cDC1_pDC”潦闲。在這個(gè)文件夾里，有一個(gè)scatter plot圖迫皱，展示的是在兩種細(xì)胞之間TF動(dòng)力學(xué)活性歉闰。y軸代表TF的不同的activity，紅色標(biāo)記的是有顯著性差異的TF(p<0.05)卓起。x軸是隨機(jī)數(shù)（for jittering purposes）：
這個(gè)圖上面的紅色點(diǎn)代表的TF是在pDC細(xì)胞里結(jié)合比較強(qiáng)的和敬，下面的紅點(diǎn)代表的TF是在cDC1細(xì)胞里結(jié)合比較強(qiáng)的。這個(gè)圖有點(diǎn)丑既绩，另外標(biāo)記也都重疊在一起了概龄，你可以自己用R重新畫個(gè)圖（因?yàn)槭鞘纠龜?shù)據(jù)，就不單獨(dú)美化這個(gè)散點(diǎn)圖了）饲握，這個(gè)點(diǎn)圖對應(yīng)的txt文件是生成的differential_statistics.txt：
這個(gè)txt打開后是這樣的：
另外蚕键，在上面文件夾截圖里顯示的“Lineplots”里，你還可以找到所有TF的曲線圖衰粹，比如挑兩個(gè)有顯著差異的：
