221-bam转bigwig的binsize问题,你遇到了吗?

刘小泽写于2020.12.09

0 前言

在ChIP-Seq的分析中,我们常常需要将bam文件转为bigwig,以便后续的可视化及分析。一般使用deeptoolsbamCoverage即可快速实现

官网链接:https://deeptools.readthedocs.io/en/develop/content/tools/bamCoverage.html

bamCoverage的作用

将比对的bam转为基因组 coverage track (bigWig or bedGraph)

coverage的计算

the number of reads per bin( where bins are short consecutive counting windows of a defined size)

重要参数

--scaleFactor 1.092883 -bs 10 --normalizeUsing RPGC --effectiveGenomeSize 2685511504 --extendReads 150
  • scaleFactor就是利用ChIPSeqSpike计算的

  • --binSize, -bs :默认是50,设置的数值越小,窗口越小,整个coverage track的分辨率越高

  • normalizeUsing :几个选项RPKM, CPM, BPM, RPGC, None(默认选择是None)

    • RPKM = Reads Per Kilobase per Million mapped reads 【计算方法 RPKM (per bin) = number of reads per bin / (number of mapped reads (in millions) * bin length (kb))】
    • CPM = Counts Per Million mapped reads, same as CPM in RNA-seq;
    • BPM = Bins Per Million mapped reads, same as TPM in RNA-seq;
    • RPGC = reads per genomic content (1x normalization); 设置这个选项就需要提供一个effective genome size
  • effective genome size :有一些默认的数值,可以看https://deeptools.readthedocs.io/en/latest/content/feature/effectiveGenomeSize.html image-20201127164957687

  • extendReads :默认情况下,read length是不允许延长的,它主要是考虑到一些spliced-read data like RNA-seq;而对于ChIP-seq,本身fragments就是连续比对的,因此可以使用。需要注意:如果reads本身已经超过了这个设定的值,那就不会被延长

最简单的运行方法

bamCoverage -b reads.bam -o coverage.bw

bamCoverage --bam reads.bam -o coverage.bdg -p 10 --outFileFormat bedgraph

其他参数

bamCoverage --bam a.bam -o a.SeqDepthNorm.bw \
    --binSize 10
    --normalizeUsing RPGC
    --effectiveGenomeSize 2685511504
    --ignoreForNormalization chrX
    --extendReads 150
    --scaleFactor 1.092883

对于多个样本

config=~/samples.txt
cat $config |while read i
do
    conf=($i)
  	sample=${conf[0]}
  	# echo $sample
  	
  	bamCoverage --bam align/CHIP_${sample}_rm.mm.bam -o bw/${sample}_scale_rmbl_rpkm.bw \
    --binSize 10 \
    --normalizeUsing RPKM \
    --extendReads 150 \
    --scaleFactor $factor \
    -p 4  \
    -bl mm10.merge.blacklist.bed
done

还可以将bw转为bed

因为后面进行call peaks时,可以使用bam或者bed。这里的bw经过了normalization,因此可以接着使用。通过bigWigToWig(ucsc) + wig2bed (BEDOPS) 工具可以进行转换

# example: H3K27me3-SL_scale_rmbl_rpkm.bigWig
ls *bigWig | while read i;do 
	bigWigToWig $i ${i%%_*}.wig
	#  --zero-indexed (-x) option to convert zero-indexed WIG data
	wig2bed -x < ${i%%_*}.wig > ${i%%_*}.bed
done

1 问题来了~

虽然deeptools会将bam转为bw,也支持自定义binSize(“窗口大小”也就是每隔多少碱基计算一次测序深度)和多种标准化计算,但是如果你看看它的结果

下图我其实是定义了binSize为100,但很显然并不是从0-100,100-200这样的

image-20201209133210570

这是因为:

(来自官网)As you can see, each row corresponds to one region. If consecutive bins have the same number of reads overlapping, they will be merged. 即连续的区间都是相同的测序深度,那么它会自动把这些bins合并在一起。作者也说这样做的目的是为了减少bw文件的空间

当然,如果只是进行可视化还好,没什么影响。但是如果后续需要用到bw进一步计算(后面的推送会提到应用场景),就会报错,说两个bw的坐标是不同的,无法继续进行计算

Screen Shot 2020-12-09 at 1.19.03 PM

2 怎么办呢?

还是想利用bamCoverage,毕竟使用简单,只是查看了它的全部参数,也没有发现有用的信息

然后看到之前也有人提到了这个问题:https://github.com/deeptools/deepTools/issues/907,并且作者Devon也给出了一个解决方案:

image-20201209132300258

https://github.com/deeptools/deepTools/blob/develop/deeptools/writeBedGraph.py#L250-L258

打开一看,其实就是让我们把这段代码的注释去掉即可

image-20201209132341294

3 尝试解决~

我之前是利用conda安装的,因此如何才能找到这个writeBedGraph.py 文件呢?

find ./ -name "writeBedGraph.py"
# /home/miniconda3/envs/epi/lib/python3.6/site-packages/deeptools/writeBedGraph.py

然后使用vim编辑它,小提示:

# 设置vim显示行号
vi ~/.vimrc
# 然后在里面编辑
:set number
# 保存退出

找到第250行和258行代码,将python的注释"""" 去掉,保存退出

之后重新运行

bamCoverage --bam reads.bam -o coverage2.bdg -p 10 --outFileFormat bedgraph

image-20201209133124250

至此,我们转换bw时又多了一个选择

Yunze Liu
Yunze Liu
Bioinformatics Sharer

Co-founder of Bioinfoplanet(生信星球)

Next
Previous

Related