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
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:
completely overlap the bin
overlap the bin on n bp, n < 10
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 -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
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?