Merge remote-tracking branch 'unstable/master'

.gitignore

@@ -0,0 +1,20 @@
+/*.bam
+/*.bai
+/*.bed
+*.idx
+*~
+/*.vcf
+/*.txt
+/*.csh
+/.*
+/*.pdf
+/*.eval
+*.ipr
+*.iws
+*.iml
+.DS_Store
+queueScatterGather
+/foo*
+/bar*
+integrationtests/
+public/testdata/onTheFlyOutputTest.vcf

public/R/plot_residualError_OtherCovariate.R → public/R/scripts/org/broadinstitute/sting/analyzecovariates/plot_residualError_OtherCovariate.R

@@ -1,5 +1,7 @@
 #!/bin/env Rscript
 
+library(tools)
+
 args <- commandArgs(TRUE)
 verbose = TRUE
 
@@ -47,6 +49,9 @@ if( is.numeric(c$Covariate) ) {
 }
 dev.off()
 
+if (exists('compactPDF')) {
+  compactPDF(outfile)
+}
 
 #
 # Plot mean quality versus the covariate
@@ -69,6 +74,10 @@ if( is.numeric(c$Covariate) ) {
 }
 dev.off()
 
+if (exists('compactPDF')) {
+  compactPDF(outfile)
+}
+
 #
 # Plot histogram of the covariate
 #
@@ -106,3 +115,7 @@ if( is.numeric(c$Covariate) ) {
 	axis(2,axTicks(2), format(axTicks(2), scientific=F))
 }
 dev.off()
+
+if (exists('compactPDF')) {
+  compactPDF(outfile)
+}

public/R/plot_residualError_QualityScoreCovariate.R → public/R/scripts/org/broadinstitute/sting/analyzecovariates/plot_residualError_QualityScoreCovariate.R

@@ -1,5 +1,7 @@
 #!/bin/env Rscript
 
+library(tools)
+
 args <- commandArgs(TRUE)
 
 input = args[1]
@@ -33,6 +35,10 @@ points(f$Qreported, f$Qempirical, type="p", col="maroon1", pch=16)
 abline(0,1, lty=2)
 dev.off()
 
+if (exists('compactPDF')) {
+  compactPDF(outfile)
+}
+
 #
 # Plot Q empirical histogram
 #
@@ -52,6 +58,10 @@ points(hst2$f.Qempirical, hst2$f.nBases, type="h", lwd=4, col="maroon1")
 axis(2,axTicks(2), format(axTicks(2), scientific=F))
 dev.off()
 
+if (exists('compactPDF')) {
+  compactPDF(outfile)
+}
+
 #
 # Plot Q reported histogram
 #
@@ -68,3 +78,7 @@ plot(hst$e.Qreported, hst$e.nBases, type="h", lwd=4, xlim=c(0,maxQ), ylim=c(0,yM
 points(hst2$f.Qreported, hst2$f.nBases, type="h", lwd=4, col="maroon1")
 axis(2,axTicks(2), format(axTicks(2), scientific=F))
 dev.off()
+
+if (exists('compactPDF')) {
+  compactPDF(outfile)
+}

public/R/plot_Tranches.R → public/R/scripts/org/broadinstitute/sting/gatk/walkers/variantrecalibration/plot_Tranches.R

@@ -1,5 +1,7 @@
 #!/bin/env Rscript
 
+library(tools)
+
 args <- commandArgs(TRUE)
 verbose = TRUE
 
@@ -85,3 +87,7 @@ if ( ! is.null(sensitivity) ) {
 }
 
 dev.off()
+
+if (exists('compactPDF')) {
+  compactPDF(outfile)
+}

public/R/queueJobReport.R → public/R/scripts/org/broadinstitute/sting/queue/util/queueJobReport.R

@@ -12,28 +12,28 @@ if ( onCMDLine ) {
   inputFileName = args[1]
   outputPDF = args[2]
 } else {
-  #inputFileName = "~/Desktop/broadLocal/GATK/unstable/report.txt"
-  inputFileName = "/humgen/gsa-hpprojects/dev/depristo/oneOffProjects/[email protected]"
+  inputFileName = "~/Desktop/broadLocal/GATK/unstable/wgs.jobreport.txt"
+  #inputFileName = "/humgen/gsa-hpprojects/dev/depristo/oneOffProjects/[email protected]"
   #inputFileName = "/humgen/gsa-hpprojects/dev/depristo/oneOffProjects/rodPerformanceGoals/history/report.082711.txt"
   outputPDF = NA
 }
 
-RUNTIME_UNITS = "(sec)"
-ORIGINAL_UNITS_TO_SECONDS = 1/1000
+RUNTIME_UNITS = "(hours)"
+ORIGINAL_UNITS_TO_RUNTIME_UNITS = 1/1000/60/60
 
 # 
 # Helper function to aggregate all of the jobs in the report across all tables
 #
 allJobsFromReport <- function(report) {
-  names <- c("jobName", "startTime", "analysisName", "doneTime", "exechosts")
+  names <- c("jobName", "startTime", "analysisName", "doneTime", "exechosts", "runtime")
   sub <- lapply(report, function(table) table[,names])
   do.call("rbind", sub)
 }
 
 #
 # Creates segmentation plots of time (x) vs. job (y) with segments for the duration of the job
 #
-plotJobsGantt <- function(gatkReport, sortOverall) {
+plotJobsGantt <- function(gatkReport, sortOverall, includeText) {
   allJobs = allJobsFromReport(gatkReport)
   if ( sortOverall ) {
     title = "All jobs, by analysis, by start time"
@@ -44,16 +44,18 @@ plotJobsGantt <- function(gatkReport, sortOverall) {
   }
   allJobs$index = 1:nrow(allJobs)
   minTime = min(allJobs$startTime)
-  allJobs$relStartTime = allJobs$startTime - minTime
+  allJobs$relStartTime = allJobs$startTime - minTime 
   allJobs$relDoneTime = allJobs$doneTime - minTime
   allJobs$ganttName = paste(allJobs$jobName, "@", allJobs$exechosts)
   maxRelTime = max(allJobs$relDoneTime)
   p <- ggplot(data=allJobs, aes(x=relStartTime, y=index, color=analysisName))
-  p <- p + geom_segment(aes(xend=relDoneTime, yend=index), size=2, arrow=arrow(length = unit(0.1, "cm")))
-  p <- p + geom_text(aes(x=relDoneTime, label=ganttName, hjust=-0.2), size=2)
+  p <- p + theme_bw()
+  p <- p + geom_segment(aes(xend=relDoneTime, yend=index), size=1, arrow=arrow(length = unit(0.1, "cm")))
+  if ( includeText )
+    p <- p + geom_text(aes(x=relDoneTime, label=ganttName, hjust=-0.2), size=2)
   p <- p + xlim(0, maxRelTime * 1.1)
-  p <- p + xlab(paste("Start time (relative to first job)", RUNTIME_UNITS))
-  p <- p + ylab("Job")
+  p <- p + xlab(paste("Start time, relative to first job", RUNTIME_UNITS))
+  p <- p + ylab("Job number")
   p <- p + opts(title=title)
   print(p)
 }
@@ -119,7 +121,7 @@ plotGroup <- function(groupTable) {
   if ( length(groupAnnotations) == 1 && dim(sub)[1] > 1 ) {
     # todo -- how do we group by annotations?
     p <- ggplot(data=sub, aes(x=runtime)) + geom_histogram()
-    p <- p + xlab("runtime in seconds") + ylab("No. of jobs")
+    p <- p + xlab(paste("runtime", RUNTIME_UNITS)) + ylab("No. of jobs")
     p <- p + opts(title=paste("Job runtime histogram for", name))
     print(p)
   }
@@ -139,9 +141,9 @@ print(paste("Project          :", inputFileName))
 
 convertUnits <- function(gatkReportData) {
   convertGroup <- function(g) {
-    g$runtime = g$runtime * ORIGINAL_UNITS_TO_SECONDS
-    g$startTime = g$startTime * ORIGINAL_UNITS_TO_SECONDS
-    g$doneTime = g$doneTime * ORIGINAL_UNITS_TO_SECONDS
+    g$runtime = g$runtime * ORIGINAL_UNITS_TO_RUNTIME_UNITS
+    g$startTime = g$startTime * ORIGINAL_UNITS_TO_RUNTIME_UNITS
+    g$doneTime = g$doneTime * ORIGINAL_UNITS_TO_RUNTIME_UNITS
     g
   }
   lapply(gatkReportData, convertGroup)
@@ -157,8 +159,8 @@ if ( ! is.na(outputPDF) ) {
   pdf(outputPDF, height=8.5, width=11)
 } 
 
-plotJobsGantt(gatkReportData, T)
-plotJobsGantt(gatkReportData, F)
+plotJobsGantt(gatkReportData, T, F)
+plotJobsGantt(gatkReportData, F, F)
 plotProgressByTime(gatkReportData)
 for ( group in gatkReportData ) {
  plotGroup(group)

public/R/src/gsalib/R/gsa.read.gatkreport.R → public/R/src/org/broadinstitute/sting/utils/R/gsalib/R/gsa.read.gatkreport.R

@@ -99,5 +99,5 @@ gsa.read.gatkreport <- function(filename) {
         .gsa.assignGATKTableToEnvironment(tableName, tableHeader, tableRows, tableEnv);
     }
 
-    gatkreport = as.list(tableEnv);
+    gatkreport = as.list(tableEnv, all.names=TRUE);
 }

public/java/src/org/broadinstitute/sting/alignment/AlignmentValidationWalker.java

@@ -25,7 +25,6 @@
 
 package org.broadinstitute.sting.alignment;
 
-import net.sf.samtools.SAMRecord;
 import org.broadinstitute.sting.alignment.bwa.BWAConfiguration;
 import org.broadinstitute.sting.alignment.bwa.BWTFiles;
 import org.broadinstitute.sting.alignment.bwa.c.BWACAligner;
@@ -35,6 +34,7 @@
 import org.broadinstitute.sting.gatk.walkers.ReadWalker;
 import org.broadinstitute.sting.utils.BaseUtils;
 import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
+import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
 
 import java.util.Iterator;
 
@@ -72,12 +72,13 @@ public void initialize() {
 
     /**
      * Aligns a read to the given reference.
+     *
      * @param ref Reference over the read.  Read will most likely be unmapped, so ref will be null.
      * @param read Read to align.
      * @return Number of reads aligned by this map (aka 1).
      */
     @Override
-    public Integer map(ReferenceContext ref, SAMRecord read, ReadMetaDataTracker metaDataTracker) {
+    public Integer map(ReferenceContext ref, GATKSAMRecord read, ReadMetaDataTracker metaDataTracker) {
         //logger.info(String.format("examining read %s", read.getReadName()));
 
         byte[] bases = read.getReadBases();

public/java/src/org/broadinstitute/sting/alignment/AlignmentWalker.java

@@ -39,6 +39,7 @@
 import org.broadinstitute.sting.gatk.refdata.ReadMetaDataTracker;
 import org.broadinstitute.sting.gatk.walkers.ReadWalker;
 import org.broadinstitute.sting.gatk.walkers.WalkerName;
+import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
 
 import java.io.File;
 
@@ -92,12 +93,13 @@ public void initialize() {
 
     /**
      * Aligns a read to the given reference.
+     *
      * @param ref Reference over the read.  Read will most likely be unmapped, so ref will be null.
      * @param read Read to align.
      * @return Number of alignments found for this read.
      */
     @Override
-    public Integer map(ReferenceContext ref, SAMRecord read, ReadMetaDataTracker metaDataTracker) {
+    public Integer map(ReferenceContext ref, GATKSAMRecord read, ReadMetaDataTracker metaDataTracker) {
         SAMRecord alignedRead = aligner.align(read,header);
         out.addAlignment(alignedRead);
         return 1;

public/java/src/org/broadinstitute/sting/alignment/CountBestAlignmentsWalker.java

@@ -25,7 +25,6 @@
 
 package org.broadinstitute.sting.alignment;
 
-import net.sf.samtools.SAMRecord;
 import org.broadinstitute.sting.alignment.bwa.BWAConfiguration;
 import org.broadinstitute.sting.alignment.bwa.BWTFiles;
 import org.broadinstitute.sting.alignment.bwa.c.BWACAligner;
@@ -34,6 +33,7 @@
 import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
 import org.broadinstitute.sting.gatk.refdata.ReadMetaDataTracker;
 import org.broadinstitute.sting.gatk.walkers.ReadWalker;
+import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
 
 import java.io.PrintStream;
 import java.util.Iterator;
@@ -79,12 +79,13 @@ public void initialize() {
 
     /**
      * Aligns a read to the given reference.
+     *
      * @param ref Reference over the read.  Read will most likely be unmapped, so ref will be null.
      * @param read Read to align.
      * @return Number of alignments found for this read.
      */
     @Override
-    public Integer map(ReferenceContext ref, SAMRecord read, ReadMetaDataTracker metaDataTracker) {
+    public Integer map(ReferenceContext ref, GATKSAMRecord read, ReadMetaDataTracker metaDataTracker) {
         Iterator<Alignment[]> alignmentIterator = aligner.getAllAlignments(read.getReadBases()).iterator();
         if(alignmentIterator.hasNext()) {
             int numAlignments = alignmentIterator.next().length;