Merge remote-tracking branch 'unstable/master'

build.xml

@@ -76,7 +76,8 @@
     <property name="staging.dir" value="staging" />
     <property name="default.executable" value="none" />
 
-    <!-- Javadoc/Scaladoc directories -->
+    <!-- GATKDocs/Javadoc/Scaladoc directories -->
+    <property name="gatkdocs.dir" value="gatkdocs" />
     <property name="javadoc.dir" value="javadoc" />
     <property name="scaladoc.dir" value="scaladoc" />
 
@@ -92,12 +93,10 @@
     <!-- To disable for test targets, run with -Duse.contracts=false -->
     <!-- To enable for non-test targets, run with -Duse.contracts=true -->
     <property name="java.contracts.dir" value="${build.dir}/java/contracts" />
-    <property name="contracts.version" value="1.0-20110609" />
+    <property name="contracts.version" value="1.0-r139" />
     <property name="cofoja.jar" value="${lib.dir}/cofoja-${contracts.version}.jar"/>
     <property name="contract.dump.dir" value="dump" />
 
-    <property name="gatkdocs.dir" value="gatkdocs" />
-
     <!-- do we want to halt on failure of a unit test? default to yes (Bamboo uses 'no') -->
     <property name="halt" value="yes" />
 
@@ -208,19 +207,19 @@
         <include name="**/*.java" />
     </fileset>
 
-    <pathconvert property="external.build.dir">
-        <path path="${java.classes}"/>
-    </pathconvert>
+    <path id="external.build.dir">
+        <path path="${java.classes}" />
+    </path>
 
-    <pathconvert property="external.dist.dir">
+    <path id="external.dist.dir">
         <path path="${dist.dir}" />
-    </pathconvert>
+    </path>
 
     <!-- GATK dependencies consist of 3rd party plugins plus compiled GATK classes -->
-    <pathconvert property="external.gatk.classpath">
+    <path id="external.gatk.classpath">
         <path path="${java.classes}"/>
         <path refid="external.dependencies" />
-    </pathconvert>
+    </path>
 
     <!-- the path for resources that need to go into the GATK jar;
          any additional resources should go into this set. -->
@@ -430,9 +429,9 @@
 
     <target name="gatk.compile.external.source" depends="gatk.compile.internal.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}" />
-            <property name="gatk.classpath" value="${external.gatk.classpath}" />
+            <property name="build.dir" refid="external.build.dir" />
+            <property name="dist.dir" refid="external.dist.dir" />
+            <property name="gatk.classpath" refid="external.gatk.classpath" />
             <fileset dir="${external.dir}" includes="*/build.xml" erroronmissingdir="false" />
         </subant>
     </target>
@@ -680,9 +679,9 @@
         </jar>
 
         <subant target="dist" genericantfile="build.xml">
-            <property name="build.dir" value="${external.build.dir}" />
-            <property name="dist.dir" value="${external.dist.dir}" />
-            <property name="gatk.classpath" value="${external.gatk.classpath}" />
+            <property name="build.dir" refid="external.build.dir" />
+            <property name="dist.dir" refid="external.dist.dir" />
+            <property name="gatk.classpath" refid="external.gatk.classpath" />
             <fileset dir="${external.dir}" includes="*/build.xml" erroronmissingdir="false" />
         </subant>
     </target>

ivy.xml

@@ -52,7 +52,7 @@
         <dependency org="gov.nist" name="Jama" rev="1.0.2"/>
 
         <!-- Dependencies for the graph aligner -->
-        <dependency org="org.jgrapht" name="jgrapht-jdk1.5" rev="0.7.3"/>
+        <dependency org="net.sf.jgrapht" name="jgrapht" rev="0.8.3"/>
 
         <!-- Dependencies for the html walker documention -->
         <dependency org="org.freemarker" name="freemarker" rev="2.3.18"/>
@@ -87,7 +87,7 @@
         <dependency org="com.google.code.caliper" name="caliper" rev="1.0-SNAPSHOT" conf="test"/>
 
         <!-- Contracts for Java and dependencies -->
-        <dependency org="com.google.code.cofoja" name="cofoja" rev="1.0-20110609"/>
+        <dependency org="com.google.code.cofoja" name="cofoja" rev="1.0-r139"/>
         <dependency org="asm" name="asm-all" rev="3.3.1"/>
 
         <!-- POI, for reading pipeline files -->

public/R/scripts/org/broadinstitute/sting/gatk/walkers/bqsr/BQSR.R → public/R/scripts/org/broadinstitute/sting/utils/recalibration/BQSR.R

@@ -1,8 +1,18 @@
 library("ggplot2")
+library(gplots)
+library("reshape")
+library("grid")
 library("tools") #For compactPDF in R 2.13+
+library(gsalib)
 
-args <- commandArgs(TRUE)
+
+if ( interactive() ) {
+  args <- c("NA12878.6.1.dedup.realign.recal.bqsr.grp.csv", "NA12878.6.1.dedup.realign.recal.bqsr.grp", NA)
+} else {
+  args <- commandArgs(TRUE)
+} 
 data <- read.csv(args[1])
+gsa.report <- gsa.read.gatkreport(args[2])
 data <- within(data, EventType <- factor(EventType, levels = rev(levels(EventType))))
 
 numRG = length(unique(data$ReadGroup))
@@ -82,20 +92,45 @@ for(cov in levels(data$CovariateName)) {    # for each covariate in turn
 
     p <- ggplot(d, aes(x=CovariateValue)) +
       xlab(paste(cov,"Covariate")) +
-      ylab("Number of Observations") +
+      ylab("No. of Observations (area normalized)") +
       blankTheme
-    d <- p + geom_histogram(aes(fill=Recalibration,weight=Observations),alpha=0.6,binwidth=1,position="identity") + scale_fill_manual(values=c("maroon1","blue")) + facet_grid(.~EventType) +     
-      scale_y_continuous(formatter="comma")
-
+    d <- p + geom_histogram(aes(fill=Recalibration,weight=Observations,y=..ndensity..),alpha=0.6,binwidth=1,position="identity")
+    d <- d + scale_fill_manual(values=c("maroon1","blue"))
+    d <- d + facet_grid(.~EventType) 
+#    d <- d + scale_y_continuous(formatter="comma")
   }
 }
 
-pdf(args[2],height=9,width=15)
+if ( ! is.na(args[3]) )
+  pdf(args[3],height=9,width=15)
+
+#frame()
+textplot(gsa.report$Arguments, show.rownames=F)
+title(
+  main="GATK BaseRecalibration report",
+  sub=date())
+
 distributeGraphRows(list(a,b,c), c(1,1,1))
 distributeGraphRows(list(d,e,f), c(1,1,1))
-dev.off()
 
+# format the overall information
+rt0 <- data.frame(
+  ReadGroup = gsa.report$RecalTable0$ReadGroup,
+  EventType = gsa.report$RecalTable0$EventType,
+  EmpiricalQuality = sprintf("%.1f", gsa.report$RecalTable0$EmpiricalQuality),
+  EstimatedQReported = sprintf("%.1f", gsa.report$RecalTable0$EstimatedQReported),
+  Observations = sprintf("%.2e", gsa.report$RecalTable0$Observations),
+  Errors = sprintf("%.2e", gsa.report$RecalTable0$Errors))  
+textplot(t(rt0), show.colnames=F)
+title("Overall error rates by event type")
+
+# plot per quality score recalibration table
+textplot(gsa.report$RecalTable1, show.rownames=F)
+title("Error rates by event type and initial quality score")
 
-if (exists('compactPDF')) {
-  compactPDF(args[2])
+if ( ! is.na(args[3]) ) {
+  dev.off()
+  if (exists('compactPDF')) {
+    compactPDF(args[2])
+  }
 }

public/R/src/org/broadinstitute/sting/utils/R/gsalib/R/gsa.variantqc.utils.R

@@ -207,7 +207,7 @@ plotVariantQC <- function(metrics, measures, requestedStrat = "Sample",
 
   if ( requestedStrat == "Sample" ) {
     perSampleGraph <- perSampleGraph + geom_text(aes(label=strat), size=1.5) + geom_blank() # don't display a scale
-    perSampleGraph <- perSampleGraph + scale_x_discrete("Sample (ordered by nSNPs)", formatter=function(x) "")
+    perSampleGraph <- perSampleGraph + scale_x_discrete("Sample (ordered by nSNPs)")
   } else { # by AlleleCount
     perSampleGraph <- perSampleGraph + geom_point(aes(size=log10(nobs))) #+ geom_smooth(aes(weight=log10(nobs)))
     perSampleGraph <- perSampleGraph + scale_x_log10("AlleleCount")

public/java/src/net/sf/picard/reference/FastaSequenceIndexBuilder.java

@@ -208,6 +208,7 @@ else if (basesThisLine != basesPerLine || bytesPerLine != bytesRead - endOfLastL
                         break;
                 }
             }
+            in.close();
             return sequenceIndex;
         }
         catch (IOException e) {

public/java/src/org/broadinstitute/sting/alignment/reference/bwt/Bases.java

@@ -12,10 +12,10 @@
  */
 public class Bases implements Iterable<Byte>
 {
-    public static byte A = 'A';
-    public static byte C = 'C';
-    public static byte G = 'G';
-    public static byte T = 'T';
+    public static final byte A = 'A';
+    public static final byte C = 'C';
+    public static final byte G = 'G';
+    public static final byte T = 'T';
 
     public static final Bases instance = new Bases();
 

public/java/src/org/broadinstitute/sting/commandline/ArgumentDefinition.java

@@ -26,6 +26,7 @@
 package org.broadinstitute.sting.commandline;
 
 import org.broadinstitute.sting.utils.Utils;
+import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
 
 import java.lang.annotation.Annotation;
 import java.util.List;
@@ -147,6 +148,9 @@ public ArgumentDefinition( ArgumentIOType ioType,
         this.exclusiveOf = exclusiveOf;
         this.validation = validation;
         this.validOptions = validOptions;
+
+        validateName(shortName);
+        validateName(fullName);
     }
 
     /**
@@ -192,6 +196,9 @@ else if( !isFullNameProvided )
         else
             shortName = null;
 
+        validateName(shortName);
+        validateName(fullName);
+
         this.ioType = ioType;
         this.argumentType = argumentType;
         this.fullName = fullName;
@@ -277,4 +284,14 @@ public static String getValidationRegex( Annotation annotation ) {
         String validation = (String)CommandLineUtils.getValue(annotation, "validation");
         return validation.trim().length() > 0 ? validation.trim() : null;
     }
+
+    /**
+     * Make sure the argument's name is valid
+     *
+     * @param name
+     */
+    private void validateName(final String name) {
+        if ( name != null && name.startsWith("-") )
+            throw new ReviewedStingException("Invalid argument definition: " + name + " begins with a -");
+    }
 }

public/java/src/org/broadinstitute/sting/commandline/ArgumentDefinitionGroup.java

@@ -55,10 +55,8 @@ public ArgumentDefinitionGroup( String groupName, List<ArgumentDefinition> argum
      * Does the name of this argument group match the name of another?
      */
     public boolean groupNameMatches( ArgumentDefinitionGroup other ) {
-        if( this.groupName == null && other.groupName == null )
-            return true;
-        if( this.groupName == null && other.groupName != null )
-            return false;
+        if( this.groupName == null )
+            return other.groupName == null;
         return this.groupName.equals(other.groupName);
     }
 

public/java/src/org/broadinstitute/sting/commandline/ArgumentTypeDescriptor.java

@@ -53,7 +53,7 @@ public abstract class ArgumentTypeDescriptor {
     /**
      * our log, which we want to capture anything from org.broadinstitute.sting
      */
-    protected static Logger logger = Logger.getLogger(ArgumentTypeDescriptor.class);
+    protected static final Logger logger = Logger.getLogger(ArgumentTypeDescriptor.class);
 
     /**
      * Fetch the given descriptor from the descriptor repository.

public/java/src/org/broadinstitute/sting/commandline/ParsingMethod.java

@@ -120,8 +120,8 @@ public static Tags parseTags(String argument, String tagString) {
      */
     private static final String TAG_TEXT = "[\\w\\-\\.\\=]*";
 
-    public static ParsingMethod FullNameParsingMethod = new ParsingMethod(Pattern.compile(String.format("\\s*--(%1$s)(?:\\:(%2$s(?:,%2$s)*))?\\s*",ARGUMENT_TEXT,TAG_TEXT)),
+    public static final ParsingMethod FullNameParsingMethod = new ParsingMethod(Pattern.compile(String.format("\\s*--(%1$s)(?:\\:(%2$s(?:,%2$s)*))?\\s*",ARGUMENT_TEXT,TAG_TEXT)),
                                                                           ArgumentDefinitions.FullNameDefinitionMatcher) {};
-    public static ParsingMethod ShortNameParsingMethod = new ParsingMethod(Pattern.compile(String.format("\\s*-(%1$s)(?:\\:(%2$s(?:,%2$s)*))?\\s*",ARGUMENT_TEXT,TAG_TEXT)),
+    public static final ParsingMethod ShortNameParsingMethod = new ParsingMethod(Pattern.compile(String.format("\\s*-(%1$s)(?:\\:(%2$s(?:,%2$s)*))?\\s*",ARGUMENT_TEXT,TAG_TEXT)),
                                                                            ArgumentDefinitions.ShortNameDefinitionMatcher) {};
 }

public/java/src/org/broadinstitute/sting/gatk/CommandLineExecutable.java

@@ -130,8 +130,8 @@ private void authorizeGATKRun() {
              getArgumentCollection().phoneHomeType == GATKRunReport.PhoneHomeOption.STDOUT ) {
             if ( getArgumentCollection().gatkKeyFile == null ) {
                 throw new UserException("Running with the -et NO_ET or -et STDOUT option requires a GATK Key file. " +
-                                        "Please see http://www.broadinstitute.org/gsa/wiki/index.php/Phone_home " +
-                                        "for more information and instructions on how to obtain a key.");
+                                        "Please see " + GATKRunReport.PHONE_HOME_DOCS_URL +
+                                        " for more information and instructions on how to obtain a key.");
             }
             else {
                 PublicKey gatkPublicKey = CryptUtils.loadGATKDistributedPublicKey();

public/java/src/org/broadinstitute/sting/gatk/CommandLineGATK.java

@@ -130,6 +130,12 @@ private static void checkForMaskedUserErrors(final Throwable t) {
 
         // can't close tribble index when writing
         if ( message.indexOf("Unable to close index for") != -1 )
+            exitSystemWithUserError(new UserException(t.getCause() == null ? message : t.getCause().getMessage()));
+
+        // disk is full
+        if ( message.indexOf("No space left on device") != -1 )
+            exitSystemWithUserError(new UserException(t.getMessage()));
+        if ( t.getCause() != null && t.getCause().getMessage().indexOf("No space left on device") != -1 )
             exitSystemWithUserError(new UserException(t.getCause().getMessage()));
     }
 

public/java/src/org/broadinstitute/sting/gatk/GenomeAnalysisEngine.java

@@ -233,10 +233,6 @@ public Object execute() {
         if (args.nonDeterministicRandomSeed)
             resetRandomGenerator(System.currentTimeMillis());
 
-        // TODO -- REMOVE ME WHEN WE STOP BCF testing
-        if ( args.USE_SLOW_GENOTYPES )
-            GenotypeBuilder.MAKE_FAST_BY_DEFAULT = false;
-
         // if the use specified an input BQSR recalibration table then enable on the fly recalibration
         if (args.BQSR_RECAL_FILE != null)
             setBaseRecalibration(args.BQSR_RECAL_FILE, args.quantizationLevels, args.disableIndelQuals, args.PRESERVE_QSCORES_LESS_THAN, args.emitOriginalQuals);
@@ -797,6 +793,14 @@ private SAMDataSource createReadsDataSource(GATKArgumentCollection argCollection
         if ( getWalkerBAQApplicationTime() == BAQ.ApplicationTime.FORBIDDEN && argCollection.BAQMode != BAQ.CalculationMode.OFF)
             throw new UserException.BadArgumentValue("baq", "Walker cannot accept BAQ'd base qualities, and yet BAQ mode " + argCollection.BAQMode + " was requested.");
 
+        if (argCollection.removeProgramRecords && argCollection.keepProgramRecords)
+            throw new UserException.BadArgumentValue("rpr / kpr", "Cannot enable both options");
+
+        boolean removeProgramRecords = argCollection.removeProgramRecords || walker.getClass().isAnnotationPresent(RemoveProgramRecords.class);
+
+        if (argCollection.keepProgramRecords)
+            removeProgramRecords = false;
+
         return new SAMDataSource(
                 samReaderIDs,
                 threadAllocation,
@@ -813,7 +817,8 @@ private SAMDataSource createReadsDataSource(GATKArgumentCollection argCollection
                 getWalkerBAQQualityMode(),
                 refReader,
                 getBaseRecalibration(),
-                argCollection.defaultBaseQualities);
+                argCollection.defaultBaseQualities,
+                removeProgramRecords);
     }
 
     /**
@@ -840,20 +845,9 @@ private List<ReferenceOrderedDataSource> getReferenceOrderedDataSources(Collecti
                                                                             SAMSequenceDictionary sequenceDictionary,
                                                                             GenomeLocParser genomeLocParser,
                                                                             ValidationExclusion.TYPE validationExclusionType) {
-        VCFHeader header = null;
-        if ( getArguments().repairVCFHeader != null ) {
-            try {
-                final PositionalBufferedStream pbs = new PositionalBufferedStream(new FileInputStream(getArguments().repairVCFHeader));
-                header = (VCFHeader)new VCFCodec().readHeader(pbs).getHeaderValue();
-                pbs.close();
-            } catch ( IOException e ) {
-                throw new UserException.CouldNotReadInputFile(getArguments().repairVCFHeader, e);
-            }
-        }
-
-        RMDTrackBuilder builder = new RMDTrackBuilder(sequenceDictionary,genomeLocParser, header, validationExclusionType);
+        final RMDTrackBuilder builder = new RMDTrackBuilder(sequenceDictionary,genomeLocParser, validationExclusionType);
 
-        List<ReferenceOrderedDataSource> dataSources = new ArrayList<ReferenceOrderedDataSource>();
+        final List<ReferenceOrderedDataSource> dataSources = new ArrayList<ReferenceOrderedDataSource>();
         for (RMDTriplet fileDescriptor : referenceMetaDataFiles)
             dataSources.add(new ReferenceOrderedDataSource(fileDescriptor,
                                                            builder,