diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java index f55d0a4a3..10cffa09c 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java @@ -6,6 +6,7 @@ import org.broadinstitute.sting.gatk.refdata.rodDbSNP; import org.broadinstitute.sting.gatk.walkers.LocusWalker; import org.broadinstitute.sting.playground.utils.AlleleFrequencyEstimate; import org.broadinstitute.sting.utils.cmdLine.Argument; +import org.broadinstitute.sting.utils.*; import net.sf.samtools.SAMRecord; import java.util.List; @@ -15,7 +16,6 @@ import java.io.PrintStream; public class AlleleFrequencyWalker extends LocusWalker { - @Argument public int N; @Argument public int DOWNSAMPLE; @Argument public String GFF_OUTPUT_FILE; @@ -25,7 +25,6 @@ public class AlleleFrequencyWalker extends LocusWalker rodData, char ref, LocusContext context) { - // Convert context data into bases and 4-base quals // Convert context data into bases and 4-base quals String bases = getBases(context); double quals[][] = getOneBaseQuals(context); @@ -354,16 +353,67 @@ public class AlleleFrequencyWalker extends LocusWalker= 5) || (alleleFreq.lodVsRef <= -5)) { this.output.print(alleleFreq.asGFFString()); } + + if (inside_confident_ref_interval && (alleleFreq.lodVsRef > -5.0)) + { + // No longer hom-ref, so output a ref line. + String[] tokens; + tokens = confident_ref_interval_start.split(":"); + + String contig = tokens[0]; + int start = Integer.parseInt(tokens[1]); + + tokens = alleleFreq.location.split(":"); + int end = Integer.parseInt(tokens[1])-1; + + double lod = confident_ref_interval_LOD_sum / confident_ref_interval_length; + + output.format("%s\tCALLER\tREFERENCE\t%d\t%d\t%f\t.\t.\tLENGTH %d\n", + contig, + start, + end, + lod, + (int)(confident_ref_interval_length)); + + inside_confident_ref_interval = false; + } + else if (inside_confident_ref_interval && (alleleFreq.lodVsRef <= -5.0)) + { + // Still hom-ref so increment the counters. + confident_ref_interval_LOD_sum += alleleFreq.lodVsRef; + confident_ref_interval_length += 1; + } + else if ((!inside_confident_ref_interval) && (alleleFreq.lodVsRef > -5.0)) + { + // do nothing. + } + else if ((!inside_confident_ref_interval) && (alleleFreq.lodVsRef <= -5.0)) + { + // We moved into a hom-ref region so start a new interval. + confident_ref_interval_start = alleleFreq.location; + confident_ref_interval_LOD_sum = alleleFreq.lodVsRef; + confident_ref_interval_length = 1; + inside_confident_ref_interval = true; + } + + if (alleleFreq.lodVsRef >= 5) { this.output.print(alleleFreq.asGFFString()); } return ""; }