deeptools工具系列 |
您所在的位置：网站首页 › bw几点到比较好 › deeptools工具系列 |
|
目录
前言 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 fragments1.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的区域，默认不跳过，这些区域设为02.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的方法选择有关。 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |