這里是佳奧！
獲得了peaks之后我們導(dǎo)入IGV檢查膀哲。
1 deeptools初探
##安裝
conda install -c bioconda/label/cf201901 deeptools
bamCoverage的基本用法
source activate chipseq
bamCoverage -e 170 -bs 10 -b ap2_chip_rep1_2_sorted.bam -o ap2_chip_rep1_2.bw
# ap2_chip_rep1_2_sorted.bam是前期比對(duì)得到的BAM文件
得到的bw文件就可以送去IGV/Jbrowse進(jìn)行可視化往产。 這里的參數(shù)僅使用了-e/--extendReads和-bs/--binSize即拓展了原來(lái)的read長(zhǎng)度，且設(shè)置分箱的大小某宪。其他參數(shù)還有
-
--filterRNAstrand {forward, reverse}: 僅統(tǒng)計(jì)指定正鏈或負(fù)鏈 -
--region/-r CHR:START:END: 選取某個(gè)區(qū)域統(tǒng)計(jì) -
--smoothLength: 通過(guò)使用分箱附近的read對(duì)分箱進(jìn)行平滑化
如果為了其他結(jié)果進(jìn)行比較仿村，還需要進(jìn)行標(biāo)準(zhǔn)化，deeptools提供了如下參數(shù)：
-
--scaleFactor: 縮放系數(shù) - `--normalizeUsingRPKMReads``: Per Kilobase per Million mapped reads (RPKM)標(biāo)準(zhǔn)化
-
--normalizeTo1x: 按照1x測(cè)序深度(reads per genome coverage, RPGC)進(jìn)行標(biāo)準(zhǔn)化 -
--ignoreForNormalization： 指定那些染色體不需要經(jīng)過(guò)標(biāo)準(zhǔn)化
2.1 首先把bam文件轉(zhuǎn)為bw文件
##合并的bam文件
cd mergeBam
ls *.bam |xargs -i samtools index {}
ls *.bam |while read id;do
bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.bw &
done
##去除PCR重復(fù)的bam
cd dup
ls *.bam |xargs -i samtools index {}
ls *.bam |while read id;do
bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.rm.bw &
done
好像原本chipseq環(huán)境是python2的兴喂，deeptools有些不兼容蔼囊，我新建了一個(gè)python3的chipseq2環(huán)境試一試。
QQ截圖20220812113829.png
這一次跑起來(lái)了衣迷。
##合并的bam文件
(chipseq2) root 11:48:09 /home/kaoku/chipseq/mouse_project/mergeBam
$ ls *.bw
Control.bw H2Aub1.bw H3K36me3.bw Ring1B.bw RNAPII_8WG16.bw RNAPII_S2P.bw RNAPII_S5P.bw RNAPII_S5PRepeat.bw RNAPII_S7P.bw
##去除PCR重復(fù)的bam
(chipseq) root 14:11:18 /home/kaoku/chipseq/mouse_project/dup
$ ls *.bw
Control.rm.bw H3K36me3.rm.bw RNAPII_8WG16.rm.bw RNAPII_S5PRepeat.rm.bw RNAPII_S7P.rm.bw
H2Aub1.rm.bw Ring1B.rm.bw RNAPII_S2P.rm.bw RNAPII_S5P.rm.bw
2 IGV可視化
以H2Aub1這個(gè)樣本（去除PCR重復(fù)）為例
-rw-r--r-- 1 root root 926M 8月 11 22:20 H2Aub1_rmdup_control_lambda.bdg
-rw-r--r-- 1 root root 77K 8月 11 22:16 H2Aub1_rmdup_model.r
-rw-r--r-- 1 root root 86K 8月 11 22:20 H2Aub1_rmdup_peaks.narrowPeak
-rw-r--r-- 1 root root 98K 8月 11 22:20 H2Aub1_rmdup_peaks.xls
-rw-r--r-- 1 root root 59K 8月 11 22:20 H2Aub1_rmdup_summits.bed
-rw-r--r-- 1 root root 1.1G 8月 11 22:20 H2Aub1_rmdup_treat_pileup.bdg
我們把下列文件導(dǎo)入IGV：File——Load from File
QQ截圖20220812122910.png
QQ截圖20220812123016.png
QQ截圖20220812132114.png
QQ截圖20220812132122.png
我們?cè)賮?lái)對(duì)比一下去除PCR重復(fù)的和沒(méi)去出的H2Aub1有什么區(qū)別
參考基因組選擇mm10
QQ截圖20220812135444.png
選擇基因BRCA1
QQ截圖20220812140838.png
可以看到這里去除PCR重復(fù)有效果
QQ截圖20220812140927.png
重新導(dǎo)入一批
QQ截圖20220812141709.png
2號(hào)染色體中間的peaks是老鼠染色體特有的压真。
回到Linux看一下跑完了的peaks目錄，這些是軟件認(rèn)為的peaks蘑险。
(chipseq) root 14:24:26 /home/kaoku/chipseq/mouse_project/peaks
$ head H2Aub1_rmdup_summits.bed
chr1 4492669 4492670 H2Aub1_rmdup_peak_1 3.99494
chr1 4493158 4493159 H2Aub1_rmdup_peak_2 8.99601
chr1 4493469 4493470 H2Aub1_rmdup_peak_3 3.99494
chr1 4496552 4496553 H2Aub1_rmdup_peak_4 3.60380
chr1 4497092 4497093 H2Aub1_rmdup_peak_5 5.04333
chr1 4497408 4497409 H2Aub1_rmdup_peak_6 3.99494
chr1 5916554 5916555 H2Aub1_rmdup_peak_7 5.04333
chr1 5917020 5917021 H2Aub1_rmdup_peak_8 5.94642
chr1 12991930 12991931 H2Aub1_rmdup_peak_9 3.13457
chr1 14506458 14506459 H2Aub1_rmdup_peak_10 2.69568
##IGV可識(shí)別的格式
$ awk '{print $1":"$2"-"$3}' H2Aub1_rmdup_summits.bed | head
chr1:4492669-4492670
chr1:4493158-4493159
chr1:4493469-4493470
chr1:4496552-4496553
chr1:4497092-4497093
chr1:4497408-4497409
chr1:5916554-5916555
chr1:5917020-5917021
chr1:12991930-12991931
chr1:14506458-14506459
可以清晰看到第一個(gè)peaks：chr1:4492669-4492670
QQ截圖20220812143007.png
單擊Refseq Genes一欄可以看到窗口
QQ截圖20220812143303.png
SOX17
總結(jié)：
H2Aub1_rmdup_summits.bed生成peaks的位置信息滴肿。
-
把文件導(dǎo)入IGV查看。
H2Aub1.bw Control.bw H2Aub1.merge.bam H2Aub1.merge.rmdup.bam H2Aub1_rmdup_summits.bed IGV內(nèi)選擇小鼠參考基因組mm10 結(jié)合背景Control.bw與H2Aub1.bw比較佃迄。
判斷是否是真的peaks泼差。
下一篇我們繼續(xù)學(xué)習(xí)deeptools在可視化的應(yīng)用。
我們下一篇再見(jiàn)呵俏！
