(cross posted on Biostar: http://www.biostars.org/post/show/52698 )
I want to get the coverage-depth at a given position : (chr1-12450)
using picard:SamLocusIterator
I can clearly see my reads overlapping that position using samtools tview.
12441 12451 12461 12471 12481 12491 12501
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
CCAGTGGATTGGCCTAGGTGGGATCTCTGAGCTCAACAAGCCCTCTCTGGGTGGTAGGTGCAGAGACGGGAGGGGCAGAG
ccagtggattggcctaggtgggatctctgagctcaacaagccctctctgggtggtaggtgcagagacgggaggggcagag
CCAGTGGA
tggcctaggtgggatctctgagctcaacaagccctctctgggtggtaggtgcagagacgggaggggcagag
ccagtgga
GCCTAGGTGGGATCTCTGAGCTCAACAAGCCCTCTCTGGGTGGTAGGTGCAGAGACGGGAGGGGCAGAG
CCAGTGGAT
CTAGGAGGGATCTCTGAGCTCAACAAGCCCTCTATGGGGGGGAGGTGCAGAGACGGGAGGGGCAGAG
CCAGTGGATTG
taggtgggatctctgagctcaacaagccctctctgggtggtaggtgcagagacgggaggggcagag
CCAGTGGATTGGCCTAG
tgggatctctgagctcaacaagccctctctgggtggtaggtgcagagacgggaggggcagag
CCAGTGGATTGGCCTAGG
gatctctgagctcaacaagccctctctgggtggtaggtgcagagacgggaggggcagag
CCAGTGGATTGGCCTAGGTG
atctctgagctcaacaagccctctctgggtggtaggtgcagagacgggaggggcagag
CCAGTGGATTGGCCTAGGTGG CTGAGCTCAACAAGCCCTCTCTGGGTGGTAGGTGCAGAGACGGGAGGGGCAGAG
(...)
However, using the following program:
| import net.sf.picard.reference.*;
import net.sf.picard.util.*;
import net.sf.samtools.*;
import java.io.*;
import java.util.*;
public class Test
{
public static void main(String args[])
{
SAMFileReader samReader=new SAMFileReader(new File(args[0]));
String chromId="chr1";
int chromStart=12450;
int chromEnd=chromStart;
Interval interval=new Interval(chromId,chromStart,chromEnd);
IntervalList iL=new IntervalList(samReader.getFileHeader());
iL.add(interval);
SamLocusIterator sli=new SamLocusIterator(samReader,iL,true);
for(Iterator<SamLocusIterator.LocusInfo> iter=sli.iterator();
iter.hasNext();
)
{
SamLocusIterator.LocusInfo locusInfo=iter.next();
System.out.println(
"POS="+locusInfo.getPosition()+
" depth:"+locusInfo.getRecordAndPositions().size()
);
}
sli.close();
samReader.close();
}
}
|
compile and run:
| javac-cp/sam-1.77.jar:picard-1.77.jar:. Test.java
java-cp/sam-1.77.jar:picard-1.77.jar:. Test my.bam|
produces the following output (depth==0)
| POS=12450 depth:0|
What's wrong with my program ?
Another problem, my program also raises an exception when sli.close() is
invoked:
| Exception in thread"main" java.lang.NullPointerException
at net.sf.picard.util.SamRecordIntervalIterator.close(SamRecordIntervalIterator.java:162)
(...)
at net.sf.picard.util.SamLocusIterator.close(SamLocusIterator.java:219)
atTest.main(Test.java:32)|
why ?
Many thanks for your help,
Pierre
|