Harvard FAS Informatics出品的ATAC-seq测序指南

github链接：harvardinformatics/ATAC-seq

参考文献：ATAC-seq: A Method for Assaying Chromatin Accessibility Genome-Wide

ATAC-seq测序推荐双端测序，推荐数据量：15G (150PE)

分析流程

1. FastQC查看测序数据质量

ls *.fastq.gz | xargs -I {} -P 20 -n 1 echo ~/biosoft/fastqc/FastQC/fastqc {}  -o ~/disks/diskd/data/project/lungCancer/gse79209/fastq/qc -t 30

2. 去除接头序列

（1）如果知道接头序列的话，用cutadapt去除接头序列

cutadapt \
            --cores=CORES \ 
            -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
            -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \
            -o 101N.trimmed_1.fq.gz -p 101N.trimmed_2.fq.gz \
            /mnt/data/data/lncRNA/101N/101N_1.fq.gz /mnt/data/data/lncRNA/101N/101N_2.fq.gz 

（2）如果不知道接头序列的话，使用NGmerge去除接头序列

NGmerge -a -1 <sample>.R1.fastq.gz -2 <sample>.R2.fastq.gz -o <sample>

3. 序列比对

（1）构建bowtie2索引

（2）序列比对

推荐使用bowtie2进行序列比对（alignment）。常用的参数有以下几个：

• ==-X== 最大DNA片段长度（默认是500bp）
• ==--very-sensitive== 默认参数是--sensitive，使用--very-sensitive可以取得更好的序列比对效果
• ==-p== 多核运算

默认输出格式问SAM文件，包含了每个序列的比对信息。SAM文件可以压缩为BAM文件，并使用SAMtools进行排序。最佳的分析流程为：

bowtie2 --very-sensitive -x <genomeIndexName> -1 <sample>_1.fastq.gz -2 <sample>_2.fastq.gz | samtools view -u - | samtools sort - > <BAM>

（3） 比对序列调整

a. 线粒体reads

ATAC-seq数据一般含有相当比例的线粒体DNA，详见该讨论。可以使用CRISPR技术去除线粒体基因组污染。也可以参考最近发表的Omni-ATAC方法用去污剂去除线粒体DNA，这种方法对于多数研究人员可能更易于操作，但是不要按照他们的分析流程来分析数据。

不管使用何种实验方案，得到的测序数据中多多少少都会有线粒体DNA。因为线粒体基因组中没有ATAC-seq感兴趣的峰，这些数据会使后续分析变得复杂，因此，推荐在下一步分析之前去除线粒体DNA序列。可以通过下面两种方法中任一种实现：

• 在比对之前，将线粒体DNA从比对使用的参考基因组中去除。该方法的缺点使比对数会看起来比较差，因为所有的线粒体DNA序列都将被记为未匹配序列。
• 在进行比对之后去除线粒体DNA。可以使用Harvard FAS Informatics编写的python脚本去除线粒体RNA序列removeChrom。

samtools view -h <inBAM> | removeChrom - - chrM | samtools view -b - > <outBAM> #该程序在Harvard超算运行代码，可以根据自己服务器进行更改

# 可能的解决方案如下：
samtools view -h <inBAM> | python /path/to/removeChrom.py - - chrM | samtools view -b - > <outBAM> 

b. PRC重复序列

PCR重复序列是完全一致的reads。

java -jar $PICARD_TOOLS_HOME/picard.jar MarkDuplicates I=<inBAM> O=<outBAM> M=dups.txt REMOVE_DUPLICATES=true

c. 非特异性的序列

二代测序得到的短序列中，有一些序列会匹配到基因组的多个区域。一些学者在分析数据时会通过samtools view中的-q参数将这些非特异性的序列去除。对于这种非特异性的序列，bowtie2或bwa只会报告一个比对的结果，同时会给一个低匹配质量的（mapping quality, MAPQ）评分，匹配质量定义为: -10 * log10Pr{mapping position is wrong}。可以通过下面命令来去除MAPQ<10的序列。

samtools view -b  -q 10  <inBAM>  >  <outBAM>

4. Peak calling（找峰）

Model-based Analysis of ChIP-Seq (MACS2)是一款用于检测基因富集区域的软件。尽管是为ChIP-seq设计的，但也适用于ATAC-seq测序及其他有小峰的基因组富集分析。MACS2软件的主程序是callpeak。
可以将前期处理好的比对序列文件作为MACS2的输入文件。然而，一定要记住：比对上的序列只是ATAC技术产生的DNA片段的一部分。因此，必须考虑如何用MACS2来解读这些比对序列。

（1）用于分析的比对序列

对于双端测序书来说，比对序列主要分两种类型：成对匹配的和单一匹配的序列。当使用MACS2分析数据时，应事先决定哪种比对序列可用于后续分析。主要有以下三个选择：

• 只分析成对匹配的序列，单是忽略R2序列，并将R1视为单一匹配序列。这是默认的选项(-f AUTO)。

• 用-f BAMPE选项，只分析成对匹配的序列。相比第一种方法，更推荐使用这种方法来处理成对匹配的序列。

• 分析所有的序列。ATAC-seq分析流程中SAMtoBED脚本可以实现。这个脚本可以将比对的序列转换成BED区间，对于成对匹配的序列，直接转换，对于单一匹配的序列，有以下四个选项：忽略他们；保持原状；将他们延申到任意长度（与MACS2的--extsize选项类似）；或将他们延申至到完美匹配序列的平均长度（-x选项）。示例如下：

samtools view -h  <BAM>  |  SAMtoBED  -i -  -o <BED>  -x  -v

# 如果在自己服务器上运行，可以相应调整代码
samtools view -h <BAM> | python /path/to/SAMtoBED.py -i - -o <BED> -x -v

输出为BED文件，可以用MACS2 -f BEDPE命令处理。

• 注意：此处不能使用BEDTools中的bamtobed软件来处理，因为它的输出文件不是MACS2可以处理的标准文件。

（2）其他参数

除了前文中描述的一些参数，MACS2还有一些其他参数可以设置。下面是一些可以考虑的参数：

• -n <str> 样本名称。输出文件会以此为前缀命名
• -g <int> 有效基因组大小。小于实际基因组大小。默认选项是hs，对应的基因组大小事2.7e9。
• -q <fload> 找峰的最小q值（矫正p值或FDR），默认值为0.05。降低q值会减少找到的峰的个数，但是会提高可信度。
• --keep-dup <arg> 如何处理PCR重复序列。默认值是--keep-dup 1，去除所有可能的PCR重复。如果PCR重复序列已经用其他软件去除了，如Picards的MarkDuplicates选项，将该参数设置为--keep-dup all。
• --max-gap <int> 将显著区域聚类的最大间隙，默认值是50bp(v2.1.2_dev only)。
• --min-length <int> 峰的最小长度，默认值是1000bp（(v2.1.2_dev only)）。
  ==注意==：MACS2不支持多核运算，因此只能使用一个核.
  MACS2分析线虫ATAC-seq数据，使用所有序列的完整命令大概长这样：

macs2 callpeak  -t <BED>  -f BEDPE  -n NAME  -g ce  --keep-dup all

（3） 输出文件

标准的macs2 callpeak程序有3个输出文件。如果有-n NAME参数，输出文件分布为：NAME_peaks.xls，NAME_peaks.narrowPeak,和NAME_summits.bed。最有用的文件是NAME_peaks.narrowPeak，它是一个bed文档，包含了所有峰的基因组坐标和不同的统计值（倍数改变、p值和q值等）。

5. 后续步骤

一旦完成一组样本的MACS2分析之后，可根据实验设计进行一些后续分析。

（1）可视化

可以将peak文件在基因组中可视化。模式生物的ATAC-seq数据，peak文件（NAME_peaks.narrowPeak）可以直接上传到UCSC genome browser中进行查看。如果peak文件没有表头的话，可以将下面的内容添加至首行：track type=narrowPeak。
另外，可以使用IGV进行可视化。peak文件可以直接通过File -->Load from File选项上传。如果要对BAM文件可视化，需要用samtools对BAM文件进行排序并构建索引。

（2）比较峰文件

可以使用BEDTools来比较一组peak文件的异同。例如：bedtools intersect可以比较两个peak文件中相同的峰；寻找两个峰文件中差异的峰，如实验组和对照组，可以通过bedtools subtract命令实现。

（3）注释

ChIPseeker最早是用于注释ChIP-seq数据的，但是也适用于ATAC-seq数据的注释。ChIPseeker可以对基因的位置和其他特征进行注释，详细使用方法见ChIPseeker的使用指南。

（4）寻找motif

HOMER可用于寻找motif。可以将peak文件作为HOMER文件的输入，来检测已知的motif和新的motif。

参考文献：

Andrews S. (2010). FastQC: a quality control tool for high throughput sequence data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc

Buenrostro JD, Giresi PG, Zaba LC, Chang HY, Greenleaf WJ. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. Nat Methods. 2013 Dec;10(12):1213-8.

Buenrostro JD, Wu B, Chang HY, Greenleaf WJ. ATAC-seq: A Method for Assaying Chromatin Accessibility Genome-Wide. Curr Protoc Mol Biol. 2015 Jan 5;109:21.29.1-9.

Corces MR, Trevino AE, Hamilton EG, Greenside PG, Sinnott-Armstrong NA, Vesuna S, Satpathy AT, Rubin AJ, Montine KS, Wu B, Kathiria A, Cho SW, Mumbach MR, Carter AC, Kasowski M, Orloff LA, Risca VI, Kundaje A, Khavari PA, Montine TJ, Greenleaf WJ, Chang HY. An improved ATAC-seq protocol reduces background and enables interrogation of frozen tissues. Nat Methods. 2017 Oct;14(10):959-962.

Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell. 2010 May 28;38(4):576-89.

Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods. 2012 Mar 4;9(4):357-9.

Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009 Aug 15;25(16):2078-9.

Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal. 2011;17:10-2.

Montefiori L, Hernandez L, Zhang Z, Gilad Y, Ober C, Crawford G, Nobrega M, Jo Sakabe N. Reducing mitochondrial reads in ATAC-seq using CRISPR/Cas9. Sci Rep. 2017 May 26;7(1):2451.

Quinlan AR. BEDTools: The Swiss-Army Tool for Genome Feature Analysis. Curr Protoc Bioinformatics. 2014 Sep 8;47:11.12.1-34.

Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics. 2015 Jul 15;31(14):2382-3.