Merge branch 'master' of ssh://nickel.broadinstitute.org/humgen/gsa-scr1/gsa-engineering/git/unstable

build.xml

@@ -281,6 +281,10 @@
         <equals arg1="${gatk.target}" arg2="private" casesensitive="false" />
       </condition>
 
+      <condition property="include.external">
+          <available file="${external.dir}"/>
+      </condition>
+
       <condition property="include.contracts">
         <equals arg1="${use.contracts}" arg2="true" />
       </condition>
@@ -331,7 +335,7 @@
         </javac>
     </target>
 
-    <target name="gatk.compile.external.source" depends="gatk.compile.public.source,gatk.compile.private.source">
+    <target name="gatk.compile.external.source" depends="gatk.compile.public.source,gatk.compile.private.source" if="include.external">
         <subant target="compile" genericantfile="build.xml">
             <property name="build.dir" value="${external.build.dir}" />
             <property name="dist.dir" value="${external.dist.dir}" />
@@ -761,6 +765,7 @@
     <property name="java.test.classes" value="${build.dir}/java/testclasses"/>
     <property name="java.public.test.classes" value="${java.test.classes}/public"/>
     <property name="java.private.test.classes" value="${java.test.classes}/private"/>
+    <property name="java.external.test.classes" value="${java.test.classes}/external"/>
     <property name="java.public.test.sources" value="${public.dir}/java/test"/>
     <property name="java.private.test.sources" value="${private.dir}/java/test"/>
     <property name="scala.test.classes" value="${build.dir}/scala/testclasses"/>
@@ -811,7 +816,23 @@
         </javac>
     </target>
 
-    <target name="test.java.compile" depends="test.java.public.compile, test.java.private.compile"/>
+    <target name="test.java.external.compile" depends="dist,test.init.compile,test.java.public.compile" if="include.external">
+        <mkdir dir="${java.external.test.classes}"/>
+        <echo message="Sting: Compiling external test cases!"/>
+        <javac fork="true" memoryMaximumSize="512m" destdir="${java.external.test.classes}" debug="true" optimize="on" tempdir="${java.io.tmpdir}" srcdir="${external.dir}">
+            <include name="*/test/**/*.java"/>
+            <classpath>
+                <path refid="external.dependencies" />
+                <pathelement location="${java.public.test.classes}"/>
+                <pathelement location="${java.classes}"/>
+                <pathelement location="${java.contracts}"/>
+                <pathelement location="${testng.jar}"/>
+            </classpath>
+            <compilerarg value="-proc:none"/>
+        </javac>
+    </target>
+
+    <target name="test.java.compile" depends="test.java.public.compile, test.java.private.compile, test.java.external.compile"/>
 
     <target name="test.scala.public.compile" depends="test.java.compile,scala.compile" if="scala.include">
         <mkdir dir="${scala.public.test.classes}"/>
@@ -852,6 +873,7 @@
         <pathelement location="${java.contracts}" />
         <pathelement location="${java.public.test.classes}" />
         <pathelement location="${java.private.test.classes}" />
+        <pathelement location="${java.external.test.classes}" />
         <pathelement location="${scala.public.test.classes}" />
         <pathelement location="${scala.private.test.classes}" />
         <pathelement location="${R.tar.dir}" />
@@ -934,6 +956,9 @@
                 <classfileset dir="${java.private.test.classes}" erroronmissingdir="false">
                     <include name="**/@{testtype}.class" if="include.private"/>
                 </classfileset>
+                <classfileset dir="${java.external.test.classes}" erroronmissingdir="false">
+                    <include name="**/@{testtype}.class" if="include.external"/>
+                </classfileset>
                 <classfileset dir="${scala.public.test.classes}" erroronmissingdir="false">
                     <include name="**/@{testtype}*.class" if="scala.include"/>
                 </classfileset>

public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/FisherStrand.java

@@ -54,7 +54,7 @@ public class FisherStrand extends InfoFieldAnnotation implements StandardAnnotat
     private static final double MIN_PVALUE = 1E-320;
 
     public Map<String, Object> annotate(RefMetaDataTracker tracker, AnnotatorCompatibleWalker walker, ReferenceContext ref, Map<String, AlignmentContext> stratifiedContexts, VariantContext vc) {
-        if ( ! vc.isVariant() || vc.isFiltered() )
+        if ( !vc.isVariant() )
             return null;
 
         int[][] table;

public/java/src/org/broadinstitute/sting/gatk/walkers/filters/VariantFiltrationWalker.java

@@ -139,6 +139,12 @@ public class VariantFiltrationWalker extends RodWalker<Integer, Integer> {
     @Argument(fullName="missingValuesInExpressionsShouldEvaluateAsFailing", doc="When evaluating the JEXL expressions, missing values should be considered failing the expression", required=false)
     protected Boolean FAIL_MISSING_VALUES = false;
 
+    /**
+     * Invalidate previous filters applied to the VariantContext, applying only the filters here
+     */
+    @Argument(fullName="invalidatePreviousFilters",doc="Remove previous filters applied to the VCF",required=false)
+    boolean invalidatePrevious = false;
+
     // JEXL expressions for the filters
     List<VariantContextUtils.JexlVCMatchExp> filterExps;
     List<VariantContextUtils.JexlVCMatchExp> genotypeFilterExps;
@@ -215,6 +221,9 @@ public class VariantFiltrationWalker extends RodWalker<Integer, Integer> {
 
         for ( VariantContext vc : VCs ) {
 
+            if ( invalidatePrevious ) {
+                vc = (new VariantContextBuilder(vc)).filters(new HashSet<String>()).make();
+            }
             // filter based on previous mask position
             if ( previousMaskPosition != null &&                                       // we saw a previous mask site
                  previousMaskPosition.getContig().equals(vc.getChr()) &&               // it's on the same contig

public/java/src/org/broadinstitute/sting/gatk/walkers/recalibration/CountCovariatesWalker.java

@@ -155,7 +155,9 @@ public class CountCovariatesWalker extends LocusWalker<CountCovariatesWalker.Cou
     /////////////////////////////
     /**
      * This algorithm treats every reference mismatch as an indication of error. However, real genetic variation is expected to mismatch the reference,
-     * so it is critical that a database of known polymorphic sites is given to the tool in order to skip over those sites.
+     * so it is critical that a database of known polymorphic sites is given to the tool in order to skip over those sites. This tool accepts any number of RodBindings (VCF, Bed, etc.)
+     * for use as this database. For users wishing to exclude an interval list of known variation simply use -XL my.interval.list to skip over processing those sites.
+     * Please note however that the statistics reported by the tool will not accurately reflected those sites skipped by the -XL argument.
      */
     @Input(fullName="knownSites", shortName = "knownSites", doc="A database of known polymorphic sites to skip over in the recalibration algorithm", required=false)
     public List<RodBinding<Feature>> knownSites = Collections.emptyList();

public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibratorEngine.java

@@ -27,7 +27,6 @@ package org.broadinstitute.sting.gatk.walkers.variantrecalibration;
 
 import org.apache.log4j.Logger;
 import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
-import org.broadinstitute.sting.utils.exceptions.UserException;
 
 import java.util.List;
 

public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/CombineVariants.java

@@ -105,7 +105,7 @@ public class CombineVariants extends RodWalker<Integer, Integer> {
      * and each named argument will be labeled as such in the output (i.e., set=name rather than
      * set=variants2).  The order of arguments does not matter unless except for the naming, so
      * if you provide an rod priority list and no explicit names than variants, variants2, etc
-     * are techincally order dependent.  It is strongly recommended to provide explicit names when
+     * are technically order dependent.  It is strongly recommended to provide explicit names when
      * a rod priority list is provided.
      */
     @Input(fullName="variant", shortName = "V", doc="Input VCF file", required=true)

public/java/src/org/broadinstitute/sting/utils/MathUtils.java

@@ -1569,4 +1569,48 @@ public class MathUtils {
         }
         return shuffled;
     }
+
+    /**
+     * Vector operations
+     */
+    public static double[] vectorSum(double v1[], double v2[]) {
+        if (v1.length != v2.length)
+            throw new UserException("BUG: vectors v1, v2 of different size in vectorSum()");
+
+        double result[] = new double[v1.length];
+        for (int k=0; k < v1.length; k++)
+            result[k] = v1[k]+v2[k];
+
+        return result;
+    }
+
+    public static double[] scalarTimesIntVector(double a, int[] v1) {
+
+        double result[] = new double[v1.length];
+        for (int k=0; k < v1.length; k++)
+            result[k] = a*v1[k];
+
+        return result;
+    }
+
+    public static double dotProduct(double v1[], double v2[]) {
+        if (v1.length != v2.length)
+            throw new UserException("BUG: vectors v1, v2 of different size in vectorSum()");
+
+        double result = 0.0;
+        for (int k=0; k < v1.length; k++)
+            result += v1[k]*v2[k];
+
+        return result;
+
+    }
+
+    public static double[] vectorLog10(double v1[]) {
+        double result[] = new double[v1.length];
+        for (int k=0; k < v1.length; k++)
+            result[k] = Math.log10(v1[k]);
+
+        return result;
+
+    }
 }

public/java/src/org/broadinstitute/sting/utils/sam/ReadUtils.java

@@ -608,7 +608,7 @@ public class ReadUtils {
      *    Example: Locus => {read1, read2, ..., readN}
      *
      *  - Second: readToLocusMap -- a HashMap that describes for each read what loci it contributes to the coverage.
-     *    Note: Locus is a boolean array, indexed from 0 (= startLocation) to N (= stopLocation), with true meaning it contributes to the coverage.
+     *    Note: Locus is a boolean array, indexed from 0 (= startLocation) to N (= stopLocation), with value==true meaning it contributes to the coverage.
      *    Example: Read => {true, true, false, ... false}
      *
      * @param readList the list of reads to generate the association mappings

public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java

@@ -464,7 +464,11 @@ public class VariantContextUtils {
         /**
          * Requires all records present at site to be unfiltered. VCF files that don't contain the record don't influence this.
          */
-        KEEP_IF_ALL_UNFILTERED
+        KEEP_IF_ALL_UNFILTERED,
+        /**
+         * If any record is present at this site (regardless of possibly being filtered), then all such records are kept and the filters are reset.
+         */
+        KEEP_UNCONDITIONAL
     }
 
     /**
@@ -635,7 +639,7 @@ public class VariantContextUtils {
         }
 
         // if at least one record was unfiltered and we want a union, clear all of the filters
-        if ( filteredRecordMergeType == FilteredRecordMergeType.KEEP_IF_ANY_UNFILTERED && nFiltered != VCs.size() )
+        if ( (filteredRecordMergeType == FilteredRecordMergeType.KEEP_IF_ANY_UNFILTERED && nFiltered != VCs.size()) || filteredRecordMergeType == FilteredRecordMergeType.KEEP_UNCONDITIONAL )
             filters.clear();
 
 

public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIntegrationTest.java

@@ -44,7 +44,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
     public void testWithAllelesPassedIn2() {
         WalkerTest.WalkerTestSpec spec2 = new WalkerTest.WalkerTestSpec(
                 baseCommand + " --output_mode EMIT_ALL_SITES --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles " + validationDataLocation + "allelesForUG.vcf -I " + validationDataLocation + "pilot2_daughters.chr20.10k-11k.bam -o %s -L 20:10,000,000-10,025,000", 1,
-                Arrays.asList("43e7a17d95b1a0cf72e669657794d802"));
+                Arrays.asList("1899bdb956c62bbcbf160b18cd3aea60"));
         executeTest("test MultiSample Pilot2 with alleles passed in and emitting all sites", spec2);
     }
 
@@ -275,7 +275,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
                 baseCommandIndels + " --output_mode EMIT_ALL_SITES --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles "
                         + validationDataLocation + "indelAllelesForUG.vcf -I " + validationDataLocation +
                         "pilot2_daughters.chr20.10k-11k.bam -o %s -L 20:10,000,000-10,100,000", 1,
-                Arrays.asList("36ce53ae4319718ad9c8ae391deebc8c"));
+                Arrays.asList("320f61c87253aba77d6dc782242cba8b"));
         executeTest("test MultiSample Pilot2 indels with alleles passed in and emitting all sites", spec2);
     }
 

public/scala/src/org/broadinstitute/sting/queue/library/ipf/vcf/VCFExtractIntervals.scala

@@ -26,9 +26,6 @@ class VCFExtractIntervals(inVCF: File, outList: File, useFilterSites: Boolean) e
     var cur : String = null
     if ( elems.hasNext ) {
       cur = elems.next
-    } else {
-      out.printf("%s%n",prev)
-    }
       while ( elems.hasNext ) {
         out.printf("%s%n",prev)
         while ( cur.equals(prev) && elems.hasNext && !cur.equals("") ) {
@@ -39,10 +36,13 @@ class VCFExtractIntervals(inVCF: File, outList: File, useFilterSites: Boolean) e
           if ( elems.hasNext ) {
             prev = cur
             cur = elems.next
-        } else {
+          }
+        }
+      }
+      out.printf("%s%n",prev)
       out.printf("%s%n",cur)
-        }
+    } else {
-      }
+      out.printf("%s%n",prev)
     }
 
     out.close