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!
2 answers
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).
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.
Log in to answer this question.