zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

No more hunting down R "resources". As a tradeoff Rscript cannot be specified on the commandline and will be found in the environment path. 

Other minor cleanup.
This commit is contained in:
Khalid Shakir 2011-10-26 23:05:41 -04:00
parent 8c4dbce6d8
commit b80d407dc7
22 changed files with 407 additions and 223 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
build.xml
View File

@ -38,6 +38,7 @@
<property name="java.private.source.dir" value="${private.dir}/java/src" />
<property name="java.classes" value="${build.dir}/java/classes" />
<property name="R.public.scripts.dir" value="${public.dir}/R/scripts" />
<property name="R.private.scripts.dir" value="${private.dir}/R/scripts" />
<property name="R.public.src.dir" value="${public.dir}/R/src" />
<!-- Legacy: Installing libraries back into the source directory
instead of the build or dist directory... intentionally avoids ant clean?? -->
@ -571,6 +572,9 @@
<fileset dir="${R.public.scripts.dir}">
<include name="**/utils/**/*.R"/>
</fileset>
<fileset dir="${R.private.scripts.dir}" erroronmissingdir="false">
<include name="**/utils/**/*.R"/>
</fileset>
<manifest>
<attribute name="Premain-Class" value="org.broadinstitute.sting.utils.instrumentation.Sizeof" />
</manifest>
@ -603,6 +607,10 @@
<include name="**/gatk/**/*.R"/>
<include name="**/alignment/**/*.R"/>
</fileset>
<fileset dir="${R.private.scripts.dir}" erroronmissingdir="false">
<include name="**/gatk/**/*.R"/>
<include name="**/alignment/**/*.R"/>
</fileset>
<manifest>
<attribute name="Main-Class" value="org.broadinstitute.sting.gatk.CommandLineGATK" />
</manifest>
@ -621,6 +629,10 @@
<include name="**/analyzecovariates/**/*.R"/>
<include name="**/gatk/walkers/recalibration/**/*.R"/>
</fileset>
<fileset dir="${R.private.scripts.dir}" erroronmissingdir="false">
<include name="**/analyzecovariates/**/*.R"/>
<include name="**/gatk/walkers/recalibration/**/*.R"/>
</fileset>
<manifest>
<attribute name="Main-Class" value="org.broadinstitute.sting.analyzecovariates.AnalyzeCovariates" />
</manifest>
@ -653,6 +665,9 @@
<fileset dir="${R.public.scripts.dir}">
<include name="org/broadinstitute/sting/queue/**/*.R"/>
</fileset>
<fileset dir="${R.private.scripts.dir}" erroronmissingdir="false">
<include name="org/broadinstitute/sting/queue/**/*.R"/>
</fileset>
<manifest>
<attribute name="Main-Class" value="org.broadinstitute.sting.queue.QCommandLine" />
</manifest>
@ -842,6 +857,7 @@
<pathelement location="${scala.private.test.classes}" />
<pathelement location="${R.tar.dir}" />
<pathelement location="${R.public.scripts.dir}" />
<pathelement location="${R.private.scripts.dir}" />
</path>
<path id="testng.gatk.releasetest.classpath">

0
public/R/plot_residualError_OtherCovariate.R → public/R/scripts/org/broadinstitute/sting/analyzecovariates/plot_residualError_OtherCovariate.R
View File

0
public/R/plot_residualError_QualityScoreCovariate.R → public/R/scripts/org/broadinstitute/sting/analyzecovariates/plot_residualError_QualityScoreCovariate.R
View File

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

126
public/java/src/org/broadinstitute/sting/analyzecovariates/AnalyzeCovariates.java
View File

@ -25,6 +25,9 @@
package org.broadinstitute.sting.analyzecovariates;
import org.apache.commons.io.FileUtils;
import org.apache.commons.io.IOUtils;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.commandline.Argument;
import org.broadinstitute.sting.commandline.Hidden;
import org.broadinstitute.sting.commandline.CommandLineProgram;
@ -33,14 +36,16 @@ import org.broadinstitute.sting.gatk.walkers.recalibration.Covariate;
import org.broadinstitute.sting.gatk.walkers.recalibration.RecalDatum;
import org.broadinstitute.sting.gatk.walkers.recalibration.RecalibrationArgumentCollection;
import org.broadinstitute.sting.utils.R.RScriptExecutor;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.classloader.PluginManager;
import org.broadinstitute.sting.utils.exceptions.DynamicClassResolutionException;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.help.DocumentedGATKFeature;
import org.broadinstitute.sting.utils.io.Resource;
import org.broadinstitute.sting.utils.text.XReadLines;
import java.io.*;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.Map;
import java.util.regex.Pattern;
@ -71,15 +76,13 @@ import java.util.regex.Pattern;
* </ul>
*
* <p>
* NOTE: For those running this tool externally from the Broad, it is crucial to note that both the -Rscript and -resources options
* must be changed from the default. -Rscript needs to point to your installation of Rscript (this is the scripting version of R,
* not the interactive version) while -resources needs to point to the folder holding the R scripts that are used. For those using
* this tool as part of the Binary Distribution the -resources should point to the resources folder that is part of the tarball.
* For those using this tool by building from the git repository the -resources should point to the R/ subdirectory of the Sting checkout.
* NOTE: Rscript needs to be in your environment PATH (this is the scripting version of R, not the interactive version).
* See <a target="r-project" href="http://www.r-project.org">http://www.r-project.org</a> for more info on how to download and install R.
*
* <p>
* See the GATK wiki for a tutorial and example recalibration accuracy plots.
* http://www.broadinstitute.org/gsa/wiki/index.php/Base_quality_score_recalibration
* <a target="gatkwiki" href="http://www.broadinstitute.org/gsa/wiki/index.php/Base_quality_score_recalibration"
* >http://www.broadinstitute.org/gsa/wiki/index.php/Base_quality_score_recalibration</a>
*
* <h2>Input</h2>
* <p>
@ -91,7 +94,6 @@ import java.util.regex.Pattern;
* java -Xmx4g -jar AnalyzeCovariates.jar \
* -recalFile /path/to/recal.table.csv \
* -outputDir /path/to/output_dir/ \
* -resources resources/ \
* -ignoreQ 5
* </pre>
*
@ -101,6 +103,11 @@ import java.util.regex.Pattern;
groupName = "AnalyzeCovariates",
summary = "Package to plot residual accuracy versus error covariates for the base quality score recalibrator")
public class AnalyzeCovariates extends CommandLineProgram {
final private static Logger logger = Logger.getLogger(AnalyzeCovariates.class);
private static final String PLOT_RESDIUAL_ERROR_QUALITY_SCORE_COVARIATE = "plot_residualError_QualityScoreCovariate.R";
private static final String PLOT_RESDIUAL_ERROR_OTHER_COVARIATE = "plot_residualError_OtherCovariate.R";
private static final String PLOT_INDEL_QUALITY_RSCRIPT = "plot_indelQuality.R";
/////////////////////////////
// Command Line Arguments
@ -114,11 +121,7 @@ public class AnalyzeCovariates extends CommandLineProgram {
@Input(fullName = "recal_file", shortName = "recalFile", doc = "The input recal csv file to analyze", required = false)
private String RECAL_FILE = "output.recal_data.csv";
@Argument(fullName = "output_dir", shortName = "outputDir", doc = "The directory in which to output all the plots and intermediate data files", required = false)
private String OUTPUT_DIR = "analyzeCovariates/";
@Argument(fullName = "path_to_Rscript", shortName = "Rscript", doc = "The path to your implementation of Rscript. For Broad users this is maybe /broad/software/free/Linux/redhat_5_x86_64/pkgs/r_2.12.0/bin/Rscript", required = false)
private String PATH_TO_RSCRIPT = "Rscript";
@Argument(fullName = "path_to_resources", shortName = "resources", doc = "Path to resources folder holding the Sting R scripts.", required = false)
private String PATH_TO_RESOURCES = "public/R/";
private File OUTPUT_DIR = new File("analyzeCovariates");
@Argument(fullName = "ignoreQ", shortName = "ignoreQ", doc = "Ignore bases with reported quality less than this number.", required = false)
private int IGNORE_QSCORES_LESS_THAN = 5;
@Argument(fullName = "numRG", shortName = "numRG", doc = "Only process N read groups. Default value: -1 (process all read groups)", required = false)
@ -154,29 +157,26 @@ public class AnalyzeCovariates extends CommandLineProgram {
protected int execute() {
// create the output directory where all the data tables and plots will go
try {
Process p = Runtime.getRuntime().exec("mkdir " + OUTPUT_DIR);
} catch (IOException e) {
System.out.println("Couldn't create directory: " + OUTPUT_DIR);
System.out.println("User is responsible for making sure the output directory exists.");
}
if( !OUTPUT_DIR.endsWith("/") ) { OUTPUT_DIR = OUTPUT_DIR + "/"; }
if( !PATH_TO_RESOURCES.endsWith("/") ) { PATH_TO_RESOURCES = PATH_TO_RESOURCES + "/"; }
if (!OUTPUT_DIR.exists() && !OUTPUT_DIR.mkdirs())
throw new UserException.BadArgumentValue("--output_dir/-outDir", "Unable to create output directory: " + OUTPUT_DIR);
if (!RScriptExecutor.RSCRIPT_EXISTS)
Utils.warnUser(logger, "Rscript not found in environment path. Plots will not be generated.");
// initialize all the data from the csv file and allocate the list of covariates
System.out.println("Reading in input csv file...");
logger.info("Reading in input csv file...");
initializeData();
System.out.println("...Done!");
logger.info("...Done!");
// output data tables for Rscript to read in
System.out.println("Writing out intermediate tables for R...");
logger.info("Writing out intermediate tables for R...");
writeDataTables();
System.out.println("...Done!");
logger.info("...Done!");
// perform the analysis using Rscript and output the plots
System.out.println("Calling analysis R scripts and writing out figures...");
logger.info("Calling analysis R scripts and writing out figures...");
callRScripts();
System.out.println("...Done!");
logger.info("...Done!");
return 0;
}
@ -287,37 +287,40 @@ public class AnalyzeCovariates extends CommandLineProgram {
if(NUM_READ_GROUPS_TO_PROCESS == -1 || ++numReadGroups <= NUM_READ_GROUPS_TO_PROCESS) {
String readGroup = readGroupKey.toString();
RecalDatum readGroupDatum = (RecalDatum) dataManager.getCollapsedTable(0).data.get(readGroupKey);
System.out.print("Writing out data tables for read group: " + readGroup + "\twith " + readGroupDatum.getNumObservations() + " observations" );
System.out.println("\tand aggregate residual error = " + String.format("%.3f", readGroupDatum.empiricalQualDouble(0, MAX_QUALITY_SCORE) - readGroupDatum.getEstimatedQReported()));
logger.info(String.format(
"Writing out data tables for read group: %s\twith %s observations\tand aggregate residual error = %.3f",
readGroup, readGroupDatum.getNumObservations(),
readGroupDatum.empiricalQualDouble(0, MAX_QUALITY_SCORE) - readGroupDatum.getEstimatedQReported()));
// for each covariate
for( int iii = 1; iii < requestedCovariates.size(); iii++ ) {
Covariate cov = requestedCovariates.get(iii);
// Create a PrintStream
PrintStream output = null;
File outputFile = new File(OUTPUT_DIR, readGroup + "." + cov.getClass().getSimpleName()+ ".dat");
PrintStream output;
try {
output = new PrintStream(new FileOutputStream(OUTPUT_DIR + readGroup + "." + cov.getClass().getSimpleName()+ ".dat"));
} catch (FileNotFoundException e) {
System.err.println("Can't create file: " + OUTPUT_DIR + readGroup + "." + cov.getClass().getSimpleName()+ ".dat");
System.exit(-1);
output = new PrintStream(FileUtils.openOutputStream(outputFile));
} catch (IOException e) {
throw new UserException.CouldNotCreateOutputFile(outputFile, e);
}
// Output the header
output.println("Covariate\tQreported\tQempirical\tnMismatches\tnBases");
try {
// Output the header
output.println("Covariate\tQreported\tQempirical\tnMismatches\tnBases");
for( Object covariateKey : ((Map)dataManager.getCollapsedTable(iii).data.get(readGroupKey)).keySet()) {
output.print( covariateKey.toString() + "\t" ); // Covariate
RecalDatum thisDatum = (RecalDatum)((Map)dataManager.getCollapsedTable(iii).data.get(readGroupKey)).get(covariateKey);
output.print( String.format("%.3f", thisDatum.getEstimatedQReported()) + "\t" ); // Qreported
output.print( String.format("%.3f", thisDatum.empiricalQualDouble(0, MAX_QUALITY_SCORE)) + "\t" ); // Qempirical
output.print( thisDatum.getNumMismatches() + "\t" ); // nMismatches
output.println( thisDatum.getNumObservations() ); // nBases
for( Object covariateKey : ((Map)dataManager.getCollapsedTable(iii).data.get(readGroupKey)).keySet()) {
output.print( covariateKey.toString() + "\t" ); // Covariate
RecalDatum thisDatum = (RecalDatum)((Map)dataManager.getCollapsedTable(iii).data.get(readGroupKey)).get(covariateKey);
output.print( String.format("%.3f", thisDatum.getEstimatedQReported()) + "\t" ); // Qreported
output.print( String.format("%.3f", thisDatum.empiricalQualDouble(0, MAX_QUALITY_SCORE)) + "\t" ); // Qempirical
output.print( thisDatum.getNumMismatches() + "\t" ); // nMismatches
output.println( thisDatum.getNumObservations() ); // nBases
}
} finally {
// Close the PrintStream
IOUtils.closeQuietly(output);
}
// Close the PrintStream
output.close();
}
} else {
break;
@ -327,10 +330,6 @@ public class AnalyzeCovariates extends CommandLineProgram {
}
private void callRScripts() {
RScriptExecutor.RScriptArgumentCollection argumentCollection =
new RScriptExecutor.RScriptArgumentCollection(PATH_TO_RSCRIPT, Arrays.asList(PATH_TO_RESOURCES));
RScriptExecutor executor = new RScriptExecutor(argumentCollection, true);
int numReadGroups = 0;
// for each read group
@ -338,23 +337,32 @@ public class AnalyzeCovariates extends CommandLineProgram {
if(++numReadGroups <= NUM_READ_GROUPS_TO_PROCESS || NUM_READ_GROUPS_TO_PROCESS == -1) {
String readGroup = readGroupKey.toString();
System.out.println("Analyzing read group: " + readGroup);
logger.info("Analyzing read group: " + readGroup);
// for each covariate
for( int iii = 1; iii < requestedCovariates.size(); iii++ ) {
Covariate cov = requestedCovariates.get(iii);
final String outputFilename = OUTPUT_DIR + readGroup + "." + cov.getClass().getSimpleName()+ ".dat";
final File outputFile = new File(OUTPUT_DIR, readGroup + "." + cov.getClass().getSimpleName()+ ".dat");
if (DO_INDEL_QUALITY) {
executor.callRScripts("plot_indelQuality.R", outputFilename,
cov.getClass().getSimpleName().split("Covariate")[0]); // The third argument is the name of the covariate in order to make the plots look nice
RScriptExecutor executor = new RScriptExecutor();
executor.addScript(new Resource(PLOT_INDEL_QUALITY_RSCRIPT, AnalyzeCovariates.class));
// The second argument is the name of the covariate in order to make the plots look nice
executor.addArgs(outputFile, cov.getClass().getSimpleName().split("Covariate")[0]);
executor.exec();
} else {
if( iii == 1 ) {
// Analyze reported quality
executor.callRScripts("plot_residualError_QualityScoreCovariate.R", outputFilename,
IGNORE_QSCORES_LESS_THAN, MAX_QUALITY_SCORE, MAX_HISTOGRAM_VALUE); // The third argument is the Q scores that should be turned pink in the plot because they were ignored
RScriptExecutor executor = new RScriptExecutor();
executor.addScript(new Resource(PLOT_RESDIUAL_ERROR_QUALITY_SCORE_COVARIATE, AnalyzeCovariates.class));
// The second argument is the Q scores that should be turned pink in the plot because they were ignored
executor.addArgs(outputFile, IGNORE_QSCORES_LESS_THAN, MAX_QUALITY_SCORE, MAX_HISTOGRAM_VALUE);
executor.exec();
} else { // Analyze all other covariates
executor.callRScripts("plot_residualError_OtherCovariate.R", outputFilename,
cov.getClass().getSimpleName().split("Covariate")[0]); // The third argument is the name of the covariate in order to make the plots look nice
RScriptExecutor executor = new RScriptExecutor();
executor.addScript(new Resource(PLOT_RESDIUAL_ERROR_OTHER_COVARIATE, AnalyzeCovariates.class));
// The second argument is the name of the covariate in order to make the plots look nice
executor.addArgs(outputFile, cov.getClass().getSimpleName().split("Covariate")[0]);
executor.exec();
}
}
}

60
public/java/src/org/broadinstitute/sting/gatk/walkers/variantrecalibration/VariantRecalibrator.java
View File

@ -35,9 +35,11 @@ import org.broadinstitute.sting.gatk.walkers.RodWalker;
import org.broadinstitute.sting.gatk.walkers.TreeReducible;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.QualityUtils;
import org.broadinstitute.sting.utils.R.RScriptExecutor;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.collections.ExpandingArrayList;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.io.Resource;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.variantcontext.VariantContextUtils;
@ -101,6 +103,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
public static final String VQS_LOD_KEY = "VQSLOD"; // Log odds ratio of being a true variant versus being false under the trained gaussian mixture model
public static final String CULPRIT_KEY = "culprit"; // The annotation which was the worst performing in the Gaussian mixture model, likely the reason why the variant was filtered out
private static final String PLOT_TRANCHES_RSCRIPT = "plot_Tranches.R";
@ArgumentCollection private VariantRecalibratorArgumentCollection VRAC = new VariantRecalibratorArgumentCollection();
@ -158,12 +161,8 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
private double[] TS_TRANCHES = new double[] {100.0, 99.9, 99.0, 90.0};
@Argument(fullName="ignore_filter", shortName="ignoreFilter", doc="If specified the variant recalibrator will use variants even if the specified filter name is marked in the input VCF file", required=false)
private String[] IGNORE_INPUT_FILTERS = null;
@Argument(fullName="path_to_Rscript", shortName = "Rscript", doc = "The path to your implementation of Rscript. For Broad users this is maybe /broad/software/free/Linux/redhat_5_x86_64/pkgs/r_2.12.0/bin/Rscript", required=false)
private String PATH_TO_RSCRIPT = "Rscript";
@Argument(fullName="rscript_file", shortName="rscriptFile", doc="The output rscript file generated by the VQSR to aid in visualization of the input data and learned model", required=false)
private String RSCRIPT_FILE = null;
@Argument(fullName = "path_to_resources", shortName = "resources", doc = "Path to resources folder holding the Sting R scripts.", required=false)
private String PATH_TO_RESOURCES = "public/R/";
@Output(fullName="rscript_file", shortName="rscriptFile", doc="The output rscript file generated by the VQSR to aid in visualization of the input data and learned model", required=false)
private File RSCRIPT_FILE = null;
@Argument(fullName="ts_filter_level", shortName="ts_filter_level", doc="The truth sensitivity level at which to start filtering, used here to indicate filtered variants in the model reporting plots", required=false)
private double TS_FILTER_LEVEL = 99.0;
@ -192,9 +191,13 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
//---------------------------------------------------------------------------------------------------------------
public void initialize() {
if( !PATH_TO_RESOURCES.endsWith("/") ) { PATH_TO_RESOURCES = PATH_TO_RESOURCES + "/"; }
dataManager = new VariantDataManager( new ArrayList<String>(Arrays.asList(USE_ANNOTATIONS)), VRAC );
if (RSCRIPT_FILE != null && !RScriptExecutor.RSCRIPT_EXISTS)
Utils.warnUser(logger, String.format(
"Rscript not found in environment path. %s will be generated but PDF plots will not.",
RSCRIPT_FILE));
if( IGNORE_INPUT_FILTERS != null ) {
ignoreInputFilterSet.addAll( Arrays.asList(IGNORE_INPUT_FILTERS) );
}
@ -327,20 +330,13 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
createVisualizationScript( dataManager.getRandomDataForPlotting( 6000 ), goodModel, badModel, lodCutoff );
}
// Execute Rscript command to create the tranche plot
// Print out the command line to make it clear to the user what is being executed and how one might modify it
final String rScriptTranchesCommandLine = PATH_TO_RSCRIPT + " " + PATH_TO_RESOURCES + "plot_Tranches.R" + " " + TRANCHES_FILE.getAbsolutePath() + " " + TARGET_TITV;
logger.info( "Executing: " + rScriptTranchesCommandLine );
// Execute the RScript command to plot the table of truth values
try {
Process p;
p = Runtime.getRuntime().exec( rScriptTranchesCommandLine );
p.waitFor();
} catch ( Exception e ) {
Utils.warnUser("Unable to execute the RScript command. While not critical to the calculations themselves, the script outputs a report that is extremely useful for confirming that the recalibration proceded as expected. We highly recommend trying to rerun the script manually if possible.");
}
RScriptExecutor executor = new RScriptExecutor();
executor.addScript(new Resource(PLOT_TRANCHES_RSCRIPT, VariantRecalibrator.class));
executor.addArgs(TRANCHES_FILE.getAbsoluteFile(), TARGET_TITV);
// Print out the command line to make it clear to the user what is being executed and how one might modify it
logger.info("Executing: " + executor.getApproximateCommandLine());
executor.exec();
}
private void createVisualizationScript( final ExpandingArrayList<VariantDatum> randomData, final GaussianMixtureModel goodModel, final GaussianMixtureModel badModel, final double lodCutoff ) {
@ -348,15 +344,18 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
try {
stream = new PrintStream(RSCRIPT_FILE);
} catch( FileNotFoundException e ) {
throw new UserException.CouldNotCreateOutputFile(RSCRIPT_FILE, "", e);
throw new UserException.CouldNotCreateOutputFile(RSCRIPT_FILE, e);
}
// We make extensive use of the ggplot2 R library: http://had.co.nz/ggplot2/
stream.println("library(ggplot2)");
// For compactPDF in R 2.13+
stream.println("library(tools)");
createArrangeFunction( stream );
stream.println("pdf(\"" + RSCRIPT_FILE + ".pdf\")"); // Unfortunately this is a huge pdf file, BUGBUG: need to work on reducing the file size
stream.println("outputPDF <- \"" + RSCRIPT_FILE + ".pdf\"");
stream.println("pdf(outputPDF)"); // Unfortunately this is a huge pdf file, BUGBUG: need to work on reducing the file size
for(int iii = 0; iii < USE_ANNOTATIONS.length; iii++) {
for( int jjj = iii + 1; jjj < USE_ANNOTATIONS.length; jjj++) {
@ -431,18 +430,17 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
}
stream.println("dev.off()");
stream.println("if (exists(\"compactPDF\")) {");
stream.println("compactPDF(ouputPDF)");
stream.println("}");
stream.close();
// Execute Rscript command to generate the clustering plots
final String rScriptTranchesCommandLine = PATH_TO_RSCRIPT + " " + RSCRIPT_FILE;
logger.info( "Executing: " + rScriptTranchesCommandLine );
try {
Process p;
p = Runtime.getRuntime().exec( rScriptTranchesCommandLine );
p.waitFor();
} catch ( Exception e ) {
Utils.warnUser("Unable to execute the RScript command. While not critical to the calculations themselves, the script outputs a report that is extremely useful for visualizing the recalibration results. We highly recommend trying to rerun the script manually if possible.");
}
RScriptExecutor executor = new RScriptExecutor();
executor.addScript(RSCRIPT_FILE);
logger.info("Executing: " + executor.getApproximateCommandLine());
executor.exec();
}
// The Arrange function is how we place the 4 model plots on one page

119
public/java/src/org/broadinstitute/sting/utils/R/RScriptExecutor.java
View File

@ -27,8 +27,6 @@ package org.broadinstitute.sting.utils.R;
import org.apache.commons.io.FileUtils;
import org.apache.commons.lang.StringUtils;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.commandline.Advanced;
import org.broadinstitute.sting.commandline.Argument;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.exceptions.StingException;
import org.broadinstitute.sting.utils.exceptions.UserException;
@ -36,48 +34,33 @@ import org.broadinstitute.sting.utils.io.IOUtils;
import org.broadinstitute.sting.utils.io.Resource;
import org.broadinstitute.sting.utils.runtime.ProcessController;
import org.broadinstitute.sting.utils.runtime.ProcessSettings;
import org.broadinstitute.sting.utils.runtime.RuntimeUtils;
import java.io.File;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
/**
* Generic service for executing RScripts
*/
public class RScriptExecutor {
private static final String RSCRIPT_BINARY = "Rscript";
private static final File RSCRIPT_PATH = RuntimeUtils.which(RSCRIPT_BINARY);
public static final boolean RSCRIPT_EXISTS = (RSCRIPT_PATH != null);
private static final String RSCRIPT_MISSING_MESSAGE = "Please add the Rscript directory to your environment ${PATH}";
/**
* our log
*/
private static Logger logger = Logger.getLogger(RScriptExecutor.class);
public static class RScriptArgumentCollection {
@Advanced
@Argument(fullName = "path_to_Rscript", shortName = "Rscript", doc = "The path to your implementation of Rscript. Defaults Rscript meaning to use the first available on the environment PATH. For Broad users should 'use R-2.12' or later.", required = false)
public String PATH_TO_RSCRIPT = "Rscript";
@Advanced
@Argument(fullName = "path_to_Rresources", shortName = "Rresources", doc = "Path to resources folder holding the Sting R scripts.", required = false)
public List<String> PATH_TO_RESOURCES = Arrays.asList("public/R/", "private/R/");
public RScriptArgumentCollection() {}
/* For testing and convenience */
public RScriptArgumentCollection(final String PATH_TO_RSCRIPT, final List<String> PATH_TO_RESOURCES) {
this.PATH_TO_RSCRIPT = PATH_TO_RSCRIPT;
this.PATH_TO_RESOURCES = PATH_TO_RESOURCES;
}
}
private final RScriptArgumentCollection myArgs;
private final boolean exceptOnError;
private boolean exceptOnError = false;
private final List<RScriptLibrary> libraries = new ArrayList<RScriptLibrary>();
private final List<Resource> scriptResources = new ArrayList<Resource>();
private final List<File> scriptFiles = new ArrayList<File>();
private final List<String> args = new ArrayList<String>();
public RScriptExecutor(final RScriptArgumentCollection myArgs, final boolean exceptOnError) {
this.myArgs = myArgs;
public void setExceptOnError(boolean exceptOnError) {
this.exceptOnError = exceptOnError;
}
@ -103,7 +86,27 @@ public class RScriptExecutor {
this.args.add(arg.toString());
}
public void exec() {
public String getApproximateCommandLine() {
StringBuilder command = new StringBuilder("Rscript");
for (Resource script: this.scriptResources)
command.append(" (resource)").append(script.getFullPath());
for (File script: this.scriptFiles)
command.append(" ").append(script.getAbsolutePath());
for (String arg: this.args)
command.append(" ").append(arg);
return command.toString();
}
public boolean exec() {
if (!RSCRIPT_EXISTS) {
if (exceptOnError) {
throw new UserException.CannotExecuteRScript(RSCRIPT_MISSING_MESSAGE);
} else {
logger.warn("Skipping: " + getApproximateCommandLine());
return false;
}
}
List<File> tempFiles = new ArrayList<File>();
try {
File tempLibDir = IOUtils.tempDir("R.", ".lib");
@ -126,7 +129,7 @@ public class RScriptExecutor {
expression.append(");");
for (RScriptLibrary library: this.libraries) {
expression.append("require('").append(library.getLibraryName()).append("', lib.loc=tempLibDir);");
expression.append("library('").append(library.getLibraryName()).append("', lib.loc=tempLibDir);");
}
}
@ -142,7 +145,7 @@ public class RScriptExecutor {
String[] cmd = new String[this.args.size() + 3];
int i = 0;
cmd[i++] = myArgs.PATH_TO_RSCRIPT;
cmd[i++] = RSCRIPT_BINARY;
cmd[i++] = "-e";
cmd[i++] = expression.toString();
for (String arg: this.args)
@ -156,52 +159,30 @@ public class RScriptExecutor {
ProcessController controller = ProcessController.getThreadLocal();
logger.debug("Executing: " + Utils.join(" ", cmd));
logger.debug("Result: " + controller.exec(processSettings).getExitValue());
if (logger.isDebugEnabled()) {
logger.debug("Executing:");
for (String arg: cmd)
logger.debug(" " + arg);
}
int exitValue = controller.exec(processSettings).getExitValue();
logger.debug("Result: " + exitValue);
if (exitValue != 0)
throw new RScriptExecutorException(
"RScript exited with " + exitValue +
(logger.isDebugEnabled() ? "" : ". Run with -l DEBUG for more info."));
return true;
} catch (StingException e) {
generateException(e);
if (exceptOnError) {
throw e;
} else {
logger.warn(e.getMessage());
return false;
}
} finally {
for (File temp: tempFiles)
FileUtils.deleteQuietly(temp);
}
}
public void callRScripts(String scriptName, Object... scriptArgs) {
final File pathToScript = findScript(scriptName);
if (pathToScript == null) return; // we failed but shouldn't exception out
addScript(pathToScript);
addArgs(scriptArgs);
exec();
}
public File findScript(final String scriptName) {
for ( String pathToResource : myArgs.PATH_TO_RESOURCES ) {
final File f = new File(pathToResource + "/" + scriptName);
if ( f.exists() ) {
if ( f.canRead() )
return f;
else
generateException("Script exists but couldn't be read: " + scriptName);
}
}
generateException("Couldn't find script: " + scriptName + " in " + myArgs.PATH_TO_RESOURCES);
return null;
}
private void generateException(String msg) {
generateException(msg, null);
}
private void generateException(Throwable e) {
generateException("", e);
}
private void generateException(String msg, Throwable e) {
if ( exceptOnError )
throw new UserException(msg, e);
else
logger.warn(msg + (e == null ? "" : ":" + e.getMessage()));
}
}

33
public/java/src/org/broadinstitute/sting/utils/R/RScriptExecutorException.java 100644
View File

@ -0,0 +1,33 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.utils.R;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
public class RScriptExecutorException extends ReviewedStingException {
public RScriptExecutorException(String msg) {
super(msg);
}
}

10
public/java/src/org/broadinstitute/sting/utils/Utils.java
View File

@ -80,20 +80,24 @@ public class Utils {
}
public static void warnUser(final String msg) {
warnUser(logger, msg);
}
public static void warnUser(final Logger logger, final String msg) {
logger.warn(String.format("********************************************************************************"));
logger.warn(String.format("* WARNING:"));
logger.warn(String.format("*"));
prettyPrintWarningMessage(msg);
prettyPrintWarningMessage(logger, msg);
logger.warn(String.format("********************************************************************************"));
}
/**
* pretty print the warning message supplied
*
* @param logger logger for the message
* @param message the message
*/
private static void prettyPrintWarningMessage(String message) {
private static void prettyPrintWarningMessage(Logger logger, String message) {
StringBuilder builder = new StringBuilder(message);
while (builder.length() > 70) {
int space = builder.lastIndexOf(" ", 70);

3
public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java
View File

@ -237,6 +237,9 @@ public class UserException extends ReviewedStingException {
}
public static class CannotExecuteRScript extends UserException {
public CannotExecuteRScript(String message) {
super(String.format("Unable to execute RScript command: " + message));
}
public CannotExecuteRScript(String message, Exception e) {
super(String.format("Unable to execute RScript command: " + message), e);
}

12
public/java/src/org/broadinstitute/sting/utils/io/IOUtils.java
View File

@ -90,6 +90,18 @@ public class IOUtils {
}
}
/**
* Writes content to a temp file and returns the path to the temporary file.
*
* @param content to write.
* @param prefix Prefix for the temp file.
* @param suffix Suffix for the temp file.
* @return the path to the temp file.
*/
public static File writeTempFile(String content, String prefix, String suffix) {
return writeTempFile(content, prefix, suffix, null);
}
/**
* Writes content to a temp file and returns the path to the temporary file.
*

13
public/java/src/org/broadinstitute/sting/utils/io/Resource.java
View File

@ -24,6 +24,8 @@
package org.broadinstitute.sting.utils.io;
import java.io.File;
/**
* Stores a resource by path and a relative class.
*/
@ -50,4 +52,15 @@ public class Resource {
public String getPath() {
return path;
}
public String getFullPath() {
if (relativeClass == null)
return path;
if (new File(path).isAbsolute())
return path;
return String.format("%s%s%s",
relativeClass.getPackage().getName().replace('.', File.separatorChar),
File.separator,
path);
}
}

58
public/java/src/org/broadinstitute/sting/utils/runtime/RuntimeUtils.java 100644
View File

@ -0,0 +1,58 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.utils.runtime;
import org.apache.commons.lang.StringUtils;
import java.io.File;
public class RuntimeUtils {
public static final String[] PATHS;
static {
String path = System.getenv("PATH");
if (path == null)
path = System.getenv("path");
if (path == null) {
PATHS = new String[0];
} else {
PATHS = StringUtils.split(path, File.pathSeparatorChar);
}
}
/**
* Returns the path to an executable or null if it doesn't exist.
* @param executable Relative path
* @return The absolute file path.
*/
public static File which(String executable) {
for (String path: PATHS) {
File file = new File(path, executable);
if (file.exists())
return file.getAbsoluteFile();
}
return null;
}
}

114
public/java/test/org/broadinstitute/sting/utils/R/RScriptExecutorUnitTest.java
View File

@ -26,60 +26,84 @@ package org.broadinstitute.sting.utils.R;
import org.apache.commons.io.FileUtils;
import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.io.IOUtils;
import org.testng.Assert;
import org.testng.annotations.Test;
import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
/**
* Basic unit test for RScriptExecutor in reduced reads
*/
public class RScriptExecutorUnitTest extends BaseTest {
final static String testrscript = "print(\"hello, world\")\n";
final static String publicRScript = "plot_Tranches.R";
// --------------------------------------------------------------------------------
//
// Difference testing routines
//
// --------------------------------------------------------------------------------
private void testOne(String script, String pathToRscript, String anotherSearchPath, boolean exceptOnError) {
RScriptExecutor.RScriptArgumentCollection collection =
new RScriptExecutor.RScriptArgumentCollection();
if ( pathToRscript != null )
collection.PATH_TO_RSCRIPT = pathToRscript;
if ( anotherSearchPath != null ) {
List<String> x = new ArrayList<String>(collection.PATH_TO_RESOURCES);
x.add(anotherSearchPath);
collection.PATH_TO_RESOURCES = x;
}
RScriptExecutor executor = new RScriptExecutor(collection, exceptOnError);
executor.callRScripts(script);
}
private static final String HELLO_WORLD_SCRIPT = "print('hello, world')";
private static final String GSALIB_LOADED_SCRIPT = "if (!'package:gsalib' %in% search()) stop('gsalib not loaded')";
@Test
public void testPublic() { testOne(publicRScript, null, null, true); }
@Test(expectedExceptions = UserException.class)
public void testNonExistantScriptException() { testOne("does_not_exist.R", null, null, true); }
@Test()
public void testNonExistantScriptNoException() { testOne("does_not_exist.R", null, null, false); }
@Test(expectedExceptions = UserException.class)
public void testNonExistantRScriptException() { testOne(publicRScript, "badRScriptValue", null, true); }
@Test()
public void testNonExistantRScriptNoException() { testOne(publicRScript, "badRScriptValue", null, false); }
@Test()
public void testScriptInNewPath() throws IOException {
File t = createTempFile("myTestScript", ".R");
FileUtils.writeStringToFile(t, testrscript);
testOne(t.getName(), null, t.getParent(), true);
public void testRscriptExists() {
Assert.assertTrue(RScriptExecutor.RSCRIPT_EXISTS, "Rscript not found in environment ${PATH}");
}
}
@Test(dependsOnMethods = "testRscriptExists")
public void testExistingScript() {
File script = writeScript(HELLO_WORLD_SCRIPT);
try {
RScriptExecutor executor = new RScriptExecutor();
executor.addScript(script);
executor.setExceptOnError(true);
Assert.assertTrue(executor.exec(), "Exec failed");
} finally {
FileUtils.deleteQuietly(script);
}
}
@Test(dependsOnMethods = "testRscriptExists", expectedExceptions = RScriptExecutorException.class)
public void testNonExistantScriptException() {
RScriptExecutor executor = new RScriptExecutor();
executor.setExceptOnError(true);
executor.addScript(new File("does_not_exists.R"));
executor.exec();
}
@Test(dependsOnMethods = "testRscriptExists")
public void testNonExistantScriptNoException() {
logger.warn("Testing that warning is printed an no exception thrown for missing script.");
RScriptExecutor executor = new RScriptExecutor();
executor.setExceptOnError(false);
executor.addScript(new File("does_not_exists.R"));
Assert.assertFalse(executor.exec(), "Exec should have returned false when the job failed");
}
@Test(dependsOnMethods = "testRscriptExists")
public void testLibrary() {
File script = writeScript(GSALIB_LOADED_SCRIPT);
try {
RScriptExecutor executor = new RScriptExecutor();
executor.addScript(script);
executor.addLibrary(RScriptLibrary.GSALIB);
executor.setExceptOnError(true);
Assert.assertTrue(executor.exec(), "Exec failed");
} finally {
FileUtils.deleteQuietly(script);
}
}
@Test(dependsOnMethods = "testRscriptExists", expectedExceptions = RScriptExecutorException.class)
public void testLibraryMissing() {
File script = writeScript(GSALIB_LOADED_SCRIPT);
try {
RScriptExecutor executor = new RScriptExecutor();
executor.addScript(script);
// GSALIB is not added nor imported in the script
executor.setExceptOnError(true);
executor.exec();
} finally {
FileUtils.deleteQuietly(script);
}
}
private File writeScript(String content) {
return IOUtils.writeTempFile(content, "myTestScript", ".R");
}
}

2
public/java/test/org/broadinstitute/sting/utils/runtime/ProcessControllerUnitTest.java
View File

@ -168,7 +168,7 @@ public class ProcessControllerUnitTest extends BaseTest {
File input = null;
try {
String fileText = "Hello from file";
input = IOUtils.writeTempFile(fileText, "stdin.", ".txt", null);
input = IOUtils.writeTempFile(fileText, "stdin.", ".txt");
ProcessSettings job = new ProcessSettings(new String[] {"cat"});
job.getStdoutSettings().setBufferSize(-1);

41
public/java/test/org/broadinstitute/sting/utils/runtime/RuntimeUtilsUnitTest.java 100644
View File

@ -0,0 +1,41 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.utils.runtime;
import org.broadinstitute.sting.BaseTest;
import org.testng.Assert;
import org.testng.annotations.Test;
public class RuntimeUtilsUnitTest extends BaseTest {
@Test
public void testWhichExists() {
Assert.assertNotNull(RuntimeUtils.which("ls"), "Unable to locate ls");
}
@Test
public void testWhichNotExists() {
Assert.assertNull(RuntimeUtils.which("does_not_exist"), "Found nonexistent binary: does_not_exist");
}
}

7
public/packages/AnalyzeCovariates.xml
View File

@ -7,11 +7,8 @@
<!-- Recalibration analysis script -->
<class name="org.broadinstitute.sting.analyzecovariates.AnalyzeCovariates" />
<package name="org.broadinstitute.sting.gatk.walkers.recalibration" />
<!-- Supplemental scripts for graph generation, etc. -->
<dir name="org/broadinstitute/sting/analyzecovariates" includes="**/*.R" />
</dependencies>
</executable>
<resources>
<!-- Supplemental scripts for graph generation, etc. -->
<file name="public/R/plot_residualError_OtherCovariate.R" />
<file name="public/R/plot_residualError_QualityScoreCovariate.R" />
</resources>
</package>

4
public/packages/GATKEngine.xml
View File

@ -33,7 +33,9 @@
<!-- Workaround - depend on the logger impl required by JEXL -->
<package name="org.apache.commons.logging.impl" />
<!-- R packages -->
<dir name="org/broadinstitute/sting/utils/R" includes="*.tar.gz" />
<dir name="org/broadinstitute/sting/utils/R" includes="**/*.tar.gz" />
<!-- All R scripts in org.broadinstitute.sting -->
<dir name="org/broadinstitute/sting" includes="**/*.R" />
</dependencies>
</executable>
<resources>

1
public/packages/GenomeAnalysisTK.xml
View File

@ -19,7 +19,6 @@
<file name="public/java/src/org/broadinstitute/sting/gatk/walkers/qc/CountLociWalker.java" />
<file name="public/java/src/org/broadinstitute/sting/gatk/walkers/qc/CountReadsWalker.java" />
<file name="public/java/src/org/broadinstitute/sting/gatk/walkers/qc/ValidatingPileupWalker.java" />
<file name="public/R/plot_Tranches.R" />
<file name="public/R/titvFPEst.R" />
</resources>
<release>

2
public/scala/src/org/broadinstitute/sting/queue/QCommandLine.scala
View File

@ -131,7 +131,7 @@ class QCommandLine extends CommandLineProgram with Logging {
val pdfFile = new File(jobStringName + ".pdf")
logger.info("Plotting JobLogging GATKReport to file " + pdfFile)
QJobReport.plotReport(settings.rScriptArgs, reportFile, pdfFile)
QJobReport.plotReport(reportFile, pdfFile)
}
}
}

4
public/scala/src/org/broadinstitute/sting/queue/engine/QGraphSettings.scala
View File

@ -28,7 +28,6 @@ import java.io.File
import org.broadinstitute.sting.queue.QSettings
import org.broadinstitute.sting.queue.util.SystemUtils
import org.broadinstitute.sting.commandline.{Advanced, ArgumentCollection, Argument}
import org.broadinstitute.sting.utils.R.RScriptExecutor
/**
* Command line options for a QGraph.
@ -77,9 +76,6 @@ class QGraphSettings {
@Argument(fullName="disableJobReport", shortName="disabpleJobReport", doc="If provided, we will not create a job report", required=false)
var disableJobReport: Boolean = false
@ArgumentCollection
var rScriptArgs = new RScriptExecutor.RScriptArgumentCollection
@ArgumentCollection
val qSettings = new QSettings
}

5
public/scala/src/org/broadinstitute/sting/queue/util/QJobReport.scala
View File

@ -29,7 +29,6 @@ import org.broadinstitute.sting.gatk.report.{GATKReportTable, GATKReport}
import org.broadinstitute.sting.utils.exceptions.UserException
import org.broadinstitute.sting.queue.engine.JobRunInfo
import java.io.{FileOutputStream, PrintStream, File}
import org.broadinstitute.sting.utils.R.RScriptExecutor.RScriptArgumentCollection
import org.broadinstitute.sting.utils.R.{RScriptLibrary, RScriptExecutor}
import org.broadinstitute.sting.utils.io.Resource
@ -104,8 +103,8 @@ object QJobReport {
stream.close()
}
def plotReport(args: RScriptArgumentCollection, reportFile: File, pdfFile: File) {
val executor = new RScriptExecutor(args, false) // don't except on error
def plotReport(reportFile: File, pdfFile: File) {
val executor = new RScriptExecutor
executor.addLibrary(RScriptLibrary.GSALIB)
executor.addScript(new Resource(JOB_REPORT_QUEUE_SCRIPT, classOf[QJobReport]))
executor.addArgs(reportFile.getAbsolutePath, pdfFile.getAbsolutePath)