【原】九月學徒ChIP-seq學習成果展（6萬字總結(jié)）（上篇）

健明 2021-07-14

展開全文

學徒第4周是ChIP-seq數(shù)據(jù)分析實戰(zhàn)訓練，講義大綱文末的閱讀原文，配套視頻在B站：
九月學徒已經(jīng)結(jié)業(yè)，表現(xiàn)還不錯，學了幾個NGS組學數(shù)據(jù)處理加上部分單細胞，隨機安排的文獻數(shù)據(jù)處理圖表復現(xiàn)也完成的還不賴，昨天在生信技能樹的WGCNA代碼就是他寫的；重復一篇WGCNA分析的文章（代碼版）
因為公眾號排版真的是力氣活，同樣我也是逼著學徒自己排版的，反正這輩子就苦這么一回！“作詞作曲”都是他自己！
因為學徒提交作業(yè)實在是太長，超過了微信公眾號發(fā)完限制，所以分成上下文！

首先是上游Chip-Seq

1.文章數(shù)據(jù)選取

這里緊跟jimmy大神的腳本，選擇同樣的文章數(shù)據(jù)，來自于《CARM1 Methylates Chromatin Remodeling Factor BAF155 to Enhance Tumor Progression and Metastasis》

作者首先實驗證明了用small haripin RNA來knockout CARM1 只能達到90%的敲除效果，有趣的是，對CARM1的功能影響非常小，說明只需要極少量的CARM1就可以發(fā)揮很好的作用。

因此作者通過zinc finger nuclease這種基因組編輯技術(shù)設計了100%敲除CARM1的實驗材料，這樣就能比較CARM1有無時各種蛋白被催化狀態(tài)了。

其中SWI/SNF(BAF) chromatin remodeling complex 染色質(zhì)重構(gòu)復合物的一個亞基 BAF155，非常明顯的只有在CARM1這個基因完好無損的細胞系里面才能被正常的甲基化。作者證明了BAF155是CARM1這個基因非常好的一個底物，而且通過巧妙的實驗設計，證明了BAF155這個蛋白的第1064位氨基酸(R) 是 CARM1的作用位點。

因為早就有各種文獻說明了SWI/SNF(BAF) chromatin remodeling complex 染色質(zhì)重構(gòu)復合物在癌癥的重要作用，所以作者也很自然想探究BAF155在癌癥的功能詳情，這里作者選擇的是ChIP-seq技術(shù)。BAF155是作為SWI/SNF(BAF) chromatin remodeling complex 染色質(zhì)重構(gòu)復合物的一個組分，必然可以直接或者間接的結(jié)合DNA。

而ChIP-seq技術(shù)最適合來探究能直接或者間接結(jié)合DNA的蛋白的功能，所以作者構(gòu)造了一種細胞系（MCF7），它的BAF155蛋白的第1064位氨基酸(R) 突變而無法被CARM1這個基因催化而甲基化，然后比較突變的細胞系和野生型細胞系的BAF155的兩個ChIP-seq結(jié)果，這樣就可以研究BAF155是否必須要被CARM1這個基因催化而甲基化后才能行使生物學功能。

作者用me-BAF155特異性抗體+western bloting 證明了正常的野生型MCF7細胞系里面有~74%的BAF155被甲基化。

有一個細胞系SKOV3，可以正常表達除了BAF155之外的其余14種SWI/SNF(BAF) chromatin remodeling complex 染色質(zhì)重構(gòu)復合物，而不管是把突變的細胞系和野生型細胞系的BAF155混在里面都可以促進染色質(zhì)重構(gòu)復合物的組裝，所以甲基化與否并不影響這個染色質(zhì)重構(gòu)復合物的組裝，重點應該研究的是甲基化會影響B(tài)AF155在基因組其它地方結(jié)合。

2.數(shù)據(jù)下載

考慮到原始文件過大（主要是練習的服務器配置不夠），所以只下載了4個文件，就是WT和MUT的input和IP的數(shù)據(jù)。

mkdir -p /project/home/lyang/ChIP-seq_test && cd /project/home/lyang/ChIP-seq_test
mkdir rawData && cd rawData
prefetch SRR1042595 -O ./
prefetch SRR1042596 -O ./
prefetch SRR1042597 -O ./
prefetch SRR1042598 -O ./

3.格式轉(zhuǎn)換——單端測序結(jié)果

mkdir 1.trans_fastq && cd 1.trans_fastq
ln -s /project/home/lyang/ChIP-seq_test/rawData/*sra ./
fastq-dump --gzip -O ~/ChIP-seq_test/1.trans_fastq SRR1042595.sra
fastq-dump --gzip -O ~/ChIP-seq_test/1.trans_fastq SRR1042596.sra
fastq-dump --gzip -O ~/ChIP-seq_test/1.trans_fastq SRR1042597.sra
fastq-dump --gzip -O ~/ChIP-seq_test/1.trans_fastq SRR1042598.sra

4.質(zhì)控fastqc

mkdir /project/home/lyang/ChIP-seq_test/2.fasqc && cd /project/home/lyang/ChIP-seq_test/2.fasqc
ln -s ../1.trans_fastq/*.fastq.gz ./
fastqc -t 2 *

得到html報告可以自行查看，或者實驗multiqc整合報告

5.用multiqc整合報告

cd /project/home/lyang/ChIP-seq_test/2.fasqc
mkdir multiqc
multiqc *zip -o multiqc/

6.過濾低質(zhì)量的fq文件

mkdir /project/home/lyang/ChIP-seq_test/3.trim_galore && cd /project/home/lyang/ChIP-seq_test/3.trim_galore
mkdir multiqc
ln -s ../1.trans_fastq/*.fastq.gz ./

nohup trim_galore --phred33 -q 5 --length 36 --stringency 3 --fastqc -o ./ SRR1042595.fastq.gz &
nohup trim_galore --phred33 -q 5 --length 36 --stringency 3 --fastqc -o ./ SRR1042596.fastq.gz &
nohup trim_galore --phred33 -q 5 --length 36 --stringency 3 --fastqc -o ./ SRR1042597.fastq.gz &
nohup trim_galore --phred33 -q 5 --length 36 --stringency 3 --fastqc -o ./ SRR1042598.fastq.gz &

multiqc *zip -o multiqc/

因為我們只有4個樣本，所以懶得寫循環(huán)，也懶得寫控制腳本

可以看到與修剪之前相比，%Dups和%GC比例都降低了。

7.bowtie2比對

7.1索引的構(gòu)建

文中用到的參考基因組是hg19，所以去下載hg19參考基因組，然后用bowtie2構(gòu)建索引。

mkdir -p /project/home/lyang/refdata/bowtie2/ && cd /project/home/lyang/refdata/bowtie2/
wget -c http://hgdownload.cse./goldenPath/hg19/bigZips/hg19.fa.gz
gunzip hg19.fa.gz
##構(gòu)建索引
nohup bowtie2-build --threads 2 -f hg19.fa hg19 &

7.2進行比對得到bam文件

mkdir /project/home/lyang/ChIP-seq_test/4.mapping && cd /project/home/lyang/ChIP-seq_test/4.mapping
ln -s ../3.trim_galore/*trimmed.fq.gz ./
###可以做成一個腳本文件后用nohup bash xxx.bash &掛到后臺運行
index=/project/home/lyang/refdata/bowtie2/hg19
ls *_trimmed.fq.gz|while read id
do
bowtie2 -x $index -U $id -S ${id%%_trimmed*}.sam
samtools view -bhS -q 30 ${id%%_trimmed*}.sam > ${id%%_trimmed*}.bam
samtools sort ${id%%_trimmed*}.bam -@ 10 -o ${id%%_trimmed*}.sorted.bam
samtools index -@ 10 ${id%%_trimmed*}.sorted.bam
done

上面的samtools步驟很多，可以管道起來讓它簡潔一點，不過需要注意的是服務器內(nèi)存要好！

對bam文件進行QC

ls  *sorted.bam  |xargs -i samtools index {} 
ls  *sorted.bam  | while read id ;do (samtools flagstat $id > $(basename $id ".sorted.bam").stat);done

查看成功率：

cat *stat*|grep %

結(jié)果似乎不錯：

11810228 + 0 mapped (100.00% : N/A)
47955599 + 0 mapped (100.00% : N/A)
17790926 + 0 mapped (100.00% : N/A)
42749895 + 0 mapped (100.00% : N/A)

這個百分百比對成功看起來有點不真實。

8.如果有必要，則合并bam文件（這里不需要）

有些數(shù)據(jù)是同一個樣本的測序數(shù)據(jù)由于被分到了不同的lane中進行測序，所以最后得到的sra文件也分成了數(shù)個。

如果存在這種情況則需要對不同的文件進行合并，但是前提是檢驗了不同bam文件的相關(guān)性非常好！

這里的數(shù)據(jù)單個文件就是一個完整的樣本，所以無需合并，但是還是先把代碼記錄下來，以免日后需要使用：

ls *.bam|sed 's/_[0-9]_trimmed.bam//g' |sort -u |while read id;do samtools merge ../mergeBam/$id.merge.bam $id*.bam ;done

#解釋：
#ls *.bam|sed 's/_[0-9]_trimmed.bam//g' |sort -u 是為了得到所以的樣本名，去重以后的樣本名
#samtools merge ../mergeBam/$id.merge.bam $id*.bam         samtools merge后接上merge后的文件名，再加上所有需要merge的$id*.bam文件

9.關(guān)于PCR重復的去除

經(jīng)過bowtie2比對并將bam文件排序，構(gòu)建索引后得到的文件：

557M Oct  4 10:32 SRR1042595.sorted.bam
4.1M Oct  4 15:55 SRR1042595.sorted.bam.bai
2.6G Oct  4 12:49 SRR1042596.sorted.bam
2.3M Oct  4 15:56 SRR1042596.sorted.bam.bai
889M Oct  4 13:38 SRR1042597.sorted.bam
3.2M Oct  4 15:56 SRR1042597.sorted.bam.bai
2.4G Oct  4 15:36 SRR1042598.sorted.bam
2.3M Oct  4 15:57 SRR1042598.sorted.bam.bai

去除PCR重復并建立索引

mkdir /project/home/lyang/ChIP-seq_test/5.RemovePCRduplicated/ && cd /project/home/lyang/ChIP-seq_test/5.RemovePCRduplicated/
cp ../4.mapping/*sort* ./

ls  *sorted.bam  | while read id ;do (nohup samtools markdup -r $id $(basename $id ".sorted.bam").rmdup.bam & );done
ls  *.rmdup.bam  |xargs -i samtools index {} 

##對rmdup.bam文件進行QC
ls  *.rmdup.bam  | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done

QC質(zhì)控結(jié)果：

cat *stat*|grep %

4159545 + 0 mapped (100.00% : N/A)
46290036 + 0 mapped (100.00% : N/A)
9519997 + 0 mapped (100.00% : N/A)
40969756 + 0 mapped (100.00% : N/A)

去除重復前后比較（主要看文件大?。?/p>ls -lh *bam|cut -d" " -f5-10 230M Oct 4 16:29 SRR1042595.rmdup.bam 557M Oct 4 10:32 SRR1042595.sorted.bam 2.6G Oct 4 16:33 SRR1042596.rmdup.bam 2.6G Oct 4 12:49 SRR1042596.sorted.bam 527M Oct 4 16:30 SRR1042597.rmdup.bam 889M Oct 4 13:38 SRR1042597.sorted.bam 2.3G Oct 4 16:33 SRR1042598.rmdup.bam 2.4G Oct 4 15:36 SRR1042598.sorted.bam

可以看到至少去除了一半的序列。

10.利用macs2來callpeak

10.1 首先了解關(guān)于macs2這個軟件：

macs2包含一系列的子命令，其中最主要的就是callpeak，官方提供了使用實例

macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01

一般而言，我們照葫蘆畫瓢，按照這個實例替換對應部分就行了，介紹一下各個參數(shù)的意義

-t: 實驗組的輸出結(jié)果
-c: 對照組的輸出結(jié)果
-f: -t和-c提供文件的格式，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一個。如果不提供這項，就是自動檢測選擇。
-g: 基因組大小，默認提供了hs, mm, ce, dm選項，不在其中的話，比如說擬南芥，就需要自己提供了。
-n: 輸出文件的前綴名
-B: 會保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
-q: q值，也就是最小的PDR閾值，默認是0.05。q值是根據(jù)p值利用BH計算，也就是多重試驗矯正后的結(jié)果。
-p：這個是p值，指定p值后MACS2就不會用q值了。
-m: 和MFOLD有關(guān)，而MFOLD和MACS預構(gòu)建模型有關(guān)，默認是5：50，MACS會先尋找100多個peak區(qū)構(gòu)建模型，一般不用改，因為我們很大概率上不會懂。

10.2 應用到自己的實例中：

首先確定哪些樣本是當作control來消除噪音，哪些是treatment來真正尋找peaks的。

可以看到SRR104259{4，6，8}，SRR1042600都是作為input的control組。

而SRR104259{3，5，7，9}都是作為treatment組。

但是不知道這里control和treatment組是不是存在一一對應的關(guān)系，就當作是一一對應的關(guān)系往下做吧

mkdir /project/home/lyang/ChIP-seq_test/6.macs2-callpeak && cd /project/home/lyang/ChIP-seq_test/6.macs2-callpeak
ln -s ../5.RemovePCRduplicated/*.bam ./


nohup macs2 callpeak -c SRR1042596.rmdup.bam -t SRR1042595.rmdup.bam -f BAM -g hs -n Xu_MUT_rep2_rmdup &
nohup macs2 callpeak -c SRR1042598.rmdup.bam -t SRR1042597.rmdup.bam -f BAM -g hs -n Xu_WT_rep1_rmdup &

nohup macs2 callpeak -c SRR1042596.sorted.bam -t SRR1042595.sorted.bam -f BAM -g hs -n Xu_MUT_rep2_sort &
nohup macs2 callpeak -c SRR1042598.sorted.bam -t SRR1042597.sorted.bam -f BAM -g hs -n Xu_WT_rep1_sort &

#-c 后面加上control組樣本名
#-t 后面加上treatment組樣本名
#-f {AUTO,BAM,SAM,BED,ELAND,ELANDMULTI,ELANDEXPORT,BOWTIE,BAMPE,BEDPE}
#            后面加上format 輸出文件的格式
#-q QVALUE 后面加上人為設定的q值。默認為0.05
#-B, --bdg    加上這個選項則表示需要輸出extended fragment pileup和local lambda tracks (two files)
#-n NAME  指定輸出文件名
#-g GSIZE    Effective genome size，可以是具體數(shù)字，也可以用簡寫：
#        'hs' for human (2.7e9), 'mm' for mouse (1.87e9), 'ce' for C. elegans (9e7) 
#        'dm' for fruitfly (1.2e8), Default:hs

10.3 macs2結(jié)果輸出

macs2的輸出文件解讀

得到的文件有：

10.3.1 NAME_peaks.xls

雖然后綴名是xls，但實際上就是一個普通的文本文件。包含peak信息的tab分割的文件，前幾行會顯示callpeak時的命令。輸出信息包含：

染色體號
peak起始位點
peak結(jié)束位點
peak區(qū)域長度
peak的峰值位點（summit position）
peak 峰值的屬性（包括pileup峰高和可信度）（pileup height at peak summit, -log10(pvalue) for the peak summit）都是值越高越好
peak的富集倍數(shù)（相對于random Poisson distribution with local lambda）

10.3.2 NAME_peaks.narrowPeak

narrowPeak文件是BED6+4格式，可以上傳到UCSC瀏覽。輸出文件每列信息分別包含：

1；染色體號
2：peak起始位點
3：結(jié)束位點
4：peak name
5：int(-10*log10qvalue)可信度，值越大越好
6 ：正負鏈
7：fold change
8：-log10pvalue可信度，值越大越好
9：-log10qvalue可信度，值越大越好
10：relative summit position to peak start相對于起始位點來說peaks峰值的位置

bed格式的文件和其他Peak文件唯一不同的是——Peak文件中，由于以1為基，而bed文件是以0為基，所以peak的起始位置需要減1才是bed格式的文件。

10.3.3 NAME_summits.bed

BED格式的文件，包含peak的summits位置，第5列是-log10pvalue。如果想找motif，推薦使用此文件。

能夠直接載入UCSC browser，用其他軟件分析時需要去掉第一行。？？？這里有點奇怪

10.3.4 NAME_model.r

NAME_model.r能通過 $ Rscript NAME_model.r作圖，得到是基于你提供數(shù)據(jù)的peak模型。

可以看到非常貼心的幫我們寫好了腳本，直接運行即可出pdf圖。

10.4 macs2背后的統(tǒng)計學原理

參考Chip-seq處理流程

我們對ChIP-seq測序數(shù)據(jù)尋找peaks的本質(zhì)就是得到所有測序數(shù)據(jù)比對在全基因組之后在基因組上測序深度里面尋找比較突出的部分。比如對WES數(shù)據(jù)來說，各個外顯子，或者外顯子的5端到3端，理論上測序深度應該是一致的，都是50X~200X，畫一個測序深度曲線，應該是近似于一條直線。但對我們的ChIP-seq測序數(shù)據(jù)來說，在所捕獲的區(qū)域上面，理論上測序深度是絕對不一樣的，應該是近似于一個山峰。

而那些覆蓋度高的地方，山頂，就是我們的IP所結(jié)合的熱點，也就是我們想要找的peaks。我們通常說的ChIP-seq測序的IP，可以是各個組蛋白上各修飾位點對應的抗體，或者是各種轉(zhuǎn)錄因子的抗體等等。

熱點是這樣一些位置，這些位置多次被測得的read所覆蓋（我們測的是一個細胞群體，read出現(xiàn)次數(shù)多，說明該位置被TF結(jié)合的幾率大）。那么，read數(shù)達到多少才叫多？這就要用到統(tǒng)計檢驗來分析。假設TF在基因組上的分布是沒有任何規(guī)律的，那么，測序得到的read在基因組上的分布也必然是隨機的，某個堿基上覆蓋的read的數(shù)目應該服從二項分布。

TF在基因組上的結(jié)合其實是一個隨機過程，基因組的每個位置其實都有機會結(jié)合某個TF，只是概率不一樣，說白了，peak出現(xiàn)的位置，是TF結(jié)合的熱點，而peak-calling就是為了找到這些熱點。當n很大，p很小時，二項分布可以近似用泊松分布替代，在這里：

n是測序得到的read總數(shù)目，l是單個read的長度，s是基因組的大小。

我們可以算出在某個置信概率（如0.00001）下，隨機情況下，某個堿基上可以覆蓋的read的數(shù)目的最小值，當實際觀察到的read數(shù)目超過這個值（單側(cè)檢驗）時，我們認為該堿基是TF的一個結(jié)合熱點。反過來，針對每一個read數(shù)目，我們也可以算出對應的置信概率P。

但是，這只是一個簡化的模型，實際情況要復雜好多。比如，由于測序、mapping過程內(nèi)在的偏好性，以及不同染色質(zhì)間的差異性，相比全基因組，某些堿基可能內(nèi)在地會被更多的read所覆蓋，這種情況得到的很多peak可能都是假的。

MACS考慮到了這一點，當對某個堿基進行假設檢驗時，MACS只考慮該堿基附近的染色質(zhì)區(qū)段（如10k），此時，上述公式中n表示附近10k區(qū)間內(nèi)的read數(shù)目，s被置為10k。當有對照組實驗（Control，相比實驗組，沒有用抗體捕獲TF，或用了一個通用抗體）存在時，利用Control組的數(shù)據(jù)構(gòu)建泊松分布，當沒有Control時，利用實驗組，稍大一點的局部區(qū)間（比如50k）的數(shù)據(jù)構(gòu)建泊松分布。

read只是跟隨著TF一起沉淀下來的DNA fragment的末端，read的位置并不是真實的TF結(jié)合的位置。所以在peak-calling之前，延伸read是必須的。不同TF大小不一樣，對read延伸的長度也理應不同。我們知道，測得的read最終其實會近似地平均分配到正負鏈上，這樣，對于一個TF結(jié)合熱點而言，read在附近正負鏈上會近似地形成“雙峰”。

MACS會以某個window size掃描基因組，統(tǒng)計每個window里面read的富集程度，然后抽取（比如1000個）合適的（read富集程度適中，過少，無法建立模型，過大，可能反映的只是某種偏好性）window作樣本，建立“雙峰模型”。最后，兩個峰之間的距離就被認為是TF的長度D，每個read將延伸D/2的長度。

當有對照組實驗存在時，MACS會進行兩次peak calling。第一次以實驗組（Treatment）為實驗組，對照組為對照組，第二次顛倒，以實驗組為對照組，對照組為實驗組。之后，MACS對每一個P計算了相應的FDR（False Discovery Rate）值：

Control-p表示第二次peak calling（顛倒的）得到的置信概率小于P的peak的個數(shù)。

Treatment-p表示第一次peak calling得到的置信概率小于P的peak的個數(shù)。

FDR綜合利用了實驗組和對照組的信息，顯然，FDR越小越好。

這一段摘抄自plob的教程。

10.5 自己得到的輸出結(jié)果如下：

##sorted.bam文件得到的結(jié)果
3.5K Xu_MUT_rep2_sort_peaks.narrowPeak
5.2K Xu_MUT_rep2_sort_peaks.xls
2.4K Xu_MUT_rep2_sort_summits.bed

4.7K Xu_WT_rep1_sort_peaks.narrowPeak
6.5K Xu_WT_rep1_sort_peaks.xls
3.2K Xu_WT_rep1_sort_summits.bed


##rmdup.bam文件得到的結(jié)果
3.6K Xu_MUT_rep2_rmdup_peaks.narrowPeak
5.2K Xu_MUT_rep2_rmdup_peaks.xls
2.5K Xu_MUT_rep2_rmdup_summits.bed

4.8K Xu_WT_rep1_rmdup_peaks.narrowPeak
6.6K Xu_WT_rep1_rmdup_peaks.xls
3.3K Xu_WT_rep1_rmdup_summits.bed


###統(tǒng)計去除PCR重復和沒有去重的bam文件找到的peaks數(shù)目差別：結(jié)果居然也是沒有差別
wc -l *bed

43     Xu_MUT_rep2_rmdup_summits.bed
43     Xu_MUT_rep2_sort_summits.bed
58     Xu_WT_rep1_rmdup_summits.bed
58     Xu_WT_rep1_sort_summits.bed

BED包含有3個必須的字段和9個可選字段：

必須字段：

1 chrom - 染色體名字
2 chromStart - 染色體起始位點
3 chromEnd - 染色體終止位點

可選字段：

4 name - 名字
5 score - 分值(0-1000), 用于genome browser展示時上色。
6 strand - 正負鏈，對于ChIPseq數(shù)據(jù)來說，一般沒有正負鏈信息。
7 thickStart - 畫矩形的起點
8 thickEnd - 畫矩形的終點
9 itemRgb - RGB值
10 blockCount - 子元件（比如外顯子）的數(shù)目
11 blockSizes - 子元件的大小
12 blockStarts - 子元件的起始位點

而peak calling的MACS輸出bed文件默認是前5列，其中第五列指的是片段堆積的峰高：

chr1    121485178       121485179       Xu_MUT_rep2_rmdup_peak_1        169.64561
chr10   42385559        42385560        Xu_MUT_rep2_rmdup_peak_2        28.26733
chr10   42387315        42387316        Xu_MUT_rep2_rmdup_peak_3        21.20116
chr10   42392730        42392731        Xu_MUT_rep2_rmdup_peak_4        24.90387
chr10   42599794        42599795        Xu_MUT_rep2_rmdup_peak_5        13.26007
chr11   18127886        18127887        Xu_MUT_rep2_rmdup_peak_6        5.11958
chr12   107104197       107104198       Xu_MUT_rep2_rmdup_peak_7        4.89626
chr16   46392027        46392028        Xu_MUT_rep2_rmdup_peak_8        18.40010
#~~~內(nèi)容過多，無法全部展示~~~

所以在ChIPseeker是畫peak coverage的函數(shù)covplot要有個weightCol的參數(shù)。

11.使用deeptools是進行可視化

deeptools工具介紹

BAM神器--Deeptools使用指南

對于deeptools里的任意子命令，都支持

--help看幫助文檔

--numberOfProcessors/-p設置多核處理

--region/-r CHR:START:END處理部分區(qū)域。

還有一些過濾用參數(shù)部分子命令可用，如ignoreDuplicates,minMappingQuality,samFlagInclude,samFlagExclue.

官方文檔見http://deeptools./en/latest/index.html, 下面按照用法引入不同的工具。

BAM轉(zhuǎn)換為bigWig或bedGraph

BAM文件是SAM的二進制轉(zhuǎn)換版，bigWig是wig或bedGraph的二進制版，存放區(qū)間的坐標軸信息和相關(guān)計分(score)，主要用于在基因組瀏覽器上查看數(shù)據(jù)的連續(xù)密度圖，可用wigToBigWig從wiggle進行轉(zhuǎn)換。

為什么要用bigWig? 主要是因為BAM文件比較大，直接用于展示時對服務器要求較大。因此在GEO上僅會提供bw,即bigWig下載，便于下載和查看。如果真的感興趣，則可以下載原始數(shù)據(jù)進行后續(xù)分析。

deeptools提供bamCoverage和bamCompare進行格式轉(zhuǎn)換。

為了能夠比較不同的樣本，需要對先將基因組分成等寬分箱(bin)，統(tǒng)計每個分箱的read數(shù)，最后得到描述性統(tǒng)計值。

bamCoverage對于單個樣本，可以用SES方法進行標準化。

bamCompare 對于兩個樣本，描述性統(tǒng)計值可以是兩個樣本的比率，或是比率的log2值，或者是差值。

bamCoverage的基本用法

bamCoverage -e 170 -bs 10 -b ap2_chip_rep1_2_sorted.bam -o ap2_chip_rep1_2.bw

得到的bw文件就可以送去IGV/Jbrowse進行可視化。這里的參數(shù)僅使用了-e/--extendReads和-bs/--binSize即拓展了原來的read長度，且設置分箱的大小。其他參數(shù)還有

--outFileFormat|-of {bigwig,bedgraph}：指定輸出文件的格式(default: bigwig)
-b/--bam BAM file：輸入文件格式為bam
-o/--outFileName FILENAME：輸出文件名
-e/--extendReads int：拓展了原來的read長度
-bs/--binSize int：設置分箱的大小
-p/--numberOfProcessors int：設置線程數(shù)
-r/--region CHR:START:END: 選取某個區(qū)域統(tǒng)計
--smoothLength: 通過使用分箱附近的read對分箱進行平滑化

如果為了其他結(jié)果進行比較，還需要進行標準化，deeptools提供了如下參數(shù)：

--scaleFactor: 縮放系數(shù)
--normalizeUsing {RPKM,CPM,BPM,RPGC,None}: 選擇標準化的方法

如果需要以100為分箱，并且標準化到1x，且僅統(tǒng)計某一條染色體區(qū)域的正鏈，輸出格式為bedgraph,那么命令行可以這樣寫：

bamCoverage -e 170 -bs 100 -of bedgraph -r Chr4:12985884:12997458 -b 02-read-alignment/ap2_chip_rep1_1_sorted.bam -o chip.bedgraph

bamCompare和bamCoverage類似，只不過需要提供兩個樣本，并且采用SES方法進行標準化，于是多了--ratio參數(shù)。

把bam文件轉(zhuǎn)為bw文件

cd  /project/home/lyang/ChIP-seq_test/5.RemovePCRduplicated

#ls  *.bam  |xargs -i samtools index {} 進行轉(zhuǎn)換前需要先建立index
#由于電腦配置問題，無法使用nohup command &這種形式循環(huán)運行，所以可以將下面的代碼寫入一個腳本中，再掛到后臺運行。
cat >bam2bw.command
ls *.sorted.bam |while read id;do
bamCoverage -p 10 --normalizeUsing CPM -b $id -o ${id%%.*}.bw
done 
ls *.rmdup.bam |while read id;do
bamCoverage -p 10 --normalizeUsing CPM -b $id -o ${id%%.*}.rmdup.bw
done 
ctrl+c

nohup bash bam2bw.command &
mkdir bwfile
mv *bw bwfile

載入IGV中進行查看：

查看bam文件

找一個樣本SRR1042595.sorted.bam（未去重）和SRR1042595.rmdup.bam（去除PCR重復后）導入IGV中進行查看去除PCR重復前后的變化：

可以看到這個樣本沒有去除重復前有很多PCR重復。

因為這些片段都是經(jīng)過PCR重復擴增而來的，所以是一模一樣的。因為都是一條母鏈復制而來的。

查看同一個樣本的不同格式

bam，bw，bed文件

可以看到在rmdup后，bam文件上變化挺大的，但是奇怪的是在最后的callpeaks的bed文件居然差別不大。且bw文件也相差不大。

查看bed文件

可以先對bed文件進行處理，這樣就可以把結(jié)果直接復制到IGV中，得到peaks的坐標地址了。

awk '{print $1":"$2"-"$3}' Xu_WT_rep1_sort_summits.bed

根據(jù)上面Xu_WT_rep1_sort_summits.bed文件中的結(jié)果去IGV中查看，看到在Xu_WT_rep1樣本中有一個peaks。

再找一個Xu_MUT_rep2_rmdup_summits.bed文件中Xu_MUT_rep2樣本的peaks。

如果igv載入非常慢的時候，可以考慮點擊reload session這個選項：

查看TSS附件信號強度：

背景介紹

在真核cDNA序列詳解有詳細介紹，關(guān)于UTR，CDS，內(nèi)含子，外顯子以及TSS之間的關(guān)系：

啟動子（promoter）：與RNA聚合酶結(jié)合并能起始mRNA合成的序列。

轉(zhuǎn)錄起始點（TSS）：轉(zhuǎn)錄時，mRNA鏈第一個核苷酸相對應DNA鏈上的堿基，通常為一個嘌呤。

UTR（Untranslated Regions)：即非翻譯區(qū)，是信使RNA（mRNA）分子兩端的非編碼片段。

5'-UTR從mRNA起點的甲基化鳥嘌呤核苷酸帽延伸至AUG起始密碼子，3'-UTR從編碼區(qū)末端的終止密碼子延伸至多聚A尾巴（Poly-A）的末端。

下載TSS.bed文件

方法一：[https:///projects/seqminer/files/Reference%20coordinate/](https:///projects/seqminer/files/Reference coordinate/)（結(jié)果證明這個文件似乎有問題，后面用computeMatrix reference-point調(diào)用時總是報錯）

方法二：去UCSC官網(wǎng)（推薦在這里下載）

進入table browser工具后進行選擇：

group這個選項我不太確定，但是大家都選的這個我就沒改了。

track里可以選擇NCBI RefSeq或者UCSC gene，這里選了NCBI

group一般不變

table里如果track選擇NCBI RefSeq，這里就選擇RefSeq；如果track選擇UCSC gene，這里就選knownGene

output format根據(jù)自己的需求選擇

file type returned這里選gzip compressed，這樣就可以下載到壓縮包格式的輸出文件。選text則下載文本格式。

output file一定要寫上一個文件名字，如果為空則后面無法下載，而只能在瀏覽器上查看。

最后點擊get output即可。

這里我選擇了它默認的whole gene。點擊get bed即可開始下載文件。

TSS附件信號強度畫圖

在tss前后2kb的位置批量畫圖

mkdir /project/home/lyang/ChIP-seq_test/hg19.tss.bed && cd /project/home/lyang/ChIP-seq_test/hg19.tss.bed
bed=/project/home/lyang/ChIP-seq_test/hg19.tss.bed/hg19.tss.bed
ln -s ../5.RemovePCRduplicated/bwfile ./

### 在tss前后2kb的位置畫圖
bed=/trainee1/gz15/tss/hg19.tss.bed
for id in bwfile/*bw;
do
file=$(basename $id )
sample=${file%%.bw}
echo $sample

computeMatrix reference-point  --referencePoint TSS  -p 15  \
-b 2000 -a 2000    \
-R  $bed \
-S $id  \
--skipZeros  -o matrix1_${sample}_TSS_2K.gz  \
--outFileSortedRegions regions1_${sample}_TSS_2K.bed

###     both plotHeatmap and plotProfile will use the output from   computeMatrix
plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.png
plotHeatmap -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Heatmap_2K.pdf --plotFileFormat pdf  --dpi 720  
plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.png
plotProfile -m matrix1_${sample}_TSS_2K.gz  -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720 

done 



#computeMatrix reference-point函數(shù)的參數(shù)
#--referencePoint {TSS,TES,center}    The reference point for the plotting could be either the region start (TSS), the region end (TES) or the center of the region
#-b/--beforeRegionStartLength INT bp：區(qū)域前多少bp
#-a/--afterRegionStartLength INT bp：區(qū)域后多少bp
#-R/--regionsFileName File：指定bed文件
#-S/--scoreFileName File：指定bw文件
#--skipZeros：默認為False
#-o/-out/--outFileName OUTFILENAME：指定輸出文件名
#--outFileSortedRegions BED file

#plotHeatmap和plotProfile函數(shù)參數(shù)
#-m/--matrixFile MATRIXFILE：指定matrix名字（computeMatrix的輸出文件名）
#-o/-out/--outFileName OUTFILENAME：指定輸出文件名，根據(jù)后綴自動決定圖像的格式
#--plotFileFormat {"png", "eps", "pdf", "plotly"，"svg"}：指定輸出文件格式
#--dpi DPI：指定dpi
#--perGroup：默認是畫一個樣本的所有區(qū)域群的圖，使用這個選項后則按區(qū)域群畫所有樣本的情況

還可以給多個bed文件來繪圖，還可以畫genebody的圖，因為原理一樣，我就不做過多介紹啦。

關(guān)于genebody，等回去后有空了看看這篇文獻：Gene body methylation can alter gene expression and is a therapeutic target in cancer.

上面的批量代碼其實就是為了統(tǒng)計全基因組范圍的peak在基因特征的分布情況，也就是需要用到computeMatrix計算，用plotHeatmap以熱圖的方式對覆蓋進行可視化，用plotProfile以折線圖的方式展示覆蓋情況。

computeMatrix具有兩個模式:scale-region和reference-point。前者用來信號在一個區(qū)域內(nèi)分布，后者查看信號相對于某一個點的分布情況。無論是那個模式，都有有兩個參數(shù)是必須的，-S是提供bigwig文件，-R是提供基因的注釋信息。還有更多個性化的可視化選項。