ChIP-seq数据分析实战训练（三）
六、合并bam文件
在进行peaks calling的时候，需要把bam进行合并。
## 如果不用循环：
## samtools merge control.merge.bam Control_1_trimmed.bam Control_2_trimmed.bam
## 通常我们用批处理。
cd /public/workspace/fangwen/learn/chip-seq/
mkdir mergeBam
cd ./align
ls *.bam|sed 's/_[0-9]_trimmed.bam//g' |sort -u |while read id;do samtools merge ../mergeBam/$id.merge.bam $id*.bam ;done
得到全新的bam文件如下：
(chip) [fangwen@server mergeBam]$ ls -lh *merge.bam|cut -d " " -f 5-100
844M Jul 15 15:18 Control.merge.bam
1.3G Jul 15 15:22 H2Aub1.merge.bam
1.6G Jul 15 15:27 H3K36me3.merge.bam
1.3G Jul 15 15:31 Ring1B.merge.bam
1.1G Jul 15 15:34 RNAPII_8WG16.merge.bam
1.5G Jul 15 15:39 RNAPII_S2P.merge.bam
1.4G Jul 15 15:43 RNAPII_S5P.merge.bam
216M Jul 15 15:44 RNAPII_S5PRepeat.merge.bam
718M Jul 15 15:46 RNAPII_S7P.merge.bam
假如需要去除PCR重复
cd /public/workspace/fangwen/learn/chip-seq/mergeBam
ls *merge.bam | while read id ;do (nohup samtools markdup -r $id $(basename $id ".bam").rmdup.bam & );done
ls *.rmdup.bam |xargs -i samtools index {}
ls *.rmdup.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
去除PCR重复前后比较：
(chip) [fangwen@server mergeBam]$ ls -lh *.bam|cut -d " " -f 5-100
844M Jul 15 15:18 Control.merge.bam
750M Jul 15 15:50 Control.merge.rmdup.bam
1.3G Jul 15 15:22 H2Aub1.merge.bam
1.1G Jul 15 15:51 H2Aub1.merge.rmdup.bam
1.6G Jul 15 15:27 H3K36me3.merge.bam
1.5G Jul 15 15:52 H3K36me3.merge.rmdup.bam
1.3G Jul 15 15:31 Ring1B.merge.bam
1007M Jul 15 15:51 Ring1B.merge.rmdup.bam
1.1G Jul 15 15:34 RNAPII_8WG16.merge.bam
983M Jul 15 15:51 RNAPII_8WG16.merge.rmdup.bam
1.5G Jul 15 15:39 RNAPII_S2P.merge.bam
1.2G Jul 15 15:51 RNAPII_S2P.merge.rmdup.bam
1.4G Jul 15 15:43 RNAPII_S5P.merge.bam
778M Jul 15 15:50 RNAPII_S5P.merge.rmdup.bam
216M Jul 15 15:44 RNAPII_S5PRepeat.merge.bam
211M Jul 15 15:48 RNAPII_S5PRepeat.merge.rmdup.bam
718M Jul 15 15:46 RNAPII_S7P.merge.bam
614M Jul 15 15:50 RNAPII_S7P.merge.rmdup.bam
(chip) [fangwen@server mergeBam]$ grep "%" *.stat
Control.merge.rmdup.stat:12330969 + 0 mapped (85.16% : N/A)
H2Aub1.merge.rmdup.stat:17516222 + 0 mapped (96.82% : N/A)
H3K36me3.merge.rmdup.stat:22685679 + 0 mapped (98.51% : N/A)
Ring1B.merge.rmdup.stat:24901367 + 0 mapped (93.46% : N/A)
RNAPII_8WG16.merge.rmdup.stat:23397509 + 0 mapped (94.84% : N/A)
RNAPII_S2P.merge.rmdup.stat:26655659 + 0 mapped (95.36% : N/A)
RNAPII_S5P.merge.rmdup.stat:13680963 + 0 mapped (90.78% : N/A)
RNAPII_S5PRepeat.merge.rmdup.stat:3997567 + 0 mapped (82.22% : N/A)
RNAPII_S7P.merge.rmdup.stat:9759486 + 0 mapped (77.96% : N/A)
samtools flagstat命令简介：
统计输入文件的相关数据并将这些数据输出至屏幕显示。每一项统计数据都由两部分组成，分别是QC pass和QC failed，表示通过QC的reads数据量和未通过QC的reads数量。以“PASS + FAILED”格式显示。
命令格式：
samtools flagstat <in.bam> |<in.sam> | <in.cram>
七、 使用macs2进行找peaks
macs2包含一系列的子命令，其中最主要的就是callpeak， 官方提供了使用实例
macs2 callpeak -t ChIP.bam -c Control.bam -f BAM -g hs -n test -B -q 0.01
各个参数的意义：
- -t: 实验组的输出结果
- -c: 对照组的输出结果
- -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值是根据p值利用BH计算，也就是多重试验矫正后的结果。
- -p： 这个是p值，指定p值后MACS2就不会用q值了。
- -m: 和MFOLD有关，而MFOLD和MACS预构建模型有关，默认是5：50，MACS会先寻找100多个peak区构建模型，一般不用改，因为你很大概率上不会懂。
##方法一
cd /public/workspace/fangwen/learn/chip-seq/mergeBam
ls *merge.bam |cut -d"." -f 1 |while read id;do (nohup macs2 callpeak -c Control.merge.bam -t $id.merge.bam -f BAM -g mm -n $id 2 > $id.log & )；done
##方法二
cd /public/workspace/fangwen/learn/chip-seq/mergeBam
ls *merge.bam |cut -d"." -f 1 |while read id;
do
if [ ! -s ${id}_summits.bed ];
then
echo $id
nohup macs2 callpeak -c Control.merge.bam -t $id.merge.bam -f BAM -B -g mm -n $id --outdir ../peaks 2> $id.log &
fi
done
#将.merge.rmdup.bam分开放
mkdir dup
mv *rmdup* dup/
cd dup/
ls *.merge.rmdup.bam |cut -d"." -f 1 |while read id;
do
if [ ! -s ${id}_rmdup_summits.bed ];
then
echo $id
nohup macs2 callpeak -c Control.merge.rmdup.bam -t $id.merge.rmdup.bam -f BAM -B -g mm -n ${id}_rmdup --outdir ../peaks 2> $id.log &
fi
done
其实上面的-B 参数意义也不大，得到的bedgraph文件没啥用。
得到的bed格式的peaks文件的行数如下：
(chip) [fangwen@server peaks]$ wc -l *bed
0 Control_summits.bed
1115 H2Aub1_summits.bed
40830 H3K36me3_summits.bed
26053 Ring1B_summits.bed
41864 RNAPII_8WG16_summits.bed
20042 RNAPII_S2P_summits.bed
38663 RNAPII_S5PRepeat_summits.bed
62765 RNAPII_S5P_summits.bed
72640 RNAPII_S7P_summits.bed
(chip) [fangwen@server peaks]$ wc -l *bed
0 Control_rmdup_summits.bed
1115 H2Aub1_rmdup_summits.bed
40830 H3K36me3_rmdup_summits.bed
26053 Ring1B_rmdup_summits.bed
41841 RNAPII_8WG16_rmdup_summits.bed
20020 RNAPII_S2P_rmdup_summits.bed
38663 RNAPII_S5PRepeat_rmdup_summits.bed
62765 RNAPII_S5P_rmdup_summits.bed
72577 RNAPII_S7P_rmdup_summits.bed
303864 total
因为MockIP是control，所以它自己跟自己比较，肯定是没有peaks的。
去重前后peaks数几乎相同：macs2会自动进行mark duplicate？
八、使用deeptools进行可视化
deeptools提供bamCoverage和bamCompare进行格式转换，为了能够比较不同的样本，需要对先将基因组分成等宽分箱(bin)，统计每个分箱的read数，最后得到描述性统计值。对于两个样本，描述性统计值可以是两个样本的比率，或是比率的log2值，或者是差值。如果是单个样本，可以用SES方法进行标准化。
#`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文件
##建索引
ls *bam |xargs -i samtools index {}
##首先把bam文件转为bw文件
ls *.bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.bw &
done
##去除PCR重复的bam文件转为bw文件
cd dup
ls *.bam |while read id;do
nohup bamCoverage --normalizeUsing CPM -b $id -o ${id%%.*}.rm.bw &
done
得到的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： 指定那些染色体不需要经过标准化
如果需要以100为分箱，并且标准化到1x，且仅统计某一条染色体区域的正链，输出格式为bedgraph,那么命令行可以这样写
bamCoverage -e 170 -bs 100 -of bedgraph -r Chr4:12985884:12997458 --normalizeTo1x 100000000 -b 02-read-alignment/ap2_chip_rep1_1_sorted.bam -o chip.bedgraph
bamCompare和bamCoverage类似，只不过需要提供两个样本，并且采用SES方法进行标准化，于是多了--ratio参数。
查看TSS附件信号强度：
首先对单一样本绘图：
## both -R and -S can accept multiple files
mkdir -p /public/workspace/fangwen/learn/chip-seq/mergeBam/tss
cd /public/workspace/fangwen/learn/chip-seq/mergeBam/tss
computeMatrix reference-point --referencePoint TSS -p 5 \
-b 10000 -a 10000 \
-R /public/workspace/fangwen/learn/chip-seq/referece/ref/mm10.tss.bed \
-S /public/workspace/fangwen/learn/chip-seq/mergeBam/H2Aub1.bw \
--skipZeros -o matrix1_H2Aub1_TSS.gz \
--outFileSortedRegions regions1_H2Aub1_genes.bed
## both plotHeatmap and plotProfile will use the output from computeMatrix
plotHeatmap -m matrix1_H2Aub1_TSS.gz -out H2Aub1_Heatmap.png
plotHeatmap -m matrix1_H2Aub1_TSS.gz -out H2Aub1_Heatmap.pdf --plotFileFormat pdf --dpi 720
plotProfile -m matrix1_H2Aub1_TSS.gz -out H2Aub1_Profile.png
plotProfile -m matrix1_H2Aub1_TSS.gz -out H2Aub1_Profile.pdf --plotFileFormat pdf --perGroup --dpi 720
image.png
批处理
首先画10K附近
bed=/public/workspace/fangwen/learn/chip-seq/referece/ref/mm10.tss.bed
for id in /public/workspace/fangwen/learn/chip-seq/mergeBam/*bw ;
do
echo $id
file=$(basename $id )
sample=${file%%.*}
echo $sample
computeMatrix reference-point --referencePoint TSS -p 5 \
-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
使用命令批量提交：nohup bash 10k.sh 1>10k.log &（cat >10k.sh）
然后画2k附近
bed=/public/workspace/fangwen/learn/chip-seq/referece/ref/mm10.tss.bed
for id in /public/workspace/fangwen/learn/chip-seq/mergeBam/*bw ;
do
echo $id
file=$(basename $id )
sample=${file%%.*}
echo $sample
computeMatrix reference-point --referencePoint TSS -p 5 \
-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
使用命令批量提交：nohup bash 2k.sh 1>2k.log &（cat >2k.sh）
(base) [fangwen@server tss1]$ ls -lh |cut -d " " -f 5-100
14M Jul 17 21:18 Control_Heatmap_10K.pdf
1.4M Jul 17 21:17 Control_Heatmap_10K.png
6.6M Jul 17 20:42 Control_Heatmap_2K.pdf
1020K Jul 17 20:41 Control_Heatmap_2K.png
387K Jul 17 21:18 Control_Profile_10K.pdf
59K Jul 17 21:18 Control_Profile_10K.png
378K Jul 17 20:42 Control_Profile_2K.pdf
44K Jul 17 20:42 Control_Profile_2K.png
6.2M Jul 17 20:45 H2Aub1_Heatmap_2K.pdf
925K Jul 17 20:45 H2Aub1_Heatmap_2K.png
378K Jul 17 20:45 H2Aub1_Profile_2K.pdf
36K Jul 17 20:45 H2Aub1_Profile_2K.png
4.8M Jul 17 20:49 H3K36me3_Heatmap_2K.pdf
722K Jul 17 20:48 H3K36me3_Heatmap_2K.png
378K Jul 17 20:49 H3K36me3_Profile_2K.pdf
38K Jul 17 20:49 H3K36me3_Profile_2K.png
8.8M Jul 17 21:17 matrix1_Control_TSS_10K.gz
2.7M Jul 17 20:41 matrix1_Control_TSS_2K.gz
上面的批量代码其实就是为了统计全基因组范围的peak在基因特征的分布情况，也就是需要用到computeMatrix计算，用plotHeatmap以热图的方式对覆盖进行可视化，用plotProfile以折线图的方式展示覆盖情况。
computeMatrix具有两个模式:scale-region和reference-point。前者用来信号在一个区域内分布，后者查看信号相对于某一个点的分布情况。无论是那个模式，都有有两个参数是必须的，-S是提供bigwig文件，-R是提供基因的注释信息。还有更多个性化的可视化选项。
image.png
scale-regions模式
computeMatrix scale-regions \ # 选择模式
-b 3000 -a 5000 \ # 感兴趣的区域，-b上游，-a下游
-R ~/reference/gtf/TAIR10/TAIR10_GFF3_genes.bed \
-S 03-read-coverage/ap2_chip_rep1_1.bw \
--skipZeros \
--outFileNameMatrix 03-read-coverage/matrix1_ap2_chip_rep1_1_scaled.tab \ # 输出为文件用于plotHeatmap, plotProfile
--outFileSortedRegions 03-read-coverage/regions1_ap2_chip_re1_1_genes.bed
reference-point模式
computeMatrix reference-point \ # 选择模式
--referencePoint TSS \ # 选择参考点: TES, center
-b 3000 -a 5000 \ # 感兴趣的区域，-b上游，-a下游
-R ~/reference/gtf/TAIR10/TAIR10_GFF3_genes.bed \
-S 03-read-coverage/ap2_chip_rep1_1.bw \
--skipZeros \
-out 03-read-coverage/matrix1_ap2_chip_rep1_1_TSS.gz \ # 输出为文件用于plotHeatmap, plotProfile
--outFileSortedRegions 03-read-coverage/ons1regions1_ap2_chip_re1_1_genes.bed
结果可视化
可视化的方法有两种，一种是轮廓图，一种是热图。两则都提供了足够多的参数对结果进行细节上的修改。
plotProfile -m matrix1_ap2_chip_rep1_1_TSS.gz \
-out ExampleProfile1.png \
--numPlotsPerRow 2 \
--plotTitle "Test data profile"
plotHeatmap -m matrix1_ap2_chip_rep1_1_TSS.gz \
-out ExampleHeatmap1.png \
image.png
参考资料：
https://blog.csdn.net/u013553061/article/details/53402232
https://vip.biotrainee.com/d/226
bam文件转为bw文件：http://www.bio-info-trainee.com/1815.html
IGV:https://www.nature.com/articles/nbt.1754
https://www.jianshu.com/p/b494426ecd32
网友评论