使用HINT-ATAC进行ATAC-Seq的footprinting分析

关于ATAC-seq分析，在网上看到两篇关于同一片综述的翻译，写的很好：
ATAC-seq数据分析工具的比较和推荐（Genome Biology综述）
ATAC-seq生信分析综述翻译(Genome Biology)

这一篇文献是对现有的ATAC-sea分析给出的一个总结性的介绍。里面提到了一个分析步骤，称为footprinting。

什么是转录因子的footprinting？

真核生物基因表达是一个复杂而有序的过程，它是众多反式作用因子和顺式作用元件之间相互作用的结果。反式作用因子是指能直接或间接识别和结合在顺式作用元件上，调控靶基因表达的蛋白质因子，一般也称为转录因子（transcriptional factor,TF），转录因子结合位点((Transcription factor binding site, TFBS)是与转录因子结合的DNA 序列。确定TFBS是理解转录调控机制, 建立转录调控网络的关键问题。
（参考：http://www.ebiotrade.com/custom/Genomics/180613/index.html）

ATAC-seq footprints可以帮助我们查看转录因子在全基因组上结合的状态。ATAC-seq中的足迹是指活性TF与DNA结合，从而阻止Tn5在结合位点切割。所以客观上形成一个被保护的区域，一个低coverage区域。比如上面这张图。在转录因子结合的位置是一个“低谷”的形状。（ATAC-seq数据分析）

目前有哪些软件可以进行footprinting分析？

在本文开头提到的综述里有一个表，列出了目前所有可以进行footprinting分析的软件：

在这个表里的倒数第二列，标出了这些软件是否可以用来分析ATAC-seq数据。而在这些软件里，HINT-ATAC考虑了ATAC-seq特有的偏向性(对链特异性的Tn5切割偏好进行校正)，而且在低质量和发现新的motif方面有优势。

所以综合看来HINT-ATAC是个不错的选择。那么就来学习一下如何利用HINT-ATAC做footprinting分析。

HINT-ATAC是附属在RGT软件里的一个功能。RGT全称Regulatory Genomics Toolbox(网站：https://www.regulatory-genomics.org/)，顾名思义它是一个工具箱，里面很多工具，而HINT-ATAC就是其中之一。RGT是一个开源的python库，所以你需要确保你的电脑里安装了python。需要注意的是，这个软件目前只支持python2，不支持python3!

HINT-ATAC的官方教程：
http://www.regulatory-genomics.org/hint/introduction/
http://www.regulatory-genomics.org/hint/tutorial/
在本教程中，将展示如何使用HINT-ATAC来比较活化的转录因子的足迹变化。使用在HINT-ATAC论文中提供的数据，对比两个树突细胞的footprint。从小鼠骨髓中提取并培养经典的树突细胞1型(cDC1)和浆细胞样树突细胞 (pDC)，进行了Omni ATAC-seq实验（原始fastq文件：here）。

这个教程里可供下载的文件是已经进行了比对和peak calling之后的文件。

NOTE:这个软件的最后一步分析，一定要用服务器运行，普通电脑是根本带不起来的。

（一）安装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

（二）下载示例数据

比较大，一共十几个G：

$ wget --no-check-certificate https://costalab.ukaachen.de/open_data/hint/tutorial/DendriticCells.tar.gz
$ tar xvfz DendriticCells.tar.gz
$ cd DendriticCells

（三）处理基因组

参考官网步骤：http://www.regulatory-genomics.org/rgt/rgt-data-folder/，这个软件支持的基因组文件有：hg19, hg38, mm9, mm10, zv9, and zv10。

这里使用的基因组不是我们平时下载的基因组的fa或者gtf文件，是你在上面安装RGT的时候，里面自带的基因组信息，然后需要用python进行配置（configuration）：

$ cd ~/rgtdata
$ python setupGenomicData.py --mm9　#这一步会运行几分钟

（四）对两种细胞进行call footprint

NOTE：从这一步你需要注意的是，如果你用的是linux系统，你先需要保证上面提到的rgtdata这个文件夹放在home目录里。因为这个软件是高度依赖python的，我的conda是安装在了home里。如果你的rgtdata这个文件夹没有拷贝到home里，这里会报错，例如：

FileNotFoundError: [Errno 2] No such file or directory: '/home/yf/rgtdata/data.config'

这一步的时候会比较费时（如果你的电脑配置不是很高）：

$ 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]＃运行后会显示一共花费了多少时间，我的电脑是i7,8G内存

$ 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需要一个bam文件(indexed and sorted)和narrowpeak文件作为输入文件。这个软件只考虑peak区域进行footprinting，以加快分析速度。每个命令运行后，将在当前文件夹中生成一个输出文件，其中包含footprinting的基因组位置。它还生成一个以“.info”结尾的文件，其中包含来自库的统计信息，即总读取次数等等。当指定paired-end时，会显著提高footprinting的预测准确度。另外你还需要指定带有标记-有机体的基因组/有机体。

上面两行代码得到４个文件：

打开一个info文件，里面是这样的信息：

Number of reads: 

（五）生成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数量，这个数量是在信号标准化和偏差校正后，由HINT-ATAC估计的在每个基因组位置的read数。要注意的是，上面的代码目前只支持偏差校正和ATAC-seq的移位。这比简单地查看bam文件的覆盖率概要更精确。

用IGV打开我们上面生成的bed文件(footprint)和bw文件，需要注意的是你要选择mm9作为参考基因组，因为实验是基于小鼠基因组mm9的。比如查看基因Zbtb46周围的peak情况，这个基因在cDC里表达，在pDC不表达。但是在IGV里两个细胞的bed文件在这个基因附近都有peak，不同的是cDC1细胞的bed文件里有一个特定的区域有footprint，说明这个footprint有TF的结合，有可能促进这个基因的表达:

（六）预测细胞特异性的转录因子

HINT-ATAC的另一个重要的应用是根据细胞特异性寻找转录因子。RGT还提供了一个工具，利用JASPAR, UNIPROBE或者HOCOMOCO数据库里的motifs寻找motif结合位点（MPBS）。

比如下面的代码，使用JASPAR数据库中的motifs作为默认值进行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]

在上面的第四步里，生成了两个细胞的bed文件，这一步的input文件就是bed文件。根据你的bed文件里footprint数量的多少和motif数据库的大小，这一步需要的时间不同。运行后会生成一个新文件夹，名字与bed文件相匹配，并且产生一个新的bed文件，包含每一个footprint区域匹配上的motif。

其中第四列是motif的名字，第五列是motif匹配的bit-score。你可以阅读这个网页了解更多的信息：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结合位点周围的平均ATAC-seq profile。这个分析使得我们检查每个特定TF的染色质可接近性。此外，通过比较两个ATAC-seq文件(cDC1细胞和pDC细胞)，我们可以了解两个细胞中TF结合的变化。

下面这行代码将读取上面第四步里生成的motif 匹配文件和原始的bam文件。其中代码里的-bc的意思是使用bias矫正信号（目前只支持ATAC-seq数据）。-nc是cpu数量：

#这一行代码是教程里的原代码
$ 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

我在服务器里的使用的代码：

#你可以看到这里我调用了128G的内存，之前也试过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的内存来运行，还要1个小时。。。

运行后，产生一个新的文件夹名为“cDC1_pDC”。在这个文件夹里，有一个scatter plot图，展示的是在两种细胞之间TF动力学活性。y轴代表TF的不同的activity，红色标记的是有显著性差异的TF(p<0.05)。x轴是随机数（for jittering purposes）：

这个图上面的红色点代表的TF是在pDC细胞里结合比较强的，下面的红点代表的TF是在cDC1细胞里结合比较强的。这个图有点丑，另外标记也都重叠在一起了，你可以自己用R重新画个图（因为是示例数据，就不单独美化这个散点图了），这个点图对应的txt文件是生成的differential_statistics.txt：

这个txt打开后是这样的：

另外，在上面文件夹截图里显示的“Lineplots”里，你还可以找到所有TF的曲线图，比如挑两个有显著差异的：