Counting Reads in BAM file by Region

[CodeHippo]

Greetings!

I am a neophyte Bioinformaticist working with SAMtools for the very first time. To make matters worse, I have to work in Java, which is not a language I am particularly familiar with. To make matters even worse, it seems that most individuals run SAMtools functions through command line, while I am working towards compilable code in an IDE instead.

My goal is to create java code that will count the number of reads in a region of a BAM file specified by the user. Could anyone start me off in the right direction by pointing out what specific SAMtools functions I should be looking at to accomplish this? I would be eternally grateful for any help.

[GenoMax]

Generally this is done using bedtools or featurecounts.

http://bedtools.readthedocs.org/en/l.../coverage.html (https://github.com/arq5x/bedtools2)

http://www.ncbi.nlm.nih.gov/pubmed/24227677 (http://subread.sourceforge.net)

You will need samtools to generate the necessary BAM.

[Brian Bushnell]

To make matters worse, I have to work in Java

Nonono, Java makes things BETTER!

But have you considered Picard? It has most or all of the functionality of samtools, and is written in Java, so you'll be able to access the bam-related functions directly.

[CodeHippo]

Nonono, Java makes things BETTER!

But have you considered Picard? It has most or all of the functionality of samtools, and is written in Java, so you'll be able to access the bam-related functions directly.

Picard is actually exactly what I'm using! I've found a very nice piece of Java code that uses Picard to count the total number of reads in a BAM file right here (courtesy of Andrew E. Bruno):

http://left.subtree.org/2012/04/13/counting-the-number-of-reads-in-a-bam-file/

My main issue now I suppose is figuring out how to take it a step further and use Picard to count only the number of reads within a specified region, and not the whole BAM file. I guess the real problem here isn't the Java, but the SAM/BAM terminology and interfaces that I'm not familiar with.

I am searching through the javadoc for Picard and the BAMIndexMetaData Class looks like it has what I'm looking for, but I'm still not entirely sure how it works. Oh well, I guess it's back to flipping burgers for a living.

[mozack]

The methods queryOverlapping() or queryContained() in SAMFileReader should give you want you want:

http://picard.sourceforge.net/javadoc/net/sf/samtools/SAMFileReader.html