寫在前面
- 其實很多作者已經探索的差不多了闪檬，相關腳本、教程也很多氛悬，但我還是自己系統(tǒng)的總結一下吧则剃。有時候感覺是浪費時間，但是在寫的時候能發(fā)現并學習到很多新的知識如捅，不斷記錄反思棍现，也算是對自己這段時間的學習成果有個交代。但由于我是初學者镜遣，還望各位不吝賜教己肮。
ATAC-seq的首次提出：【文章原文】【scATAC原文】
ATAC-seq植物界的里程碑：【文章原文】【中文解讀】
ATAC-seq針對不同植物物種進行研究：【文章原文】【中文解讀】
ATAC-seq針對不同植物組織進行研究：【文章原文】【中文解讀】
ATAC-seq針對不同處理條件進行研究：【文章原文】【中文解讀】
ATAC-seq生信分析方法綜述：【文章原文】【中文解讀1】【中文解讀2】
ATAC-seq基礎原理一文了解：【一文了解ATAC-seq】
一、基礎知識
1. 解密表觀遺傳學的三個方向與測序方法
-
1. 探索染色質的開放性（chromatin accessibility）
ATAC-seq: Assay of Transposase Accessible Chromatin sequencing
DNase-seq: DNase I hypersensitive sites sequencing
FAIRE-seq: Formaldehyde-Assisted Isolation of Regulatory Elements sequencing -
2. 探索轉錄因子的綁定與組蛋白修飾（TF binding and histone modifications）
ChIP-seq: Chromatin Immuno-Precipitation sequencing -
3. 探索核小體占位（nucleosome positioning and occupancy）
MNase-seq: Micrococcal Nuclease sequencing
2. 名詞解釋
-
epigenetics：表觀遺傳學悲关。表觀遺傳學修飾在不改變DNA序列的情況下控制著基因的表達谎僻，包括染色質重塑、組蛋白修飾寓辱、DNA甲基化和microRNA通路等 -
peaks: 峰艘绍。常用來表示染色質的開放程度，因為是測序的reads落在了染色質的開放區(qū)秫筏，堆疊后被可視化的一種豐度的體現诱鞠。 -
THSs: 轉座酶超敏感位點（transposase hypersensitive sites）挎挖。 -
CREs: 順式調控元件（cis-regulatory elements）。即DNA分子中具有轉錄調節(jié)功能的特異DNA序列航夺。按功能特性蕉朵，真核基因順式作用元件分為啟動子、增強子及沉默子阳掐。 -
ACRs: 染色質開放區(qū)域（accessible chromatin regions）始衅。即正常或核小體被酶切裸露出來的DNA片段所在的區(qū)域锚烦。這篇【NP | 2019】根據ACRs距離最近基因的距離將ACRs分為三種類型：genic (gACRs; overlapping a gene), proximal (pACRs; within 2?kb of a gene) or distal (dACRs; >2?kb from a gene)，分別是跨越基因的帝雇，近端的涮俄，遠端的染色質開放區(qū)。 -
transposon: 轉座子尸闸。一段可以從原位上單獨復制或斷裂下來彻亲，環(huán)化后插入另一位點，并對其后的基因起調控作用的DNA序列吮廉。戳看轉座子科普 -
promoter: 啟動子苞尝。啟動子是位于結構基因5'端上游的DNA序列，能活化RNA聚合酶宦芦，使之與模板DNA準確的結合并具有轉錄起始的特異性宙址。每個啟動子包括至少一個轉錄起始點以及一個以上的功能組件（典型的如TATA盒子） -
proximal promoters: 近端啟動子。是DNA上位于基因開始之前的一個區(qū)域调卑，在那里蛋白質和其他分子結合在一起準備讀取該基因抡砂。 -
enhancer: 增強子。增強子是遠離轉錄起始點恬涧、決定基因的時間注益、空間特異性表達、增強啟動子轉錄活性的DNA序列溯捆，其發(fā)揮作用的方式通常與方向丑搔、距離無關，可位于轉錄起始點的上游或下游提揍。從功能上講啤月，沒有增強子存在，啟動子通常不能表現活性劳跃；沒有啟動子時顽冶，增強子也無法發(fā)揮作用。 -
TFs: 轉錄因子（transcription factors）是保證目的基因以特定的強度在特定的時間與空間表達的蛋白質分子售碳。與RNA聚合酶Ⅱ形成轉錄起始復合體强重，共同參與轉錄起始的過程绞呈。 -
TSS: 轉錄起始位點（transcription start site）。在一個典型的基因內部间景，排列順序為轉錄起始位點(TSS佃声，一個堿基)-起始密碼子編碼序列 (ATG)-終止密碼子編碼序列(TGA)-轉錄終止位點 (TTS)，即TSS-ATG-TGA-TTS -
histone：組蛋白倘要。通常含有H1圾亏，H2A，H2B封拧，H3志鹃，H4等5種成分，其中H1與H3極度富含賴氨酸（lysine）泽西，H1不保守曹铃，其他組蛋白的基因非常保守。除H1外捧杉，其他4種組蛋白均分別以二聚體(共八聚體）相結合陕见，形成核小體核心。DNA便纏繞在核小體的核心上味抖。而H1則與核小體間的DNA結合 -
nucleosome: 核小體评甜。是由DNA和組蛋白形成的染色質基本結構單位。每個核小體由146bp的DNA纏繞組蛋白八聚體1.75圈形成仔涩。核小體核心顆粒之間通過50bp左右的連接DNA相連忍坷，暴露在核小體表面的DNA能被特定的核酸酶接近并切割 -
H3K4me1: 組蛋白H3上的第4位賴氨酸發(fā)生單甲基化（mono-methylation of H3 at lysine 4） （NG | H3K4me1與增強子的關系） -
H3K9ac：組蛋白H3上的第9位賴氨酸發(fā)生乙酰化（H3 acetylation marks at lysine 9）關于H3K4me1與H3K9ac的不同熔脂，依舊可以參考這篇文獻【NP | 2019】或【NC | 2019】 -
組蛋白甲基化：甲基化可發(fā)生在組蛋白的賴氨酸和精氨酸殘基上承匣，而且賴氨酸殘基能夠發(fā)生單、雙锤悄、三甲基化韧骗，而精氨酸殘基能夠單、雙甲基化零聚，這些不同程度的甲基化極大地增加了組蛋白修飾和調節(jié)基因表達的復雜性 -
組蛋白乙跖郾化：四種類型的組蛋白相互作用，將細胞核里的DNA緊緊地包裝起來隶症。這樣的緊密包裝能夠有效阻止酶讀取DNA上的遺傳信息政模。然而，乙趼旎幔基連到組蛋白上能削弱它們對DNA的占據淋样。因此局部乙酰化能暴露出相應的基因胁住，讓它們更容易激活
真核生物細胞核內染色質結構
二趁猴、基礎分析
1. 分析前的準備
1.1 下載ATAC測序數據（或者用你自己的數據也行）
所需軟件：Aspera Connect
軟件安裝與使用：參考【秘籍 | 驚了刊咳，使用aspera下載SRA數據速度高達 291Mb/s】
數據下載：項目ID【PRJNA480717】【文章原文】，SRR序列號如下圖儡司，下載方法和上文相同娱挨。由于數據量大，測試浪費時間捕犬，所以我在每個文件中取了100w條reads跷坝，由于百度網盤太慢了，你可以去扣扣群559758504里面隨便下載一個用于測試碉碉。
# 提取100w條reads代碼
ls *gz|while read id;do zcat $id|head -4000000|gzip > ./100/${id};done
1.2 下載tair10基因組與gff3注釋文件
wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas
wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes_transposons.gff
# 改名
mv TAIR10_chr_all.fas TAIR10_genome.fa
mv TAIR10_GFF3_genes_transposons.gff TAIR10.gff3
1.3 安裝miniconda
- 由于后續(xù)分析大多數軟件都是來源于conda柴钻，因此你必須知道怎么使用。參考【conda的安裝與使用】
1.4 所需軟件安裝（分析中所用到的軟件都匯集于此垢粮，下文不再介紹如何下載）
conda install -c bioconda -y fastqc
conda install -c bioconda -y trimmomatic
conda install -c bioconda -y bwa
conda install -c bioconda -y samtools
conda install -c bioconda -y picard
conda install -c bioconda -y macs2
conda install -c bioconda -y bedtools
conda install -c bioconda -y deeptools
conda install -c bioconda -y homer
2. 測序質量檢測
fastqc -t 8 -o ./ *.fastq.gz
-t: 線程
-o: 存放路徑斩跌，不用指定前綴，默認為.fastq.gz前面的字段
*.fastq.gz: fastqc可以同時接受多個fastq.gz文件捞慌，因此采用正則表達式【*】表示全部
- 結果解讀
由于開始測序時不準確耀鸦，因此前15bp左右的ATGC含量有波動，需要在后面截掉啸澡。而最長reads為76bp袖订，則暗示著我最后至少要保留多長的序列。
一般我的計算思路是：保留的序列長度 =（最長的的reads長度-前面波動的堿基長度）× 40%
2. 測序質量檢控制
- 所需軟件：Trimmomatic【官方網站】【中文介紹】
其實這批數據已經做過質控了嗅虏，如何判斷的呢洛姑？QC的為一個范圍且短于提供給NCBI上的數據。因此我們沒必要質控皮服，但為了流程的完整性楞艾，還是假裝質控一下。質控之后龄广，還要進行質檢硫眯，如此循環(huán)，直到合格择同。
# 定義一些變量两入，后面修改起來方便
filename=SRR7512044
trimmomatic PE -threads 10 \
./${filename}_1.fastq.gz ./${filename}_2.fastq.gz \
./${filename}_1.clean.fq.gz ./${filename}_1.drop.fq.gz \
./${filename}_2.clean.fq.gz ./${filename}_2.drop.fq.gz \
HEADCROP:15 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:24
PE: 雙端模式。需要兩個輸入文件敲才，正向測序序列和反向測序序列：sample_R1.fastq.gz和sample_R2.fastq.gz裹纳；
以及四個輸出文件：sample_1.clean.fq.gz和sample_1.drop.fq.gz择葡；sample_2.clean.fq.gz和sample_2.drop.fq.gz
-threads: 線程數
HEADCROP: 從 reads 的開頭切掉指定數量的堿基，即從fastqc中看的
LEADING: 從 reads 的開頭切除質量值低于閾值的堿基
TRAILING: 從 reads 的末尾開始切除質量值低于閾值的堿基
SLIDINGWINDOW: 從 reads 的 5' 端開始痊夭，進行滑窗質量過濾刁岸，切掉堿基質量平均值低于閾值的滑窗
MINLEN: 如果經過剪切后 reads 的長度低于閾值則丟棄這條 reads，即上面的計算方法：（76-15）×40%≈24
3. 建立索引與序列比對
- 所需軟件：bwa【官方網站】【中文介紹】與samtools【官方網站】【中文介紹】
- 建立基因組索引她我。需要的時間挺長的虹曙，建議放后臺。建好之后以后就可以跳過這一步了番舆，且在經典的比對軟件STAR酝碳、hisat2、bowtie2和bwa-mem中我更喜歡使用STAR和bwa-mem恨狈。本次使用bwa-mem操作演示
bwa index -p /路徑/bwaIndexTair /路徑/TAIR10_genome.fa
-p: 索引前綴疏哗，后綴自動補充
TAIR10_genome.fa: 基因組文件
- 序列比對
由于sam文件較大，因此我們直接跳過sam禾怠，使用samtools轉換為排序后的bam文件返奉，注意【|】不要漏掉，且排序是依據reads比對到基因組的位置排的
# 定義一些變量吗氏，后面用
filename=SRR7512044
index=/路徑/bwaIndexTair
bwa mem -M -t 8 \
-R "@RG\tID:${filename}\tSM:${filename}\tLB:WXS\tPL:Illumina" \
$index \
${filename}_1.clean.fq.gz \
${filename}_2.clean.fq.gz \
| samtools sort -O bam -@ 10 -o ./${filename}.raw.bam
#無論干啥芽偏，生成bam文件后就要建立bam索引，并檢查比對情況
samtools index -@ 10 -b ./${filename}.raw.bam
samtools flagstat ./${filename}.raw.bam > ./${filename}.raw.stat
cat ${filename}.raw.stat
-M: 將 shorter split hits 標記為次優(yōu)弦讽，以兼容 Picard’s markDuplicates 軟件
-R: STR 完整的read group的頭部污尉，可以用 '\t' 作為分隔符， 在輸出的SAM文件中被解釋為制表符TAB. read group 的ID往产，會被添加到輸出文件的每一個read的頭部被碗，強推去看關于-R的解釋【-R的使用】
-O: 表示輸出的bam文件
-o: 輸出的bam文件名
-t與@: 線程數
- 關于stat中內容的解釋如下
14608455 + 0 in total (QC-passed reads + QC-failed reads) reads總數
37967 + 0 secondary 出現比對到參考基因組多個位置的reads數
0 + 0 supplementary 可能存在嵌合的reads數
0 + 0 duplicates 重復的reads數
14590894 + 0 mapped (99.88% : N/A) 比對到參考基因組上的reads數
14570488 + 0 paired in sequencing 屬于PE read的reads總數。
7285244 + 0 read1 PE read中Read 1 的reads 總數仿村。
7285244 + 0 read2 PE read中Read 2 的reads 總數锐朴。
14507068 + 0 properly paired (99.56% : N/A) 完美比對的reads總數。PE兩端reads比對到同一條序列蔼囊，且根據比對結果推斷的插入片段大小符合設置的閾值焚志。
14551500 + 0 with itself and mate mapped PE兩端reads都比對上參考序列的reads總數。
1427 + 0 singletons (0.01% : N/A) PE兩端reads压真，其中一端比上娩嚼，另一端沒比上的reads總數。
26260 + 0 with mate mapped to a different chr PE read中滴肿，兩端分別比對到兩條不同的序列的reads總數岳悟。
17346 + 0 with mate mapped to a different chr (mapQ>=5) PE read中，兩端分別比對到兩條不同
4. 去除PCR重復并再次進行質控
- 所需軟件：picard【官方網站】【中文介紹】與samtools【官方網站】【中文介紹】
不建議使用sambamba，雖然快贵少，但是毛病多呵俏，尤其是在你不知情的情況下，軟件卡死或者在最后輸出文件時提示輸出文件夾過多等滔灶。
picard MarkDuplicates \
REMOVE_DUPLICATES=true \
I=${filename}.raw.bam \
O=${filename}.rmdup.bam \
M=${filename}.rmdup.log
# 實時監(jiān)測我們的數據發(fā)生了什么變化
samtools index ${filename}.rmdup.bam
samtools flagstat ${filename}.rmdup.bam > ./${filename}.rmdup.stat
REMOVE_DUPLICATES=true: 表示將檢測出的PCR duplication直接去除普碎；
I: 或者INPUT指定輸入的bam文件；
O: 或者OUTPUT指定輸出去除PCR duplication的bam文件录平；
M: 或者METRICS_FILE指定輸出的matrix log文件麻车。
- 下面是再次進行質控
samtools view -h -f 2 -q 30 ${filename}.rmdup.bam \
| grep -v -e "mitochondria" -e "*" -e "chloroplast" \
| samtools sort -O bam -@ 10 -o - > ${filename}.last.bam
# 實時監(jiān)測我們的數據發(fā)生了什么變化
samtools index ${filename}.last.bam
samtools flagstat ${filename}.last.bam > ./${filename}.last.stat
-h: 輸出的sam文件帶header，默認不帶
-f: 輸出含有所有flag的reads
-q 比對的最低質量值斗这，一般20/30都可以
其他參數介紹看 samtools 使用簡述动猬，需要說明的是由于比對到線粒體和葉綠體的reads不是我們關注的重點，因此使用grep剔除掉表箭，當然了赁咙，做動物的話需要剔除什么視情況而定了。
- 好了免钻，下面來看看我們的數據從
到底發(fā)生了那些變化
- 然后彼水，從IGV看看發(fā)生了什么變化
關于使用IGV查看bam文件中reads呈現不同顏色的說明，請參考：【玩轉基因組瀏覽器之IGV展示bam文件】
5. 統(tǒng)計過濾后的結果中Insertsize長度分布
- 所需軟件：picard【官方網站】【中文介紹】
Insertsize =（read2的比對位點 - read1的比對位點）+ read2的長度极舔。也可以隨便點擊上圖IGV中任一一條reads查看插入片段的大小
picard CollectInsertSizeMetrics \
I=./${filename}.last.bam \
O=./${filename}.last.insertsize \
H=./${filename}.last.hist.pdf
有關下圖的介紹凤覆，請看【給你bam文件，你會畫插入片段長度分布圖嗎姆怪？】
5. 
與FRiP_score/SPOT_score的計算
- 所需軟件：macs2【官方網站】【中文介紹】與床上用品bedtools 【官方網站】【中文介紹】
因為macs2需要用到基因組的有效大小叛赚，而軟件本身就提供了4個物種的有效基因組大小澡绩，遠遠不夠用盎摇！因此如何計算其他物種基因組的有效大蟹士ā溪掀？請看我這篇文章【秘籍 | 計算基因組的有效大小（effective genome size）】
effective_genome_size=119481543
bedtools bamtobed -i ${filename}.last.bam > ${filename}.last.bed
macs2 callpeak \
-t ${filename}.last.bed \
-g ${effective_genome_size} \
--nomodel --shift -100 --extsize 200 \
-n ${filename}.last \
-q 0.01 --outdir ./peaks
- 其實步鉴，這里面討論最多的是
--nomodel --shift -100 --extsize 200這些參數如何選擇揪胃，下面的圖很形象的展示了參數的作用。當然氛琢，我也是查閱了很多資料與文獻喊递，一般默認在ATAC-seq，DNase-seq阳似，FAIRE-seq的時候將shift設置為extsize的一半骚勘，且參數固定為：--nomodel --shift -100 --extsize 200。而在MNase-seq的時候苍姜，參數固定為：--nomodel --shift 37 --extsize 73续扔。在ChiP-seq的時候不用移峰，所以只使用-nomodel踱阿，當做組蛋白修飾的時候泽疆，由于peak并不典型户矢，所以使用–nomodel –shift 73參數。另外如何使用macs2 predictd預測雙峰模型殉疼，這里就不展開介紹了梯浪，可以私聊我，我把腳本發(fā)你瓢娜。
- 下面是我見過的圖中解釋移峰現象最清楚的了，可以參考
- 最后恋腕，就是輸出文件的解讀
- ${filename}_peaks.xls
前幾行是callpeak時的命令抹锄。后面是：
1.染色體號
2.peak起始位點
3.peak結束位點
4.peak區(qū)域長度
5.peak的峰值位點（summit position）
6.peak 峰值的高度（pileup height at peak summit, -log10(pvalue) for the peak summit）
7.peak的富集倍數（相對于random Poisson distribution with local lambda）
- ${filename}_peaks.narrowPeak
narrowPeak文件是BED6+4格式，可以上傳到IGV查看
1.染色體號
2.peak起始位點
3.peak結束位點
4.peak name（樣本名+_peak_1的后綴之類的）
5.-10*log10qvalue
6.正負鏈
7.fold change
8.-log10pvalue
9.-log10qvalue
10.相對于峰的峰頂位置（relative summit position to peak start）
- 注意：
- ${filename}_summits.bed
BED格式的文件伙单，包含peak的summits（峰頂）位置
1.染色體號
2.峰的峰頂起始位點（注意是峰頂，不是峰哈肖，他只是一個點）
3.峰的峰頂起始位點（因此二者就差一個堿基的位置）
4.峰頂（樣本名+_peak_1的后綴之類的）
5.-log10pvalue
注意：
另外，需要告訴你的是淤井，bed文件中第2列中的峰頂起始位點是如何計算的布疼？
bed文件中第2列（峰的峰頂起始位點） = xls第二列 + xls第十列 -1”液荩或者↓↓↓↓
bed文件中第2列（峰的峰頂起始位點） = narrowPeak第二列 + narrowPeak第十列-
總覽比較：
bed游两，narrowPeak，bdg漩绵，xls四種類型輸出文件的比較 FRiP_score的計算【官方介紹】
簡單來說贱案，FRiP_score說的就是統(tǒng)計落在Peaks里面的reads數量與所有比對到基因組上reads數量的比值≈雇拢可以作為判定樣本測序質量的一個指標宝踪。
Reads=$(bedtools intersect -a ${filename}.last.bed -b ${filename}_peaks.narrowPeak |wc -l|awk '{print $1}')
totalReads=$(wc -l ${filename}.last.bed|awk '{print $1}')
echo ${Reads} ${totalReads} ${filename} 'FRiP_score:' $(echo "scale=2;100*${Reads}/${totalReads}"|bc)'%'
- 官網上說ATAC-Seq的FRiP score should be >30%（https://www.encodeproject.org/atac-seq/）我猜是如果大多數reads都不在peaks區(qū)域，說明測序質量或者實驗效果不是太好碍扔，因為酶切切開的就是peaks區(qū)域瘩燥，太少不是就偏離實驗目的了嘛
- SPOT_score的計算【官方介紹】【軟件安裝】
簡單來說，SPOT_score就是在FRiP_score的基礎上增加了測序深度的條件限制不同，比如FRiP是統(tǒng)計全局的reads厉膀，而SPOT只統(tǒng)計測序深度為30×的reads，大于30×的reads不再參與統(tǒng)計。目前我還在探索中如何計算站蝠。
二汰具、高級分析
所謂的高級并不高大上，做了有可能也不能讓你的文章提升檔次菱魔，望周知留荔！
1. IGV可視化
前面的bam文件，narrowPeak文件澜倦，bed文件都有了聚蝶，現在還差一個bw文件，因此我們先跑下面的代碼藻治，然后在IGV里面統(tǒng)一看一下
bamCoverage --bam ${filename}.last.bam -o ${filename}.last.bw \
--binSize 10 \
--normalizeUsing RPGC \
--effectiveGenomeSize ${effective_genome_size} \
--ignoreForNormalization chrX \
--extendReads
注：如果發(fā)現沒有這個包pysam.libchtslib碘勉，直接安裝pip install numpydoc
-o: 生成bw文件
binSize: 自行設定覆蓋度計算的窗口大小(bin)
normalizeUsing: 數據準化方法，有三種RPGC桩卵，RPKM验靡，CPM，BPM
ignoreForNormalization: 設置你想忽略的染色體名稱
extendReads: 讓reads擴展至fragment的大小
- IGV一覽
2. TSS/TES 可視化
- 所需軟件：deeptools【官方網站】【中文介紹】
主要是如何生成Refseq文件辞州，雖然官網上提供了一些，但總歸是不能覆蓋所有物種寥粹，因此变过，我們可以從gff/gtf文件中自己創(chuàng)造。
附1：gff格式說明：http://www.reibang.com/p/b26c285bd027
附2：Refseq格式說明：http://www.bio-info-trainee.com/2494.html - 先提取基因的位置信息（不是cds或者exon）涝涤，根據不同的gff進行修改媚狰，但格式始終是
chrom/chromStart/chromEnd/name/score/strand
awk -vFS='[\t=;]' -vOFS="\t" \
'{if ($3=="gene") print $1,$4,$5,$10,$6,$7}' \
TAIR10.gff|sed 's/Chr//g' > refseq.bed
#由于tair官網上提供的基因組文件中染色人命名沒有Chr，為保持一致妄痪，這里去掉
- 然后進行TSS可視化
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 2000 -a 2000 \
-R ./refseq.bed \
-S *.bw \
-o TSS.gz \
--skipZeros \
--outFileSortedRegions Heatmap1sortedRegions.bed
plotHeatmap -m TSS.gz -out TSS.Heatmap.pdf --plotFileFormat pdf
plotProfile -m TSS.gz -out TSS.Profile.pdf --plotFileFormat pdf --perGroup
- 然后進行TES 可視化
computeMatrix scale-regions -p 15 \
-b 2000 -a 2000 \
-R ./refseq.bed \
-S *.bw \
--skipZeros -o body.gz
plotHeatmap -m body.gz -out body.Heatmap.pdf --plotFileFormat pdf
plotProfile -m body.gz -out body.Profile.pdf --plotFileFormat pdf --perGroup
- 另外哈雏，參考下面這個圖楞件，看每個參數控制的哪里衫生，如何同時繪制多個等。這是我好久之前總結的土浸，清晰度堪憂
3. 相關性可視化（依舊是使用deeptools）
- 因為我測試樣本少罪针，所以就使用官網的圖了哈，但代碼是能跑通的黄伊，放心使用泪酱。
- 先使用 [bw] 文件生成 [npz] 文件。所有樣本做相關性，用 [*.last.bw]墓阀，多個單列出來也行
multiBigwigSummary bins -b *.last.bw -o number.of.bins.npz
- 畫相關性熱圖
plotCorrelation -in number.of.bins.npz \
--corMethod pearson \
--skipZeros \
--plotTitle "Pearson Correlation of Average Scores Per Transcript" \
--plotFileFormat pdf \
--whatToPlot heatmap --colorMap RdYlBu --plotNumbers \
-o heatmap_PearsonCorr_bigwigScores.pdf \
--outFileCorMatrix PearsonCorr_bigwigScores.tab
- 畫相關性點圖
plotCorrelation -in number.of.bins.npz \
--corMethod pearson \
--skipZeros \
--plotTitle "Pearson Correlation of Average Scores Per Transcript" \
--plotFileFormat pdf \
--whatToPlot scatterplot \
-o scatterplot_PearsonCorr_bigwigScores.pdf \
--outFileCorMatrix PearsonCorr_bigwigScores.tab
5. ChIPseeker對peaks進行注釋和可視化
- 對peak的注釋分為兩個部分——
會將peak所落在基因組上的區(qū)域結構注釋出來毡惜，比如說啟動子區(qū)域，UTR區(qū)域斯撮，內含子區(qū)域等等经伙。同時，也會將與peak最臨近的基因注釋出來勿锅，非常好用帕膜。
其實并不是對peak進行注釋的，而是對peak最臨近的基因注釋的溢十，因此這部分就和傳統(tǒng)的GO/KEGG等分析如出一轍啦垮刹。
5.1 結構注釋
5.1.1 準備工作：導入包，沒有的安裝吧张弛。另外需要使用gff/gtf構建一個TxDb荒典。因為官網上TxDb也不是每一個物種都有（目前共43個），因此我們自己手動構建吞鸭。
#if (!requireNamespace("BiocManager", quietly = TRUE))
# install.packages("BiocManager")
#BiocManager::install("clusterProfiler")
library(clusterProfiler)
library(ChIPseeker)
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("public/bioinfoYu/genome/tair10/TAIR10.gff",
format="gff") #也可以使用gtf
keytypes(txdb) #感興趣的話种蝶，可以用以下方法探索txdb都包含了什么內容
keys(txdb)
5.1.2 單個文件的結構注釋（不推薦在這一步加OrgDB的內容，沒用）
#讀入單個summits文件
peaks <- readPeakFile("SRR7512044_summits.bed")
#結構注釋
peakAnno <- annotatePeak(peaks,
TxDb=txdb,
tssRegion=c(-1000, 1000))
上面的tssRegion取值有說法的瞒大，因為判定"Promoter", "Exon", "Intron", "Intergenic"等結構的依據就是根據范圍來的螃征。比如定義1000以內的，則在1kbp之內規(guī)定為"Promoter (<=1kb)"透敌，而1k-2k內定義為"Downstream (1-2kb)"盯滚。但你定義2000以內呢，則在1kbp之內規(guī)定為"Promoter (<=1kb)"酗电，而1k-2k內定義為"Promoter (1-2kb)"魄藕。雖然，ChIPseeker會標記出來區(qū)域撵术，但還是三思后再取范圍背率。
#注釋完，進行可視化嫩与，多種圖可供選擇
plotAnnoBar(peakAnno)
plotDistToTSS(peakAnno)
vennpie(peakAnno)
plotAnnoPie(peakAnno)
#install.packages("ggupset")
library(ggupset)
upsetplot(peakAnno)
#install.packages("ggimage")
library(ggimage)
upsetplot(peakAnno, vennpie=TRUE)
#最后將我們的注釋結果轉為數據框寝姿，便于查看
df <- as.data.frame(peakAnno)
#將注釋到的基因提取出來（第14列），用于后續(xù)功能分析
gene <- df[,14]
5.1.3 多個文件的結構注釋
#一次也可以讀入多個summits文件划滋，使用list存儲饵筑，然后使用lapply注釋
files = list(peak1 = "peaks.bed", peak2 = ("peak.bed"))
peakAnnoList <- lapply(files,
annotatePeak,
TxDb=txdb,
tssRegion=c(-2000, 2000))
plotAnnoBar(peakAnnoList)
plotDistToTSS(peakAnnoList)
5.2 功能注釋
5.2.1 可選支線
由于這里需要物種的OrgDb注釋庫，但并不是所有的物種都有按ζ骸根资！怎么辦架专？戳看【構建自己的OrgDb庫】。如果你足夠幸運玄帕，研究的物種正好在官網僅提供的20個物種中存在部脚，那就拿來用吧。除了以上這些都好說裤纹，還有是睛低，后面再說。
5.2.2 OrgDb下載
本次測試是用擬南芥的樣本服傍，在官網中可以直接將所需要的的OrgDb庫下載到服務器上钱雷，這樣就不用在R里面重復在網上獲取浪費時間了。比如說擬南芥吹零，就可以下載罩抗。當然也可以直接install了，看你喜歡哪種了灿椅。
wget http://www.bioconductor.org/packages/release/data/annotation/src/contrib/org.At.tair.db_3.13.0.tar.gz
tar -zvxf org.At.tair.db_3.13.0.tar.gz
5.2.3 導入OrgDb
#從本地下載OrgDb庫（如果你直接下載到了R的lib里面套蒂，就直接library也可）
install.packages("/public/bioinfoYu/genome/tair10/org.At.tair.db", repos=NULL)
#另外一種方法
#if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
#BiocManager::install("org.At.tair.db")
#載入庫
library(org.At.tair.db)
#感興趣的話，看看org.At.tair.db包含的都是些啥
#查看包含的列名
keytypes(org.At.tair.db)
#查看包含的基因ID
keys(org.At.tair.db)
#隨便選擇一個基因和一些列名茫蛹，看一看到第是什么東東操刀，不過因為導入了其他包，下面命令報錯婴洼，我至今還沒有解決
select(org.At.tair.db,keys="AT1G01010",columns = c("ENTREZID","GO","SYMBOL"))
5.2.4 ID轉換（非必須）
上面5.1.2已經得到peak附近的基因集【gene】了骨坑，如果需要ID轉換，就是用Y叔的clusterProfiler中的bitr函數
#首先柬采，去除版本號欢唾，如AT1G50040.1后面的.1，當然我們的測試集沒有粉捻，因此不用做
#gene_rm <- gsub(".[0-9]+$","",gene)
#下面進行ID轉換
keytypes(org.At.tair.db)
gene_last <- bitr(geneID = gene,
fromType = "TAIR",
toType = c("ENTREZID", "SYMBOL", "GENENAME"),
OrgDb = org.At.tair.db)
# geneID 需要轉換的ID
# fromType 當前ID類型
# toType 轉換成什么ID礁遣，使用keytypes()查看有哪些類型
# OrgDb 注釋數據庫
5.2.5 GO功能注釋
# 因為keyType不接受數據框類型的數據，因此我們轉換為字符型
genelist <- gene.symbol[,2] #注意我這里把ENTREZID單獨拿出來了肩刃，所以下面自然是keyType = "ENTREZID",
go_result <- enrichGO(genelist,
keyType = "ENTREZID",
OrgDb = org.At.tair.db,
ont = "BP", #BP,MF,CC,ALL可選
pAdjustMethod = "BH",
qvalueCutoff = 0.05,
readable = FALSE)
#可視化祟霍，showCategory=20調整顯示的條目多少
dotplot(go_result, showCategory=20)
#將結果轉為數據框，可以選擇導出
goresult <- as.data.frame(go_result)
5.2.6 KEGG通路注釋
#geneClusters必須是數據框類型盈包，ID必須是ENTREZID類型
class(gene.symbol)
compKEGG <- compareCluster(geneClusters=gene.symbol,
fun = "enrichKEGG",
organism = "ath", #在 http://www.genome.jp/kegg/catalog/org_list.html 上查看物種簡寫
pvalueCutoff = 0.05,
pAdjustMethod = "BH")
dotplot(compKEGG, showCategory = 20, title = "KEGG Pathway Enrichment Analysis")
6. homer尋找motif
# 提取peaks的位置信息文件
awk '{print $4"\t"$1"\t"$2"\t"$3"\t+"}' \
SRR7512044_peaks.narrowPeak > SRR7512044.homer_peaks.tmp
# 尋找motif
findMotifsGenome.pl SRR7512044.homer_peaks.tmp \
./TAIR10_genome.fa \ #指定參考基因組
./SRR7512044.motif \ #輸出文件夾的名字
-size 200 -len 8,10,12
- 參數與結果的解讀參考這篇文章【motif分析——從實戰(zhàn)到原理（Homer篇）】
7. DiffBind找差異peak
首先根據我自己的理解將找差異peak的原理說一下：既然是找差異的peak沸呐，那首先就要保證不同樣本間peak所在的區(qū)間基本一致才可以。因此第一步就是找到多個樣本相同的peak區(qū)間续语。第二步垂谢，基于這些相同區(qū)域的peak開始找差異，這個地方的差異和傳統(tǒng)的RNA-Seq找差異基因的原理一模一樣疮茄，都是對落在峰區(qū)間的reads進行計數滥朱，數量不同的話差異就由此而來。但軟件都幫我們做了力试，就不用考慮那么多了
-
我廢話真多徙邻，其實最麻煩的就是準備輸入數據，下面來看怎么做
一般情況相下需要包括以下幾列："SampleID"畸裳， "Tissue"缰犁， "Factor"， "Condition" 怖糊，"Treatment"帅容，"Replicate"， "bamReads" 伍伤，"ControlID" 并徘，"bamControl" ，"Peaks"和"PeakCaller"
-
但我們做ATAC-Seq的哪里有什么
control扰魂，controlID麦乞，bamControl，因此刪掉劝评，而"Tissue"姐直， "Factor"， "Condition" 蒋畜，"Treatment"声畏，"Replicate"，表示的都是分組的意思姻成，因此我保留了兩個"Factor"砰识，"Replicate"，最后csv文件的樣子如下：
SampleID：不能有重復
Factor：進行分組佣渴，后面要做差異分析的依據
Replicate：重復的編號
bamReads：最后生成的bam文件辫狼。全路徑+bam文件名
Peaks：macs2生成的narrowPeak文件。全路徑+narrowPeak文件名
PeakCaller：說明peak類型辛润，使用narrowPeak文件則為narrow膨处。使用xls文件，則用macs砂竖。使用bed文件真椿，則用bed
- 開始正題
# 安裝包
if (!requireNamespace("BiocManager", quietly = TRUE))install.packages("BiocManager")
BiocManager::install("DiffBind")
library(DiffBind)
#導入數據，dba函數會將會將bam文件一同導入乎澄，因此csv中的路徑一定要準確
tamoxifen <- dba(sampleSheet="atDiff.csv")
#查看樣本間的相關性
plot(tamoxifen)
#計算每個樣本中的reads數
tamoxifen <- dba.count(tamoxifen)
#簡單查看計算結果
info <- dba.show(tamoxifen)
libsizes <- cbind(LibReads=info$Reads, FRiP=info$FRiP, PeakReads=round(info$Reads * info$FRiP))
rownames(libsizes) <- info$ID
libsizes
#統(tǒng)計結果
LibReads FRiP PeakReads
cmt2_1 6720462 0.29 1948934
cmt2_2 7747491 0.33 2556672
wild_1 7078309 0.36 2548191
wild_2 8030752 0.34 2730456
ddcc_1 6370232 0.30 1911069
ddcc_2 6956429 0.35 2434750
#默認基于測序深度對數據進行標椎化
tamoxifen <- dba.normalize(tamoxifen)
norm <- dba.normalize(tamoxifen, bRetrieve=TRUE)
norm
#查看進行標準化的文庫大小(library-size)突硝，其實就是各樣本文庫大小總和的平均值，你可以計算驗證一下
normlibs <- cbind(FullLibSize=norm$lib.sizes, NormFacs=norm$norm.factors, NormLibSize=round(norm$lib.sizes/norm$norm.factors))
rownames(normlibs) <- info$ID
normlibs
#標椎化結果
FullLibSize NormFacs NormLibSize
cmt2_1 6720462 0.9398443 7150612
cmt2_2 7747491 1.0834724 7150612
wild_1 7078309 0.9898885 7150612
wild_2 8030752 1.1230859 7150612
ddcc_1 6370232 0.8908652 7150612
ddcc_2 6956429 0.9728438 7150612
#分組置济，格式是表頭在最前面解恰，要分的組依次寫在后面锋八，只能兩兩比較，因此后面只能寫兩組护盈，但可以多執(zhí)行幾次挟纱，都會追加到tamoxifen 中
#分組1，后面使用contrast=1單獨查看
tamoxifen <- dba.contrast(tamoxifen, contrast=c("Factor","ddcc","wild"))
#分組2腐宋，后面使用contrast=2單獨查看
tamoxifen <- dba.contrast(tamoxifen, contrast=c("Factor","cmt2","wild"))
#當然還可有分組3,4,5等紊服，均可以使用contrast=number單獨查看
tamoxifen
#按照分組分別進行差異分析，默認使用DESeq2進行計算胸竞，可以選擇method = DBA_EDGER(edgR)欺嗤，或者兩個都要method = DBA_ALL_METHODS
tamoxifen <- dba.analyze(tamoxifen)
dba.show(tamoxifen, bContrasts=TRUE)
#查看差異分析的結果與導出為csv文件
#查看第1組差異分析的結果
tamoxifen.DB_1 <- dba.report(tamoxifen,contrast=1)
tamoxifen.DB_1
write.csv(tamoxifen.DB_1, file="diffBind_result.csv")
#查看第2組差異分析的結果
tamoxifen.DB_2 <- dba.report(tamoxifen,contrast=2)
tamoxifen.DB_2
write.csv(tamoxifen.DB_2, file="diffBind_result.csv")
##傳統(tǒng)的上下調在找差異peak中稱為“Gain”或“Loss”
#查看第1組的差異peak數量“Fold>0或<0”控制是“Gain”或“Loss”
sum(tamoxifen.DB$Fold>0,contrast=1)
sum(tamoxifen.DB$Fold<0,contrast=1)
#查看第2組的差異peak數量“Fold>0或<0”控制是“Gain”或“Loss”
sum(tamoxifen.DB$Fold>0,contrast=2)
sum(tamoxifen.DB$Fold<0,contrast=2)
- 由于我是使用兩組進行測試的，所以可以使用contrast=1或contrast=2分別查看卫枝，由于都是相同的代碼煎饼，下面進行可視化的時候，只選擇contrast=1進行可視化
#韋恩圖可視化
dba.plotVenn(tamoxifen, contrast=1, bDB=TRUE, bGain=TRUE, bLoss=TRUE, bAll=FALSE)
#PCA
dba.plotPCA(tamoxifen) #這里可以使用所有樣本進行PCA
dba.plotPCA(tamoxifen, contrast=1, label=DBA_FACTOR)#單獨對分組1進行PCA
#MA plots
dba.plotMA(tamoxifen, contrast=1)
#Volcano plots
dba.plotVolcano(tamoxifen, contrast=1)
#Box plots
dba.plotBox(tamoxifen, contrast=1)
#Heatmap
hmap <- colorRampPalette(c("red", "black", "green"))(n = 13)
#對所有樣本的所有的差異peak畫熱圖
dba.plotHeatmap(tamoxifen,correlations=FALSE,scale="row",colScheme = hmap)
#單獨對分組1中所有的差異peak畫熱圖
dba.plotHeatmap(tamoxifen, contrast=1, correlations=FALSE, scale="row", colScheme = hmap)
dba.plotHeatmap(tamoxifen,correlations=FALSE,scale="row",colScheme = hmap)
- 想畫下面的圖剃盾，可是一直提示
腺占，因此還沒有成功，用官方提供的圖展示一下吧痒谴，大家能畫出來可以和我說一聲哦衰伯。這是【官方文檔】可以參考
#Profiling and Profile Heatmaps
dba.plotProfile(tamoxifen)
8. footprint分析
目前我還沒有找到使用于所有物種的軟件，HINT-ATAC用來做footprint分析的話积蔚，應該只支持不幾個物種意鲸，比較雞肋。ATACseqQC應該是一款不錯的軟件尽爆，最近還在探索怎顾。各位如有建議，歡迎提出漱贱。
三槐雾、ATAC-seq與多組學數據聯(lián)合分析
轉錄因子ChIP-seq：由于大部分轉錄因子結合的是染色質開放區(qū)域，所以ATAC-seq的Peak可能和轉錄因子ChIP-seq的Peak存在部分重疊的情況幅狮，而且ATAC-seq得到的Peak長度往往更長募强，因此ATAC-seq數據和轉錄因子ChIP-seq數據可以相互驗證。轉錄因子在ChIP-seq中獨有的Peak暗示這個轉錄因子可能是結合在異染色質區(qū)域的驅動型轉錄因子（Pioneer TFs）崇摄，驅動型轉錄因子隨后招募染色質重塑復合體以及其它轉錄因子開始轉錄擎值。另外，聯(lián)合分析已經報道的ChIP-seq數據可以更準確地分析轉錄因子的足跡逐抑。
組蛋白修飾ChIP-seq：ATAC-seq數據同樣可以和組蛋白修飾ChIP-seq數據進行聯(lián)合分析鸠儿，其中轉錄激活性修飾（H3K4me3，H3K4me1和H3K27ac等）與染色質開放程度呈正相關厕氨，轉錄抑制性修飾（H3K27me3）與染色質開放程度呈負相關进每。聯(lián)合已知的增強子和啟動子之間的相互作用數據也可以幫助構建調控網絡汹粤。
RNA-seq：ATAC-seq數據可以通過聯(lián)合分析RNA-seq數據來發(fā)現哪些差異表達的基因是受染色質可及性調控的，進一步可以推測這些差異表達的基因哪些是受開放染色質中具有motif和footprint的轉錄因子調控的品追，因此ATAC-seq與RNA-seq的聯(lián)合分析有助于破譯基因調控網絡和細胞異質性玄括。
-
前前后后用了1周的時間冯丙，將近段所學匯總于此肉瓦，如有錯誤望各位多多諒解并指出，感謝胃惜！
