背景： 染色質和染色體的結構和功能
每一條染色單體由單個線性DNA分子組成漆诽。細胞核中的DNA是經過高度有序的包裝加叁，否則就是一團亂麻哥力，不利于DNA復制和表達調控蔗怠。這種有序的狀態(tài)才能保證基因組的復制和表達調控能準確和高效進行。
包裝分為多個水平吩跋，核小體核心顆粒(nucleosome core particle)寞射、染色小體(chromatosome)、 30 nm水平染色質纖絲(30 nm fibre)和高于30 nm水平的染色體包裝锌钮。在細胞周期的不同時期桥温，DNA的濃縮程度不同，間期表現為染色質具有轉錄活性轧粟，而中期染色體是轉錄惰性策治。細胞主要處于分裂間期，所以DNA大部分時間都是染色質而不是染色體兰吟，只不過大家喜歡用染色體泛指染色質和染色體通惫。
很久之前大家喜歡研究中期的染色體，原因是光學顯微鏡只能看的到這種高度濃縮狀態(tài)的DNA結構混蔼。不過中期染色體在轉錄上是惰性的履腋，沒有研究間期染色體的意義大。后來技術發(fā)展了，大家就開始通過熒光蛋白標記技術以及顯微鏡技術研究間期染色質的三維結構和動態(tài)遵湖。比如說悔政，間期染色體其實并非隨機地彌漫在細胞核中，不同的染色體占據相對獨立的空間延旧，染色體在細胞核所占的空間稱之為染色體領地(chromosome territory, CT)谋国。研究發(fā)現，貧基因(gene-porr)的染色體領域一般傾向于靠近核膜迁沫，而富含基因(gene-rich)的染色體領地通常位于細胞核內部芦瘾。這也反應了人類社會的情況，富人處于核心區(qū)集畅，窮人在邊緣地帶近弟。
除了染色體細胞核內的三維結構外，還需要談談和轉錄調控相關的染色質的核小體挺智。用內切核糖酶--微球菌核酸酶(micrococcal nuclease, MNase, MN酶)處理染色質可以得到單個核小體祷愉。核小體是染色質的基本結構，由DNA赦颇、蛋白質和RNA組成的一種致密結構二鳄。組蛋白是由2個H3-H4二聚體，2個H2A-H2B二聚體形成的八聚體沐扳，直徑約為10 nm泥从， 組蛋白八聚體和DNA結合在一起形成的核心顆粒包含146bp DNA句占。DNA暴露在核小體表面使得其能被特定的核酸酶接近并切割沪摄。
染色質結構改變會發(fā)生在與轉錄起始相關或與DNA的某種結構特征相關的特定位點。當染色質用DNA酶I(DNase)消化時纱烘，第一個效果就是在雙鏈體中特定的超敏位點(hypersenitive site)引入缺口杨拐，這種敏感性可以反應染色質中DNA的可及性(accessible)，也就是說這些是染色質中DNA由于未組裝成通常核小體結構而特別暴露出的結構擂啥。
許多超敏位點與基因表達有關哄陶。每個活性基因在啟動子區(qū)域都存在一個超敏位點。大部分超敏位點僅存在于相關基因正在被表達的或正在準備表達的細胞染色中哺壶；基因表達不活躍時他們則不出現屋吨。
染色質開放區(qū)域和ATAC-Seq
背景已經談到，超敏位點和基因表達有關山宾，并且超敏位點反應了染色質的可及性至扰。也就可以反推出“可及性”的染色質結構區(qū)域可能與基因表達調控相關。于是2015年的一篇文章Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNDNA-binding proteins and nucleosome position就使用了超敏Tn5轉座酶切割染色質的開放區(qū)域资锰，并且加上接頭(adapter)進行高通量測序敢课。
那篇文章通過ATAC-Seq得到了如下結論：
- ATAC-seq insert sizes disclose nucleosome positions
- ATAC-seq reveals patterns of nucleosome-TF spacing
- ATAC-seq footprints infer factor occupancy genome wide
- ATAC-seq enables epigenomic analysis on clinical timescales
也就是說ATAC-Seq能幫助你從全基因組范圍內推測可能的轉錄因子，還能通過比較不同時間的染色質開放區(qū)域解答發(fā)育問題。
數據分析概要
在前面的鋪墊工作中直秆，一共提到了三種酶,能切割出單個核小體的MNase, 能識別超敏位點的DNase 和ATAC-Seq所需要的Tn5 transposase濒募，這三種酶的異同如下圖：
分析ATAC-Seq從本質上來看和分析ChIP-Seq沒啥區(qū)別，都是peak-calling圾结，也就是從比對得到BAM文件中找出reads覆蓋區(qū)瑰剃，也就是那個峰。(尷尬的是筝野，這句話對于老司機而言是廢話培他，對于新手而言則是他們連ChIP-Seq都不知道)那么問題集中在如何找到peak，peak的定義是啥?
“Peak不就是HOMER/MACS2/ZINBA這些peak-finder工具找到的結果嗎遗座？”
找Peak就好像找美女舀凛，你覺得美女要手如柔荑，膚如凝脂途蒋，領如蝤蠐猛遍，齒如瓠犀，螓首蛾眉号坡，巧笑倩兮懊烤，美目盼兮。但實際情況下宽堆，是先給你看一個長相平平的人或者有點缺陷的人腌紧，然后再把那個人PS一下，你就覺得是一個美女了畜隶。理想情況下壁肋， peak應該是一個對稱的等腰三角形，并且底角要足夠的大籽慢。實際情況下是稍微不那么平坦似乎就行了浸遗。
假設目前已經找到了peak，這是不是意味著我們找到轉錄因子了箱亿？不好意思跛锌，這不存在的，因為ATAC-Seq只是找到了全基因組范圍的開放區(qū)域届惋，而這些開放區(qū)域的產生未必是轉錄因子引起髓帽，所以需要一些預測性工作。
數據分析實戰(zhàn)
基本信息
以目前預發(fā)表在bioRxiv的文章“Chromatin accessibility changes between Arabidopsis stem cells and mesophyll cells illuminate cell type-specific transcription factor networks” 為例脑豹，介紹ATAC-Seq數據分析的套路郑藏。
GEO編號：GSE101940，一共6個樣本晨缴，SRR為SRR5874657~SRR587462
實驗設計：用INTACT方法提取植物干細胞(stem cell)和葉肉細胞(mesophyll cells)的細胞核译秦，然后通過ATAC-Seq比較兩者在轉錄因子上的差別。 INTACT技術提供數據的利用率，因為ATAC-Seq主要是分析染色體的開放區(qū)域筑悴，所以比對到細胞器的序列在后期分析中會被丟棄们拙。
分析流程：分為數據預處理,Peak-calling和后續(xù)分析三步。
數據預處理
數據預處理步驟分為：質量控制阁吝，原始序列比對砚婆，比對后去除重復序列和細胞器序列。當然在這之前突勇，先得做一下準備工作装盯，創(chuàng)建工作環(huán)境，從SRA下載數據并進行數據解壓甲馋。
# 創(chuàng)建項目文件下
mkdir -p ATAC-Seq/{data/raw_data,analysis,script,ref}
# 使用sra-tool prefetch下載數據, 數據保存在~/ncbi/public/sra
for i in `seq 57 62`;
do
prefetch SRR5874${i} &
done
# 數據下載完埂奈，用fastq-dump解壓
for i in `seq 57 62`;
do
fastq-dump --split-3 --defline-qual '+' --defline-seq '@$ac-$si/$ri length=$rl' --gzip SRR5874${i} -O data/raw_data &
done
質量控制：在數據分析之前先要大致了解手頭數據的質量，目前基本就用fastqc了
mkdir -p analysis/fastqc
for i in `seq 57 62`; do fastqc data/raw_data/SRR58746${i}_{1,2}.fastq.gz -o analysis/qc & done
# multiqc匯總
multiqc analysis/qc/ -o analysis/qc/
FastQC結果大部分都過關定躏，除了在read的各位置堿基含量圖上fail账磺。具體原因我還不知道，文章中并沒有提到要對原始數據進行預處理痊远。
序列比對: 目前在Peak-calling這里分析的流程中垮抗，最常用的比對軟件就是Bowtie, 分為Bowtie1和Bowtie2，前者適合25~50bp碧聪，后者適用于 > 50bp的情況冒版。此處分析的read長度為50 bp，因此選擇bowtie1逞姿。
使用之前建議先看看兩個工具的手冊辞嗡，了解一下參數說明。我也是通過Bowtie2的手冊才發(fā)現Bowtie2和Bowite其實是兩個用途不同的工具哼凯，而不是說bowtie2是用來替代bowtie1.
# 下載或者鏈接已有的參考基因組到ref文件下
ln -s ~/db/Genomes/Athalina/TAIR10/Sequence/TAIR10.fa ref/
# 建立BOWTIE2 index欲间，或者下載已有的index, 或軟連接已有的索引
bowtie-build --threads 8 ref/TAIR10.fa ref/TAIR10
# 注意不要用bowtie2-build
# 序列比對
for i in `seq 57 62`;
do
bowtie -p 10 -S -m 1 -X 2000 --sam-RG "ID:sample_${i}" \
--sam-RG "PL:illumina" --sam-RG "SM:SRR58746${i}" \
ref/TAIR10 \
-1 <(zcat data/raw_data/SRR58746${i}_1.fastq.gz) \
-2 <(zcat data/raw_data/SRR58746${i}_2.fastq.gz) | \
samtools sort -@ 6 -m 1G -o analysis/BAM/SRR58746${i}_sorted.bam ;
done &
# 建索引
ls analysis/BAM/*sorted.bam | while read id; do sambamba index -t 6 $id ;done
這里使用Bowtie比對到擬南芥參考基因組TAIR10楚里，參數為-X2000, 允許長達2 Kb的片段断部， -m1僅保留唯一聯配。
初步比對后班缎，可以統(tǒng)計下比對到organellar genomesh和nuclear genome的read數量蝴光。這個工具可以在shell腳本中用samtools處理，也可以用python的pysam模塊.Pysam封裝了htslib C-API达址，提供了SAM/BAM/VCF/BCF/BED/GFF/GTF/FASTA/FASTQ的操作蔑祟，最新版本已經支持Python3,強烈推薦學習。
在script目錄下新建alignment_stat.py沉唠， 添加如下內容
import pysam
import glob
import os
import numpy as np
from pandas import Series, DataFrame
# 利用glob讀取文件路徑
bam_files = glob.glob('../analysis/BAM/*sorted.bam')
threads = "10"
# 統(tǒng)計每條染色體的reads數
Chrs = ['Chr1','Chr2','Chr3','Chr4','Chr5','ChrM','ChrCh']
reads_count = DataFrame(np.ones((len(Chrs),len(bam_files))),index=Chrs,columns=bam_files)
# 需要索引
for bam in bam_files:
reads_count[bam] = [int(pysam.view("-@", threads, "-c", bam, Chr).strip()) for Chr in Chrs]
sam.close()
reads_count = reads_count.T
reads_count['nuc_ratio'] = np.sum(reads_count.iloc[:,:5],axis=1) / np.sum(reads_count.iloc[:,:],axis=1)
reads_count['org_ratio'] = np.sum(reads_count.iloc[:,5:],axis=1) / np.sum(reads_count.iloc[:,:],axis=1)
print(reads_count)
運行python alignment.py,得到如下結果：
比對后去(標記)重復和細胞器reads：由于細胞器DNA蛋白結合少疆虚，所以顯然更容易被Tn5 轉座酶切割，普通的ATAC-Seq的read就會有大量是細胞器的DNA，這就是為啥需要用INTACT技術径簿。此外如果不是PCR-free的建庫方法罢屈，會有大量重復的read，也就需要標記或去除重復篇亭。
bigwig定量文件: 使用deepTools進行標準化和可視化, 一般以RPKM做標準化缠捌，默認bin為1
# 項目文件下
# 標記重復
mkdir -p analysis/dupbam
filenames=$(ls analysis/BAM/*sorted.bam)
for fn in ${filenames}
do
fnn=$(basename $fn)
output=${fnn%%.bam}
gatk-launch MarkDuplicates -I ${fn} \
-O analysis/dupbam/${output}_nuc_markdup.bam \
-M analysis/dupbam/${output}_nuc_markdup_matrix.txt
done
# bigwig, 必須要有索引
mkdir -p analysis/bigwig
markdup_files=$(ls analysis/dupbam/*_markdup.bam)
for markdup in ${markdup_files}
do
fn=$(basename $markdup)
output=$(fn%%.bam)
sambamba index -t 6 $markdup
bamCoverage -b $deup --ignoreDuplicates \
--skipNonCoveredRegions \
--normalizeUsingRPKM \
--binSize 1 -p 5 -o analysis/bigwig/${output}.bw
done
得到的BW文件可以用IGV進行可視化
