目录
前言
- bamCompare 工具
1.1 工具原理
1.2 使用说明
1.3 使用实例 - bigwigCompare 工具
2.1 使用说明
2.2 使用实例
前言
今天主要介绍下，如果我有2个 bw/bam 文件，我想要比较这2个文件信号之间的差异时，要如何进行分析。
这个听上去是不是很像差异分析？对的，其实就是差异分析的简化版。
而这个的实际用途，可以是：分别有实验组和对照组ATAC-seq的 bam/bw 文件，此时我们可以直接通过 deeptools 工具，对某些感兴趣区域的染色质开放程度进行可视化展示。
1. bamCompare 工具
1.1 工具原理
该工具的算法主要分成2步：
Normalization：通过使用不同的方法对测序深度进行normalization，让不同样本之间可比。
以bin为单位，对两个不同数据之间的信号进行比较差异。
关于normalization的方法，大致又可以分成2类：
对每个样本计算scale factor，代表方法：SES（Signal Extraction Scaling）【PMID: 22499706】
样本内部normalization，代表方法：RPKM，BPM，CPM等
在该方法中，如果是PE数据，则数据会被当作两个SE数据独立地处理。
1.2 使用说明
usage: bamCompare -b1 treatment.bam -b2 control.bam -o log2ratio.bw
# 必需参数
--bamfile1, -b1 BAM1 # 排序后的bam文件,一般是实验组样本
--bamfile2, -b2 BAM2 # 排序后的bam文件,一般是对照组样本
--outFileName, -o FILENAME # 输出文件名
# 可选参数
--outFileFormat, -of {bigwig,bedgraph} # 输出文件格式,bigwig或bedgraph
--binSize, -bs INT bp # bin的size,默认50
--region, -r CHR:START:END # 指定某个区域,一般用于测试代码
--blackListFileName, -bl BED # 指定blackList
--numberOfProcessors, -p INT # 指定线程
# scale factor相关参数
--scaleFactorsMethod {readCount,SES,None} # 计算scale factor的方法,如果设置为None,则--normalizeUsing参数生效. 默认readCount
--sampleLength, -l SAMPLELENGTH # 当选择SES计算scale factor时,随机抽取SAMPLELENGTH长度基因组去计算. 默认1000
--numberOfSamples, -n NUMBEROFSAMPLES # 当选择SES计算scale factor时,随机抽取NUMBEROFSAMPLES个基因组长度去计算. 默认100000
--scaleFactors SCALEFACTORS # 手动指定scale factor,例如0.7:1,则会将BAM1信号全部乘以0.7
--operation {log2,ratio,subtract,add,mean,reciprocal_ratio,first,second} # 输出结果,默认输出log2(BAM1/BAM2), 即log2
--pseudocount PSEUDOCOUNT # 添加一个较小的数值避免分母为0，默认是1
--skipZeroOverZero # 跳过2个bam文件都是0的bin
# 内部normalization相关参数
--effectiveGenomeSize EFFECTIVEGENOMESIZE # 基因组长度,具体见下
--normalizeUsing {RPKM,CPM,BPM,RPGC,None} # normalization方法选择,具体见下
--exactScaling # 利用全部数据进行计算精确的scaling factor. 默认False
--ignoreForNormalization, -ignore IGNOREFORNORMALIZATION # 计算normalization时不考虑的染色体号，不同染色体号之间用空格间隔。对于ChIP-seq数据可以设置chrX chrM
--skipNonCoveredRegions, --skipNAs # 是否跳过没有mapping的区域，默认不跳过，这些区域设为0
# Read相关处理
--extendReads, -e [INT bp] # 对于SE数据,根据实验过程中打断的片段大小则直接填写该数值.对于PE数据无需指定
--ignoreDuplicates # 忽略重复序列
--minMappingQuality INT # 至少达到最低mapping质量得分的read才被考虑
--centerReads # 相对于片段长度，reads处于中心位置。
--samFlagInclude INT # 根据bam文件Flag进行过滤
--samFlagExclude INT # 根据bam文件Flag进行过滤
--minFragmentLength # Fragment最小长度,主要用于ATACseq中过滤mono-/di-nucleosome fragments
--maxFragmentLength # Fragment最大长度,主要用于ATACseq中过滤mono-/di-nucleosome fragments
1.3 使用实例
bamCompare -b1 A.bam -b2 B.bam \
--scaleFactorsMethod readCount --operation log2 \
--outFileName log2ratio.bw \
--binSize 50 \
--region chr10 -p 16 \
-ignore chrX
为了对结果有更好的理解，我们把 bw 文件转为可以查看的 bedgraph 文件：
$ bigWigToBedGraph log2ratio.bw log2ratio.bedgraph
$ head log2ratio.bedgraph
chr10 50 49350 0
chr10 49350 49450 -0.802724
chr10 49450 49750 0
chr10 49750 49800 1
chr10 49800 50050 0
chr10 50050 50100 -0.802724
chr10 50100 50750 0
chr10 50750 50800 -0.802724
chr10 50800 50850 1
chr10 50850 51600 0
一般来说，做到这里，我们接下来可以对这个 bw 进行可视化，也就是画成热图来看：
computeMatrix reference-point \
--referencePoint TSS \
--scoreFileName log2ratiw.bw \
--regionsFileName gencode.v40.annotation.gtf \
--outFileName log2ratio.tss.gz \
--samplesLabel log2ratio \
--binSize 50 \
-a 5000 -b 5000 \
--averageTypeBins mean --numberOfProcessors 16 \
--skipZeros
plotProfile -m log2ratio.tss.gz -out log2ratio.png
这里的话，我们就可以根据这个结果，得到 “A样本的信号在TSS比B样本信号有显著下调” 的结论。
当然了，我们也同样画出A和B样本的信号分布图，用来验证：
# get bw file
for id in A B
do
bamCoverage --bam ${id}.bam --outFileName ${id}.bw \
--outFileFormat bigwig --binSize 50 \
--normalizeUsing CPM \
--region chr10 \
--numberOfProcessors 20
done
# get plot
for id in A B
do
computeMatrix reference-point \
--referencePoint TSS \
--scoreFileName ${id}.bw \
--regionsFileName gencode.v40.annotation.gtf \
--outFileName ${id}.tss.gz \
--samplesLabel ${id} \
--binSize 50 \
-a 5000 -b 5000 \
--averageTypeBins mean --numberOfProcessors 16 \
--skipZeros
plotProfile -m ${id}.tss.gz -out ${id}.png &
done
5d7cf687ea3293439fba8cd991795044.png
可以看到，B的信号整体比A更高，且在TSS初，因为B达到了最高信号，所以A/B的信号达到最低，和前面的结果是一致的。
并且该结果是在使用另外一种normalization方式的前提下完成的，更加说明了结果的一致性。
2. bigwigCompare 工具
2.1 使用说明
usage: bigwigCompare --bigwig1 Bigwig1 --bigwig2 Bigwig2
# 必需参数
--bigwig1, -b1 Bigwig1 # 排序后的bam文件,一般是实验组样本
--bigwig1, -b2 Bigwig2 # 排序后的bam文件,一般是对照组样本
--outFileName, -o FILENAME # 输出文件名
# 可选参数
--outFileFormat, -of {bigwig,bedgraph} # 输出文件格式,bigwig或bedgraph
--binSize, -bs INT bp # bin的size,默认50
--region, -r CHR:START:END # 指定某个区域,一般用于测试代码
--blackListFileName, -bl BED # 指定blackList
--numberOfProcessors, -p INT # 指定线程
# scale factor相关参数
--scaleFactors SCALEFACTORS # 手动指定scale factor,例如0.7:1,则会将BAM1信号全部乘以0.7
--operation {log2,ratio,subtract,add,mean,reciprocal_ratio,first,second} # 输出结果,默认输出log2(BAM1/BAM2), 即log2
--pseudocount PSEUDOCOUNT # 添加一个较小的数值避免分母为0，默认是1
--skipZeroOverZero # 跳过2个bam文件都是0的bin
--skipNonCoveredRegions, --skipNAs # 是否跳过没有mapping的区域，默认不跳过，这些区域设为0
2.2 使用实例
bigwigCompare -b1 A.bw -b2 B.bw \
--operation log2 \
--outFileName log2ratio_bigwigCompare.bw \
--binSize 50 \
--region chr10 -p 16
为了对结果有更好的理解，我们把 bw 文件转为可以查看的 bedgraph 文件：
$ bigWigToBedGraph log2ratio_bigwigCompare.bw log2ratio_bigwigCompare.bedgraph
$ head log2ratio_bigwigCompare.bedgraph
chr10 50 49350 0
chr10 49350 49450 -0.802724
chr10 49450 49750 0
chr10 49750 49800 1
chr10 49800 50050 0
chr10 50050 50100 -0.802724
chr10 50100 50750 0
chr10 50750 50800 -0.802724
chr10 50800 50850 1
chr10 50850 51600 0
$ head log2ratio.bedgraph
chr10 50 49350 0
chr10 49350 49450 -0.151532
chr10 49450 49750 0
chr10 49750 49800 0.184749
chr10 49800 50050 0
chr10 50050 50100 -0.151532
chr10 50100 50750 0
chr10 50750 50800 -0.151532
chr10 50800 50850 0.18
可以看到，这里的结果虽然有一定的差异，但是整体是一致的。
为了更好的展示，我还是将其画成图：
computeMatrix reference-point \
--referencePoint TSS \
--scoreFileName log2ratio_bigwigCompare.bw \
--regionsFileName gencode.v40.annotation.gtf \
--outFileName log2ratio_bigwigCompare.tss.gz \
--samplesLabel log2ratio \
--binSize 50 \
-a 5000 -b 5000 \
--averageTypeBins mean --numberOfProcessors 16 \
--skipZeros
plotProfile -m log2ratio_bigwigCompare.tss.gz -out log2ratio_bigwigCompare.png
720af44981ae500e71128bae7abe5235.png
同样的，我把两张图放在一起看：
5b0a4e3b5293ea7bc7eac7a4530a4c93.png
可以看到，两者的pattern是一模一样的，只是在数值上有一定的差异。这可能与normalization的方法选择有关。
