I start with a sorted and indexed bam file ("mapped.bam") representing the mapping of small reads on a reference genome, and a bed file ("genes.bed") containing the coordinates of a set of features of interest (let's say they are genes), for which I want to compute an average profile using programs from deeptools. I would like to understand the steps involved to be sure of what the vertical axis of the final profile represents.

First step: making a bigwig file

I create a bigwig file ("mapped.bw") from the bam file using bamCoverage as follows:

bamCoverage -b mapped.bam -bs 10 -of=bigwig -o mapped.bw

The help of bamCoverage says:

The coverage is calculated as the number of reads per bin, where bins are short consecutive counting windows of a defined size.

In my case, the bins are 10 bp long. My reads are longer than that.

For a given bin, a given read can:

1. completely overlap the bin

2. overlap the bin on n bp, n < 10

3. not overlap the bin at all

Please correct me if I'm wrong: My guess it that the read is counted as 1 in cases 1. and 2., and 0 otherwise, and I also suppose that a read can be counted for several successive bins if it is long enough.

Second step: averaging over genes and plotting

I compute a "meta profile matrix" ("mapped_on_genes.gz") using computeMatrix scale-regions as follows:

computeMatrix scale-regions \
    -S mapped.bw \
    -R genes.bed \
    --upstream 300 \
    --unscaled5prime 500 \
    --regionBodyLength 2000 \
    --unscaled3prime 500 \
    --downstream 300 \
    -out mapped_on_genes.gz

(There is sa -bs parameter which default value is 10 according to the help of the command.)

I use this to plot a profile using plotProfile:

plotProfile -m mapped_on_genes.gz \
    -out mapped_on_genes_meta_profile.pdf

I obtain a profile in with values on the y axis. In what units are these values?

My guess is the following:

For the upstream (300 bp) and internal 5-prime (500 bp), since the bin size was the same in bamCoverage and computeMatrix, each point on the x axis probably represents a 10 bp window, and its y coordinate is the average over the regions present in the bed file of the corresponding bins in the bigwig file, so it is an average number of reads overlapping a 10 bp bin.

Same thing at the 3-prime and downstream side.

For the central 100 bp portion, before averaging over regions some shrinking or spreading of the bins must have been performed, I guess by averaging between neighbouring bins. So the final unit is still a number of reads overlapping a 10 bp bin.

And if I use larger bins, I should end up with proportionally higher values.

Am I correct?