diff --git a/java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/TiTvVariantEvaluator.java b/java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/TiTvVariantEvaluator.java new file mode 100755 index 000000000..f72fc40eb --- /dev/null +++ b/java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/TiTvVariantEvaluator.java @@ -0,0 +1,65 @@ +package org.broadinstitute.sting.oneoffprojects.variantcontext.varianteval2; + +import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker; +import org.broadinstitute.sting.gatk.contexts.AlignmentContext; +import org.broadinstitute.sting.gatk.contexts.ReferenceContext; +import org.broadinstitute.sting.oneoffprojects.variantcontext.VariantContext; +import org.broadinstitute.sting.oneoffprojects.variantcontext.Genotype; +import org.broadinstitute.sting.utils.StingException; + +import java.util.List; +import java.util.Arrays; + +public class TiTvVariantEvaluator extends VariantEvaluator { + long nTi = 0, nTv = 0; + long nTiInStd = 0, nTvInStd = 0; + + private double ratio(long num, long denom) { + return ((double)num) / (Math.max(denom, 1)); + } + + public String getName() { + return "titv"; + } + + public int getComparisonOrder() { + return 2; // we only need to see each eval track + } + + public void updateTiTv(VariantContext vc, boolean updateStandard) { + if ( vc != null && vc.isSNP() && vc.isBiallelic() ) { + if ( vc.isTransition() ) { + if ( updateStandard ) nTiInStd++; else nTi++; + } else { + if ( updateStandard ) nTvInStd++; else nTv++; + } + } + } + + public void update2(VariantContext vc1, VariantContext vc2, RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) { + if ( vc1 != null ) updateTiTv(vc1, false); + if ( vc2 != null ) updateTiTv(vc2, true); + } + + public String toString() { + return getName() + ": " + summaryLine(); + } + + private String summaryLine() { + return String.format("%d %d %.2f %d %d %.2f", + nTi, nTv, ratio(nTi, nTv), + nTiInStd, nTvInStd, ratio(nTiInStd, nTvInStd)); + } + + private static List HEADER = + Arrays.asList("nTi", "nTv", "TiTvRatio", "nTiStandard", "nTvStandard", "TiTvRatioStandard"); + + // making it a table + public List getTableHeader() { + return HEADER; + } + + public List> getTableRows() { + return Arrays.asList(Arrays.asList(summaryLine().split(" "))); + } +} \ No newline at end of file diff --git a/java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/VariantEval2Walker.java b/java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/VariantEval2Walker.java index 41ef3df4a..ed83fbac5 100755 --- a/java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/VariantEval2Walker.java +++ b/java/src/org/broadinstitute/sting/oneoffprojects/variantcontext/varianteval2/VariantEval2Walker.java @@ -69,13 +69,15 @@ public class VariantEval2Walker extends RodWalker { for ( VariantContextUtils.MatchExp e : selectExps ) { addNewContext(d.getName(), d.getName() + "." + e.name, e); } - addNewContext(d.getName(), d.getName(), null); + addNewContext(d.getName(), d.getName() + ".all", null); } else if ( d.getName().startsWith("dbsnp") || d.getName().startsWith("hapmap") || d.getName().startsWith("comp") ) { compNames.add(d.getName()); } else { logger.info("Not evaluating ROD binding " + d.getName()); } } + + determineContextNamePartSizes(); } private void determineAllEvalations() { @@ -144,7 +146,7 @@ public class VariantEval2Walker extends RodWalker { // now call the single or paired update function switch ( evaluation.getComparisonOrder() ) { case 1: - if ( vc != null ) evaluation.update1(vc, tracker, ref, context); + if ( evalWantsVC && vc != null ) evaluation.update1(vc, tracker, ref, context); break; case 2: for ( VariantContext comp : comps.values() ) { @@ -249,16 +251,59 @@ public class VariantEval2Walker extends RodWalker { return null; } - public > List sorted(Collection c ) { + private > List sorted(Collection c ) { List l = new ArrayList(c); Collections.sort(l); return l; } + private final static String CONTEXT_HEADER = "track.subset.novelty.filter"; + private final static int N_CONTEXT_NAME_PARTS = CONTEXT_HEADER.split("\\.").length; + private static int[] nameSizes = new int[N_CONTEXT_NAME_PARTS]; + static { + int i = 0; + for ( String elt : CONTEXT_HEADER.split("\\.") ) + nameSizes[i++] = elt.length(); + } + + private void determineContextNamePartSizes() { + for ( String contextName : sorted(contexts.keySet()) ) { + EvaluationContext group = contexts.get(contextName); + for ( String evalSubgroupName : sorted(group.keySet()) ) { + String keyWord = contextName + "." + evalSubgroupName; + String[] parts = keyWord.split("\\."); + if ( parts.length != N_CONTEXT_NAME_PARTS ) { + throw new StingException("Unexpected number of eval name parts " + keyWord + " length = " + parts.length + ", expected " + N_CONTEXT_NAME_PARTS); + } else { + for ( int i = 0; i < parts.length; i++ ) + nameSizes[i] = Math.max(nameSizes[i], parts[i].length()); + } + } + } + } + + private String formatKeyword(String keyWord) { + //System.out.printf("keyword %s%n", keyWord); + + StringBuilder s = new StringBuilder(); + int i = 0; + for ( String part : keyWord.split("\\.") ) { + //System.out.printf("part %s %d%n", part, nameSizes[i]); + s.append(String.format("%"+nameSizes[i]+"s ", part)); + i++; + } + + return s.toString(); + } + public void onTraversalDone(Integer result) { + // todo -- this really needs to be pretty printed; use some kind of table organization + // todo -- so that we can load up all of the data in one place, analyze the widths of the columns + // todo -- and pretty print it for ( String evalName : variantEvaluationNames ) { boolean first = true; out.printf("%n%n"); + // todo -- show that comp is dbsnp, etc. is columns for ( String contextName : sorted(contexts.keySet()) ) { EvaluationContext group = contexts.get(contextName); @@ -267,13 +312,14 @@ public class VariantEval2Walker extends RodWalker { Set evalSet = group.get(evalSubgroupName); VariantEvaluator eval = getEvalByName(evalName, evalSet); String keyWord = contextName + "." + evalSubgroupName; + if ( first ) { - out.printf("%20s\t%40s\t%s%n", evalName, "context", Utils.join("\t\t", eval.getTableHeader())); + out.printf("%20s %s %s%n", evalName, formatKeyword(CONTEXT_HEADER), Utils.join("\t\t", eval.getTableHeader())); first = false; } for ( List row : eval.getTableRows() ) - out.printf("%20s\t%40s\t%s%n", evalName, keyWord, Utils.join("\t\t", row)); + out.printf("%20s %s %s%n", evalName, formatKeyword(keyWord), Utils.join("\t\t", row)); } } }