counting reads per user-defined region. Bedtools intersect?

[mcnach]

I have the results of my first few experiments, and aligned them to my genome etc, but now I would like to count the number of reads I get in a user-defined set of regions.
I am reading about Bedtools and it seems to me that I should be able to use bedtools intersect to obtain this information, does this seem the right way to go?

For example, I have my genome alignment file (as a bam file or summarised into bed/bedgraph), and I want to divide the whole genome into fragments flanked by the motif GATC. I will generate a bed file with the coordinates of every "inter-GATC" fragment, and then I want to ask how many reads I get per fragment. Is bedtools intersect the way to do it? Or are there better ways?

I'm completely new to this, as you may have guessed... just trying to get my head around all the tools available.

Thanks for your help.

[crazyhottommy]

I have the results of my first few experiments, and aligned them to my genome etc, but now I would like to count the number of reads I get in a user-defined set of regions.
I am reading about Bedtools and it seems to me that I should be able to use bedtools intersect to obtain this information, does this seem the right way to go?

For example, I have my genome alignment file (as a bam file or summarised into bed/bedgraph), and I want to divide the whole genome into fragments flanked by the motif GATC. I will generate a bed file with the coordinates of every "inter-GATC" fragment, and then I want to ask how many reads I get per fragment. Is bedtools intersect the way to do it? Or are there better ways?

I'm completely new to this, as you may have guessed... just trying to get my head around all the tools available.

Thanks for your help.

I think you need to use coveragebed http://bedtools.readthedocs.org/en/l.../coverage.html

[mcnach]

I think you need to use coveragebed http://bedtools.readthedocs.org/en/l.../coverage.html

Yes, thank you, that looks more appropriate.

Now I have to figure out how to make it work

[mcnach]

Indeed...

coverageBed -abam <bamfile.bam> -b <ref.bed>

works beautifully.

I tested it with some made up bed files first and it's just what I was after.

Thank you!

[crazyhottommy]

Indeed...

coverageBed -abam <bamfile.bam> -b <ref.bed>

works beautifully.

I tested it with some made up bed files first and it's just what I was after.

Thank you!

Not a problem and good luck!

[Gordon Smyth]

I have the results of my first few experiments, and aligned them to my genome etc, but now I would like to count the number of reads I get in a user-defined set of regions.
...
I'm completely new to this, as you may have guessed... just trying to get my head around all the tools available.

Here's a new tool that you could consider, with some comparisons:

http://www.ncbi.nlm.nih.gov/pubmed/24227677

[mcnach]

Here's a new tool that you could consider, with some comparisons:

http://www.ncbi.nlm.nih.gov/pubmed/24227677

Thanks, Gordon. I'll check it out. I've been a fan of Limma for a long time, for microarray analyses, so I'm very interested to see what you came up with here

Jose