Fixes for VariantEval for genotyping mode

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1659 348d0f76-0448-11de-a6fe-93d51630548a

java/src/org/broadinstitute/sting/gatk/refdata/RodGeliText.java

@@ -68,7 +68,12 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
             refBase = parts[2].charAt(0);
             depth = Integer.valueOf(parts[3]);
             maxMappingQuality = Integer.valueOf(parts[4]);
-            bestGenotype = parts[5];
+
+            System.out.printf("%s%n", parts[5]);
+            char[] x = parts[5].toUpperCase().toCharArray();
+            Arrays.sort(x);
+            bestGenotype = new String(x);
+
             lodBtr = Double.valueOf(parts[6]);
             lodBtnb = Double.valueOf(parts[7]);
 
@@ -349,7 +354,7 @@ public class RodGeliText extends BasicReferenceOrderedDatum implements Variation
      */
     @Override
     public Genotype getGenotype(DiploidGenotype x) {
-        if (x.toString() != this.getAltBasesFWD()) throw new IllegalStateException("We don't contain that genotype!"); 
+        if (x.toString() != this.getAltBasesFWD()) throw new IllegalStateException("We don't contain genotype " + x); 
         return new BasicGenotype(getLocation(), x.toString(), refBase, lodBtnb);
     }
 

java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/CallableBasesAnalysis.java

@@ -69,14 +69,14 @@ public class CallableBasesAnalysis extends BasicVariantAnalysis implements Genot
         // For every threshold, updated discoverable and callable
         for (int i = 0; i < thresholds.length; i++) {
             double threshold = thresholds[i];
-            DiploidGenotype g = DiploidGenotype.createHomGenotype(ref);
-            Genotype genotype = ((VariantBackedByGenotype) eval).getGenotype(g);
+            Genotype genotype = ((VariantBackedByGenotype)eval).getGenotypes().get(0);
+
             // update discoverable
-            if (eval.isSNP() && eval.getNegLog10PError() >= threshold)
+            if ( eval.isSNP() && eval.getNegLog10PError() >= threshold)
                 discoverable_bases[i]++;
-            if (!eval.isSNP() && genotype.getNegLog10PError() >= threshold)
+            if ( eval.isReference() && genotype.getNegLog10PError() >= threshold)
                 discoverable_bases[i]++;
-            if (genotype.getNegLog10PError() >= threshold)
+            if ( genotype.getNegLog10PError() >= threshold)
                 genotypable_bases[i]++;
         }
 

java/src/org/broadinstitute/sting/playground/gatk/walkers/varianteval/VariantCounter.java

@@ -5,6 +5,7 @@ import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
 import org.broadinstitute.sting.utils.genotype.DiploidGenotype;
 import org.broadinstitute.sting.utils.genotype.VariantBackedByGenotype;
 import org.broadinstitute.sting.utils.genotype.Variation;
+import org.broadinstitute.sting.utils.genotype.Genotype;
 
 import java.util.ArrayList;
 import java.util.List;
@@ -31,8 +32,12 @@ public class VariantCounter extends BasicVariantAnalysis implements GenotypeAnal
     public String update(Variation eval, RefMetaDataTracker tracker, char ref, AlignmentContext context) {
         nSNPs += eval == null ? 0 : 1;
 
-        if ( this.getMaster().evalContainsGenotypes && eval != null && eval.isSNP() && ((VariantBackedByGenotype)eval).getGenotype( DiploidGenotype.valueOf(eval.getAlternateBase())).isHet() )
-            nHets++;
+        if ( this.getMaster().evalContainsGenotypes && eval != null ) {
+            List<Genotype> genotypes = ((VariantBackedByGenotype)eval).getGenotypes();
+            if ( eval.isSNP() && eval.isBiallelic() && genotypes.get(0).isHet() ) {
+                nHets++;
+            }
+        }
 
         return null;
     }

python/MergeBAMBatch.py

@@ -37,6 +37,9 @@ if __name__ == "__main__":
     parser.add_option("-n", "--naIDPops", dest="NAIDS2POP",
                         type="string", default=None,
                         help="Path to file contains NAID POP names.  If provided, input files will be merged by population")
+    parser.add_option("-p", "--pop", dest="onlyPop",
+                        type="string", default=None,
+                        help="If provided, only this population will be processed")
                         
     (OPTIONS, args) = parser.parse_args()
     if len(args) != 1:
@@ -58,6 +61,9 @@ if __name__ == "__main__":
             allSources = reduce( operator.__add__, map( glob.glob, s[1:] ), [] )
             print 'Merging info:'
             for spec in splitSourcesByPopulation(allSources, merged_filename_base, NAID2Pop):
+                if OPTIONS.onlyPop <> None and spec.group() <> OPTIONS.onlyPop:
+                    continue
+            
                 spec.setPath(directory)
                 spec.pprint()
                 

python/MergeBAMsUtils.py

@@ -68,7 +68,8 @@ class MergeFilesSpec:
         print ' Population:         ', self.group()
         print ' Merged filename:    ', self.getMergedBAM()
         print ' N sources:          ', len(self.sources())
-        print ' Sources:            ', self.sources()
+        for source in sorted(self.sources()):
+            print ' Source:             ', source
         print ' Sizes:              ', self.sourceSizes(humanReadable=True)
         print ' Est. merged size:   ', greek(reduce(operator.__add__, self.sourceSizes(), 0))