Samlocusiterator Skips Every Second Locus

Question

[Johan]

I'm writing some Scala code which uses the SamLocusIterator after a excellent tip here: A: Java fasta reference iterator

However I'm getting some pretty weird output, which makes it seem like the iterator is skipping every other locus. Here's an example:

pos: chr1:99984
pos: chr1:99986
pos: chr1:99988
pos: chr1:99990
pos: chr1:99992
pos: chr1:99994
pos: chr1:99996
pos: chr1:99998
pos: chr1:100000

I can replicate this both in Scala and in Java with the following example programs:

The java version

import java.io.File;
import java.util.Iterator;

import net.sf.picard.util.SamLocusIterator;
import net.sf.picard.util.SamLocusIterator.LocusInfo;
import net.sf.samtools.SAMFileReader;


public class JavaMinimalIteratorTest {
    /**
     * @param args
     */
    public static void main(String[] args) {
        File testFile = new File("exampleBAM.bam");
        SAMFileReader fileReader = new SAMFileReader(testFile);

        SamLocusIterator samLocusIterator = new SamLocusIterator(fileReader);
        Iterator<LocusInfo> iterator = samLocusIterator.iterator();

        while (iterator.hasNext()) {
            LocusInfo locus = iterator.next();
            System.out.println("pos: " + locus.getSequenceName() + ":" + locus.getPosition());
        }
    }
}

The Scala version

import net.sf.samtools.SAMFileReader
import net.sf.picard.util.SamLocusIterator
import java.io.File

object MinimalIteratorExample extends App {

    val testFile = new File("exampleBAM.bam")
    val fileReader = new SAMFileReader(testFile)

    val samLocusIterator: SamLocusIterator = new SamLocusIterator(fileReader)
    val iterator = samLocusIterator.iterator()

    while (iterator.hasNext()) {
        val locus = iterator.next()
        println("pos: " + locus.getSequenceName() + ":" + locus.getPosition())
    }
}

Am I'm using the iterator in the wrong way, or is this some bug in the Picard code? I've been looking debugging and looking into the Picard classes, but this far I've not been able to make sense of it. Any help figuring this out would be much appreciated.

[Pierre Lindenbaum]

I'm afraid 'hasNext()' calls 'next()' ...

[Johan]

Well, I see. That's not great, now, is it? I'll post it as an issue to the picard/samtools-dev list.

[Pierre Lindenbaum]

Thank you. I just had a glance at github: the only place I saw this API used was a test in GATK (a simple loop to get the time for scanning a bam.... )

[Pierre Lindenbaum]

did you post your question on the samtools mailing list ? I cannot see the post. http://sourceforge.net/mailarchive/forum.php?forum_name=samtools-devel

[Johan]

I did, but only realized this morning that I had to be a member on the list to post. I posted it now though.

[Sangwoo Kim]

I think it is not the case. I just tried the similar code above replacing the "hasNext()" condition to true. But got the same result.

[Pierre Lindenbaum]

I've got the same problem here. Thanks I've never noticed this problem. you'd better ask the picard/samtools-dev mailing list.

[Fabian Bull]

In case you did not know already: You can import scala.collection.JavaConversions.asScalaIterator and you the SamLocusIterator as every iterator. e.g.: for( entry <- iterator) { do something meaningful with entry } Saves some LOC and looks much cleaner.

Answer 1

[Sangwoo Kim]

I just found that this occurs in regions where no reads have been aligned. Because these regions are not of our interest, we should skip them by setting setEmitUncoveredLoci()

Try this.

SamLocusIterator samLocusIterator = new SamLocusIterator(fileReader);
samLocusIterator.setEmitUncoveredLoci(false);
Iterator<LocusInfo> iterator = samLocusIterator.iterator();

[Johan]

You're right, this only happens in uncovered regions. :) Still somewhat unexpected, though.

[Pierre Lindenbaum]

Nice to know ! Thank you.

Answer 2

[agolicz]

This bug is fixed now, as of picard tools v1.88. In case somebody also needs non-covered regions.