Comparing input read subtracted ChIP-read profiles using Deeptools v3.4.3

bamCompare --bamfile1 merged_reps_sampleA.bam --bamfile2 input.bam --scaleFactorsMethod readCount --operation subtract --binSize 50 --blackListFileName blacklist_file.bed -p 8 --effectiveGenomeSize 2900000000 --ignoreDuplicates -o subtracted_input_sampleA.bw

plotProfile --matrixFile SampleA_vs_SampleB_matrix.gz --plotType lines --numPlotsPerRow 1 --perGroup --SampleA_vs_SampleB_matrix_geneProfile_TSS_TES_5k_matrix.tab -o SampleA_vs_SampleB_matrix.gz_geneProfile_TSS_TES_5k.png --regionsLabel "HG19 UCSC canonical genes"

4 months ago by

Devon Ryan ♦ 98k

Freiburg, Germany

Devon Ryan ♦ 98k wrote:

You seem to be using deepTools 3.x rather than deepTools2.

The guts of read count normalization is the following bit of code:

mapped_reads = [bam1_mapped, bam2_mapped]

# new scale_factors (relative to min of two bams)
scale_factors = float(min(bam1_mapped, bam2_mapped)) / np.array(mapped_reads)

bam1_mapped and bam2_mapped contain an estimate of the number of mapped reads in each file after filtering, so the scale factors are 1.0 for the smaller file and something <1 for the larger file. For a given genomic bin, with the number of reads stored in tileCoverage[0] and tileCoverage[1]:

value1 = scale_factors[0] * tileCoverage[0]
value2 = scale_factors[1] * tileCoverage[1]
bin_value = value1 - value2

What the resulting bigWig contains and what you have plotted is bin_value, which is the "depth normalized difference in coverage" in this case. It would seem that sample B had much lower normalized depth in this region (if this is unexpected, then perhaps there are genomic regions that should be excluded (blacklisted)).

ADD COMMENT • link written 4 months ago by Devon Ryan ♦ 98k

Thank you Devon and for pointing out that I should be referring to it as Deeptools2, old habits die hard. Sample B was expected to be lower. Just to clarify, the Y-axis is explained as "the depth normalized difference in coverage", which is the relative difference in coverage after normalisation? Is my use of bamCompare correct for what I am trying to achieve, i.e. showing the difference in histone binding profiles? Does "subtract" reduce the number of normalised ChIP reads (per bin) by the number of normalised input reads (per bin)?

ADD REPLY • link written 4 months ago by Ian ♦ 5.7k

Yes to all of the questions in your comment.

ADD REPLY • link written 4 months ago by Devon Ryan ♦ 98k

Similar posts • Search »