determining read depth for every position in a sequence

Question

[blur]

I am trying to create a table for a specific sequence, about 1900bp long and have for each nucleotide the read depth achieved by sequencing. I saw you can check read depth by:

samtools view file.bam chr 1:100-500 | wc

but that would require 1900 commands (not good...) I am sure that there is a smarter way to do this - can anyone offer any help? Thanks!

Answer 1

[Medhat]

you should use genomeCoverageBed from bedtools

bedtools genomecov computes histograms (default), per-base reports (-d) and BEDGRAPH (-bg) summaries of feature coverage (e.g., aligned sequences) for a given genome.

that is I think exactly what you need

Answer 2

[Ron]

From samtools:

samtools depth input.bam >input.coverage

Refer to this post:

How to plot coverage and depth statistics of a bam file

[blur]

super! worked like a charm!