From 91030e9afa344b44f6551551f1f5687b48d4e51c Mon Sep 17 00:00:00 2001
From: Chris Hartl <chartl@broadinstitute.org>
Date: Fri, 18 Jan 2013 09:49:48 -0500
Subject: [PATCH] Bugfix: records that get paired up during the resolution of
 multiple-records-per-site were not going into genotype-level filtering.
 Caught via testing.

Testing for moltenized output, and for genotype-level filtering. This tool is now fully functional. There are three todo items:

1) Docs
2) An additional output table that gives concordance proportions normalized by records in both eval and comp (not just total in eval or total in comp)
3) Code cleanup for table creation (putting a table together the way I do takes -way- too many lines of code)
---
 .../variantutils/GenotypeConcordance.java     | 15 +++----
 .../GenotypeConcordanceIntegrationTest.java   | 44 +++++++++++++++++++
 2 files changed, 50 insertions(+), 9 deletions(-)

diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordance.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordance.java
index 2acff956c..e8965dfc8 100644
--- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordance.java
+++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordance.java
@@ -94,9 +94,8 @@ public class GenotypeConcordance extends RodWalker<List<Pair<VariantContext,Vari
     private List<VariantContextUtils.JexlVCMatchExp> evalJexls = null;
     private List<VariantContextUtils.JexlVCMatchExp> compJexls = null;
 
-    // todo -- table with "proportion of overlapping sites" (not just eval/comp margins)
-    // todo -- moltenize
-
+    // todo -- table with "proportion of overlapping sites" (not just eval/comp margins) [e.g. drop no-calls]
+    //  (this will break all the integration tests of course, due to new formatting)
 
     public void initialize() {
         evalJexls = initializeJexl(genotypeFilterExpressionsEval);
@@ -201,7 +200,7 @@ public class GenotypeConcordance extends RodWalker<List<Pair<VariantContext,Vari
             }
             if ( pairedComp != null ) {
                 compList.remove(pairedComp);
-                resolvedPairs.add(new Pair<VariantContext, VariantContext>(eval,pairedComp));
+                resolvedPairs.add(new Pair<VariantContext, VariantContext>(filterGenotypes(eval,ignoreFilters,evalJexls),filterGenotypes(pairedComp,ignoreFilters,compJexls)));
                 pairedEval.add(eval);
                 if ( compList.size() < 1 )
                     break;
@@ -209,11 +208,11 @@ public class GenotypeConcordance extends RodWalker<List<Pair<VariantContext,Vari
         }
         evalList.removeAll(pairedEval);
         for ( VariantContext unpairedEval : evalList ) {
-            resolvedPairs.add(new Pair<VariantContext, VariantContext>(unpairedEval,createEmptyContext(unpairedEval,compSamples)));
+            resolvedPairs.add(new Pair<VariantContext, VariantContext>(filterGenotypes(unpairedEval,ignoreFilters,evalJexls),createEmptyContext(unpairedEval,compSamples)));
         }
 
         for ( VariantContext unpairedComp : compList ) {
-            resolvedPairs.add(new Pair<VariantContext, VariantContext>(createEmptyContext(unpairedComp,evalSamples),unpairedComp));
+            resolvedPairs.add(new Pair<VariantContext, VariantContext>(createEmptyContext(unpairedComp,evalSamples),filterGenotypes(unpairedComp,ignoreFilters,compJexls)));
         }
 
         return resolvedPairs;
@@ -233,6 +232,7 @@ public class GenotypeConcordance extends RodWalker<List<Pair<VariantContext,Vari
     }
 
     public void onTraversalDone(ConcordanceMetrics metrics) {
+        // todo -- this is over 200 lines of code just to format the output and could use some serious cleanup
         GATKReport report = new GATKReport();
         GATKReportTable concordanceCounts = new GATKReportTable("GenotypeConcordance_Counts","Per-sample concordance tables: comparison counts",2+GenotypeType.values().length*GenotypeType.values().length);
         GATKReportTable concordanceEvalProportions = new GATKReportTable("GenotypeConcordance_EvalProportions", "Per-sample concordance tables: proportions of genotypes called in eval",2+GenotypeType.values().length*GenotypeType.values().length);
@@ -460,9 +460,6 @@ public class GenotypeConcordance extends RodWalker<List<Pair<VariantContext,Vari
     }
 
     public VariantContext filterGenotypes(VariantContext context, boolean ignoreSiteFilter, List<VariantContextUtils.JexlVCMatchExp> exps) {
-        // placeholder method for genotype-level filtering. However if the site itself is filtered,
-        // and such filters are not ignored, the genotype-level data should be altered to reflect this
-
         if ( ! context.isFiltered() || ignoreSiteFilter ) {
             List<Genotype> filteredGenotypes = new ArrayList<Genotype>(context.getNSamples());
             for ( Genotype g : context.getGenotypes() ) {
diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordanceIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordanceIntegrationTest.java
index 113f098e3..117032ac9 100644
--- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordanceIntegrationTest.java
+++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/variantutils/GenotypeConcordanceIntegrationTest.java
@@ -61,6 +61,17 @@ public class GenotypeConcordanceIntegrationTest extends WalkerTest {
         executeTest("test non-overlapping samples", spec);
     }
 
+    @Test
+    public void testNonoverlappingSamplesMoltenized() {
+        WalkerTestSpec spec = new WalkerTestSpec(
+                baseTestString("GenotypeConcordanceNonOverlapTest_Eval.vcf", "GenotypeConcordanceNonOverlapTest_Comp.vcf"),
+                0,
+                Arrays.asList("")
+        );
+
+        executeTest("Test moltenized output",spec);
+    }
+
     @Test
     public void testMultipleRecordsPerSite() {
         WalkerTestSpec spec = new WalkerTestSpec(
@@ -71,4 +82,37 @@ public class GenotypeConcordanceIntegrationTest extends WalkerTest {
 
         executeTest("test multiple records per site",spec);
     }
+
+    @Test
+    public void testGQFilteringEval() {
+        WalkerTestSpec spec = new WalkerTestSpec(
+                baseTestString("genotypeConcordanceFilterTest.vcf","genotypeConcordanceFilterTest.vcf") + " -gfe 'GQ<30'",
+                0,
+                Arrays.asList("b7b495ccfa6d50a6be3e095d3f6d3c52")
+        );
+
+        executeTest("Test filtering on the EVAL rod",spec);
+    }
+
+    @Test
+    public void testFloatFilteringComp() {
+        WalkerTestSpec spec = new WalkerTestSpec(
+                baseTestString("genotypeConcordanceFilterTest.vcf","genotypeConcordanceFilterTest.vcf") + " -gfc 'LX<0.50'",
+                0,
+                Arrays.asList("6406b16cde7960b8943edf594303afd6")
+        );
+
+        executeTest("Test filtering on the COMP rod", spec);
+    }
+
+    @Test
+    public void testCombinedFilters() {
+        WalkerTestSpec spec = new WalkerTestSpec(
+                baseTestString("genotypeConcordanceFilterTest.vcf","genotypeConcordanceFilterTest.vcf") + " -gfc 'LX<0.52' -gfe 'DP<5' -gfe 'GQ<37'",
+                0,
+                Arrays.asList("26ffd06215b6177acce0ea9f35d73d31")
+        );
+
+        executeTest("Test filtering on both rods",spec);
+    }
 }