Genomecov -scale flag usage

Question

[john.gagnon.42]

I am trying to use bedrolls genomecov in order to scale my RNA seq data to reads per million.

I have tried

bedtools genomecov -ibam mybamfile.bam -bg -scale RPM -g mm10.chrom.sizes > myBedGraphfile.BedGraph

which generates a bed graph file with 0 reads across all positions.

I have tried to find documentation on how to use the scale flag but it only seems to provide info on how to scale by a factor such as multiplying all reads by 10.

I'd really appreciate any help with this!

Answer 1

[ATpoint]

The -scale takes an integer or a float, so you need to do the calculation of an appropriate factor yourself. It is nothing more than a factor that all bedGraph values are multiplied by. A simple workaround would be to use a scaling factor that is directly derived by the bam file you used:

TmpScale=$(bc <<< "scale=6;1000000/$(samtools view -f 0 -c mybamfile.bam)")
bedtools genomecov -ibam mybamfile.bam -bg -scale $TmpScale -g mm10.chrom.sizes > myBedGraphfile.BedGraph

The TmpScale counts the number of reads in the BAM and multiplies it by 1mio, rounded to 6 digits. This is then the factor that all values in the bedGraph are multiplied by. Depending on your needs, you can change the flag in SAMtools view, lets say only count the forward read in case of paired-end data (that would then be -f 1).

Answer 2

[Kevin Blighe]

You may also try bamCoverage (http://deeptools.readthedocs.io/en/latest/content/tools/bamCoverage.html#usage-examples-for-rna-seq), which produces RPKM values from a BAM File in bedGraph format. It is designed for both ChIP- and RNA-seq. I have used it to good effect for ChIP-seq, but not yet tried it for RNA-seq.