What unit do I get on the y-axis of a metagene profile plot?

Question

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 a -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?

Answer 1

Feel free to @ me in deepTools questions, since I'm the primary developer.

For a given bin, the count assigned to it is the number of reads overlapping it, regardless of whether they overlap by 1 or 10 bases. So a read overlapping only partially and one overlapping completely are treated the same.

Since your bigWig file is in units of "alignments" (i.e., it's not 1x normalized), the resulting profile will also be in units of "alignments" (i.e., profiles and heatmaps are in whatever units the input files are in).

Upstream/downstream regions and unscaled regions are also 10 base bins. Note that these are then the average of the per-base value, since the bins here may not perfectly correspond to the bins in the bigWig files. The line in the profile plot is indeed the average (by default, you can choose median, max, min, etc.) of the underlying regions for each bin.

Regarding the scaled section in the middle, the the number of genomic bases per bin is changed such that the region will have "length"/(regionBodyLength/binSize) bases each. As above, the per-base value is then averaged (or whatever you specify) to derive the per-bin value. The length here is decreased if you have unscaled regions, since otherwise bases would get counted twice.