Question: Genomecov -scale flag usage
0
gravatar for john.gagnon.42
3.0 years ago by
john.gagnon.420 wrote:

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!

genomecov bedtools • 2.7k views
ADD COMMENTlink modified 3.0 years ago by Kevin Blighe65k • written 3.0 years ago by john.gagnon.420
2
gravatar for ATpoint
3.0 years ago by
ATpoint38k
Germany
ATpoint38k wrote:

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).

ADD COMMENTlink modified 3.0 years ago • written 3.0 years ago by ATpoint38k
0
gravatar for Kevin Blighe
3.0 years ago by
Kevin Blighe65k
Kevin Blighe65k wrote:

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.

ADD COMMENTlink written 3.0 years ago by Kevin Blighe65k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 617 users visited in the last hour