【ChIP-seq 实战】七、可视化（IGV部分）

2022-08-13 本文已影响0人佳奥

这里是佳奥！

获得了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在可视化的应用。

我们下一篇再见！