zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Now outputs four-base SAM record (read name prefixed with KIR) and bustard SAM record (prefixed with BUS) for easy debugging. 

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@285 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
kiran 2009-04-03 15:48:51 +00:00
parent 15151ac125
commit 5019971290
1 changed files with 31 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
java/src/org/broadinstitute/sting/playground/fourbasecaller/FourBaseRecaller.java
View File

@ -38,14 +38,17 @@ public class FourBaseRecaller extends CommandLineProgram {
protected int execute() { protected int execute() {
boolean isPaired = (END > 0); boolean isPaired = (END > 0);
BustardFileParser bfp = new BustardFileParser(DIR, LANE, isPaired, "BS"); BustardFileParser bfp;
BustardReadData bread = bfp.next(); BustardReadData bread;
bfp = new BustardFileParser(DIR, LANE, isPaired, "FB");
bread = bfp.next();
int cycle_offset = (END <= 1) ? 0 : bread.getIntensities().length/2; int cycle_offset = (END <= 1) ? 0 : bread.getIntensities().length/2;
BasecallingReadModel p = new BasecallingReadModel(bread.getFirstReadSequence().length()); BasecallingReadModel p = new BasecallingReadModel(bread.getFirstReadSequence().length());
int queryid; int queryid;
// learn parameters // learn initial parameters
queryid = 0; queryid = 0;
do { do {
String bases = (END <= 1) ? bread.getFirstReadSequence() : bread.getSecondReadSequence(); String bases = (END <= 1) ? bread.getFirstReadSequence() : bread.getSecondReadSequence();
@ -68,6 +71,9 @@ public class FourBaseRecaller extends CommandLineProgram {
SAMFileHeader sfh = new SAMFileHeader(); SAMFileHeader sfh = new SAMFileHeader();
SAMFileWriter sfw = new SAMFileWriterFactory().makeSAMOrBAMWriter(sfh, false, OUT); SAMFileWriter sfw = new SAMFileWriterFactory().makeSAMOrBAMWriter(sfh, false, OUT);
bfp = new BustardFileParser(DIR, LANE, isPaired, "FB");
bread = bfp.next();
queryid = 0; queryid = 0;
do { do {
String bases = (END <= 1) ? bread.getFirstReadSequence() : bread.getSecondReadSequence(); String bases = (END <= 1) ? bread.getFirstReadSequence() : bread.getSecondReadSequence();
@ -79,19 +85,21 @@ public class FourBaseRecaller extends CommandLineProgram {
byte[] nextbestqual = new byte[bases.length()]; byte[] nextbestqual = new byte[bases.length()];
for (int cycle = 0; cycle < bases.length(); cycle++) { for (int cycle = 0; cycle < bases.length(); cycle++) {
char basePrev = (cycle == 0) ? '*' : bases.charAt(cycle - 1); //char basePrev = (cycle == 0) ? '*' : bases.charAt(cycle - 1);
byte qualPrev = (cycle == 0) ? 0 : quals[cycle - 1]; //byte qualPrev = (cycle == 0) ? 0 : quals[cycle - 1];
char basePrev = (cycle == 0) ? '*' : (char) asciiseq[cycle - 1];
byte qualPrev = (cycle == 0) ? 0 : bestqual[cycle - 1];
double[] fourintensity = intensities[cycle + cycle_offset]; double[] fourintensity = intensities[cycle + cycle_offset];
FourProb fp = p.computeProbabilities(cycle, basePrev, qualPrev, fourintensity); FourProb fp = p.computeProbabilities(cycle, basePrev, qualPrev, fourintensity);
asciiseq[cycle] = (byte) fp.baseAtRank(0); asciiseq[cycle] = (byte) fp.baseAtRank(0);
bestqual[cycle] = fp.qualAtRank(0); bestqual[cycle] = fp.qualAtRank(0);
nextbestqual[cycle] = QualityUtils.qualAndProbToCompressedQuality(fp.indexAtRank(1), fp.probAtRank(1)); nextbestqual[cycle] = QualityUtils.baseAndProbToCompressedQuality(fp.indexAtRank(1), fp.probAtRank(1));
} }
SAMRecord sr = new SAMRecord(sfh); SAMRecord sr = new SAMRecord(sfh);
sr.setReadName(bread.getReadName()); sr.setReadName("KIR_" + bread.getReadName());
sr.setReadUmappedFlag(true); sr.setReadUmappedFlag(true);
sr.setReadBases(asciiseq); sr.setReadBases(asciiseq);
sr.setBaseQualities(bestqual); sr.setBaseQualities(bestqual);
@ -105,10 +113,23 @@ public class FourBaseRecaller extends CommandLineProgram {
} }
sfw.addAlignment(sr); sfw.addAlignment(sr);
SAMRecord sr2 = new SAMRecord(sfh);
sr2.setReadName("BUS_" + bread.getReadName());
sr2.setReadUmappedFlag(true);
sr2.setReadString(bases);
sr2.setBaseQualities(quals);
sr2.setReadFailsVendorQualityCheckFlag(!bread.isPf());
sr2.setReadPairedFlag(isPaired);
if (isPaired) {
sr2.setMateUnmappedFlag(true);
sr2.setFirstOfPairFlag(END <= 1);
sr2.setSecondOfPairFlag(END > 1);
}
sfw.addAlignment(sr2);
/* /*
System.out.println(sr.format()); System.out.println(sr.format());
System.out.println(sr.getReadString()); System.out.println(sr2.format());
System.out.println(bases);
System.out.println("\n"); System.out.println("\n");
*/ */