diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/DeNovoSNPWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/DeNovoSNPWalker.java index e8e4a0a79..961014981 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/DeNovoSNPWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/DeNovoSNPWalker.java @@ -5,10 +5,19 @@ import org.broadinstitute.sting.gatk.contexts.ReferenceContext; import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; import org.broadinstitute.sting.gatk.refdata.VariationRod; import org.broadinstitute.sting.gatk.walkers.*; +import org.broadinstitute.sting.gatk.walkers.genotyper.SingleSampleGenotyper; +import org.broadinstitute.sting.gatk.walkers.genotyper.SSGenotypeCall; +import org.broadinstitute.sting.gatk.filters.ZeroMappingQualityReadFilter; import org.broadinstitute.sting.utils.StingException; import org.broadinstitute.sting.utils.genotype.VariantBackedByGenotype; import org.broadinstitute.sting.utils.genotype.Variation; +import java.util.Set; +import java.util.List; +import java.util.ArrayList; + +import net.sf.samtools.SAMRecord; + /** * Created by IntelliJ IDEA. * User: andrewk @@ -17,45 +26,82 @@ import org.broadinstitute.sting.utils.genotype.Variation; * To change this template use File | Settings | File Templates. */ @By(DataSource.REFERENCE) -@Requires(value={DataSource.REFERENCE},referenceMetaData={@RMD(name="child",type= VariationRod.class), @RMD(name="parent1",type= VariationRod.class), @RMD(name="parent2",type= VariationRod.class)}) -//@Allows(value={DataSource.REFERENCE},referenceMetaData = {@RMD(name="eval",type=AllelicVariant.class), @RMD(name="dbsnp",type=AllelicVariant.class),@RMD(name="hapmap-chip",type=AllelicVariant.class), @RMD(name="interval",type=IntervalRod.class)}) -//@Allows(value={DataSource.REFERENCE},referenceMetaData = {@RMD(name="dbsnp",type=AllelicVariant.class)}) +@Requires(value={DataSource.REFERENCE, DataSource.REFERENCE_BASES, DataSource.READS},referenceMetaData={@RMD(name="child",type= VariationRod.class)}) +@Allows({DataSource.READS, DataSource.REFERENCE}) +@ReadFilters(ZeroMappingQualityReadFilter.class) +//, @RMD(name="parent1",type= VariationRod.class), @RMD(name="parent2",type= VariationRod.class)}) + public class DeNovoSNPWalker extends RefWalker{ + SingleSampleGenotyper SSG; + private List> readGroupSets; + + public void initialize() { + SSG = new SingleSampleGenotyper(); + SSG.initialize(); + + readGroupSets = getToolkit().getMergedReadGroupsByReaders(); + } + public String map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { Variation child = (Variation)tracker.lookup("child", null); - Variation parent1 = (Variation)tracker.lookup("parent1", null); - Variation parent2 = (Variation)tracker.lookup("parent2", null); Variation dbsnp = (Variation)tracker.lookup("dbSNP", null); - if (child != null && parent1 != null && parent2 != null) { - if (!(parent1 instanceof VariantBackedByGenotype) || !(parent2 instanceof VariantBackedByGenotype)) - throw new StingException("Both parents ROD tracks must be backed by genotype data. Ensure that your parent rod(s) contain genotyping information"); - if (child.isSNP() && - child.getNegLog10PError() > 5 && // BTR > 5 - parent1.isReference() && - ((VariantBackedByGenotype)parent1).getCalledGenotype().getNegLog10PError() > 5 && - parent2.isReference() && - ((VariantBackedByGenotype)parent2).getCalledGenotype().getNegLog10PError() > 5 + if (child != null) { + if (child.isSNP() && child.getNegLog10PError() > 5) { // BTR > 5 + + List reads = context.getReads(); + List offsets = context.getOffsets(); + + List parent1_reads = new ArrayList(); + List parent2_reads = new ArrayList(); + List parent1_offsets = new ArrayList(); + List parent2_offsets = new ArrayList(); + + assert( reads.size() == offsets.size() ); // should be same number or we're in trouble + int num_reads = reads.size(); + for (int i=0; i 5 && + !parent2.isVariant(parent2.getReference()) && + parent2.getNegLog10PError() > 5 ) { - double sumConfidences = 0.5 * (0.5 * child.getNegLog10PError() + - Math.min(((VariantBackedByGenotype)parent1).getCalledGenotype().getNegLog10PError(), - ((VariantBackedByGenotype)parent2).getCalledGenotype().getNegLog10PError())); + double sumConfidences = 0.5 * (0.5 * child.getNegLog10PError() + + Math.min(parent1.getNegLog10PError(), parent2.getNegLog10PError())); - out.format("%s\t", child.getLocation().getContig()); - out.format("%s\t", child.getLocation().getStart()); - out.format("%.4f\t", sumConfidences); - out.format("%.4f\t", child.getNegLog10PError()); - out.format("%.4f\t", ((VariantBackedByGenotype)parent1).getCalledGenotype().getNegLog10PError()); - out.format("%.4f\t", ((VariantBackedByGenotype)parent2).getCalledGenotype().getNegLog10PError()); - out.format("%s\t", dbsnp != null); + out.format("%s\t", child.getLocation().getContig()); + out.format("%s\t", child.getLocation().getStart()); + out.format("%.4f\t", sumConfidences); + out.format("%.4f\t", child.getNegLog10PError()); + out.format("%.4f\t", parent1.getNegLog10PError()); + out.format("%.4f\t", parent2.getNegLog10PError()); + out.format("%s\t", dbsnp != null); - out.format ("%s\t", child.toString()); - out.format ("%s\t", parent1.toString()); - out.format ("%s", parent2.toString()); - if (dbsnp != null) - out.format ("\tDBSNP\t:%s", dbsnp.toString()); - out.println(); + out.format ("%s\t", child.toString()); + out.format ("%s\t", parent1.toString()); + out.format ("%s", parent2.toString()); + if (dbsnp != null) + out.format ("\tDBSNP\t:%s", dbsnp.toString()); + out.println(); + } } } @@ -67,5 +113,5 @@ public class DeNovoSNPWalker extends RefWalker{ return 1; } - + public void onTraversalDone(Integer result) {} // Don't print the reduce result }