ATAC-Seq 分析流程

    本文分享我的 ATAC-Seq 分析流程，為了控制篇幅，關鍵的 Peak calling 步驟會詳細介紹，而 Motif 分析等一些步驟則拆分到以後文章，在未來幾天推送。雖然如此，篇幅還是挺長的，建議一次讀完。

    簡單來說，基因的轉錄需要 DNA 從纏繞在組蛋白的緊密狀態釋放，成為鬆散的開放區域，這樣轉錄複合體才能結合上去。ATAC-Seq 技術能高通量分析開放的染色質區域，揭示基因的轉錄狀態。

只有開放的染色質區域 Tn5 才能夠攜帶測序接頭結合併將接頭轉座上去，PCR 擴增後就可以進行測序。

基本的實驗要求
    一、實驗最好有生物學重複，無論什麼技術都會有一定的隨機誤差。理論上還需要一個實驗的技術對照，類似於 Chip-Seq 實驗的 Input 對照。即用普通的破碎法(如超聲)打斷 DNA 然後同樣加接頭測序，作為與轉座法的對照，確定基因組哪些區域是更難進行測序的。一般不做這個對照實驗。
    二、依據不同的物種和估計的染色質開放程度改變，一般人種要求 50 M 比對上 Reads，TF footprinting 分析要 200 M 以上 Reads. 合格的實驗唯一比對率應該在 80% 以上。

For inferring differences in open chromatin within human samples we generally use >50M mapped reads and for transcription factor foot-printing we use >200M mapped reads (Neph et al., 2012).

本教程主要的分析步驟是質控、比對、Peak calling、差異分析、注釋、Motif finding和可視化。下面一一給出步驟示範。

fastq 質控
    先用 fastqc 檢查質量，fastp 進行過濾和去接頭。有人認為除了去接頭，不額外進行其他截短之類操作，防止影響後續去重。也有人認為根據染色質開放區域實際片段大多在 150 bp 以下，如果是太長的測序，應該進行截短，如統一截短到 75bp.

典型的 ATAC-Seq 片段分布，小片段佔比高

    在我這次數據我進行了測試，統一截短到 75bp 在 fastqc 對單端檢測重複率下降 1% 左右，用 fastp 對雙端的檢測重複率不變。同時使用 Bowtie2 的比對率也沒變化或只有輕微改變(< 1%)。也有可能因為過短的插入片段因為導致 Reads 3' 端太多 G 鹼基且低質量，早被 fastp 過濾了，所以無法體現出影響。總的來說，自己的數據 2 方案都試試看看情況決定。

比對和過濾
    用 BWA, bowtie2 等常規軟體比對。這裡 -k 設置一條 reads 至多報告的匹配次數，默認是多個匹配報告最佳；-X 設置插入片段最大長度，默認 500. 同時用 samtools 排序輸出。
    哈佛教程推薦用 -k 10 是因為他們的 Genrich 支持利用多比對信息。如果用其他軟體做 Peak calling 還是默認不用 -k 操作。使用 Genrich 要求用 -n 進行排序。

  bowtie2 --very-sensitive -k 10 -X 2000 --threads 4 -x ${GRCh38} \ 
  -1 ${CleanDir}/${Sample}_R1.fq.gz -2 ${CleanDir}/${Sample}_R2.fq.gz \ 
  -S ${BamDir}/${Sample}.sam

  samtools view -b -h ${BamDir}/${Sample}.sam | samtool sort -n -O BAM \ 
  -o ${BamDir}/${Sample}_sorted.bam

    用 picard 對 PCR 重複進行標記。後續用 Genrich 的話不需要標記，Genrich 自帶移除 PCR 重複參數。

gatk MarkDuplicates -I ${BamDir}/${Sample}_sorted.bam -M ${BamDir}/${Sample}_dup_matrics.txt -O ${BamDir}/${Sample}_MD.bam

    順便可以看看插入片段分布情況。插入片段分布應該是隨著長度下降的，但是在 <100,200,400,600bp 處有峰，分別對應著 NER(nucleosome-free regions) 和 單、雙、三核小體。因為一個核小體纏繞的 DNA 長度約為 200bp.

 gatk CollectInsertSizeMetrics -H ${BamDir}/${Sample}_InsertSize.pdf -I ${BamDir}/${Sample}_MD.bam -O ${BamDir}/${Sample}_InsertSize.txt

測序插入片段形成典型的計數分布原因

一個樣本的實際片段分布案例

    線粒體 DNA 與核基因組不同，它沒有纏繞組蛋白形成核小體，因此都是開放的。一般來說會導致測序結果裡有不少線粒體信號，需要移除。同時移除基因組 Blacklist 區域。
    移除後過濾低質量的比對， Bowtie2 設置比對質量大於 30. 要注意不同比對軟體對 MAPQ 值設置不同。同樣，如果使用 Genrich 進行 Peak calling 那麼這些步驟不需要進行，他有相應參數實現。

  bedtools intersect -nonamecheck -v -a ${BamDir}/${Sample}_MD.bam -b ${Blacklist} > ${BamDir}/${Sample}_RMBL.bam
  samtools view -h -F 1804 -q 30 -O SAM ${BamDir}/${Sample}_RMBL.bam | grep -v -e "^chrM" | samtools sort -O BAM -o ${BamDir}/${Sample}_Filtered.bam 
  samtools index ${BamDir}/${Sample}_Filtered.bam 

移動 Reads 位置
    如下圖所示 Tn5 酶結合和切割的位點是開放染色質區域，所以我們分析出來的峰中心應該在這個切割位點，也就是 reads 的 5' 端。

Tn5 結合的地方就是測序讀長開始的地方

    如下面示意圖，Tn5 建庫時會插入 9bp 的序列，所以比對到正鏈（+）的 reads 移動 +4 bp, 比對負鏈（-）的 reads 移動 -5 bp.

Tn5 建庫示意圖

   這個移動可以用 deeptools 的 alignmentSieve 命令完成。

alignmentSieve --numberOfProcessors 8 --ATACshift -b ${BamDir}/${Sample}_Filtered.bam -o ${BamDir}/${Sample}_Shift.bam

    同樣，Genrich 會完成這個移動。

FRiP 計算
    Fraction of reads in peaks (FRiP) 指落入峰區域的 reads 佔比，應該要高於 0.3 較好。這個指標也不用硬性要求，所以計算也不需要那麼嚴謹，簡單化處理。

samtools view -c ${BamDir}/KO1_ATAC_ShiftSort.bam
bedtools intersect -u -a ${BamDir}/KO1_ATAC_ShiftSort.bam -b ${PeakDir}/KO1_ATAC_peaks.narrowPeak | samtools view -c

得到數值分別是 24634649, 1217509 除一下得出 0.049 的 FRiP.

Peak calling
    這裡分別用 Genrich 和 MACS2 進行 Peak calling. Genrich 更加方便，把前面那些過濾等步驟都集合在一起了，同時他可以利用生物學重複樣本，計算出顯著的開放區域。不需要手動去分開每個樣本進行 Peak calling 然後合併。MACS2 是很常用的 Chip-Seq 軟體。

Genrich
    這裡 -j 表示 ATAC 模式，一定要指定。多個生物學重複樣本用逗號分隔。-c 參數是前文介紹的實驗技術對照樣本，一般 ATAC-Seq 沒有，所以不指定。
    Genrich 集合了移除 PCR 重複、線粒體區域、Blacklist 區域等功能，-r 移除 PCR 重複，-e chrM 移除線粒體區域，-E bed 移除 bed 文件區域。

Genrich -t ${BamDir}/KO1_ATAC.bam,${BamDir}/KO2_ATAC.bam -o ${GenrichDir}/KO.narrowPeak -f ${GenrichDir}/KO_pq.bed -k ${GenrichDir}/KO_pileup_p.bed -b ${GenrichDir}/KO_reads.bed -r -e chrM -E ${Blacklist} -m 30 -j

Genrich -t ${BamDir}/WT1_ATAC.bam,${BamDir}/WT2_ATAC.bam -o ${GenrichDir}/WT.narrowPeak -f ${GenrichDir}/WT_pq.bed -k ${GenrichDir}/WT_pileup_p.bed -b ${GenrichDir}/WT_reads.bed -r -e chrM -E ${Blacklist} -m 30 -j

    參數 -m 過濾比對 MAPQ 值，Bowtie2 比對往往取 30 以上認為是好的比對，不同軟體 MAPQ 值設定不相同，不可以等同使用。

-b 生成的 fragments 文件需要用 tabix 建立索引才能給 IGV 使用。

bedtools sort -i KO_reads.bed | bgzip > KO_reads_sorted.bed.gz
tabix -p bed KO_reads_sorted.bed.gz

bedtools sort -i WT_reads.bed | bgzip > WT_reads_sorted.bed.gz
tabix -p bed WT_reads_sorted.bed.gz

    參數 -o 輸出 narrowPeak 格式結果文件，表示每個峰區域的範圍和強度。可以在開頭添加 track type=narrowPeak 後用 USUC genome browser 打開。narrowPeak 格式也是一種 bed 格式。

chr1    9946    10511   peak_0  1000    .       7480.733887     19.876038       -1      134
chr1    180045  180182  peak_1  1000    .       403.740967      6.181179        -1      56
chr1    180465  181043  peak_2  1000    .       4381.033203     16.855917       -1      330

    參數 -f 輸出的類 bedgraph 格式文件保存匯總的 p,q 值。如果有樣本重複，會給每個重複的 p 值和合併後的 p 值，以及是否顯著差異。
兩樣本重複的舉例

chr     start   end     -log(p)_0       -log(p)_1       -log(p)_comb    signif
chr1    0       9944    0.000000        0.000000        0.000000
chr1    9944    9945    0.252055        0.574477        0.363661
chr1    9945    9946    1.027748        1.432146        1.636151
chr1    9946    9947    1.394456        2.120604        2.556318        *
chr1    9947    9950    2.183553        2.826582        3.911967        *

    參數 -k 輸出類 bedgraph 格式文件，會用一行注釋開始表示是哪個樣本，背景對照是什麼。把這個文件修改成 bedgraph 格式，就可以用 genome browser 去可視化。下面是注釋行。

$ grep "experimental file" KO_pileup_p.bed
# experimental file: /Example/KO1_ATAC.bam; control file: NA
# experimental file: /Example/KO2_ATAC.bam; control file: NA

下面是開頭 5 行內容。

# experimental file: /Example/KO1_ATAC.bam; control file: NA
chr     start   end     experimental    control -log(p)
chr1    0       9944    0.000000        0.900041        0.000000
chr1    9944    9945    0.500000        0.900041        0.252055
chr1    9945    9946    2.000000        0.900041        1.027748

下面是第五列背景數值總結。

$ awk '{print$5}' KO_pileup_p.bed | sort | uniq
0.000000
0.866656
0.900041
control

$ awk '{print$5}' WT_pileup_p.bed | sort | uniq
0.000000
2.114438
2.238139
control

MACS2
    不像 Genrich 直接採用生物學重複數據合併分析，MACS2 進行一個個樣本分析，後續再合併峰區域。

    ATAC-seq 要用 --nomodel 參數，-g 要根據自己測序物種選擇。

macs2 callpeak -t ${BamDir}/${Sample}_Shift.bam -n ${Sample} --outdir ${MACS2Dir} \ 
-f BAMPE -g hs --nomodel --keep-dup all --cutoff-analysis --bdg

    MACS2 peak calling 結果裡 \*_peaks.xls 是可以用 Excel 打開的 peak 結果文件，注意這裡坐標以 1 開始而不是 0. \*_peaks.narrowPeak 是 narrowPeak 格式結果；*_summits.bed 是包含頂峰（Peak summit）的位置信息的 BED 格式文件，適合用於 motif 的分析。

    對於有重複的實驗來說，也許想要將得到的 peak 進行合併，可以用 bedtools merge 命令。下面的命令順便去除了不想要的部分結果，只保留常染色體和性染色體。

cat ${Macs2Dir}/WT1_ATAC_peaks.narrowPeak ${Macs2Dir}/WT2_ATAC_peaks.narrowPeak \ 
| awk -v FS="\t" -v OFS="\t" 'length($1)<=5' | sort -k1,1 -k2,2n > ${Macs2Dir}/WT_SortPeaks.bed
bedtools merge -i ${Macs2Dir}/WT_SortPeaks.bed  > ${Macs2Dir}/WT_MergePeaks.bed

    這樣合併等於說峰區域取併集，是沒辦法的簡單化方法，個人不太建議。不過如果使用 DiffBind 進行差異的峰分析，他也會採用這樣的合併方法。所以像 Genrich 在 peak calling 階段直接考慮實驗重複是更為合適的。

差異分析
    差異分析可以採用 Peak 結果直接用 bedtools 簡單化處理，得到多組的峰共同或特有區域。個人覺得下面的命令還是多調整 -f 的值試試，默認的也許不適合。

# 取共同交集
bedtools intersect -f 0.2 -F 0.2 -e -a ${GenrichDir}/KO.narrowPeak -b ${GenrichDir}/WT.narrowPeak > ${GenrichDir}/CommonPeak.bed
# KO 特異，但是有重疊則移除整條記錄
bedtools subtract -A -f 0.3 -a ${GenrichDir}/KO.narrowPeak -b ${GenrichDir}/WT.narrowPeak > ${GenrichDir}/KO_Specific.bed
# WT 特異，有重疊則移除整條記錄
bedtools subtract -A -f 0.3 -a ${GenrichDir}/WT.narrowPeak -b ${GenrichDir}/KO.narrowPeak > ${GenrichDir}/WT_Specific.bed

    下面分別是 CommonPeak.bed 輸出和原來 KO.narrowPeak 的峰區域數據。除了區域位置信息，其餘信息會用 -a 的文件信息。也可以把後面那些列移除掉，免得製造混亂。

# Common
chr1    9957    10511   peak_0  1000    .       7480.733887     19.876038       -1      134
chr1    180691  181041  peak_2  1000    .       4381.033203     16.855917       -1      330
chr1    199594  199987  peak_3  1000    .       487.511292      5.329252    

# KO 
chr1    9946    10511   peak_0  1000    .       7480.733887     19.876038       -1      134
chr1    180045  180182  peak_1  1000    .       403.740967      6.181179        -1      56
chr1    180465  181043  peak_2  1000    .       4381.033203     16.855917       -1      330

    直接取峰區域的交集差集未免過於簡單化處理，另一個方法是用 R 包 DiffBind. DiffBind 是借鑑了 RNA-seq 差異分析。首先將所有樣本的峰區域合併，合併方法是用 bedtools merge 所以是所有的樣本峰區域的併集。有了這個合併的峰區域，那麼可以對每個樣本計算裡面每個峰區域的 Reads 數，這就像是 RNA-seq 分析每個基因/轉錄本的 Reads 數一樣，然後調用 DESeq2/edgeR 來分析有顯著差異的峰區域。
    DiffBind 我覺得主要存在兩個問題，一是用取所有峰區域併集作為總的峰區域，二是無法確定 RNA-seq 分析軟體的假設前提是否在 ATAC-seq 成立，比如說 RNA-seq 是不會大量低豐度基因的，但是 ATAC-seq 是會大量低豐度的峰 Reads.

ATAC-Seq 與 RNA-Seq 相比存在更多低豐度數據

    總而言之方法還不算完美，只能自己拿數據去試了，結合課題看是否合適。
DiffBind 教程將在未來推文推送。

注釋
    沒有注釋就只有峰位置信息，對研究人員來說不太可用。注釋可以用 ChIPseeker 包在 R 完成，還可以順便出一些圖。ChIPseeker 的使用參照官方文檔。一般來說只看常染色體和性染色體，其餘的 contig 可以移除。

awk -v FS="\t" -v OFS="\t" 'length($1)<=5' PeakResult > PeakFilter

同時 ChIPseeker 還可以做通路富集分析。

Motif finding
    只依賴於對峰區域的注釋是不夠直接的，分析轉錄因子結合的 Motif 能夠更直接反映基因轉錄調控的變化。Motif 分析可以用 HOMMER 或 MEME-Chip 完成。這方面內容將在未來推文推送。

TF footprints
    另一種轉錄因子相關的分析叫 TF footprint 分析。轉錄因子結合到開放染色質導致 Tn5 無法剪切該區域，會在覆蓋度高的開放區域（Peak）中出現小範圍的覆蓋度下降，因此這個下降的小區域稱為 TF footprint.

可視化
    ChIPseeker 有不少總體峰區域分析的可視化圖片，這裡不再做介紹。下面介紹 2 種生成 Reads 覆蓋信號 bw 文件的方法，然後用 pyGenomeTracks 生成特定區域的 Reads 覆蓋信號圖。

手動轉換到 bw 文件
    警告：這種方法不適用於建庫深度有差異的情況，看看就好。具體原因在將在未來推文推送。

    Genrich -k 參數輸出文件是類 Bedgraph 格式的，手動轉換到 Bedgraph 就可以用 UCSC genome broswer 瀏覽了。
    我這裡實驗有重複，所以多個 Bam 文件給 Genrich 同時輸入，此時 -k 產生的文件會是 2 個 Bam 文件數據，用一個注釋行隔開。所以首先我需要把這 2 數據分開，如果是每個樣本單獨跑 Genrich 那跳過此步。

# 查看注釋行的位置
$ grep -n "#" KO_pileup_p.bed
1:# experimental file: /Example/KO1_ATAC.bam; control file: NA
54091955:# experimental file: /Example/KO2_ATAC.bam; control file: NA

# 分別用 head, tail 取需要的數據
$ head -n 54091954 KO_pileup_p.bed > KO1_pileup_p.bed
$ tail -n [phone] KO_pileup_p.bed > KO2_pileup_p.bed

# 檢查一下行數對不對
$ wc -l KO*pileup_p.bed
   54091954 KO1_pileup_p.bed
   52204177 KO2_pileup_p.bed
  106296131 KO_pileup_p.bed
  212592262 total

# WT 樣本進行同樣處理

    分割後將前 2 行即注釋行和表頭移除，移除後用 cut 命令取前 4 列，輸出到 TAB 分隔文件。順便進行排序。

cut -f 1,2,3,4 KO1_pileup_p.bed | LC_COLLATE=C sort -k1,1 -k2,2n > KO1_Cut.bed

# 其餘文件相同處理

    然後在文件開頭添加一行 "track type=bedGraph", 最後用 bedGraphToBigWig 轉換到 bigwig 文件。bedGraphToBigWig 在 http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/ 下載。

bedGraphToBigWig ${GenrichDir}/${Sample}_Cut.bed ${ChromSize} ${GenrichDir}/${Sample}.bw

    其中 ChromSize 用參考基因組和 faToTwoBit, twoBitInfo 兩個工具生成，工具也在上面連結下載。

faToTwoBit GRCh38.primary_assembly.genome.fa GRCh38.primary_assembly.2bit
twoBitInfo GRCh38.primary_assembly.2bit GRCh38.primary_assembly.tab

deeptools bamCoverage
    deeptools bamCoverage 命令將 Bam 文件轉換成 bigwig 文件，免去了自己從 peak calling 結果裡各種手動轉換文件格式過程。推薦使用這個方法。
    這裡 --binSize 默認 50 對於 ATAC-seq 我覺得需要設置小一些，--effectiveGenomeSize 的設置值可以在 https://deeptools.readthedocs.io/en/latest/content/feature/effectiveGenomeSize.html 查看。

bamCoverage --binSize 10 --blackListFileName ${BlackList} --numberOfProcessors 5 --effectiveGenomeSize 2747877777 --normalizeUsing RPGC --outFileFormat bigwig -b ${BamDir}/${Sample}_ShiftSort.bam -o ${BwDir}/${Sample}.bw

其中 RPGC 計算方法：
number of reads per bin / scaling factor for 1x average coverage
Scale factor 計算方法：
(total number of mapped reads * fragment length) / effective genome size

pyGenomeTracks
    先從 bw/bed 文件生成 pyGenomeTracks 可用的 tracks 文件。

make_tracks_file --trackFiles ${bwDir}/KO1_ATAC.bw ${bwDir}/KO2_ATAC.bw ${bwDir}/WT1_ATAC.bw ${bwDir}/WT2_ATAC.bw -o ${TrackDir}/Track1.ini

    生成的 ini 文件裡面包含了畫圖配置，所以需要修改畫圖配置的，手動修改這個文件。參數表在 https://pygenometracks.readthedocs.io/en/latest/content/possible-parameters.html 查看。
下面這個配置（不展示重複部分）

[x-axis]

[spacer]
height = 1

[KO1_ATAC]
file = /Example/Vis/KO1_ATAC.bw

title = KO1
height = 2

color = #F56F08
min_value = 0
max_value = 25
number_of_bins = 2000
nans_to_zeros = true
summary_method = mean
show_data_range = true
file_type = bigwig

[gtf]
file = /Example/GRCh38/GENCODE/gencode.v35.annotation.gtf
height = 3
title = Gencode Annotation
file_type = gtf
prefered_name = gene_name
merge_transcripts = true

生成如下圖片

一個特定區域 Reads 覆蓋

EnrichedHeatmap
    EnrichedHeatmap 用於生成一些區域的信號富集圖，這個圖用 ChIPseeker 和 deeptools 也可以生成，依據自己的喜好選擇。

圖來源於 Yin, Senlin 等論文

參考資料
https://training.galaxyproject.org/training-material/topics/epigenetics/tutorials/atac-seq/tutorial.html
Yan, Feng, et al. "From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis." Genome biology 21.1 (2020): 22.
Brouilette, Scott, et al. "A simple and novel method for RNA‐seq library preparation of single cell cDNA analysis by hyperactive Tn5 transposase." Developmental Dynamics 241.10 (2012): 1584-1590.
Gontarz, Paul, et al. "Comparison of differential accessibility analysis strategies for ATAC-seq data." Scientific reports 10.1 (2020): 1-13.
Yin, Senlin, et al. "Transcriptomic and open chromatin atlas of high-resolution anatomical regions in the rhesus macaque brain." Nature communications 11.1 (2020): 1-13.
https://informatics.fas.harvard.edu/atac-seq-guidelines.html
Buenrostro, Jason D., et al. "ATAC‐seq: a method for assaying chromatin accessibility genome‐wide." Current protocols in molecular biology 109.1 (2015): 21-29.
Yan, Feng, et al. "From reads to insight: a hitchhiker’s guide to ATAC-seq data analysis." Genome biology 21.1 (2020): 22.
Buenrostro, Jason D., et al. "Transposition of native chromatin for multimodal regulatory analysis and personal epigenomics." Nature methods 10.12 (2013): 1213.
Li, Zhijian, et al. "Identification of transcription factor binding sites using ATAC-seq." Genome biology 20.1 (2019): 45.
https://www.biostars.org/p/418363/