ATAC-seq信息分析流程主要分為以下幾個部分：數(shù)據(jù)質(zhì)控、序列比對、峰檢測、motif分析、峰注釋、富集分析，下面將對各部分內(nèi)容進(jìn)行展開講解。

一、測序數(shù)據(jù)過濾與質(zhì)量評估

下機(jī)數(shù)據(jù)經(jīng)過過濾去除接頭含量過高或低質(zhì)量的reads，得到clean reads用于后續(xù)分析。常見的trim軟件有Trimmomatic、Skewer、fastp等。fastp是一款比較新的軟件，使用時可以用--adapter_sequence/--adapter_sequence_r2參數(shù)傳入接頭序列，也可以不填這兩個參數(shù)，軟件會自動識別接頭并進(jìn)行剪切。如：

fastp \

--in1 A1_1.fq.gz \ # read1原始fq文件

--out1 A1_clean_1.fq.gz \ # read1過濾后輸出的fq文件

--in2 A1_2.fq.gz ?\ # read2原始fq文件

--out2 A1_clean_2.fq.gz \ # read2過濾后輸出的fq文件

--cut_tail ?\ #從3’端向5’端滑窗，如果窗口內(nèi)堿基的平均質(zhì)量值小于設(shè)定閾值，則剪切

--cut_tail_window_size=1 \ #窗口大小

--cut_tail_mean_quality=30 \ #cut_tail參數(shù)對應(yīng)的平均質(zhì)量閾值

--average_qual=30 \ #如果一條read的堿基平均質(zhì)量值小于該值即會被舍棄

--length_required=20 ?\ #經(jīng)過剪切后的reads長度如果小于該值會被舍棄

fastp軟件的詳細(xì)使用方法可參考：https://github.com/OpenGene/fastp。fastp軟件對于trim結(jié)果會生成網(wǎng)頁版的報告，可參考官網(wǎng)示例http://opengene.org/fastp/fastp.html和http://opengene.org/fastp/fastp.json，也可以用FastQC軟件對trim前后的數(shù)據(jù)質(zhì)量進(jìn)行評估，F(xiàn)astQC軟件會對單端的數(shù)據(jù)給出結(jié)果，如果是PE測序需要分別運(yùn)行兩次來評估read1和read2的數(shù)據(jù)質(zhì)量。

如：

fastqc A1_1.fq.gz

fastqc A1_2.fq.gz

FastQC會對reads從堿基質(zhì)量、接頭含量、N含量、高重復(fù)序列等多個方面對reads質(zhì)量進(jìn)行評估，生成詳細(xì)的網(wǎng)頁版報告，可參考官網(wǎng)示例：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/good_sequence_short_fastqc.html

二、序列比對

經(jīng)過trim得到的reads可以使用BWA、bowtie2等軟件進(jìn)行比對。首先需要確定參考基因組fa文件，對fa文件建立索引。不同的軟件有各自建立索引的命令，BWA軟件可以參考如下方式建立索引：

bwa index genome.fa

建立好索引后即可開始比對，ATAC-seq推薦使用mem算法，輸出文件經(jīng)samtools排序輸出bam：

bwa mem genome.fa ?A1_clean_1.fq.gz A1_clean_2.fq.gz

| samtools sort -O bam -T A1 > A1.bam

值得注意的是，在實(shí)驗(yàn)過程中質(zhì)體并不能完全去除，因此會有部分reads比對到質(zhì)體序列上，需要去除比對到質(zhì)體上的序列，去除質(zhì)體序列可以通過samtools提取，具體方法如下：首先將不含質(zhì)體的染色體名稱寫到一個chrlist文件中，一條染色體的名稱寫成一行，然后執(zhí)行如下命令即可得到去除質(zhì)體的bam

samtools view -b A1.bam $chrlist > A1.del_MT_PT.bam

用于后續(xù)分析的reads需要時唯一比對且去重復(fù)的，bwa比對結(jié)果可以通過MAPQ值來提取唯一比對reads，可以用picard、sambamba等軟件去除dup，最終得到唯一比對且去重復(fù)的bam文件。

三、reads在染色體上分布的可視化

比對后得到的bam文件可以轉(zhuǎn)化為bigWig（bw）格式，通過可視化軟件進(jìn)行展示。deeptools軟件可以實(shí)現(xiàn)bw格式轉(zhuǎn)化和可視化展示。首先需要在linux環(huán)境中安裝deeptools軟件，可以用以下命令實(shí)現(xiàn)bam向bw格式的轉(zhuǎn)換：

bamCoverage -b A1.bam -o A1.bw

此外，可以使用deeptools軟件展示reads在特定區(qū)域的分布，如：

computeMatrix reference-point ??\ # reference-pioint表示計算一個參照點(diǎn)附近的reads分布，與之相對的是scale-regions，計算一個區(qū)域附近的reads分布

--referencePoint TSS ??\#以輸入的bed文件的起始位置作為參照點(diǎn)

-S ?A1.bw \ #可以是一個或多個bw文件

-R ?gene.bed \ #基因組位置文件

-b 3000 ??\ #計算邊界為參考點(diǎn)上游3000bp

-a 3000 ??\ #計算邊界為參考點(diǎn)下游3000bp，與-b合起來就是繪制參考點(diǎn)上下游3000bp以內(nèi)的reads分布

-o ?A1.matrix.mat.gz \ #輸出作圖數(shù)據(jù)名稱

#圖形繪制

plotHeatmap \

-m ?new_A1.matrix.mat.gz \ #上一步生成的作圖數(shù)據(jù)

-out A1.pdf \ # 輸出圖片名稱

繪圖結(jié)果展示：

reads在TSS附近的分布

四、Peak calling

MACS2能夠檢測DNA片斷的富集區(qū)域，是ATAC-seq數(shù)據(jù)call peak的主流軟件。峰檢出的原理如下：首先將所有的reads都向3'方向延伸插入片段長度，然后將基因組進(jìn)行滑窗，計算該窗口的dynamic λ，λ的計算公式為：λlocal = λBG（λBG是指背景區(qū)域上的reads數(shù)目），然后利用泊松分布模型的公式計算該窗口的顯著性P值，最后對每一個窗口的顯著性P值進(jìn)行FDR校正。默認(rèn)校正后的P值（即qvalue）小于或者等于0.05的區(qū)域?yàn)閜eak區(qū)域。需要現(xiàn)在linux環(huán)境中安裝macs2軟件，然后執(zhí)行以下命令：

macs2 callpeak \

-t A1.uni.dedup.bam \ #bam文件

-n A1 \ # 輸出文件前綴名

--shift -100 \ #extsize的一半乘以-1

--extsize 200 \ #一般是核小體大小

--call-summits #檢測峰頂信息

注：以上參數(shù)參考文獻(xiàn)（Jie Wang，et.al.2018.“ATAC-Seq analysis reveals a widespread decrease of chromatin accessibility in age-related macular degeneration.”Nature Communications）

五、motif分析

ATAC分析得到的peak是染色質(zhì)上的開放區(qū)域，這些染色質(zhì)開放區(qū)域常常預(yù)示著轉(zhuǎn)錄因子的結(jié)合，因此對peak區(qū)域進(jìn)行motif分析很有意義。常見的motif分析軟件有homer和MEME。以homer軟件為例，首先在linux環(huán)境中安裝homer，然后用以下命令進(jìn)行motif分析：

findMotifsGenome.pl \

A1_peaks.bed \ #用于進(jìn)行motif分析的bed文件

genome.fa ?\ #參考基因組fa文件

A1 ?\ #輸出文件前綴

-size ?given \ #使用給定的bed區(qū)域位置進(jìn)行分析，如果填-size -100,50則是用給定bed中間位置的上游100bp到下游50bp的區(qū)域進(jìn)行分析

homer分析motif的原理及結(jié)果參見：http://homer.ucsd.edu/homer/motif/index.html

根據(jù)motif與已知轉(zhuǎn)錄因子的富集情況可以繪制氣泡圖，從而可以看到樣本與已知轉(zhuǎn)錄因子的富集顯著性。

六、差異分析

差異peak代表著比較組合染色質(zhì)開放性有差異的位點(diǎn)，ChIP-seq和ATAC-seq都可以用DiffBind進(jìn)行差異分析。DiffBind通過可以通過bam文件和peak的bed文件計算出peak區(qū)域標(biāo)準(zhǔn)化的readcount，可以選擇edgeR、DESeq2等模型進(jìn)行差異分析。

七、峰注釋

在科研分析中我們往往需要將peak區(qū)域與基因聯(lián)系起來，也就是通過對peak進(jìn)行注釋找到peak相關(guān)基因。常見的peak注釋軟件有ChIPseeker、homer、PeakAnnotator等。以ChIPseeker為例，需要在R中安裝ChIPseeker包和GenomicFeatures包，然后就可以進(jìn)行分析了。

library(ChIPseeker)

library(GenomicFeatures)

txdb<- makeTxDbFromGFF(‘gene.gtf’)#生成txdb對象，如果研究物種沒有已知的TxDb,可以用GenomicFeatures中的函數(shù)生成

peakfile <-readPeakFile(‘A1_peaks.narrowPeak’)#導(dǎo)入需要注釋的peak文件

peakAnno <- annotatePeak(peakfile,tssRegion=c(-2000, 2000), TxDb=txdb)

# 用peak文件和txdb進(jìn)行peak注釋，這里可以通過tssRegion定義TSS區(qū)域的區(qū)間

對于peak注釋的結(jié)果，也可以進(jìn)行可視化展示，如：

p <- plotAnnoPie(peakAnno)

八、富集分析

通過注釋得到的peak相關(guān)基因可以使用goseq、topGO等R包進(jìn)行GO富集分析，用kobas進(jìn)行kegg富集分析，也可以使用DAVID在線工具來完成富集分析?？梢酝ㄟ^挑選感興趣的GO term或pathway進(jìn)一步篩選候選基因。