counting the sequencing reads in 10kb regions on a genome-wide scale

Question

[Bogdan]

Dear all,

during a step of a CNV analysis on cancer genomes, I would like to be able to COUNT and DISPLAY along the CHROMOSOME AXIS the NUMBER of READS in 10KB REGIONS of GERMLINE and CANCER GENOMES.

I would like to ask you please for your suggestions about :

-- any fast function/algorithm that COUNTS the reads from a BAM file in 10kb windows ?

-- a package/function that DISPLAY the COUNTS along the chromosome axes ?

many thanks,

bogdan

[Bogdan]

Thank you gentlemen.

Possibly, could I do it also in R, by using Rsamtools or a related package to COUNT the reads in specific genome WINDOWS ?

[h.mon]

Yes, you could. But then you wouldn't match the speed:

In our tests, we can estimate depth across 60X genomes for 30 samples in 30 seconds.

Answer 1

[igor]

You should be able to do it with bedtools. First generate a BED file of 10kb regions with bedtools makewindows. Then calculate coverage for those regions with bedtools coverage.

Alternatively, you can try deepTools bamCoverage to generate a coverage track (you can specify the size of bins).

Of course, whatever package you are using for CNV calling should have an option to generate some sort of graphical output.

Answer 2

[h.mon]

Have a look at indexcov, but its intervals are of 16kb.