This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How to quickly count the number of reads at a specific coordinate in a bam file

Hi all,

I am trying to count the number of reads (or alignments) for specific genomic locations in a bam file. I don't want reads with skipped region from the reference. I am using samtools and awk to do it now, but I hope to do it quicker.

samtools view input.bam chr10:18000-20000  | awk '($6 !~ /[N*]/)' | wc -l

Is there a way to make this process faster?

Instead of run samtools+awk for every location, should I just go through the bam files once? Would the alternative be faster?

Thanks, Woody

rna-seq sequencing snp
• 6,992 views • link
updated by dariober dariober • written by klin • klin

3 answers

bam-readcount does this: https://github.com/genome/bam-readcount

• 0 views • link
written by Chris Miller Chris Miller

usingt htsjdk (not tested):

import java.io.File;
import htsjdk.samtools.CigarOperator;
import htsjdk.samtools.SAMRecord;
import htsjdk.samtools.SAMRecordIterator;
import htsjdk.samtools.SamReader;
import htsjdk.samtools.SamReaderFactory;

public class Biostar256909 {
public static void main(String[] args) throws Exception{
    long n=0L;
    SamReader sr= SamReaderFactory.makeDefault().open(new File(args[0]));
    SAMRecordIterator iter=sr.query("chr10",18000,45500,false);
    while(iter.hasNext())
    {
        final SAMRecord rec=iter.next();
        if(!rec.getReadUnmappedFlag() && rec.getCigar().containsOperator(CigarOperator.N))
            {
            continue;
            }
        n++;
    }
    iter.close();
    sr.close();
    System.out.println(n);
}
}

compile using picard as a library:

java -cp /path/to/picard.jar Biostar256909.java

execute:

java -cp /path/to/picard.jar:. Biostar256909 input.bam
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Not tested, but this maybe faster:

samtools view input.bam chr10:18000-20000  | cut -f 6 | grep -v 'N' | wc -l

cut+grep maybe faster than awk.

But in practice I wonder how much difference it makes, if any...

• 0 views • link
written by dariober dariober

Log in to answer this question.

More posts like this