这里是佳奥！
找到了peaks以及peaks的分布，我们继续学习deeptools软件的更多内容。
先前我们学习了利用deeptools，把bam文件生成bw文件。
1 准备ucsc.refseq.bed
##需要下载
ucsc.refseq.bed
方法参考：http://www.bio-info-trainee.com/2494.html
网址：https://genome.ucsc.edu/cgi-bin/hgTables
按照博客的做法，先下载.txt文件
QQ截图20220813114954.png
然后导入Linux，使用perl处理成bed文件。
这里定义的TSS(转录起始位点)区域上下游2.5kb，代码如下：
perl -alne '{next if /^#/;if($F[3] eq "+"){$start=$F[4]-2500;$end=$F[4]+2500}else{$start=$F[5]-2500;$end=$F[5]+2500}print join("\t",$F[2],$start,$end,$F[12],0,$F[3])}' ucsc.refseq.txt |sort -u >ucsc.refseq.bed
输出文件ucsc.refseq.bed。
-rw-r--r-- 1 root root 1.2M 8月 13 11:47 ucsc.refseq.bed
-rw-r--r-- 1 kaoku kaoku 13M 8月 13 11:43 ucsc.refseq.txt
##文件位置
/home/kaoku/reference/ucsc.refseq.bed
2 deeptools分析
2.1 查看TSS（转录起始位置）附件信号强度
##首先对单一样本绘图 10K附近
mkdir tss
$ ls
alignraw clean dup mergeBam motif peaks peaksMerge qc raw sra tss
cd tss
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 10000 -a 10000 \
-R /home/kaoku/reference/ucsc.refseq.bed \
-S /home/kaoku/chipseq/mouse_project/mergeBam/H2Aub1.bw \
--skipZeros -o matrix1_test_TSS.gz \
--outFileSortedRegions regions1_test_genes.bed &
查看进程
QQ截图20220813120640.png
##both plotHeatmap and plotProfile will use the output from computeMatrix
plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.png
plotHeatmap -m matrix1_test_TSS.gz -out test_Heatmap.pdf --plotFileFormat pdf --dpi 720
plotProfile -m matrix1_test_TSS.gz -out test_Profile.png
plotProfile -m matrix1_test_TSS.gz -out test_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720
(chipseq2) root 14:53:10 /home/kaoku/chipseq/mouse_project/tss
$ ls -lh
总用量 29M
-rw-r--r-- 1 root root 12M 8月 13 12:12 matrix1_test_TSS.gz
-rw-r--r-- 1 root root 2.6M 8月 13 12:12 regions1_test_genes.bed
-rw-r--r-- 1 root root 13M 8月 13 14:50 test_Heatmap.pdf
-rw-r--r-- 1 root root 1.2M 8月 13 14:50 test_Heatmap.png
-rw-r--r-- 1 root root 386K 8月 13 14:50 test_Profile.pdf
-rw-r--r-- 1 root root 52K 8月 13 14:50 test_Profile.png
test_Heatmap.png
2.2 批处理
##首先画10K附近
bed=/home/kaoku/reference/ucsc.refseq.bed
for id in /home/kaoku/chipseq/mouse_project/mergeBam/*bw ;
do
echo $id
file=$(basename $id )
sample=${file%%.*}
echo $sample
computeMatrix reference-point --referencePoint TSS -p 15 \
-b 10000 -a 10000 \
-R $bed \
-S $id \
--skipZeros -o matrix1_${sample}_TSS_10K.gz \
--outFileSortedRegions regions1_${sample}_TSS_10K.bed
##输出的gz为文件用于plotHeatmap, plotProfile的输入
##both plotHeatmap and plotProfile will use the output from computeMatrix
plotHeatmap -m matrix1_${sample}_TSS_10K.gz -out ${sample}_Heatmap_10K.png
plotHeatmap -m matrix1_${sample}_TSS_10K.gz -out ${sample}_Heatmap_10K.pdf --plotFileFormat pdf --dpi 720
plotProfile -m matrix1_${sample}_TSS_10K.gz -out ${sample}_Profile_10K.png
plotProfile -m matrix1_${sample}_TSS_10K.gz -out ${sample}_Profile_10K.pdf --plotFileFormat pdf --perGroup --dpi 720
done
使用命令批量提交：bash 10k.sh 1>10k.log &
##然后画2K的
bed=/home/kaoku/reference/ucsc.refseq.bed
for id in /home/kaoku/chipseq/mouse_project/mergeBam/*bw ;
do
echo $id
file=$(basename $id )
sample=${file%%.*}
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
使用命令批量提交：bash 2k.sh 1>2k.log & 我是这么跑的：bash 2k.sh &
还可以给多个bed文件来绘图，还可以画genebody的图：scale-regions
上面的批量代码其实就是为了统计全基因组范围的peak在基因特征的分布情况，也就是需要用到computeMatrix计算，用plotHeatmap以热图的方式对覆盖进行可视化，用plotProfile以折线图的方式展示覆盖情况。
这里我只跑了2K的批处理，看一下结果：
(chipseq2) root 15:59:23 /home/kaoku/chipseq/mouse_project/tss
$ ls -lh
总用量 102M
-rw-r--r-- 1 root root 832 8月 13 15:28 2k.sh
-rw-r--r-- 1 root root 6.3M 8月 13 15:33 Control_Heatmap_2K.pdf
-rw-r--r-- 1 root root 1005K 8月 13 15:33 Control_Heatmap_2K.png
-rw-r--r-- 1 root root 13K 8月 13 15:34 Control_Profile_2K.pdf
-rw-r--r-- 1 root root 43K 8月 13 15:34 Control_Profile_2K.png
-rw-r--r-- 1 root root 5.8M 8月 13 15:35 H2Aub1_Heatmap_2K.pdf
-rw-r--r-- 1 root root 911K 8月 13 15:35 H2Aub1_Heatmap_2K.png
-rw-r--r-- 1 root root 13K 8月 13 15:35 H2Aub1_Profile_2K.pdf
-rw-r--r-- 1 root root 36K 8月 13 15:35 H2Aub1_Profile_2K.png
-rw-r--r-- 1 root root 4.4M 8月 13 15:37 H3K36me3_Heatmap_2K.pdf
-rw-r--r-- 1 root root 712K 8月 13 15:36 H3K36me3_Heatmap_2K.png
-rw-r--r-- 1 root root 13K 8月 13 15:37 H3K36me3_Profile_2K.pdf
-rw-r--r-- 1 root root 38K 8月 13 15:37 H3K36me3_Profile_2K.png
略
Control_Profile_2K：
Control_Profile_2K.png
H2Aub1_Profile_2K：
H2Aub1_Profile_2K.png
H3K36me3_Profile_2K：
H3K36me3_Profile_2K.png
有点点怪，没有很明显的一个波谷。这些anti-body表现和control表现太一致了。
上述文件是从mergeBam目录下，这是没有去除PCR重复的样本，我使用去除PCR重复的dup目录再运行一次2k.sh。嗯结果相似。
下一篇我们继续使用R包注释peaks。
我们下一篇再见！
