这里是佳奥！
获得了peaks之后我们导入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是前期比对得到的BAM文件
得到的bw文件就可以送去IGV/Jbrowse进行可视化。 这里的参数仅使用了-e/--extendReads和-bs/--binSize即拓展了原来的read长度，且设置分箱的大小。其他参数还有
-
--filterRNAstrand {forward, reverse}: 仅统计指定正链或负链 -
--region/-r CHR:START:END: 选取某个区域统计 -
--smoothLength: 通过使用分箱附近的read对分箱进行平滑化
如果为了其他结果进行比较，还需要进行标准化，deeptools提供了如下参数：
-
--scaleFactor: 缩放系数 - `--normalizeUsingRPKMReads``: Per Kilobase per Million mapped reads (RPKM)标准化
-
--normalizeTo1x: 按照1x测序深度(reads per genome coverage, RPGC)进行标准化 -
--ignoreForNormalization： 指定那些染色体不需要经过标准化
2.1 首先把bam文件转为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重复的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环境是python2的，deeptools有些不兼容，我新建了一个python3的chipseq2环境试一试。
QQ截图20220812113829.png
这一次跑起来了。
##合并的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重复的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这个样本（去除PCR重复）为例
-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
我们把下列文件导入IGV：File——Load from File
QQ截图20220812122910.png
QQ截图20220812123016.png
QQ截图20220812132114.png
QQ截图20220812132122.png
我们再来对比一下去除PCR重复的和没去出的H2Aub1有什么区别
参考基因组选择mm10
QQ截图20220812135444.png
选择基因BRCA1
QQ截图20220812140838.png
可以看到这里去除PCR重复有效果
QQ截图20220812140927.png
重新导入一批
QQ截图20220812141709.png
2号染色体中间的peaks是老鼠染色体特有的。
回到Linux看一下跑完了的peaks目录，这些是软件认为的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可识别的格式
$ 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
可以清晰看到第一个peaks：chr1:4492669-4492670
QQ截图20220812143007.png
单击Refseq Genes一栏可以看到窗口
QQ截图20220812143303.png
SOX17
总结：
-
H2Aub1_rmdup_summits.bed生成peaks的位置信息。
-
把文件导入IGV查看。
H2Aub1.bw Control.bw H2Aub1.merge.bam H2Aub1.merge.rmdup.bam H2Aub1_rmdup_summits.bed IGV内选择小鼠参考基因组mm10 -
结合背景Control.bw与H2Aub1.bw比较。
-
判断是否是真的peaks。
下一篇我们继续学习deeptools在可视化的应用。
我们下一篇再见！
网友评论