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

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

This commit is contained in:
Ami Levy Moonshine 2012-10-17 22:16:02 -04:00
parent 402ce963f9 b57df6cac8
commit acc0fb2f7a
71 changed files with 2026 additions and 903 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
.gitignore vendored
View File

@ -18,10 +18,8 @@ queueScatterGather
/bar*
integrationtests/
public/testdata/onTheFlyOutputTest.vcf
private/testdata/onTheFlyOutputTest.vcf
lib
html
gatkdocs
dist
build
resources
build/
dist/
dump/
lib/
out/

42
build.xml
View File

@ -22,7 +22,9 @@
~ OTHER DEALINGS IN THE SOFTWARE.
-->
<project name="Sting" default="dist" basedir="." xmlns:ivy="antlib:org.apache.ivy.ant">
<project name="Sting" default="dist" basedir="."
xmlns:artifact="antlib:org.apache.maven.artifact.ant"
xmlns:ivy="antlib:org.apache.ivy.ant">
<description>Compile and distribute the Sting toolkit</description>
<!-- ******************************************************************************** -->
@ -250,11 +252,14 @@
<property name="ivy.jar.file" value="ivy-${ivy.install.version}.jar"/>
<property name="ivy.settings.dir" value="settings"/>
<property file="${ivy.settings.dir}/ivysettings.properties"/>
<property name="maven-ant-tasks.install.version" value="2.1.3"/>
<property name="maven-ant-tasks.jar.file" value="maven-ant-tasks-${maven-ant-tasks.install.version}.jar"/>
<mkdir dir="${lib.dir}"/>
<mkdir dir="${ivy.jar.dir}"/>
<!-- Comment out the following two lines to build the GATK without a network connection, assuming you have all of the libraries cached already -->
<!-- Comment out the following lines to build the GATK without a network connection, assuming you have all of the libraries cached already -->
<get src="http://repo1.maven.org/maven2/org/apache/ivy/ivy/${ivy.install.version}/${ivy.jar.file}"
dest="${ivy.jar.dir}/${ivy.jar.file}"
usetimestamp="true"/>
@ -262,6 +267,15 @@
uri="antlib:org.apache.ivy.ant"
classpath="${ivy.jar.dir}/${ivy.jar.file}"/>
<get src="http://repo1.maven.org/maven2/org/apache/maven/maven-ant-tasks/${maven-ant-tasks.install.version}/${maven-ant-tasks.jar.file}"
dest="${ivy.jar.dir}/${maven-ant-tasks.jar.file}"
usetimestamp="true"/>
<taskdef resource="org/apache/maven/artifact/ant/antlib.xml"
uri="antlib:antlib:org.apache.maven.artifact.ant"
classpath="${ivy.jar.dir}/${maven-ant-tasks.jar.file}"/>
<!-- End network lines -->
<ivy:settings file="${ivy.settings.dir}/ivysettings.xml"/>
<property name="init.resolve.done" value="true"/>
</target>
@ -295,7 +309,7 @@
<target name="git.rev-parse" depends="git.describe" unless="git.describe.succeeded">
<exec executable="git" outputproperty="git.rev-parse.output" resultproperty="git.rev-parse.exit.value" failonerror="false">
<arg line="rev-parse HEAD" />
<arg line="rev-parse --short HEAD" />
</exec>
<condition property="git.rev-parse.succeeded">
<equals arg1="${git.rev-parse.exit.value}" arg2="0" />
@ -942,6 +956,28 @@
</ant>
</target>
<!-- Maven install a package consisting of all supporting files. Don't call this target directly. Call one of the specific packaging targets below -->
<target name="mvninstall" depends="package" description="maven install a package into .m2/repository">
<property name="mvn.build.version" value="0.0.1" />
<!--
We should use the build version or better yet a git tag version, but tags are currently missing. Alternatively how do we then depend on the LATEST?
<property name="mvn.build.version" value="${build.version}" />
-->
<artifact:pom id="${executable}.pom" groupId="org.broadinstitute.sting" artifactId="${executable}" version="${mvn.build.version}" name="${executable}" />
<artifact:writepom pomRefId="${executable}.pom" file="${package.output.dir}/${executable}-${build.version}/${executable}.pom.xml"/>
<artifact:install file="${package.output.dir}/${executable}-${build.version}/${executable}.jar">
<artifact:pom file="${package.output.dir}/${executable}-${build.version}/${executable}.pom.xml" />
</artifact:install>
</target>
<!-- Maven install specific versions of the GATK/Queue. ALWAYS do an ant clean before invoking these! -->
<target name="mvninstall.gatk.full" depends="package.gatk.full,mvninstall" />
<target name="mvninstall.gatk.lite" depends="package.gatk.lite,mvninstall" />
<target name="mvninstall.queue.full" depends="package.queue.full,mvninstall" />
<target name="mvninstall.queue.lite" depends="package.queue.lite,mvninstall" />
<!-- ******************************************************************************** -->
<!-- Clean -->

BIN
intellij_example.tar.bz2 100644
View File

Binary file not shown.

7
ivy.xml
View File

@ -46,7 +46,8 @@
<dependency org="org.apache.bcel" name="bcel" rev="5.2"/>
<!-- Dependencies for reflections mvn repository -->
<dependency org="org.reflections" name="reflections" rev="0.9.5-RC2"/>
<dependency org="org.reflections" name="reflections" rev="0.9.8"/>
<dependency org="org.slf4j" name="slf4j-log4j12" rev="1.6.1"/>
<!-- Matrix package from math.nist.gov -->
<dependency org="gov.nist" name="Jama" rev="1.0.2"/>
@ -78,8 +79,8 @@
<dependency org="net.sf.gridscheduler" name="drmaa" rev="latest.integration"/>
<!-- Scala dependancies -->
<dependency org="org.scala-lang" name="scala-compiler" rev="2.8.1"/>
<dependency org="org.scala-lang" name="scala-library" rev="2.8.1"/>
<dependency org="org.scala-lang" name="scala-compiler" rev="2.9.2"/>
<dependency org="org.scala-lang" name="scala-library" rev="2.9.2"/>
<!-- testing and evaluation dependencies -->
<dependency org="org.testng" name="testng" rev="5.14.1" conf="test"/>

17
protected/java/src/org/broadinstitute/sting/gatk/walkers/compression/reducereads/ReduceReads.java
View File

@ -181,15 +181,6 @@ public class ReduceReads extends ReadWalker<LinkedList<GATKSAMRecord>, ReduceRea
@Argument(fullName = "minimum_del_proportion_to_trigger_variant", shortName = "mindel", doc = "", required = false)
private double minIndelProportionToTriggerVariant = 0.05;
/**
* Minimum proportion of indels in a site to trigger a variant region. Anything below this will be
* considered consensus.
*/
@Argument(fullName = "contigs", shortName = "ctg", doc = "", required = false)
private int nContigs = 2;
/**
* Downsamples the coverage of a variable region approximately (guarantees the minimum to be equal to this).
* A value of 0 turns downsampling off.
@ -197,6 +188,14 @@ public class ReduceReads extends ReadWalker<LinkedList<GATKSAMRecord>, ReduceRea
@Argument(fullName = "downsample_coverage", shortName = "ds", doc = "", required = false)
private int downsampleCoverage = 250;
/**
* Number of chromossomes in the sample (this is used for the polyploid consensus compression). Only
* tested for humans (or organisms with n=2). Use at your own risk!
*/
@Hidden
@Argument(fullName = "contigs", shortName = "ctg", doc = "", required = false)
private int nContigs = 2;
@Hidden
@Argument(fullName = "", shortName = "dl", doc = "", required = false)
private int debugLevel = 0;

56
protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcPerformanceTest.java
View File

@ -5,23 +5,22 @@ import org.apache.log4j.Logger;
import org.apache.log4j.TTCCLayout;
import org.broadinstitute.sting.gatk.report.GATKReport;
import org.broadinstitute.sting.gatk.report.GATKReportTable;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.SimpleTimer;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.fasta.CachingIndexedFastaSequenceFile;
import org.broadinstitute.sting.utils.variantcontext.Allele;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.variantcontext.VariantContextBuilder;
import java.io.FileOutputStream;
import java.io.PrintStream;
import java.io.*;
import java.util.*;
/**
* Created with IntelliJ IDEA.
* User: depristo
* Date: 10/2/12
* Time: 10:25 AM
* To change this template use File | Settings | File Templates.
* A simple GATK utility (i.e, runs from command-line) for assessing the performance of
* the exact model
*/
public class AFCalcPerformanceTest {
final static Logger logger = Logger.getLogger(AFCalcPerformanceTest.class);
@ -190,7 +189,8 @@ public class AFCalcPerformanceTest {
public enum Operation {
ANALYZE,
SINGLE
SINGLE,
EXACT_LOG
}
public static void main(final String[] args) throws Exception {
final TTCCLayout layout = new TTCCLayout();
@ -204,10 +204,49 @@ public class AFCalcPerformanceTest {
switch ( op ) {
case ANALYZE: analyze(args); break;
case SINGLE: profileBig(args); break;
case EXACT_LOG: exactLog(args); break;
default: throw new IllegalAccessException("unknown operation " + op);
}
}
private static void exactLog(final String[] args) throws Exception {
final File ref = new File(args[1]);
final File exactLogFile = new File(args[2]);
final List<Integer> startsToUse = new LinkedList<Integer>();
for ( int i = 3; i < args.length; i++ )
startsToUse.add(Integer.valueOf(args[i]));
final CachingIndexedFastaSequenceFile seq = new CachingIndexedFastaSequenceFile(ref);
final GenomeLocParser parser = new GenomeLocParser(seq);
final BufferedReader reader = new BufferedReader(new FileReader(exactLogFile));
final List<ExactCallLogger.ExactCall> loggedCalls = ExactCallLogger.readExactLog(reader, startsToUse, parser);
for ( final ExactCallLogger.ExactCall call : loggedCalls ) {
final AFCalcTestBuilder testBuilder = new AFCalcTestBuilder(call.vc.getNSamples(), 1,
AFCalcFactory.Calculation.EXACT_INDEPENDENT,
AFCalcTestBuilder.PriorType.human);
logger.info(call);
final SimpleTimer timer = new SimpleTimer().start();
final AFCalcResult result = testBuilder.makeModel().getLog10PNonRef(call.vc, testBuilder.makePriors());
final long newNanoTime = timer.getElapsedTimeNano();
if ( call.originalCall.anyPolymorphic(-1) || result.anyPolymorphic(-1) ) {
logger.info("**** ONE IS POLY");
}
logger.info("\t\t getLog10PosteriorOfAFGT0: " + call.originalCall.getLog10PosteriorOfAFGT0() + " vs " + result.getLog10PosteriorOfAFGT0());
final double speedup = call.runtime / (1.0 * newNanoTime);
logger.info("\t\t runtime: " + call.runtime + " vs " + newNanoTime + " speedup " + String.format("%.2f", speedup) + "x");
for ( final Allele a : call.originalCall.getAllelesUsedInGenotyping() ) {
if ( a.isNonReference() ) {
final String warningmeMLE = call.originalCall.getAlleleCountAtMLE(a) != result.getAlleleCountAtMLE(a) ? " DANGER-MLE-DIFFERENT" : "";
logger.info("\t\t MLE " + a + ": " + call.originalCall.getAlleleCountAtMLE(a) + " vs " + result.getAlleleCountAtMLE(a) + warningmeMLE);
final String warningmePost = call.originalCall.getLog10PosteriorOfAFGt0ForAllele(a) == 0 && result.getLog10PosteriorOfAFGt0ForAllele(a) < -10 ? " DANGER-POSTERIORS-DIFFERENT" : "";
logger.info("\t\t Posterior " + a + ": " + call.originalCall.getLog10PosteriorOfAFGt0ForAllele(a) + " vs " + result.getLog10PosteriorOfAFGt0ForAllele(a) + warningmePost);
}
}
}
}
private static void profileBig(final String[] args) throws Exception {
final int nSamples = Integer.valueOf(args[1]);
final int ac = Integer.valueOf(args[2]);
@ -234,7 +273,6 @@ public class AFCalcPerformanceTest {
final List<ModelParams> modelParams = Arrays.asList(
new ModelParams(AFCalcFactory.Calculation.EXACT_REFERENCE, 10000, 10),
// new ModelParams(AFCalcTestBuilder.ModelType.GeneralExact, 100, 10),
new ModelParams(AFCalcFactory.Calculation.EXACT_CONSTRAINED, 10000, 100),
new ModelParams(AFCalcFactory.Calculation.EXACT_INDEPENDENT, 10000, 1000));
final boolean ONLY_HUMAN_PRIORS = false;

6
protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcTestBuilder.java
View File

@ -45,12 +45,16 @@ public class AFCalcTestBuilder {
human
}
public int getNumAltAlleles() {
return numAltAlleles;
}
public int getnSamples() {
return nSamples;
}
public AFCalc makeModel() {
return AFCalcFactory.createAFCalc(modelType, nSamples, 4, 4, 2);
return AFCalcFactory.createAFCalc(modelType, nSamples, getNumAltAlleles(), getNumAltAlleles(), 2);
}
public double[] makePriors() {

4
protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/LikelihoodCalculationEngine.java
View File

@ -348,12 +348,14 @@ public class LikelihoodCalculationEngine {
}
}
}
// add all filtered reads to the NO_CALL list because they weren't given any likelihoods
for( final GATKSAMRecord read : perSampleFilteredReadList.get(sample.getKey()) ) {
// only count the read if it overlaps the event, otherwise it is not added to the output read list at all
if( callLoc.overlapsP(parser.createGenomeLoc(read)) ) {
for( final Allele a : call.getFirst().getAlleles() )
for( final Allele a : call.getFirst().getAlleles() ) {
likelihoodMap.add(read, a, 0.0);
}
}
}

8
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperGeneralPloidyIntegrationTest.java
View File

@ -60,22 +60,22 @@ public class UnifiedGenotyperGeneralPloidyIntegrationTest extends WalkerTest {
@Test(enabled = true)
public void testBOTH_GGA_Pools() {
PC_LSV_Test(String.format(" -maxAltAlleles 2 -ploidy 24 -gt_mode GENOTYPE_GIVEN_ALLELES -out_mode EMIT_ALL_SITES -alleles %s",LSV_ALLELES),"LSV_BOTH_GGA","BOTH","6d60d9f3dfe8e1580214be0d170b0fff");
PC_LSV_Test(String.format(" -maxAltAlleles 2 -ploidy 24 -gt_mode GENOTYPE_GIVEN_ALLELES -out_mode EMIT_ALL_SITES -alleles %s",LSV_ALLELES),"LSV_BOTH_GGA","BOTH","67dabdbf1e6ed8a83d2e85766558a20a");
}
@Test(enabled = true)
public void testINDEL_GGA_Pools() {
PC_LSV_Test(String.format(" -maxAltAlleles 1 -ploidy 24 -gt_mode GENOTYPE_GIVEN_ALLELES -out_mode EMIT_ALL_SITES -alleles %s",LSV_ALLELES),"LSV_INDEL_GGA","INDEL","30abf3c1868a61145edbc08fe35c8150");
PC_LSV_Test(String.format(" -maxAltAlleles 1 -ploidy 24 -gt_mode GENOTYPE_GIVEN_ALLELES -out_mode EMIT_ALL_SITES -alleles %s",LSV_ALLELES),"LSV_INDEL_GGA","INDEL","d4bfae27f1b07923f381d708d8a34cf4");
}
@Test(enabled = true)
public void testINDEL_maxAltAlleles2_ploidy3_Pools_noRef() {
PC_LSV_Test_NoRef(" -maxAltAlleles 2 -ploidy 3","LSV_INDEL_DISC_NOREF_p3","INDEL","ef99bc0513d3267f43b84cb88a324376");
PC_LSV_Test_NoRef(" -maxAltAlleles 2 -ploidy 3","LSV_INDEL_DISC_NOREF_p3","INDEL","7d6f319b9edcb1ff8c290fef150a2df8");
}
@Test(enabled = true)
public void testINDEL_maxAltAlleles2_ploidy1_Pools_noRef() {
PC_LSV_Test_NoRef(" -maxAltAlleles 2 -ploidy 1","LSV_INDEL_DISC_NOREF_p1","INDEL","8ca07270717641385fe5d2e07e530782");
PC_LSV_Test_NoRef(" -maxAltAlleles 2 -ploidy 1","LSV_INDEL_DISC_NOREF_p1","INDEL","dd02890123e07e7412a49475cb6280f1");
}
@Test(enabled = true)

87
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcPerformanceUnitTest.java 100644
View File

@ -0,0 +1,87 @@
package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc;
import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.variantcontext.Allele;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.testng.Assert;
import org.testng.annotations.DataProvider;
import org.testng.annotations.Test;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
public class AFCalcPerformanceUnitTest extends BaseTest {
@DataProvider(name = "ScalingTests")
public Object[][] makepolyTestProviderLotsOfAlleles() {
List<Object[]> tests = new ArrayList<Object[]>();
// list of all high-quality models in the system
final List<AFCalcFactory.Calculation> biAllelicModels = Arrays.asList(
AFCalcFactory.Calculation.EXACT_INDEPENDENT,
AFCalcFactory.Calculation.EXACT_REFERENCE);
final List<AFCalcFactory.Calculation> multiAllelicModels = Arrays.asList(
AFCalcFactory.Calculation.EXACT_INDEPENDENT);
// for ( final int nonTypePLs : Arrays.asList(100) ) {
// for ( final int nSamples : Arrays.asList(10000) ) {
// final List<Integer> alleleCounts = Arrays.asList(50);
// for ( final int nAltAlleles : Arrays.asList(1) ) {
for ( final int nonTypePLs : Arrays.asList(100) ) {
for ( final int nSamples : Arrays.asList(100, 1000) ) {
final List<Integer> alleleCounts = Arrays.asList(0, 1, 2, 3, 4, 5, 10, 50, 500);
for ( final int nAltAlleles : Arrays.asList(1, 2, 3) ) {
final List<AFCalcFactory.Calculation> models = nAltAlleles > 1 ? multiAllelicModels : biAllelicModels;
for ( final AFCalcFactory.Calculation model : models ) {
for ( final List<Integer> ACs : Utils.makePermutations(alleleCounts, nAltAlleles, true) ) {
if ( MathUtils.sum(ACs) < nSamples * 2 ) {
final AFCalcTestBuilder testBuilder
= new AFCalcTestBuilder(nSamples, nAltAlleles, model, AFCalcTestBuilder.PriorType.human);
tests.add(new Object[]{testBuilder, ACs, nonTypePLs});
}
}
}
}
}
}
return tests.toArray(new Object[][]{});
}
private Pair<Integer, Integer> estNumberOfEvaluations(final AFCalcTestBuilder testBuilder, final VariantContext vc, final int nonTypePL) {
final int evalOverhead = 2; // 2
final int maxEvalsPerSamplePerAC = 3;
int minEvals = 0, maxEvals = 0;
for ( final Allele alt : vc.getAlternateAlleles() ) {
final int AC = vc.getCalledChrCount(alt);
minEvals += AC + evalOverhead; // everyone is hom-var
maxEvals += AC * maxEvalsPerSamplePerAC + 10;
}
return new Pair<Integer, Integer>(minEvals, maxEvals);
}
@Test(dataProvider = "ScalingTests")
private void testScaling(final AFCalcTestBuilder testBuilder, final List<Integer> ACs, final int nonTypePL) {
final AFCalc calc = testBuilder.makeModel();
final double[] priors = testBuilder.makePriors();
final VariantContext vc = testBuilder.makeACTest(ACs, 0, nonTypePL);
final AFCalcResult result = calc.getLog10PNonRef(vc, priors);
final Pair<Integer, Integer> expectedNEvaluation = estNumberOfEvaluations(testBuilder, vc, nonTypePL);
final int minEvals = expectedNEvaluation.getFirst();
final int maxEvals = expectedNEvaluation.getSecond();
logger.warn(" min " + minEvals + " obs " + result.getnEvaluations() + " max " + maxEvals + " for test " + testBuilder + " sum(ACs)=" + (int)MathUtils.sum(ACs));
Assert.assertTrue(result.getnEvaluations() >= minEvals,
"Actual number of evaluations " + result.getnEvaluations() + " < min number of evals " + minEvals);
Assert.assertTrue(result.getnEvaluations() <= maxEvals,
"Actual number of evaluations " + result.getnEvaluations() + " > max number of evals " + minEvals);
}
}

77
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResultUnitTest.java 100644
View File

@ -0,0 +1,77 @@
package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc;
import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.variantcontext.Allele;
import org.testng.Assert;
import org.testng.annotations.DataProvider;
import org.testng.annotations.Test;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
public class AFCalcResultUnitTest extends BaseTest {
private static class MyTest {
final double[] Ls, expectedPosteriors;
private MyTest(double[] ls, double[] expectedPosteriors) {
Ls = ls;
this.expectedPosteriors = expectedPosteriors;
}
@Override
public String toString() {
return "Ls [" + Utils.join(",", Ls) + "] expectedPosteriors [" + Utils.join(",", expectedPosteriors) + "]";
}
}
@DataProvider(name = "TestComputePosteriors")
public Object[][] makeTestCombineGLs() {
List<Object[]> tests = new ArrayList<Object[]>();
tests.add(new Object[]{new MyTest(log10Even, log10Even)});
for ( double L0 = -1e9; L0 < 0.0; L0 /= 10.0 ) {
for ( double L1 = -1e2; L1 < 0.0; L1 /= 100.0 ) {
final double[] input = new double[]{L0, L1};
final double[] expected = MathUtils.normalizeFromLog10(input, true);
tests.add(new Object[]{new MyTest(input, expected)});
}
}
for ( double bigBadL = -1e50; bigBadL < -1e200; bigBadL *= 10 ) {
// test that a huge bad likelihood remains, even with a massive better result
for ( final double betterL : Arrays.asList(-1000.0, -100.0, -10.0, -1.0, -0.1, -0.01, -0.001, 0.0)) {
tests.add(new Object[]{new MyTest(new double[]{bigBadL, betterL}, new double[]{bigBadL, 0.0})});
tests.add(new Object[]{new MyTest(new double[]{betterL, bigBadL}, new double[]{0.0, bigBadL})});
}
}
// test that a modest bad likelihood with an ~0.0 value doesn't get lost
for ( final double badL : Arrays.asList(-10000.0, -1000.0, -100.0, -10.0)) {
tests.add(new Object[]{new MyTest(new double[]{badL, -1e-9}, new double[]{badL, 0.0})});
tests.add(new Object[]{new MyTest(new double[]{-1e-9, badL}, new double[]{0.0, badL})});
}
return tests.toArray(new Object[][]{});
}
final static double[] log10Even = MathUtils.normalizeFromLog10(new double[]{0.5, 0.5}, true);
final static Allele C = Allele.create("C");
final static List<Allele> alleles = Arrays.asList(Allele.create("A", true), C);
@Test(enabled = true, dataProvider = "TestComputePosteriors")
private void testComputingPosteriors(final MyTest data) {
final AFCalcResult result = new AFCalcResult(new int[]{0}, 1, alleles, data.Ls, log10Even, Collections.singletonMap(C, -1.0));
Assert.assertEquals(result.getLog10PosteriorOfAFEq0(), data.expectedPosteriors[0], 1e-3, "AF = 0 not expected");
Assert.assertEquals(result.getLog10PosteriorOfAFGT0(), data.expectedPosteriors[1], 1e-3, "AF > 0 not expected");
final double[] actualPosteriors = new double[]{result.getLog10PosteriorOfAFEq0(), result.getLog10PosteriorOfAFGT0()};
Assert.assertEquals(MathUtils.sumLog10(actualPosteriors), 1.0, 1e-3, "Posteriors don't sum to 1 with 1e-3 precision");
}
}

45
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcUnitTest.java
View File

@ -372,16 +372,14 @@ public class AFCalcUnitTest extends BaseTest {
final VariantContext vc3 = new VariantContextBuilder("x","1", 1, 1, Arrays.asList(A, C, G)).make();
final AFCalcTestBuilder.PriorType priorType = AFCalcTestBuilder.PriorType.flat;
final List<AFCalcFactory.Calculation> constrainedModel = Arrays.asList(AFCalcFactory.Calculation.EXACT_CONSTRAINED);
final double TOLERANCE = 0.5;
final List<PNonRefData> initialPNonRefData = Arrays.asList(
// bi-allelic sites
new PNonRefData(vc2, makePL(AA, 0, 10, 10), 0.1666667, TOLERANCE, true),
new PNonRefData(vc2, makePL(AA, 0, 1, 10), 0.4721084, TOLERANCE, false, constrainedModel),
new PNonRefData(vc2, makePL(AA, 0, 1, 1), 0.6136992, TOLERANCE, false, constrainedModel),
new PNonRefData(vc2, makePL(AA, 0, 5, 5), 0.3874259, TOLERANCE, false, constrainedModel),
new PNonRefData(vc2, makePL(AA, 0, 1, 10), 0.4721084, TOLERANCE, false),
new PNonRefData(vc2, makePL(AA, 0, 1, 1), 0.6136992, TOLERANCE, false),
new PNonRefData(vc2, makePL(AA, 0, 5, 5), 0.3874259, TOLERANCE, false),
new PNonRefData(vc2, makePL(AC, 10, 0, 10), 0.9166667, TOLERANCE, true),
new PNonRefData(vc2, makePL(CC, 10, 10, 0), 0.9166667, TOLERANCE, true),
@ -448,7 +446,7 @@ public class AFCalcUnitTest extends BaseTest {
return tests.toArray(new Object[][]{});
}
@Test(enabled = true & ! DEBUG_ONLY, dataProvider = "Models")
@Test(enabled = true && !DEBUG_ONLY, dataProvider = "Models")
public void testBiallelicPriors(final AFCalc model) {
for ( int REF_PL = 10; REF_PL <= 20; REF_PL += 10 ) {
@ -456,26 +454,29 @@ public class AFCalcUnitTest extends BaseTest {
for ( int log10NonRefPrior = 1; log10NonRefPrior < 10*REF_PL; log10NonRefPrior += 1 ) {
final double refPrior = 1 - QualityUtils.qualToErrorProb(log10NonRefPrior);
final double[] priors = MathUtils.normalizeFromLog10(MathUtils.toLog10(new double[]{refPrior, (1-refPrior) / 2, (1-refPrior) / 2}), true);
GetGLsTest cfg = new GetGLsTest(model, 1, Arrays.asList(AB), priors, "pNonRef" + log10NonRefPrior);
final AFCalcResult resultTracker = cfg.execute();
final int actualAC = resultTracker.getAlleleCountsOfMLE()[0];
final double nonRefPrior = (1-refPrior) / 2;
final double[] priors = MathUtils.normalizeFromLog10(MathUtils.toLog10(new double[]{refPrior, nonRefPrior, nonRefPrior}), true);
if ( ! Double.isInfinite(priors[1]) ) {
GetGLsTest cfg = new GetGLsTest(model, 1, Arrays.asList(AB), priors, "pNonRef" + log10NonRefPrior);
final AFCalcResult resultTracker = cfg.execute();
final int actualAC = resultTracker.getAlleleCountsOfMLE()[0];
final double pRefWithPrior = AB.getLikelihoods().getAsVector()[0] + priors[0];
final double pHetWithPrior = AB.getLikelihoods().getAsVector()[1] + priors[1] - Math.log10(0.5);
final double nonRefPost = Math.pow(10, pHetWithPrior) / (Math.pow(10, pRefWithPrior) + Math.pow(10, pHetWithPrior));
final double log10NonRefPost = Math.log10(nonRefPost);
final double pRefWithPrior = AB.getLikelihoods().getAsVector()[0] + priors[0];
final double pHetWithPrior = AB.getLikelihoods().getAsVector()[1] + priors[1] - Math.log10(0.5);
final double nonRefPost = Math.pow(10, pHetWithPrior) / (Math.pow(10, pRefWithPrior) + Math.pow(10, pHetWithPrior));
final double log10NonRefPost = Math.log10(nonRefPost);
if ( ! Double.isInfinite(log10NonRefPost) )
Assert.assertEquals(resultTracker.getLog10PosteriorOfAFGT0(), log10NonRefPost, 1e-2);
if ( ! Double.isInfinite(log10NonRefPost) )
Assert.assertEquals(resultTracker.getLog10PosteriorOfAFGT0(), log10NonRefPost, 1e-2);
if ( nonRefPost >= 0.9 )
Assert.assertTrue(resultTracker.isPolymorphic(C, -1));
if ( nonRefPost >= 0.9 )
Assert.assertTrue(resultTracker.isPolymorphic(C, -1));
final int expectedMLEAC = 1; // the MLE is independent of the prior
Assert.assertEquals(actualAC, expectedMLEAC,
"actual AC with priors " + log10NonRefPrior + " not expected "
+ expectedMLEAC + " priors " + Utils.join(",", priors));
final int expectedMLEAC = 1; // the MLE is independent of the prior
Assert.assertEquals(actualAC, expectedMLEAC,
"actual AC with priors " + log10NonRefPrior + " not expected "
+ expectedMLEAC + " priors " + Utils.join(",", priors));
}
}
}
}

124
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ConstrainedAFCalculationModelUnitTest.java
View File

@ -1,124 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc;
import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.utils.variantcontext.Allele;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.variantcontext.VariantContextBuilder;
import org.testng.Assert;
import org.testng.annotations.DataProvider;
import org.testng.annotations.Test;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
public class ConstrainedAFCalculationModelUnitTest extends BaseTest {
static Allele A = Allele.create("A", true);
static Allele C = Allele.create("C");
static Allele G = Allele.create("G");
protected static Genotype makePL(final List<Allele> expectedGT, int ... pls) {
return AFCalcUnitTest.makePL(expectedGT, pls);
}
@DataProvider(name = "MaxACsToVisit")
public Object[][] makeMaxACsToVisit() {
List<Object[]> tests = new ArrayList<Object[]>();
final int nSamples = 10;
for (int nNonInformative = 0; nNonInformative < nSamples - 1; nNonInformative++ ) {
final int nChrom = (nSamples - nNonInformative) * 2;
for ( int i = 0; i < nChrom; i++ ) {
// bi-allelic
tests.add(new Object[]{nSamples, Arrays.asList(i), nNonInformative, AFCalcFactory.Calculation.EXACT_CONSTRAINED});
// tri-allelic
for ( int j = 0; j < (nChrom - i); j++)
tests.add(new Object[]{nSamples, Arrays.asList(i, j), nNonInformative, AFCalcFactory.Calculation.EXACT_CONSTRAINED});
}
}
return tests.toArray(new Object[][]{});
}
@Test(enabled = true, dataProvider = "MaxACsToVisit")
public void testMaxACsToVisit(final int nSamples, final List<Integer> requestedACs, final int nNonInformative, final AFCalcFactory.Calculation modelType) {
final int nAlts = requestedACs.size();
final AFCalcTestBuilder testBuilder
= new AFCalcTestBuilder(nSamples, nAlts, modelType,
AFCalcTestBuilder.PriorType.human);
final VariantContext vc = testBuilder.makeACTest(requestedACs, nNonInformative, 100);
final int[] maxACsToVisit = ((ConstrainedDiploidExactAFCalc)testBuilder.makeModel()).computeMaxACs(vc);
testExpectedACs(vc, maxACsToVisit);
}
private void testExpectedACs(final VariantContext vc, final int[] maxACsToVisit) {
// this is necessary because cannot ensure that the tester gives us back the
// requested ACs due to rounding errors
final List<Integer> ACs = new ArrayList<Integer>();
for ( final Allele a : vc.getAlternateAlleles() )
ACs.add(vc.getCalledChrCount(a));
for ( int i = 0; i < maxACsToVisit.length; i++ ) {
Assert.assertEquals(maxACsToVisit[i], (int)ACs.get(i), "Maximum AC computed wasn't equal to the max possible in the construction for alt allele " + i);
}
}
@DataProvider(name = "MaxACsGenotypes")
public Object[][] makeMaxACsForGenotype() {
List<Object[]> tests = new ArrayList<Object[]>();
final List<Allele> AA = Arrays.asList(A, A);
final List<Allele> AC = Arrays.asList(A, C);
final List<Allele> CC = Arrays.asList(C, C);
final List<Allele> AG = Arrays.asList(A, G);
final List<Allele> GG = Arrays.asList(G, G);
final List<Allele> CG = Arrays.asList(C, G);
final VariantContext vc2 = new VariantContextBuilder("x","1", 1, 1, Arrays.asList(A, C)).make();
final VariantContext vc3 = new VariantContextBuilder("x","1", 1, 1, Arrays.asList(A, C, G)).make();
tests.add(new Object[]{vc2, makePL(AA, 0, 10, 10)});
tests.add(new Object[]{vc2, makePL(AC, 10, 0, 10)});
tests.add(new Object[]{vc2, makePL(CC, 10, 10, 0)});
// make sure non-informative => 0
tests.add(new Object[]{vc2, makePL(AA, 0, 0, 0)});
tests.add(new Object[]{vc3, makePL(AA, 0, 0, 0, 0, 0, 0)});
// multi-allelics
tests.add(new Object[]{vc3, makePL(AG, 10, 10, 10, 0, 10, 10)});
tests.add(new Object[]{vc3, makePL(CG, 10, 10, 10, 10, 0, 10)});
tests.add(new Object[]{vc3, makePL(GG, 10, 10, 10, 10, 10, 0)});
// deal with non-informatives third alleles
tests.add(new Object[]{vc3, makePL(AC, 10, 0, 10, 0, 0, 10)});
tests.add(new Object[]{vc3, makePL(AC, 10, 0, 10, 10, 0, 10)});
tests.add(new Object[]{vc3, makePL(AC, 10, 0, 10, 10, 0, 0)});
tests.add(new Object[]{vc3, makePL(AC, 10, 0, 10, 0, 0, 0)});
tests.add(new Object[]{vc3, makePL(CC, 10, 10, 0, 0, 0, 10)});
tests.add(new Object[]{vc3, makePL(CC, 10, 10, 0, 10, 0, 10)});
tests.add(new Object[]{vc3, makePL(CC, 10, 10, 0, 10, 0, 0)});
tests.add(new Object[]{vc3, makePL(CC, 10, 10, 0, 0, 0, 0)});
return tests.toArray(new Object[][]{});
}
@Test(enabled = true, dataProvider = "MaxACsGenotypes")
private void testMakeACByGenotype(final VariantContext vcRoot, final Genotype g) {
final VariantContext vc = new VariantContextBuilder(vcRoot).genotypes(g).make();
final AFCalcTestBuilder testBuilder
= new AFCalcTestBuilder(1, vc.getNAlleles()-1, AFCalcFactory.Calculation.EXACT_CONSTRAINED,
AFCalcTestBuilder.PriorType.human);
final int[] maxACsToVisit = ((ConstrainedDiploidExactAFCalc)testBuilder.makeModel()).computeMaxACs(vc);
testExpectedACs(vc, maxACsToVisit);
}
}

48
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalcUnitTest.java
View File

@ -55,48 +55,14 @@ public class IndependentAllelesDiploidExactAFCalcUnitTest extends BaseTest {
return tests.toArray(new Object[][]{});
}
@DataProvider(name = "TestCombineGLsWithDrops")
public Object[][] makeTestCombineGLsWithDrops() {
List<Object[]> tests = new ArrayList<Object[]>();
final Set<Integer> noDrops = Collections.emptySet();
final Set<Integer> drop1 = Collections.singleton(1);
final Set<Integer> drop2 = Collections.singleton(2);
// AA AB BB AC BC CC
// drop1 (B): AA AC CC
// drop2 (C): AA AB BB
tests.add(new Object[]{1, 2, makePL( 0, 1, 2, 3, 4, 5), makePL(0, 2, 5), noDrops});
tests.add(new Object[]{2, 2, makePL( 0, 1, 2, 3, 4, 5), makePL(0, 4, 9), noDrops});
tests.add(new Object[]{1, 2, makePL( 0, 1, 2, 3, 4, 5), makePL(0, 1, 2), drop2});
tests.add(new Object[]{2, 2, makePL( 0, 1, 2, 3, 4, 5), makePL(0, 3, 5), drop1});
tests.add(new Object[]{1, 2, makePL( 5, 4, 3, 2, 1, 0), makePL(0, 2, 6), noDrops});
tests.add(new Object[]{2, 2, makePL( 5, 4, 3, 2, 1, 0), makePL(1, 0, 2), noDrops});
tests.add(new Object[]{1, 2, makePL( 5, 4, 3, 2, 1, 0), makePL(2, 1, 0), drop2});
tests.add(new Object[]{2, 2, makePL( 5, 4, 3, 2, 1, 0), makePL(5, 2, 0), drop1});
tests.add(new Object[]{1, 2, makePL(10,10,10,10,10, 0), makePL( 0, 8,11), noDrops});
tests.add(new Object[]{2, 2, makePL(10,10,10,10,10, 0), makePL( 5, 7, 0), noDrops});
tests.add(new Object[]{1, 2, makePL(10,10,10,10,10, 0), makePL( 0, 0, 0), drop2});
tests.add(new Object[]{2, 2, makePL(10,10,10,10,10, 0), makePL(10,10, 0), drop1});
return tests.toArray(new Object[][]{});
}
private Genotype makePL(final int ... PLs) {
return AFCalcUnitTest.makePL(Arrays.asList(Allele.NO_CALL, Allele.NO_CALL), PLs);
}
@Test(enabled = true, dataProvider = "TestCombineGLs")
private void testCombineGLs(final int altIndex, final int nAlts, final Genotype testg, final Genotype expected) {
testCombineGLsWithDrops(altIndex, nAlts, testg, expected, Collections.<Integer>emptySet());
}
@Test(enabled = true, dataProvider = "TestCombineGLsWithDrops")
private void testCombineGLsWithDrops(final int altIndex, final int nAlts, final Genotype testg, final Genotype expected, Set<Integer> allelesToDrop) {
final IndependentAllelesDiploidExactAFCalc calc = (IndependentAllelesDiploidExactAFCalc)AFCalcFactory.createAFCalc(AFCalcFactory.Calculation.EXACT_INDEPENDENT, 1, 4);
final Genotype combined = calc.combineGLs(testg, altIndex, allelesToDrop, nAlts);
final Genotype combined = calc.combineGLs(testg, altIndex, nAlts);
Assert.assertEquals(combined.getPL(), expected.getPL(),
"Combined PLs " + Utils.join(",", combined.getPL()) + " != expected " + Utils.join(",", expected.getPL()));
@ -120,22 +86,21 @@ public class IndependentAllelesDiploidExactAFCalcUnitTest extends BaseTest {
final Genotype gACG = makePL( 0, 1, 2, 3, 4, 5);
final Genotype gAGC = makePL( 0, 4, 5, 1, 3, 2);
final Genotype gACcombined = makePL(0, 2, 5);
final Genotype gACcombined2 = makePL(0, 1, 4);
final Genotype gAGcombined = makePL(0, 4, 9);
final Genotype gACdropped = makePL(0, 1, 2);
final Genotype gAGdropped = makePL(0, 3, 5);
// biallelic
tests.add(new Object[]{vcAC.genotypes(gACcombined).make(), Arrays.asList(vcAC.genotypes(gACcombined).make())});
// tri-allelic
tests.add(new Object[]{vcACG.genotypes(gACG).make(), Arrays.asList(vcAC.genotypes(gACcombined).make(), vcAG.genotypes(gAGdropped).make())});
tests.add(new Object[]{vcAGC.genotypes(gAGC).make(), Arrays.asList(vcAG.genotypes(gAGcombined).make(), vcAC.genotypes(gACdropped).make())});
tests.add(new Object[]{vcACG.genotypes(gACG).make(), Arrays.asList(vcAC.genotypes(gACcombined).make(), vcAG.genotypes(gAGcombined).make())});
tests.add(new Object[]{vcAGC.genotypes(gAGC).make(), Arrays.asList(vcAG.genotypes(gAGcombined).make(), vcAC.genotypes(gACcombined2).make())});
return tests.toArray(new Object[][]{});
}
@Test(enabled = false, dataProvider = "TestMakeAlleleConditionalContexts")
@Test(enabled = true, dataProvider = "TestMakeAlleleConditionalContexts")
private void testMakeAlleleConditionalContexts(final VariantContext vc, final List<VariantContext> expectedVCs) {
final IndependentAllelesDiploidExactAFCalc calc = (IndependentAllelesDiploidExactAFCalc)AFCalcFactory.createAFCalc(AFCalcFactory.Calculation.EXACT_INDEPENDENT, 1, 4);
final List<VariantContext> biAllelicVCs = calc.makeAlleleConditionalContexts(vc);
@ -148,7 +113,8 @@ public class IndependentAllelesDiploidExactAFCalcUnitTest extends BaseTest {
Assert.assertEquals(actual.getAlleles(), expected.getAlleles());
for ( int j = 0; j < actual.getNSamples(); j++ )
Assert.assertEquals(actual.getGenotype(j).getPL(), expected.getGenotype(j).getPL());
Assert.assertEquals(actual.getGenotype(j).getPL(), expected.getGenotype(j).getPL(),
"expected PLs " + Utils.join(",", expected.getGenotype(j).getPL()) + " not equal to actual " + Utils.join(",", actual.getGenotype(j).getPL()));
}
}

12
protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
View File

@ -21,17 +21,17 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
@Test
public void testHaplotypeCallerMultiSample() {
HCTest(CEUTRIO_BAM, "", "a305107a5ec889152aa2efbe90b249d7");
HCTest(CEUTRIO_BAM, "", "8c52c0955099cca3215a0d78fd455894");
}
@Test
public void testHaplotypeCallerSingleSample() {
HCTest(NA12878_BAM, "", "0c2217ec81f19790a6d1f98ebf8cf70d");
HCTest(NA12878_BAM, "", "01367428c26d3eaf9297c58bf8677dd3");
}
@Test
public void testHaplotypeCallerMultiSampleGGA() {
HCTest(CEUTRIO_BAM, "--max_alternate_alleles_for_indels 3 -gt_mode GENOTYPE_GIVEN_ALLELES -alleles " + validationDataLocation + "combined.phase1.chr20.raw.indels.sites.vcf", "0396c7352ab8ab98b03dca36299a0ddf");
HCTest(CEUTRIO_BAM, "--max_alternate_alleles_for_indels 3 -gt_mode GENOTYPE_GIVEN_ALLELES -alleles " + validationDataLocation + "combined.phase1.chr20.raw.indels.sites.vcf", "61c1a0fb62d909229af6b5a91dad8b35");
}
private void HCTestComplexVariants(String bam, String args, String md5) {
@ -42,7 +42,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
@Test
public void testHaplotypeCallerMultiSampleComplex() {
HCTestComplexVariants(CEUTRIO_BAM, "", "2cfb7d830d5a7eb7bc754b5f688a27a5");
HCTestComplexVariants(CEUTRIO_BAM, "", "30598abeeb0b0ae5816ffdbf0c4044fd");
}
private void HCTestSymbolicVariants(String bam, String args, String md5) {
@ -53,7 +53,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
@Test
public void testHaplotypeCallerSingleSampleSymbolic() {
HCTestSymbolicVariants(NA12878_CHR20_BAM, "", "d511848a46083c0d0b2495f65f162c2e");
HCTestSymbolicVariants(NA12878_CHR20_BAM, "", "6eb9c1026225b38ba7bd3c4c218f8269");
}
private void HCTestIndelQualityScores(String bam, String args, String md5) {
@ -70,7 +70,7 @@ public class HaplotypeCallerIntegrationTest extends WalkerTest {
@Test
public void HCTestProblematicReadsModifiedInActiveRegions() {
final String base = String.format("-T HaplotypeCaller -R %s -I %s", REF, privateTestDir + "haplotype-problem-4.bam") + " --no_cmdline_in_header -o %s -minPruning 3";
final WalkerTestSpec spec = new WalkerTestSpec(base, Arrays.asList("7e112ea4623617f1f7f8f562f54aa2aa"));
final WalkerTestSpec spec = new WalkerTestSpec(base, Arrays.asList("fa5c5eb996e95aed12c50d70e6dd74d7"));
executeTest("HCTestProblematicReadsModifiedInActiveRegions: ", spec);
}

28
public/java/src/org/broadinstitute/sting/commandline/ArgumentMatch.java
View File

@ -46,7 +46,7 @@ public class ArgumentMatch implements Iterable<ArgumentMatch> {
/**
* Maps indices of command line arguments to values paired with that argument.
*/
public final SortedMap<ArgumentMatchSite,List<String>> sites = new TreeMap<ArgumentMatchSite,List<String>>();
public final SortedMap<ArgumentMatchSite,List<ArgumentMatchValue>> sites = new TreeMap<ArgumentMatchSite,List<ArgumentMatchValue>>();
/**
* An ordered, freeform collection of tags.
@ -90,11 +90,11 @@ public class ArgumentMatch implements Iterable<ArgumentMatch> {
* @param value Value for the argument at this position.
* @param tags ordered freeform text tags associated with this argument.
*/
private ArgumentMatch(final String label, final ArgumentDefinition definition, final ArgumentMatchSite site, final String value, final Tags tags) {
private ArgumentMatch(final String label, final ArgumentDefinition definition, final ArgumentMatchSite site, final ArgumentMatchValue value, final Tags tags) {
this.label = label;
this.definition = definition;
ArrayList<String> values = new ArrayList<String>();
ArrayList<ArgumentMatchValue> values = new ArrayList<ArgumentMatchValue>();
if( value != null )
values.add(value);
sites.put(site,values );
@ -131,11 +131,11 @@ public class ArgumentMatch implements Iterable<ArgumentMatch> {
*/
@SuppressWarnings("unchecked")
ArgumentMatch transform(Multiplexer multiplexer, Object key) {
SortedMap<ArgumentMatchSite,List<String>> newIndices = new TreeMap<ArgumentMatchSite,List<String>>();
for(Map.Entry<ArgumentMatchSite,List<String>> site: sites.entrySet()) {
List<String> newEntries = new ArrayList<String>();
for(String entry: site.getValue())
newEntries.add(multiplexer.transformArgument(key,entry));
SortedMap<ArgumentMatchSite,List<ArgumentMatchValue>> newIndices = new TreeMap<ArgumentMatchSite,List<ArgumentMatchValue>>();
for(Map.Entry<ArgumentMatchSite,List<ArgumentMatchValue>> site: sites.entrySet()) {
List<ArgumentMatchValue> newEntries = new ArrayList<ArgumentMatchValue>();
for(ArgumentMatchValue entry: site.getValue())
newEntries.add(new ArgumentMatchStringValue(multiplexer.transformArgument(key,entry.asString())));
newIndices.put(site.getKey(),newEntries);
}
ArgumentMatch newArgumentMatch = new ArgumentMatch(label,definition);
@ -165,7 +165,7 @@ public class ArgumentMatch implements Iterable<ArgumentMatch> {
/**
* Iterate over each available token.
*/
private Iterator<String> tokenIterator = null;
private Iterator<ArgumentMatchValue> tokenIterator = null;
/**
* The next site to return. Null if none remain.
@ -175,7 +175,7 @@ public class ArgumentMatch implements Iterable<ArgumentMatch> {
/**
* The next token to return. Null if none remain.
*/
String nextToken = null;
ArgumentMatchValue nextToken = null;
{
siteIterator = sites.keySet().iterator();
@ -254,9 +254,9 @@ public class ArgumentMatch implements Iterable<ArgumentMatch> {
* @param site site of the command-line argument to which this value is mated.
* @param value Text representation of value to add.
*/
public void addValue( ArgumentMatchSite site, String value ) {
public void addValue( ArgumentMatchSite site, ArgumentMatchValue value ) {
if( !sites.containsKey(site) || sites.get(site) == null )
sites.put(site, new ArrayList<String>() );
sites.put(site, new ArrayList<ArgumentMatchValue>() );
sites.get(site).add(value);
}
@ -275,8 +275,8 @@ public class ArgumentMatch implements Iterable<ArgumentMatch> {
* Return the values associated with this argument match.
* @return A collection of the string representation of these value.
*/
public List<String> values() {
List<String> values = new ArrayList<String>();
public List<ArgumentMatchValue> values() {
List<ArgumentMatchValue> values = new ArrayList<ArgumentMatchValue>();
for( ArgumentMatchSite site: sites.keySet() ) {
if( sites.get(site) != null )
values.addAll(sites.get(site));

27
public/java/src/org/broadinstitute/sting/commandline/ArgumentMatchFileValue.java 100644
View File

@ -0,0 +1,27 @@
package org.broadinstitute.sting.commandline;
import java.io.File;
/**
* Holds a reference to a file as an argument match value.
*
* This is useful when the type of the stored file may be a subclass of java.io.File,
* for example a Queue RemoteFile.
*/
public class ArgumentMatchFileValue extends ArgumentMatchValue {
private final File file;
public ArgumentMatchFileValue(File file) {
this.file = file;
}
@Override
public String asString() {
return file == null ? null : file.getAbsolutePath();
}
@Override
public File asFile() {
return file;
}
}

42
public/java/src/org/broadinstitute/sting/commandline/ArgumentMatchSource.java
View File

@ -24,38 +24,36 @@
package org.broadinstitute.sting.commandline;
import java.io.File;
/**
* Where an argument match originated, via the commandline or a file.
* Where an argument match originated, via the commandline or a custom provider.
*/
public class ArgumentMatchSource implements Comparable<ArgumentMatchSource> {
public static final ArgumentMatchSource COMMAND_LINE = new ArgumentMatchSource(ArgumentMatchSourceType.CommandLine, null);
private final ArgumentMatchSourceType type;
private final File file;
private final String description;
/**
* Creates an argument match source from the specified file.
* @param file File specifying the arguments. Must not be null.
* @param description Where the arguments originated.
*/
public ArgumentMatchSource(File file) {
this(ArgumentMatchSourceType.File, file);
public ArgumentMatchSource(String description) {
this(ArgumentMatchSourceType.Provider, description);
}
private ArgumentMatchSource(ArgumentMatchSourceType type, File file) {
if (type == ArgumentMatchSourceType.File && file == null)
throw new IllegalArgumentException("An argument match source of type File cannot have a null file.");
private ArgumentMatchSource(ArgumentMatchSourceType type, String description) {
if (type == ArgumentMatchSourceType.Provider && description == null)
throw new IllegalArgumentException("An argument match source provider cannot have a null description.");
this.type = type;
this.file = file;
this.description = description;
}
public ArgumentMatchSourceType getType() {
return type;
}
public File getFile() {
return file;
public String getDescription() {
return description;
}
@Override
@ -65,13 +63,13 @@ public class ArgumentMatchSource implements Comparable<ArgumentMatchSource> {
ArgumentMatchSource that = (ArgumentMatchSource) o;
return (type == that.type) && (file == null ? that.file == null : file.equals(that.file));
return (type == that.type) && (description == null ? that.description == null : description.equals(that.description));
}
@Override
public int hashCode() {
int result = type != null ? type.hashCode() : 0;
result = 31 * result + (file != null ? file.hashCode() : 0);
result = 31 * result + (description != null ? description.hashCode() : 0);
return result;
}
@ -84,15 +82,15 @@ public class ArgumentMatchSource implements Comparable<ArgumentMatchSource> {
if (comp != 0)
return comp;
File f1 = this.file;
File f2 = that.file;
String d1 = this.description;
String d2 = that.description;
if ((f1 == null) ^ (f2 == null)) {
// If one of the files is null and the other is not
// put the null file first
return f1 == null ? -1 : 1;
if ((d1 == null) ^ (d2 == null)) {
// If one of the descriptions is null and the other is not
// put the null description first
return d1 == null ? -1 : 1;
}
return f1 == null ? 0 : f1.compareTo(f2);
return d1 == null ? 0 : d1.compareTo(d2);
}
}

4
public/java/src/org/broadinstitute/sting/commandline/ArgumentMatchSourceType.java
View File

@ -25,8 +25,8 @@
package org.broadinstitute.sting.commandline;
/**
* Type of where an argument match originated, via the commandline or a file.
* Type of where an argument match originated, via the commandline or a some other provider.
*/
public enum ArgumentMatchSourceType {
CommandLine, File
CommandLine, Provider
}

24
public/java/src/org/broadinstitute/sting/commandline/ArgumentMatchStringValue.java 100644
View File

@ -0,0 +1,24 @@
package org.broadinstitute.sting.commandline;
import java.io.File;
/**
* Argument values that originated from a string.
*/
public class ArgumentMatchStringValue extends ArgumentMatchValue {
private final String value;
public ArgumentMatchStringValue(String value) {
this.value = value;
}
@Override
public String asString() {
return value;
}
@Override
public File asFile() {
return value == null ? null : new File(value);
}
}

18
public/java/src/org/broadinstitute/sting/commandline/ArgumentMatchValue.java 100644
View File

@ -0,0 +1,18 @@
package org.broadinstitute.sting.commandline;
import java.io.File;
/**
* Returns argument values as either strings or values.
*/
public abstract class ArgumentMatchValue {
/**
* @return the value of this argument as a String object.
*/
public abstract String asString();
/**
* @return the value of this argument as a File object.
*/
public abstract File asFile();
}

34
public/java/src/org/broadinstitute/sting/commandline/ArgumentTypeDescriptor.java
View File

@ -215,8 +215,8 @@ public abstract class ArgumentTypeDescriptor {
* @param matches The matches for the given argument.
* @return The value of the argument if available, or null if not present.
*/
protected String getArgumentValue( ArgumentDefinition definition, ArgumentMatches matches ) {
Collection<String> argumentValues = getArgumentValues( definition, matches );
protected ArgumentMatchValue getArgumentValue( ArgumentDefinition definition, ArgumentMatches matches ) {
Collection<ArgumentMatchValue> argumentValues = getArgumentValues( definition, matches );
if( argumentValues.size() > 1 )
throw new UserException.CommandLineException("Multiple values associated with given definition, but this argument expects only one: " + definition.fullName);
return argumentValues.size() > 0 ? argumentValues.iterator().next() : null;
@ -244,8 +244,8 @@ public abstract class ArgumentTypeDescriptor {
* @param matches The matches for the given argument.
* @return The value of the argument if available, or an empty collection if not present.
*/
protected Collection<String> getArgumentValues( ArgumentDefinition definition, ArgumentMatches matches ) {
Collection<String> values = new ArrayList<String>();
protected Collection<ArgumentMatchValue> getArgumentValues( ArgumentDefinition definition, ArgumentMatches matches ) {
Collection<ArgumentMatchValue> values = new ArrayList<ArgumentMatchValue>();
for( ArgumentMatch match: matches ) {
if( match.definition.equals(definition) )
values.addAll(match.values());
@ -310,7 +310,7 @@ public abstract class ArgumentTypeDescriptor {
*/
protected Object parseBinding(ArgumentSource source, Type type, ArgumentMatches matches, Tags tags) {
ArgumentDefinition defaultDefinition = createDefaultArgumentDefinition(source);
String value = getArgumentValue(defaultDefinition, matches);
ArgumentMatchValue value = getArgumentValue(defaultDefinition, matches);
@SuppressWarnings("unchecked")
Class<? extends Feature> parameterType = JVMUtils.getParameterizedTypeClass(type);
String name = defaultDefinition.fullName;
@ -328,7 +328,7 @@ public abstract class ArgumentTypeDescriptor {
* @param fieldName The name of the field that was parsed. Used for error reporting.
* @return The newly created binding object of type bindingClass.
*/
public static Object parseBinding(String value, Class<? extends Feature> parameterType, Type bindingClass,
public static Object parseBinding(ArgumentMatchValue value, Class<? extends Feature> parameterType, Type bindingClass,
String bindingName, Tags tags, String fieldName) {
try {
String tribbleType = null;
@ -337,7 +337,7 @@ public abstract class ArgumentTypeDescriptor {
throw new UserException.CommandLineException(
String.format("Unexpected number of positional tags for argument %s : %s. " +
"Rod bindings only support -X:type and -X:name,type argument styles",
value, fieldName));
value.asString(), fieldName));
} else if ( tags.getPositionalTags().size() == 2 ) {
// -X:name,type style
bindingName = tags.getPositionalTags().get(0);
@ -366,7 +366,7 @@ public abstract class ArgumentTypeDescriptor {
if ( tribbleType == null ) {
// try to determine the file type dynamically
File file = new File(value);
File file = value.asFile();
if ( file.canRead() && file.isFile() ) {
FeatureManager.FeatureDescriptor featureDescriptor = manager.getByFiletype(file);
if ( featureDescriptor != null ) {
@ -379,7 +379,7 @@ public abstract class ArgumentTypeDescriptor {
// IntervalBinding can be created from a normal String
Class rawType = (makeRawTypeIfNecessary(bindingClass));
try {
return rawType.getConstructor(String.class).newInstance(value);
return rawType.getConstructor(String.class).newInstance(value.asString());
} catch (NoSuchMethodException e) {
/* ignore */
}
@ -399,14 +399,14 @@ public abstract class ArgumentTypeDescriptor {
}
Constructor ctor = (makeRawTypeIfNecessary(bindingClass)).getConstructor(Class.class, String.class, String.class, String.class, Tags.class);
return ctor.newInstance(parameterType, bindingName, value, tribbleType, tags);
return ctor.newInstance(parameterType, bindingName, value.asString(), tribbleType, tags);
} catch (Exception e) {
if ( e instanceof UserException )
throw ((UserException)e);
else
throw new UserException.CommandLineException(
String.format("Failed to parse value %s for argument %s. Message: %s",
value, fieldName, e.getMessage()));
value.asString(), fieldName, e.getMessage()));
}
}
}
@ -517,7 +517,7 @@ class SimpleArgumentTypeDescriptor extends ArgumentTypeDescriptor {
return true;
ArgumentDefinition defaultDefinition = createDefaultArgumentDefinition(source);
String value = getArgumentValue( defaultDefinition, matches );
ArgumentMatchValue value = getArgumentValue(defaultDefinition, matches);
Object result;
Tags tags = getArgumentTags(matches);
@ -527,12 +527,12 @@ class SimpleArgumentTypeDescriptor extends ArgumentTypeDescriptor {
Method valueOf = primitiveToWrapperMap.get(type).getMethod("valueOf",String.class);
if(value == null)
throw new MissingArgumentValueException(createDefaultArgumentDefinition(source));
result = valueOf.invoke(null,value.trim());
result = valueOf.invoke(null,value.asString().trim());
} else if (type.isEnum()) {
Object[] vals = type.getEnumConstants();
Object defaultEnumeration = null; // as we look at options, record the default option if it exists
for (Object val : vals) {
if (String.valueOf(val).equalsIgnoreCase(value)) return val;
if (String.valueOf(val).equalsIgnoreCase(value.asString())) return val;
try { if (type.getField(val.toString()).isAnnotationPresent(EnumerationArgumentDefault.class)) defaultEnumeration = val; }
catch (NoSuchFieldException e) { throw new ReviewedStingException("parsing " + type.toString() + "doesn't contain the field " + val.toString()); }
}
@ -544,10 +544,12 @@ class SimpleArgumentTypeDescriptor extends ArgumentTypeDescriptor {
else if (value == null)
throw new MissingArgumentValueException(createDefaultArgumentDefinition(source));
else
throw new UnknownEnumeratedValueException(createDefaultArgumentDefinition(source),value);
throw new UnknownEnumeratedValueException(createDefaultArgumentDefinition(source),value.asString());
} else if (type.equals(File.class)) {
result = value.asFile();
} else {
Constructor ctor = type.getConstructor(String.class);
result = ctor.newInstance(value);
result = ctor.newInstance(value.asString());
}
} catch (UserException e) {
throw e;

2
public/java/src/org/broadinstitute/sting/commandline/CommandLineProgram.java
View File

@ -174,7 +174,7 @@ public abstract class CommandLineProgram {
ParsingEngine parser = clp.parser = new ParsingEngine(clp);
parser.addArgumentSource(clp.getClass());
Map<ArgumentMatchSource, List<String>> parsedArgs;
Map<ArgumentMatchSource, ParsedArgs> parsedArgs;
// process the args
if (clp.canAddArgumentsDynamically()) {

13
public/java/src/org/broadinstitute/sting/commandline/ParsedArgs.java 100644
View File

@ -0,0 +1,13 @@
package org.broadinstitute.sting.commandline;
/**
* Represents a collection of parsed arguments for an argument source.
*
* Useful for printing out help documents.
*/
public abstract class ParsedArgs {
/**
* @return A compact description of the arguments from an provider/source.
*/
public abstract String getDescription();
}

30
public/java/src/org/broadinstitute/sting/commandline/ParsedListArgs.java 100644
View File

@ -0,0 +1,30 @@
package org.broadinstitute.sting.commandline;
import org.apache.commons.lang.StringUtils;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
/**
* A list of string arguments, usually from the command line or an args list file.
*/
public class ParsedListArgs extends ParsedArgs {
private final List<String> args = new ArrayList<String>();
public ParsedListArgs() {
}
public ParsedListArgs(List<String> args) {
this.args.addAll(args);
}
public void add(String... args) {
this.args.addAll(Arrays.asList(args));
}
@Override
public String getDescription() {
return StringUtils.join(this.args, " ");
}
}

102
public/java/src/org/broadinstitute/sting/commandline/ParsingEngine.java
View File

@ -30,6 +30,7 @@ import org.apache.commons.io.FileUtils;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.classloader.JVMUtils;
import org.broadinstitute.sting.utils.classloader.PluginManager;
import org.broadinstitute.sting.utils.collections.Pair;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.exceptions.UserException;
@ -61,7 +62,7 @@ public class ParsingEngine {
* Indicates as best as possible where command-line text remains unmatched
* to existing arguments.
*/
ArgumentMatches argumentMatches = null;
private ArgumentMatches argumentMatches = null;
/**
* Techniques for parsing and for argument lookup.
@ -88,7 +89,10 @@ public class ParsingEngine {
/**
* List of tags associated with the given instantiation of the command-line argument.
*/
private final Map<Object,Tags> tags = new IdentityHashMap<Object,Tags>();
private final Map<Object,Tags> tags = new IdentityHashMap<Object,Tags>();
private PluginManager<ParsingEngineArgumentProvider> argumentProviderPluginManager =
new PluginManager<ParsingEngineArgumentProvider>(ParsingEngineArgumentProvider.class);
/**
* our log, which we want to capture anything from org.broadinstitute.sting
@ -105,7 +109,10 @@ public class ParsingEngine {
argumentTypeDescriptors.addAll(clp.getArgumentTypeDescriptors());
argumentTypeDescriptors.addAll(STANDARD_ARGUMENT_TYPE_DESCRIPTORS);
addArgumentSource(ParsingEngineArgumentFiles.class);
List<Class<? extends ParsingEngineArgumentProvider>> providers = argumentProviderPluginManager.getPlugins();
for (Class<? extends ParsingEngineArgumentProvider> provider: providers) {
addArgumentSource(provider);
}
}
/**
@ -117,6 +124,10 @@ public class ParsingEngine {
addArgumentSource(null, source);
}
public ArgumentMatches getArgumentMatches() {
return argumentMatches;
}
/**
* Add an argument source. Argument sources are expected to have
* any number of fields with an @Argument annotation attached.
@ -156,29 +167,30 @@ public class ParsingEngine {
* @param tokens Tokens passed on the command line.
* @return The parsed arguments by file.
*/
public SortedMap<ArgumentMatchSource, List<String>> parse( String[] tokens ) {
public SortedMap<ArgumentMatchSource, ParsedArgs> parse( String[] tokens ) {
argumentMatches = new ArgumentMatches();
SortedMap<ArgumentMatchSource, List<String>> parsedArgs = new TreeMap<ArgumentMatchSource, List<String>>();
SortedMap<ArgumentMatchSource, ParsedArgs> parsedArgs = new TreeMap<ArgumentMatchSource, ParsedArgs>();
List<String> cmdLineTokens = Arrays.asList(tokens);
parse(ArgumentMatchSource.COMMAND_LINE, cmdLineTokens, argumentMatches, parsedArgs);
ParsingEngineArgumentFiles argumentFiles = new ParsingEngineArgumentFiles();
List<ParsingEngineArgumentProvider> providers = argumentProviderPluginManager.createAllTypes();
// Load the arguments ONLY into the argument files.
// Validation may optionally run on the rest of the arguments.
loadArgumentsIntoObject(argumentFiles);
for (ParsingEngineArgumentProvider provider: providers) {
// Load the arguments ONLY into the provider.
// Validation may optionally run on the rest of the arguments.
loadArgumentsIntoObject(provider);
}
for (File file: argumentFiles.files) {
List<String> fileTokens = getArguments(file);
parse(new ArgumentMatchSource(file), fileTokens, argumentMatches, parsedArgs);
for (ParsingEngineArgumentProvider provider: providers) {
provider.parse(this, parsedArgs);
}
return parsedArgs;
}
private void parse(ArgumentMatchSource matchSource, List<String> tokens,
ArgumentMatches argumentMatches, SortedMap<ArgumentMatchSource, List<String>> parsedArgs) {
public void parse(ArgumentMatchSource matchSource, List<String> tokens,
ArgumentMatches argumentMatches, SortedMap<ArgumentMatchSource, ParsedArgs> parsedArgs) {
ArgumentMatchSite lastArgumentMatchSite = new ArgumentMatchSite(matchSource, -1);
int i = 0;
@ -195,19 +207,44 @@ public class ParsingEngine {
}
else {
if( argumentMatches.hasMatch(lastArgumentMatchSite) &&
!argumentMatches.getMatch(lastArgumentMatchSite).hasValueAtSite(lastArgumentMatchSite))
argumentMatches.getMatch(lastArgumentMatchSite).addValue( lastArgumentMatchSite, token );
!argumentMatches.getMatch(lastArgumentMatchSite).hasValueAtSite(lastArgumentMatchSite))
argumentMatches.getMatch(lastArgumentMatchSite).addValue( lastArgumentMatchSite, new ArgumentMatchStringValue(token) );
else
argumentMatches.MissingArgument.addValue( site, token );
argumentMatches.MissingArgument.addValue( site, new ArgumentMatchStringValue(token) );
}
i++;
}
parsedArgs.put(matchSource, tokens);
parsedArgs.put(matchSource, new ParsedListArgs(tokens));
}
private List<String> getArguments(File file) {
public void parsePairs(ArgumentMatchSource matchSource, List<Pair<String, ArgumentMatchValue>> tokens,
ArgumentMatches argumentMatches, ParsedArgs matchSourceArgs,
SortedMap<ArgumentMatchSource, ParsedArgs> parsedArgs) {
int i = 0;
for (Pair<String, ArgumentMatchValue> pair: tokens) {
ArgumentMatchSite site = new ArgumentMatchSite(matchSource, i);
List<DefinitionMatcher> matchers = Arrays.asList(ArgumentDefinitions.FullNameDefinitionMatcher, ArgumentDefinitions.ShortNameDefinitionMatcher);
ArgumentDefinition definition = null;
for (DefinitionMatcher matcher: matchers) {
definition = argumentDefinitions.findArgumentDefinition( pair.getFirst(), matcher );
if (definition != null)
break;
}
if (definition == null)
continue;
ArgumentMatch argumentMatch = new ArgumentMatch(pair.getFirst(), definition, site, new Tags());
argumentMatches.mergeInto(argumentMatch);
argumentMatch.addValue(site, pair.getSecond());
i++;
}
parsedArgs.put(matchSource, matchSourceArgs);
}
protected List<String> getArguments(File file) {
try {
if (file.getAbsolutePath().endsWith(".list")) {
return getListArguments(file);
@ -283,9 +320,9 @@ public class ParsingEngine {
// Ensure that the field contents meet the validation criteria specified by the regular expression.
for( ArgumentMatch verifiableMatch: verifiableMatches ) {
for( String value: verifiableMatch.values() ) {
if( verifiableArgument.validation != null && !value.matches(verifiableArgument.validation) )
invalidValues.add( new Pair<ArgumentDefinition,String>(verifiableArgument, value) );
for( ArgumentMatchValue value: verifiableMatch.values() ) {
if( verifiableArgument.validation != null && !value.asString().matches(verifiableArgument.validation) )
invalidValues.add( new Pair<ArgumentDefinition,String>(verifiableArgument, value.asString()) );
}
}
}
@ -629,21 +666,21 @@ class UnmatchedArgumentException extends ArgumentException {
private static String formatArguments( ArgumentMatch invalidValues ) {
StringBuilder sb = new StringBuilder();
for( ArgumentMatchSite site: invalidValues.sites.keySet() )
for( String value: invalidValues.sites.get(site) ) {
for( ArgumentMatchValue value: invalidValues.sites.get(site) ) {
switch (site.getSource().getType()) {
case CommandLine:
sb.append( String.format("%nInvalid argument value '%s' at position %d.",
value, site.getIndex()) );
value.asString(), site.getIndex()) );
break;
case File:
sb.append( String.format("%nInvalid argument value '%s' in file %s at position %d.",
value, site.getSource().getFile().getAbsolutePath(), site.getIndex()) );
case Provider:
sb.append( String.format("%nInvalid argument value '%s' in %s at position %d.",
value.asString(), site.getSource().getDescription(), site.getIndex()) );
break;
default:
throw new RuntimeException( String.format("Unexpected argument match source type: %s",
site.getSource().getType()));
}
if(value != null && Utils.dupString(' ',value.length()).equals(value))
if(value.asString() != null && Utils.dupString(' ',value.asString().length()).equals(value.asString()))
sb.append(" Please make sure any line continuation backslashes on your command line are not followed by whitespace.");
}
return sb.toString();
@ -696,12 +733,3 @@ class UnknownEnumeratedValueException extends ArgumentException {
return String.format("Invalid value %s specified for argument %s; valid options are (%s).", argumentPassed, definition.fullName, Utils.join(",",definition.validOptions));
}
}
/**
* Container class to store the list of argument files.
* The files will be parsed after the command line arguments.
*/
class ParsingEngineArgumentFiles {
@Argument(fullName = "arg_file", shortName = "args", doc = "Reads arguments from the specified file", required = false)
public List<File> files = new ArrayList<File>();
}

30
public/java/src/org/broadinstitute/sting/commandline/ParsingEngineArgumentFiles.java 100644
View File

@ -0,0 +1,30 @@
package org.broadinstitute.sting.commandline;
import java.io.File;
import java.util.ArrayList;
import java.util.List;
import java.util.SortedMap;
/**
* Container class to store the list of argument files.
* The files will be parsed after the command line arguments.
*/
public class ParsingEngineArgumentFiles extends ParsingEngineArgumentProvider {
@Argument(fullName = "arg_file", shortName = "args", doc = "Reads arguments from the specified file", required = false)
public List<File> files = new ArrayList<File>();
@Override
public void parse(ParsingEngine parsingEngine, SortedMap<ArgumentMatchSource, ParsedArgs> parsedArgs) {
ArgumentMatches argumentMatches = parsingEngine.getArgumentMatches();
for (File file: this.files) {
List<String> fileTokens = parsingEngine.getArguments(file);
parsingEngine.parse(new ArgumentMatchFileSource(file), fileTokens, argumentMatches, parsedArgs);
}
}
}
class ArgumentMatchFileSource extends ArgumentMatchSource {
ArgumentMatchFileSource(File file) {
super("file " + file.getAbsolutePath());
}
}

12
public/java/src/org/broadinstitute/sting/commandline/ParsingEngineArgumentProvider.java 100644
View File

@ -0,0 +1,12 @@
package org.broadinstitute.sting.commandline;
import java.util.List;
import java.util.SortedMap;
/**
* A class that can parse arguments for the engine
*/
public abstract class ParsingEngineArgumentProvider {
public abstract void parse(ParsingEngine parsingEngine, SortedMap<ArgumentMatchSource, ParsedArgs> parsedArgs);
}

2
public/java/src/org/broadinstitute/sting/gatk/io/stubs/OutputStreamArgumentTypeDescriptor.java
View File

@ -86,7 +86,7 @@ public class OutputStreamArgumentTypeDescriptor extends ArgumentTypeDescriptor {
@Override
public Object parse( ParsingEngine parsingEngine, ArgumentSource source, Type type, ArgumentMatches matches ) {
ArgumentDefinition definition = createDefaultArgumentDefinition(source);
String fileName = getArgumentValue( definition, matches );
String fileName = getArgumentValue( definition, matches ).asString();
// This parser has been passed a null filename and the GATK is not responsible for creating a type default for the object;
// therefore, the user must have failed to specify a type default

12
public/java/src/org/broadinstitute/sting/gatk/io/stubs/SAMFileReaderArgumentTypeDescriptor.java
View File

@ -25,15 +25,11 @@
package org.broadinstitute.sting.gatk.io.stubs;
import net.sf.samtools.SAMFileReader;
import org.broadinstitute.sting.commandline.ArgumentMatches;
import org.broadinstitute.sting.commandline.ArgumentSource;
import org.broadinstitute.sting.commandline.ArgumentTypeDescriptor;
import org.broadinstitute.sting.commandline.ParsingEngine;
import org.broadinstitute.sting.commandline.*;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.sam.SAMFileReaderBuilder;
import java.io.File;
import java.lang.reflect.Type;
/**
@ -47,7 +43,7 @@ public class SAMFileReaderArgumentTypeDescriptor extends ArgumentTypeDescriptor
/**
* Create a new SAMFileReader argument, notifying the given engine when that argument has been created.
* @param engine
* @param engine engine
*/
public SAMFileReaderArgumentTypeDescriptor( GenomeAnalysisEngine engine ) {
this.engine = engine;
@ -62,12 +58,12 @@ public class SAMFileReaderArgumentTypeDescriptor extends ArgumentTypeDescriptor
public Object parse( ParsingEngine parsingEngine, ArgumentSource source, Type type, ArgumentMatches matches ) {
SAMFileReaderBuilder builder = new SAMFileReaderBuilder();
String readerFileName = getArgumentValue( createDefaultArgumentDefinition(source), matches );
ArgumentMatchValue readerFileName = getArgumentValue( createDefaultArgumentDefinition(source), matches );
if( readerFileName == null )
throw new UserException.CommandLineException("SAM file compression was supplied, but no associated writer was supplied with it.");
builder.setSAMFile(new File(readerFileName));
builder.setSAMFile(readerFileName.asFile());
// WARNING: Skipping required side-effect because stub is impossible to generate.
engine.addInput(source, builder);

45
public/java/src/org/broadinstitute/sting/gatk/io/stubs/SAMFileWriterArgumentTypeDescriptor.java
View File

@ -31,7 +31,6 @@ import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.io.StingSAMFileWriter;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import java.io.File;
import java.io.OutputStream;
import java.lang.annotation.Annotation;
import java.lang.reflect.Type;
@ -111,10 +110,10 @@ public class SAMFileWriterArgumentTypeDescriptor extends ArgumentTypeDescriptor
public Object parse( ParsingEngine parsingEngine, ArgumentSource source, Type type, ArgumentMatches matches ) {
// Extract all possible parameters that could be passed to a BAM file writer?
ArgumentDefinition bamArgumentDefinition = createBAMArgumentDefinition(source);
String writerFileName = getArgumentValue( bamArgumentDefinition, matches );
ArgumentMatchValue writerFileName = getArgumentValue( bamArgumentDefinition, matches );
String compressionLevelText = getArgumentValue( createBAMCompressionArgumentDefinition(source), matches );
Integer compressionLevel = compressionLevelText != null ? Integer.valueOf(compressionLevelText) : null;
ArgumentMatchValue compressionLevelText = getArgumentValue( createBAMCompressionArgumentDefinition(source), matches );
Integer compressionLevel = compressionLevelText != null ? Integer.valueOf(compressionLevelText.asString()) : null;
boolean indexOnTheFly = !argumentIsPresent(disableWriteIndexArgumentDefinition(source),matches);
boolean generateMD5 = argumentIsPresent(this.enableMD5GenerationArgumentDefinition(source),matches);
@ -124,32 +123,28 @@ public class SAMFileWriterArgumentTypeDescriptor extends ArgumentTypeDescriptor
// This parser has been passed a null filename and the GATK is not responsible for creating a type default for the object;
// therefore, the user must have failed to specify a type default
if(writerFileName == null) {
if(!source.isRequired())
throw new MissingArgumentValueException(bamArgumentDefinition);
if(generateMD5)
if(writerFileName.asFile() == null && generateMD5)
throw new ArgumentException("MD5 generation specified, but no output file specified. If md5 generation is desired, please specify a BAM output file and an md5 file will be written alongside.");
}
// Create the stub and set parameters.
SAMFileWriterStub stub;
if ( writerFileName != null )
stub = new SAMFileWriterStub(engine, new File(writerFileName));
else
stub = new SAMFileWriterStub(engine, defaultOutputStream);
SAMFileWriterStub stub = null; // stub = new SAMFileWriterStub(engine, defaultOutputStream);
if ( compressionLevel != null )
stub.setCompressionLevel(compressionLevel);
if ( indexOnTheFly )
stub.setIndexOnTheFly(indexOnTheFly);
if ( generateMD5 )
stub.setGenerateMD5(generateMD5);
if ( simplifyBAM )
stub.setSimplifyBAM(simplifyBAM);
if ( writerFileName.asFile() != null ) {
stub = new SAMFileWriterStub(engine, writerFileName.asFile());
// WARNING: Side effects required by engine!
parsingEngine.addTags(stub,getArgumentTags(matches));
engine.addOutput(stub);
if ( compressionLevel != null )
stub.setCompressionLevel(compressionLevel);
if ( indexOnTheFly )
stub.setIndexOnTheFly(indexOnTheFly);
if ( generateMD5 )
stub.setGenerateMD5(generateMD5);
if ( simplifyBAM )
stub.setSimplifyBAM(simplifyBAM);
// WARNING: Side effects required by engine!
parsingEngine.addTags(stub,getArgumentTags(matches));
engine.addOutput(stub);
}
return stub;
}

6
public/java/src/org/broadinstitute/sting/gatk/io/stubs/VCFWriterArgumentTypeDescriptor.java
View File

@ -138,8 +138,8 @@ public class VCFWriterArgumentTypeDescriptor extends ArgumentTypeDescriptor {
public Object parse( ParsingEngine parsingEngine, ArgumentSource source, Type type, ArgumentMatches matches ) {
ArgumentDefinition defaultArgumentDefinition = createDefaultArgumentDefinition(source);
// Get the filename for the genotype file, if it exists. If not, we'll need to send output to out.
String writerFileName = getArgumentValue(defaultArgumentDefinition,matches);
File writerFile = writerFileName != null ? new File(writerFileName) : null;
ArgumentMatchValue writerFileName = getArgumentValue(defaultArgumentDefinition,matches);
File writerFile = writerFileName != null ? writerFileName.asFile() : null;
// This parser has been passed a null filename and the GATK is not responsible for creating a type default for the object;
// therefore, the user must have failed to specify a type default
@ -151,7 +151,7 @@ public class VCFWriterArgumentTypeDescriptor extends ArgumentTypeDescriptor {
? new VariantContextWriterStub(engine, writerFile, argumentSources)
: new VariantContextWriterStub(engine, defaultOutputStream, argumentSources);
stub.setCompressed(isCompressed(writerFileName));
stub.setCompressed(isCompressed(writerFileName.asString()));
stub.setDoNotWriteGenotypes(argumentIsPresent(createSitesOnlyArgumentDefinition(),matches));
stub.setSkipWritingCommandLineHeader(argumentIsPresent(createNoCommandLineHeaderArgumentDefinition(),matches));
stub.setForceBCF(argumentIsPresent(createBCFArgumentDefinition(),matches));

65
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalc.java
View File

@ -29,24 +29,16 @@ import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.SimpleTimer;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.variantcontext.Allele;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.GenotypesContext;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import java.io.File;
import java.io.FileNotFoundException;
import java.io.FileOutputStream;
import java.io.PrintStream;
import java.util.Arrays;
import java.util.List;
/**
* Generic interface for calculating the probability of alleles segregating given priors and genotype likelihoods
*
*/
public abstract class AFCalc implements Cloneable {
private final static Logger defaultLogger = Logger.getLogger(AFCalc.class);
@ -58,8 +50,8 @@ public abstract class AFCalc implements Cloneable {
protected Logger logger = defaultLogger;
private SimpleTimer callTimer = new SimpleTimer();
private PrintStream callReport = null;
private final AFCalcResultTracker resultTracker;
private ExactCallLogger exactCallLogger = null;
protected AFCalc(final int nSamples, final int maxAltAlleles, final int maxAltAllelesForIndels, final int ploidy) {
if ( nSamples < 0 ) throw new IllegalArgumentException("nSamples must be greater than zero " + nSamples);
@ -74,7 +66,7 @@ public abstract class AFCalc implements Cloneable {
}
public void enableProcessLog(final File exactCallsLog) {
initializeOutputFile(exactCallsLog);
exactCallLogger = new ExactCallLogger(exactCallsLog);
}
public void setLogger(Logger logger) {
@ -102,8 +94,8 @@ public abstract class AFCalc implements Cloneable {
final AFCalcResult result = computeLog10PNonRef(vcWorking, log10AlleleFrequencyPriors);
final long nanoTime = callTimer.getElapsedTimeNano();
if ( callReport != null )
printCallInfo(vcWorking, log10AlleleFrequencyPriors, nanoTime, resultTracker.getLog10PosteriorOfAFzero());
if ( exactCallLogger != null )
exactCallLogger.printCallInfo(vcWorking, log10AlleleFrequencyPriors, nanoTime, result);
return result;
}
@ -170,55 +162,8 @@ public abstract class AFCalc implements Cloneable {
return Math.max(maxAlternateAllelesToGenotype, maxAlternateAllelesForIndels);
}
// ---------------------------------------------------------------------------
//
// Print information about the call to the calls log
//
// ---------------------------------------------------------------------------
private void initializeOutputFile(final File outputFile) {
try {
if (outputFile != null) {
callReport = new PrintStream( new FileOutputStream(outputFile) );
callReport.println(Utils.join("\t", Arrays.asList("loc", "variable", "key", "value")));
}
} catch ( FileNotFoundException e ) {
throw new UserException.CouldNotCreateOutputFile(outputFile, e);
}
}
private void printCallInfo(final VariantContext vc,
final double[] log10AlleleFrequencyPriors,
final long runtimeNano,
final double log10PosteriorOfAFzero) {
printCallElement(vc, "type", "ignore", vc.getType());
int allelei = 0;
for ( final Allele a : vc.getAlleles() )
printCallElement(vc, "allele", allelei++, a.getDisplayString());
for ( final Genotype g : vc.getGenotypes() )
printCallElement(vc, "PL", g.getSampleName(), g.getLikelihoodsString());
for ( int priorI = 0; priorI < log10AlleleFrequencyPriors.length; priorI++ )
printCallElement(vc, "priorI", priorI, log10AlleleFrequencyPriors[priorI]);
printCallElement(vc, "runtime.nano", "ignore", runtimeNano);
printCallElement(vc, "log10PosteriorOfAFzero", "ignore", log10PosteriorOfAFzero);
callReport.flush();
}
private void printCallElement(final VariantContext vc,
final Object variable,
final Object key,
final Object value) {
final String loc = String.format("%s:%d", vc.getChr(), vc.getStart());
callReport.println(Utils.join("\t", Arrays.asList(loc, variable, key, value)));
}
public AFCalcResultTracker getResultTracker() {
return resultTracker;
}
}

4
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcFactory.java
View File

@ -30,10 +30,6 @@ public class AFCalcFactory {
/** reference implementation of multi-allelic EXACT model */
EXACT_REFERENCE(ReferenceDiploidExactAFCalc.class, 2, -1),
/** expt. implementation */
@Deprecated
EXACT_CONSTRAINED(ConstrainedDiploidExactAFCalc.class, 2, -1),
/** expt. implementation -- for testing only */
EXACT_INDEPENDENT(IndependentAllelesDiploidExactAFCalc.class, 2, -1),

36
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/AFCalcResult.java
View File

@ -31,10 +31,7 @@ import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.variantcontext.Allele;
import java.util.Arrays;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.*;
/**
* Describes the results of the AFCalc
@ -217,6 +214,14 @@ public class AFCalcResult {
return log10PriorsOfAC[AF1p];
}
@Override
public String toString() {
final List<String> byAllele = new LinkedList<String>();
for ( final Allele a : getAllelesUsedInGenotyping() )
if ( a.isNonReference() ) byAllele.add(String.format("%s => MLE %d / posterior %.2f", a, getAlleleCountAtMLE(a), getLog10PosteriorOfAFGt0ForAllele(a)));
return String.format("AFCalc%n\t\tlog10PosteriorOfAFGT0=%.2f%n\t\t%s", getLog10LikelihoodOfAFGT0(), Utils.join("\n\t\t", byAllele));
}
/**
* Are we sufficiently confidence in being non-ref that the site is considered polymorphic?
*
@ -233,6 +238,19 @@ public class AFCalcResult {
return getLog10PosteriorOfAFGt0ForAllele(allele) >= log10minPNonRef;
}
/**
* Are any of the alleles polymorphic w.r.t. #isPolymorphic?
*
* @param log10minPNonRef the confidence threshold, in log10 space
* @return true if any are poly, false otherwise
*/
public boolean anyPolymorphic(final double log10minPNonRef) {
for ( final Allele a : getAllelesUsedInGenotyping() )
if ( a.isNonReference() && isPolymorphic(a, log10minPNonRef) )
return true;
return false;
}
/**
* Returns the log10 probability that allele is segregating
*
@ -266,15 +284,7 @@ public class AFCalcResult {
final double[] log10UnnormalizedPosteriors = new double[log10LikelihoodsOfAC.length];
for ( int i = 0; i < log10LikelihoodsOfAC.length; i++ )
log10UnnormalizedPosteriors[i] = log10LikelihoodsOfAC[i] + log10PriorsOfAC[i];
// necessary because the posteriors may be so skewed that the log-space normalized value isn't
// good, so we have to try both log-space normalization as well as the real-space normalization if the
// result isn't good
final double[] logNormalized = MathUtils.normalizeFromLog10(log10UnnormalizedPosteriors, true, true);
if ( goodLog10ProbVector(logNormalized, logNormalized.length, true) )
return logNormalized;
else
return MathUtils.normalizeFromLog10(log10UnnormalizedPosteriors, true, false);
return MathUtils.normalizeFromLog10(log10UnnormalizedPosteriors, true, false);
}
/**

107
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ConstrainedDiploidExactAFCalc.java
View File

@ -1,107 +0,0 @@
package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc;
import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import org.broadinstitute.sting.utils.MathUtils;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.GenotypeLikelihoods;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
@Deprecated
public class ConstrainedDiploidExactAFCalc extends DiploidExactAFCalc {
protected ConstrainedDiploidExactAFCalc(int nSamples, int maxAltAlleles, int maxAltAllelesForIndels, final int ploidy) {
super(nSamples, maxAltAlleles, maxAltAllelesForIndels, ploidy);
}
protected StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResultTracker resultTracker) {
final int[] maxACsToConsider = computeMaxACs(vc);
resultTracker.setAClimits(maxACsToConsider);
return new StateTracker(maxACsToConsider);
}
/**
* Computes the maximum ACs we need to consider for each alt allele
*
* Walks over the genotypes in VC, and computes for each alt allele the maximum
* AC we need to consider in that alt allele dimension. Does the calculation
* based on the PLs in each genotype g, choosing to update the max AC for the
* alt alleles corresponding to that PL. Only takes the first lowest PL,
* if there are multiple genotype configurations with the same PL value. It
* takes values in the order of the alt alleles.
*
* @param vc the variant context we will compute max alt alleles for
* @return a vector of max alt alleles, indexed by alt allele, so result[0] is the AC of the
* first alt allele.
*/
@Ensures("result != null")
protected final int[] computeMaxACs(final VariantContext vc) {
final int[] maxACs = new int[vc.getNAlleles()-1];
for ( final Genotype g : vc.getGenotypes() )
updateMaxACs(g, maxACs);
return maxACs;
}
/**
* Update the maximum achievable allele counts in maxAC according to the PLs in g
*
* Selects the maximum genotype configuration from the PLs in g, and updates
* the maxAC for this configure. For example, if the lowest PL is for 0/1, updates
* the maxAC for the alt allele 1 by 1. If it's 1/1, update is 2. Works for
* many number of alt alleles (determined by length of maxACs).
*
* If the max PL occurs at 0/0, updates nothing
* Note that this function greedily takes the first min PL, so that if 0/1 and 1/1 have
* the same PL value, then updates the first one.
*
* Also, only will update 1 alt allele, so if 0/1 and 0/2 both have the same PL,
* then only first one (1) will be updated
*
* @param g the genotype to update
* @param maxACs the max allele count vector for alt alleles (starting at 0 => first alt allele)
*/
@Requires({
"g != null",
"maxACs != null",
"goodMaxACs(maxACs)"})
private void updateMaxACs(final Genotype g, final int[] maxACs) {
final int[] PLs = g.getLikelihoods().getAsPLs();
int minPLi = 0;
int minPL = PLs[0];
for ( int i = 0; i < PLs.length; i++ ) {
if ( PLs[i] < minPL ) {
minPL = PLs[i];
minPLi = i;
}
}
final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair pair = GenotypeLikelihoods.getAllelePair(minPLi);
updateMaxACs(maxACs, pair.alleleIndex1);
updateMaxACs(maxACs, pair.alleleIndex2);
}
/**
* Simple helper. Update max alt alleles maxACs according to the allele index (where 0 == ref)
*
* If alleleI == 0 => doesn't update anything
* else maxACs[alleleI - 1]++
*
* @param maxACs array of max alt allele ACs
* @param alleleI the index (relative to 0) to update a count of 1 in max alt alleles.
*/
@Requires({
"alleleI >= 0",
"(alleleI - 1) < maxACs.length",
"goodMaxACs(maxACs)"})
private void updateMaxACs(final int[] maxACs, final int alleleI) {
if ( alleleI > 0 )
maxACs[alleleI-1]++;
}
private static boolean goodMaxACs(final int[] maxACs) {
return MathUtils.sum(maxACs) >= 0;
}
}

4
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/DiploidExactAFCalc.java
View File

@ -36,7 +36,9 @@ public abstract class DiploidExactAFCalc extends ExactAFCalc {
if ( ploidy != 2 ) throw new IllegalArgumentException("ploidy must be two for DiploidExactAFCalc and subclasses but saw " + ploidy);
}
protected abstract StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResultTracker resultTracker);
protected StateTracker makeMaxLikelihood(VariantContext vc, AFCalcResultTracker resultTracker) {
return new StateTracker();
}
@Override
protected AFCalcResult computeLog10PNonRef(final VariantContext vc,

179
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ExactCallLogger.java 100644
View File

@ -0,0 +1,179 @@
package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc;
import com.google.java.contract.Requires;
import org.apache.commons.lang.ArrayUtils;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.variantcontext.*;
import java.io.*;
import java.util.*;
/**
* Allows us to write out and read in information about exact calls (site, alleles, PLs, etc) in tabular format
*
* Once opened, calls can be writen to disk with printCallInfo
*/
public class ExactCallLogger implements Cloneable {
private PrintStream callReport = null;
/**
* Create a new ExactCallLogger writing it's output to outputFile
*
* @param outputFile
*/
public ExactCallLogger(final File outputFile) {
try {
callReport = new PrintStream(new BufferedOutputStream(new FileOutputStream(outputFile), 10000000));
callReport.println(Utils.join("\t", Arrays.asList("loc", "variable", "key", "value")));
} catch (FileNotFoundException e) {
throw new UserException.CouldNotCreateOutputFile(outputFile, e);
}
}
/**
* Summarizes information about an exact call that happened
*/
public static class ExactCall {
final VariantContext vc;
final long runtime;
final AFCalcResult originalCall;
public ExactCall(VariantContext vc, final long runtime, final AFCalcResult originalCall) {
this.vc = vc;
this.runtime = runtime;
this.originalCall = originalCall;
}
@Override
public String toString() {
return String.format("ExactCall %s:%d alleles=%s nSamples=%s orig.pNonRef=%.2f orig.runtime=%s",
vc.getChr(), vc.getStart(), vc.getAlleles(), vc.getNSamples(),
originalCall.getLog10PosteriorOfAFGT0(),
new AutoFormattingTime(runtime / 1e9).toString());
}
}
protected final void printCallInfo(final VariantContext vc,
final double[] log10AlleleFrequencyPriors,
final long runtimeNano,
final AFCalcResult result) {
printCallElement(vc, "type", "ignore", vc.getType());
int allelei = 0;
for (final Allele a : vc.getAlleles())
printCallElement(vc, "allele", allelei++, a.getDisplayString());
for (final Genotype g : vc.getGenotypes())
printCallElement(vc, "PL", g.getSampleName(), g.getLikelihoodsString());
for (int priorI = 0; priorI < log10AlleleFrequencyPriors.length; priorI++)
printCallElement(vc, "priorI", priorI, log10AlleleFrequencyPriors[priorI]);
printCallElement(vc, "runtime.nano", "ignore", runtimeNano);
printCallElement(vc, "log10PosteriorOfAFEq0", "ignore", result.getLog10PosteriorOfAFEq0());
printCallElement(vc, "log10PosteriorOfAFGt0", "ignore", result.getLog10PosteriorOfAFGT0());
for ( final Allele allele : result.getAllelesUsedInGenotyping() ) {
if ( allele.isNonReference() ) {
printCallElement(vc, "MLE", allele, result.getAlleleCountAtMLE(allele));
printCallElement(vc, "pNonRefByAllele", allele, result.getLog10PosteriorOfAFGt0ForAllele(allele));
}
}
callReport.flush();
}
@Requires({"vc != null", "variable != null", "key != null", "value != null", "callReport != null"})
private void printCallElement(final VariantContext vc,
final Object variable,
final Object key,
final Object value) {
final String loc = String.format("%s:%d", vc.getChr(), vc.getStart());
callReport.println(Utils.join("\t", Arrays.asList(loc, variable, key, value)));
}
/**
* Read in a list of ExactCall objects from reader, keeping only those
* with starts in startsToKeep or all sites (if this is empty)
*
* @param reader a just-opened reader sitting at the start of the file
* @param startsToKeep a list of start position of the calls to keep, or empty if all calls should be kept
* @param parser a genome loc parser to create genome locs
* @return a list of ExactCall objects in reader
* @throws IOException
*/
public static List<ExactCall> readExactLog(final BufferedReader reader, final List<Integer> startsToKeep, GenomeLocParser parser) throws IOException {
if ( reader == null ) throw new IllegalArgumentException("reader cannot be null");
if ( startsToKeep == null ) throw new IllegalArgumentException("startsToKeep cannot be null");
if ( parser == null ) throw new IllegalArgumentException("GenomeLocParser cannot be null");
List<ExactCall> calls = new LinkedList<ExactCall>();
// skip the header line
reader.readLine();
// skip the first "type" line
reader.readLine();
while (true) {
final VariantContextBuilder builder = new VariantContextBuilder();
final List<Allele> alleles = new ArrayList<Allele>();
final List<Genotype> genotypes = new ArrayList<Genotype>();
final double[] posteriors = new double[2];
final double[] priors = MathUtils.normalizeFromLog10(new double[]{0.5, 0.5}, true);
final List<Integer> mle = new ArrayList<Integer>();
final Map<Allele, Double> log10pNonRefByAllele = new HashMap<Allele, Double>();
long runtimeNano = -1;
GenomeLoc currentLoc = null;
while (true) {
final String line = reader.readLine();
if (line == null)
return calls;
final String[] parts = line.split("\t");
final GenomeLoc lineLoc = parser.parseGenomeLoc(parts[0]);
final String variable = parts[1];
final String key = parts[2];
final String value = parts[3];
if (currentLoc == null)
currentLoc = lineLoc;
if (variable.equals("type")) {
if (startsToKeep.isEmpty() || startsToKeep.contains(currentLoc.getStart())) {
builder.alleles(alleles);
final int stop = currentLoc.getStart() + alleles.get(0).length() - 1;
builder.chr(currentLoc.getContig()).start(currentLoc.getStart()).stop(stop);
builder.genotypes(genotypes);
final int[] mleInts = ArrayUtils.toPrimitive(mle.toArray(new Integer[]{}));
final AFCalcResult result = new AFCalcResult(mleInts, 1, alleles, posteriors, priors, log10pNonRefByAllele);
calls.add(new ExactCall(builder.make(), runtimeNano, result));
}
break;
} else if (variable.equals("allele")) {
final boolean isRef = key.equals("0");
alleles.add(Allele.create(value, isRef));
} else if (variable.equals("PL")) {
final GenotypeBuilder gb = new GenotypeBuilder(key);
gb.PL(GenotypeLikelihoods.fromPLField(value).getAsPLs());
genotypes.add(gb.make());
} else if (variable.equals("log10PosteriorOfAFEq0")) {
posteriors[0] = Double.valueOf(value);
} else if (variable.equals("log10PosteriorOfAFGt0")) {
posteriors[1] = Double.valueOf(value);
} else if (variable.equals("MLE")) {
mle.add(Integer.valueOf(value));
} else if (variable.equals("pNonRefByAllele")) {
final Allele a = Allele.create(key);
log10pNonRefByAllele.put(a, Double.valueOf(value));
} else if (variable.equals("runtime.nano")) {
runtimeNano = Long.valueOf(value);
} else {
// nothing to do
}
}
}
}
}

273
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/IndependentAllelesDiploidExactAFCalc.java
View File

@ -32,31 +32,97 @@ import org.broadinstitute.sting.utils.variantcontext.*;
import java.util.*;
public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
/**
* Computes the conditional bi-allelic exact results
*
* Suppose vc contains 2 alt allele: A* with C and T. This function first computes:
*
* (1) P(D | AF_c > 0 && AF_t == *) [i.e., T can be anything]
*
* it then computes the conditional probability on AF_c == 0:
*
* (2) P(D | AF_t > 0 && AF_c == 0)
*
* Thinking about this visually, we have the following likelihood matrix where each cell is
* the P(D | AF_c == i && AF_t == j):
*
* 0 AF_c > 0
* -----------------
* 0 | |
* |--|-------------
* a | |
* f | |
* _ | |
* t | |
* > | |
* 0 | |
*
* What we really want to know how
*
* (3) P(D | AF_c == 0 & AF_t == 0)
*
* compares with
*
* (4) P(D | AF_c > 0 || AF_t > 0)
*
* This is effectively asking for the value in the upper left vs. the sum of all cells.
*
* This class implements the conditional likelihoods summation for any number of alt
* alleles, where each alt allele has its EXACT probability of segregating calculated by
* reducing each alt B into the case XB and computing P(D | AF_b > 0 ) as follows:
*
* Suppose we have for a A/B/C site the following GLs:
*
* AA AB BB AC BC CC
*
* and we want to get the bi-allelic GLs for X/B, where X is everything not B
*
* XX = AA + AC + CC (since X = A or C)
* XB = AB + BC
* BB = BB
*
* After each allele has its probability calculated we compute the joint posterior
* as P(D | AF_* == 0) = prod_i P (D | AF_i == 0), after applying the theta^i
* prior for the ith least likely allele.
*/
public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
/**
* The min. confidence of an allele to be included in the joint posterior.
*/
private final static double MIN_LOG10_CONFIDENCE_TO_INCLUDE_ALLELE_IN_POSTERIOR = Math.log10(1e-20);
private final static int[] BIALLELIC_NON_INFORMATIVE_PLS = new int[]{0,0,0};
private final static List<Allele> BIALLELIC_NOCALL = Arrays.asList(Allele.NO_CALL, Allele.NO_CALL);
/**
* Sorts AFCalcResults by their posteriors of AF > 0, so the
*/
private final static class CompareAFCalcResultsByPNonRef implements Comparator<AFCalcResult> {
@Override
public int compare(AFCalcResult o1, AFCalcResult o2) {
return Double.compare(o1.getLog10LikelihoodOfAFGT0(), o2.getLog10LikelihoodOfAFGT0());
return Double.compare(o1.getLog10PosteriorOfAFGT0(), o2.getLog10PosteriorOfAFGT0());
}
}
private final static CompareAFCalcResultsByPNonRef compareAFCalcResultsByPNonRef = new CompareAFCalcResultsByPNonRef();
final ReferenceDiploidExactAFCalc refModel;
/**
* The AFCalc model we are using to do the bi-allelic computation
*/
final AFCalc biAlleleExactModel;
protected IndependentAllelesDiploidExactAFCalc(int nSamples, int maxAltAlleles, int maxAltAllelesForIndels, final int ploidy) {
super(nSamples, maxAltAlleles, maxAltAllelesForIndels, ploidy);
refModel = new ReferenceDiploidExactAFCalc(nSamples, 1, 1, ploidy);
}
@Override
protected StateTracker makeMaxLikelihood(VariantContext vc, AFCalcResultTracker resultTracker) {
return refModel.makeMaxLikelihood(vc, resultTracker);
biAlleleExactModel = new ReferenceDiploidExactAFCalc(nSamples, 1, 1, ploidy);
}
/**
* Trivial subclass that helps with debugging by keeping track of the supporting information for this joint call
*/
private static class MyAFCalcResult extends AFCalcResult {
/**
* List of the supporting bi-allelic AFCalcResults that went into making this multi-allelic joint call
*/
final List<AFCalcResult> supporting;
private MyAFCalcResult(int[] alleleCountsOfMLE, int nEvaluations, List<Allele> allelesUsedInGenotyping, double[] log10LikelihoodsOfAC, double[] log10PriorsOfAC, Map<Allele, Double> log10pNonRefByAllele, List<AFCalcResult> supporting) {
@ -68,121 +134,89 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
@Override
public AFCalcResult computeLog10PNonRef(final VariantContext vc,
final double[] log10AlleleFrequencyPriors) {
final double log10LikelihoodOfRef = computelog10LikelihoodOfRef(vc);
final List<AFCalcResult> independentResultTrackers = computeAlleleConditionalExact(vc, log10AlleleFrequencyPriors);
final List<AFCalcResult> independentResultTrackers = computeAlleleIndependentExact(vc, log10AlleleFrequencyPriors);
final List<AFCalcResult> withMultiAllelicPriors = applyMultiAllelicPriors(independentResultTrackers);
return combineIndependentPNonRefs(vc, log10LikelihoodOfRef, withMultiAllelicPriors);
return combineIndependentPNonRefs(vc, withMultiAllelicPriors);
}
protected final double computelog10LikelihoodOfRef(final VariantContext vc) {
// this value just the likelihood of AF == 0 in the special constrained multi-allelic calculation
final List<double[]> allGLs = getGLs(vc.getGenotypes(), false);
double log10LikelihoodOfHomRef = 0.0;
// TODO -- can be easily optimized (currently looks at all GLs via getGLs)
for ( int i = 0; i < allGLs.size(); i++ ) {
final double[] GLs = allGLs.get(i);
log10LikelihoodOfHomRef += GLs[0];
//log10LikelihoodOfHomRef += MathUtils.normalizeFromLog10(GLs, true)[0];
}
return log10LikelihoodOfHomRef;
}
/**
* Computes the conditional bi-allelic exact results
* Compute the conditional exact AFCalcResult for each allele in vc independently, returning
* the result of each, in order of the alt alleles in VC
*
* Suppose vc contains 2 alt allele: A* with C and T. This function first computes:
*
* (1) P(D | AF_c > 0 && AF_t == *) [i.e., T can be anything]
*
* it then computes the conditional probability on AF_c == 0:
*
* (2) P(D | AF_t > 0 && AF_c == 0)
*
* Thinking about this visually, we have the following likelihood matrix where each cell is
* the P(D | AF_c == i && AF_t == j):
*
* 0 AF_c > 0
* -----------------
* 0 | |
* |--|-------------
* a | |
* f | |
* _ | |
* t | |
* > | |
* 0 | |
*
* What we really want to know how
*
* (3) P(D | AF_c == 0 & AF_t == 0)
*
* compares with
*
* (4) P(D | AF_c > 0 || AF_t > 0)
*
* This is effectively asking for the value in the upper left vs. the sum of all cells.
*
* The quantity (1) is the same of all cells except those with AF_c == 0, while (2) is the
* band at the top where AF_t > 0 and AF_c == 0
*
* So (4) is actually (1) + (2).
*
* (3) is the direct inverse of the (1) and (2), as we are simultaneously calculating
*
* (1*) P(D | AF_c == 0 && AF_t == *) [i.e., T can be anything]
* (2*) P(D | AF_t == 0 && AF_c == 0) [TODO -- note this value looks like the thing we are supposed to use]
*
* This function implements the conditional likelihoods summation for any number of alt
* alleles (not just the tri-allelic case), where each subsequent variant context is
* further constrained such that each already considered allele x has AF_x == 0 in the
* compute.
*
* @param vc
* @param log10AlleleFrequencyPriors
* @return
* @param vc the VariantContext we want to analyze
* @param log10AlleleFrequencyPriors the priors
* @return a list of the AFCalcResults for each bi-allelic sub context of vc
*/
protected List<AFCalcResult> computeAlleleConditionalExact(final VariantContext vc,
final double[] log10AlleleFrequencyPriors) {
@Requires({"vc != null", "log10AlleleFrequencyPriors != null"})
@Ensures("goodIndependentResult(vc, result)")
protected final List<AFCalcResult> computeAlleleIndependentExact(final VariantContext vc,
final double[] log10AlleleFrequencyPriors) {
final List<AFCalcResult> results = new LinkedList<AFCalcResult>();
for ( final VariantContext subvc : makeAlleleConditionalContexts(vc) ) {
final AFCalcResult resultTracker = refModel.getLog10PNonRef(subvc, log10AlleleFrequencyPriors);
final AFCalcResult resultTracker = biAlleleExactModel.getLog10PNonRef(subvc, log10AlleleFrequencyPriors);
results.add(resultTracker);
}
return results;
}
protected List<VariantContext> makeAlleleConditionalContexts(final VariantContext vc) {
/**
* Helper function to ensure that the computeAlleleIndependentExact is returning reasonable results
*/
private static boolean goodIndependentResult(final VariantContext vc, final List<AFCalcResult> results) {
if ( results.size() != vc.getNAlleles() - 1) return false;
for ( int i = 0; i < results.size(); i++ ) {
if ( results.get(i).getAllelesUsedInGenotyping().size() != 2 )
return false;
if ( ! results.get(i).getAllelesUsedInGenotyping().contains(vc.getAlternateAllele(i)) )
return false;
}
return true;
}
/**
* Returns the bi-allelic variant context for each alt allele in vc with bi-allelic likelihoods, in order
*
* @param vc the variant context to split. Must have n.alt.alleles > 1
* @return a bi-allelic variant context for each alt allele in vc
*/
@Requires({"vc != null", "vc.getNAlleles() > 1"})
@Ensures("result.size() == vc.getNAlleles() - 1")
protected final List<VariantContext> makeAlleleConditionalContexts(final VariantContext vc) {
final int nAltAlleles = vc.getNAlleles() - 1;
final List<VariantContext> vcs = new LinkedList<VariantContext>();
final List<Allele> afZeroAlleles = new LinkedList<Allele>();
for ( int altI = 0; altI < nAltAlleles; altI++ ) {
final Allele altAllele = vc.getAlternateAllele(altI);
final List<Allele> biallelic = Arrays.asList(vc.getReference(), altAllele);
vcs.add(biallelicCombinedGLs(vc, biallelic, afZeroAlleles, altI + 1));
//afZeroAlleles.add(altAllele);
vcs.add(biallelicCombinedGLs(vc, altI + 1));
}
return vcs;
}
protected VariantContext biallelicCombinedGLs(final VariantContext rootVC, final List<Allele> biallelic, final List<Allele> afZeroAlleles, final int allele2) {
/**
* Create a single bi-allelic variant context from rootVC with alt allele with index altAlleleIndex
*
* @param rootVC the root (potentially multi-allelic) variant context
* @param altAlleleIndex index of the alt allele, from 0 == first alt allele
* @return a bi-allelic variant context based on rootVC
*/
@Requires({"rootVC.getNAlleles() > 1", "altAlleleIndex < rootVC.getNAlleles()"})
@Ensures({"result.isBiallelic()"})
protected final VariantContext biallelicCombinedGLs(final VariantContext rootVC, final int altAlleleIndex) {
if ( rootVC.isBiallelic() ) {
if ( ! afZeroAlleles.isEmpty() ) throw new IllegalArgumentException("Root VariantContext is biallelic but afZeroAlleles wasn't empty: " + afZeroAlleles);
return rootVC;
} else {
final Set<Integer> allelesToDiscard = new HashSet<Integer>(rootVC.getAlleleIndices(afZeroAlleles));
final int nAlts = rootVC.getNAlleles() - 1;
final List<Genotype> biallelicGenotypes = new ArrayList<Genotype>(rootVC.getNSamples());
for ( final Genotype g : rootVC.getGenotypes() )
biallelicGenotypes.add(combineGLs(g, allele2, allelesToDiscard, nAlts));
biallelicGenotypes.add(combineGLs(g, altAlleleIndex, nAlts));
final VariantContextBuilder vcb = new VariantContextBuilder(rootVC);
vcb.alleles(biallelic);
final Allele altAllele = rootVC.getAlternateAllele(altAlleleIndex - 1);
vcb.alleles(Arrays.asList(rootVC.getReference(), altAllele));
vcb.genotypes(biallelicGenotypes);
return vcb.make();
}
@ -203,30 +237,16 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
* XB = AB + BC
* BB = BB
*
* Supports the additional mode of simply dropping GLs whose allele index occurs in allelesToDiscard. This is
* useful in the case where you want to drop alleles (not combine them), such as above:
*
* AA AB BB AC BC CC
*
* and we want to get the bi-allelic GLs for X/B, where X is everything not B, but dropping C (index 2)
*
* XX = AA (since X = A and C is dropped)
* XB = AB
* BB = BB
*
* This allows us to recover partial GLs the correspond to any allele in allelesToDiscard having strictly
* AF == 0.
*
* @param original the original multi-allelic genotype
* @param altIndex the index of the alt allele we wish to keep in the bialleic case -- with ref == 0
* @param nAlts the total number of alt alleles
* @return a new biallelic genotype with appropriate PLs
*/
@Requires({"original.hasLikelihoods()", "! allelesToDiscard.contains(altIndex)"})
@Requires({"original.hasLikelihoods()"}) // TODO -- add ploidy == 2 test "original.getPLs() == null || original.getPLs().length == 3"})
@Ensures({"result.hasLikelihoods()", "result.getPL().length == 3"})
protected Genotype combineGLs(final Genotype original, final int altIndex, final Set<Integer> allelesToDiscard, final int nAlts ) {
protected Genotype combineGLs(final Genotype original, final int altIndex, final int nAlts ) {
if ( original.isNonInformative() )
return new GenotypeBuilder(original).PL(new int[]{0,0,0}).alleles(BIALLELIC_NOCALL).make();
return new GenotypeBuilder(original).PL(BIALLELIC_NON_INFORMATIVE_PLS).alleles(BIALLELIC_NOCALL).make();
if ( altIndex < 1 || altIndex > nAlts ) throw new IllegalStateException("altIndex must be between 1 and nAlts " + nAlts);
@ -236,10 +256,6 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
for ( int index = 0; index < normalizedPr.length; index++ ) {
final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair pair = GenotypeLikelihoods.getAllelePair(index);
// just continue if we shouldn't include the pair because it's in the discard set
if ( discardAllelePair(pair, allelesToDiscard) )
continue;
if ( pair.alleleIndex1 == altIndex ) {
if ( pair.alleleIndex2 == altIndex )
// hom-alt case
@ -263,11 +279,7 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
return new GenotypeBuilder(original).PL(GLs).alleles(BIALLELIC_NOCALL).make();
}
protected boolean discardAllelePair(final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair pair, Set<Integer> allelesToDiscard) {
return allelesToDiscard.contains(pair.alleleIndex1) || allelesToDiscard.contains(pair.alleleIndex2);
}
protected List<AFCalcResult> applyMultiAllelicPriors(final List<AFCalcResult> conditionalPNonRefResults) {
protected final List<AFCalcResult> applyMultiAllelicPriors(final List<AFCalcResult> conditionalPNonRefResults) {
final ArrayList<AFCalcResult> sorted = new ArrayList<AFCalcResult>(conditionalPNonRefResults);
// sort the results, so the most likely allele is first
@ -291,10 +303,11 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
/**
* Take the independent estimates of pNonRef for each alt allele and combine them into a single result
*
* TODO -- add more docs
*
* @param sortedResultsWithThetaNPriors the pNonRef result for each allele independently
*/
protected AFCalcResult combineIndependentPNonRefs(final VariantContext vc,
final double log10LikelihoodsOfACEq0,
final List<AFCalcResult> sortedResultsWithThetaNPriors) {
int nEvaluations = 0;
final int nAltAlleles = sortedResultsWithThetaNPriors.size();
@ -302,8 +315,8 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
final double[] log10PriorsOfAC = new double[2];
final Map<Allele, Double> log10pNonRefByAllele = new HashMap<Allele, Double>(nAltAlleles);
// this value is a sum in real space so we need to store values to sum up later
final double[] log10LikelihoodsOfACGt0 = new double[nAltAlleles];
// this value is a sum in log space
double log10PosteriorOfACEq0Sum = 0.0;
for ( final AFCalcResult sortedResultWithThetaNPriors : sortedResultsWithThetaNPriors ) {
final Allele altAllele = sortedResultWithThetaNPriors.getAllelesUsedInGenotyping().get(1);
@ -316,7 +329,8 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
log10PriorsOfAC[1] += sortedResultWithThetaNPriors.getLog10PriorOfAFGT0();
// the AF > 0 case requires us to store the normalized likelihood for later summation
log10LikelihoodsOfACGt0[altI] = sortedResultWithThetaNPriors.getLog10LikelihoodOfAFGT0();
if ( sortedResultWithThetaNPriors.getLog10PosteriorOfAFGT0() > MIN_LOG10_CONFIDENCE_TO_INCLUDE_ALLELE_IN_POSTERIOR )
log10PosteriorOfACEq0Sum += sortedResultWithThetaNPriors.getLog10PosteriorOfAFEq0();
// bind pNonRef for allele to the posterior value of the AF > 0 with the new adjusted prior
log10pNonRefByAllele.put(altAllele, sortedResultWithThetaNPriors.getLog10PosteriorOfAFGT0());
@ -325,14 +339,19 @@ public class IndependentAllelesDiploidExactAFCalc extends DiploidExactAFCalc {
nEvaluations += sortedResultWithThetaNPriors.nEvaluations;
}
// the log10 likelihoods are the sum of the log10 likelihoods across all alt alleles
final double[] log10LikelihoodsOfAC = new double[]{
log10LikelihoodsOfACEq0,
MathUtils.log10sumLog10(log10LikelihoodsOfACGt0)};
// In principle, if B_p = x and C_p = y are the probabilities of being poly for alleles B and C,
// the probability of being poly is (1 - B_p) * (1 - C_p) = (1 - x) * (1 - y). We want to estimate confidently
// log10((1 - x) * (1 - y)) which is log10(1 - x) + log10(1 - y). This sum is log10PosteriorOfACEq0
final double log10PosteriorOfACGt0 = Math.max(Math.log10(1 - Math.pow(10, log10PosteriorOfACEq0Sum)), MathUtils.LOG10_P_OF_ZERO);
final double[] log10LikelihoodsOfAC = new double[] {
// L + prior = posterior => L = poster - prior
log10PosteriorOfACEq0Sum - log10PriorsOfAC[0],
log10PosteriorOfACGt0 - log10PriorsOfAC[1]
};
return new MyAFCalcResult(alleleCountsOfMLE, nEvaluations, vc.getAlleles(),
MathUtils.normalizeFromLog10(log10LikelihoodsOfAC, true, true), // necessary to ensure all values < 0
MathUtils.normalizeFromLog10(log10PriorsOfAC, true), // priors incorporate multiple alt alleles, must be normalized
MathUtils.normalizeFromLog10(log10LikelihoodsOfAC, true), // necessary to ensure all values < 0
MathUtils.normalizeFromLog10(log10PriorsOfAC, true), // priors incorporate multiple alt alleles, must be normalized
log10pNonRefByAllele, sortedResultsWithThetaNPriors);
}
}

6
public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/afcalc/ReferenceDiploidExactAFCalc.java
View File

@ -1,13 +1,7 @@
package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
public class ReferenceDiploidExactAFCalc extends DiploidExactAFCalc {
protected ReferenceDiploidExactAFCalc(int nSamples, int maxAltAlleles, int maxAltAllelesForIndels, final int ploidy) {
super(nSamples, maxAltAlleles, maxAltAllelesForIndels, ploidy);
}
protected StateTracker makeMaxLikelihood(final VariantContext vc, final AFCalcResultTracker resultTracker) {
return new StateTracker();
}
}

4
public/java/src/org/broadinstitute/sting/gatk/walkers/indels/IndelRealigner.java
View File

@ -370,8 +370,6 @@ public class IndelRealigner extends ReadWalker<Integer, Integer> {
currentInterval = intervals.hasNext() ? intervals.next() : null;
writerToUse = writer;
if ( N_WAY_OUT != null ) {
boolean createIndex = true;
@ -383,9 +381,9 @@ public class IndelRealigner extends ReadWalker<Integer, Integer> {
createIndex, generateMD5s,createProgramRecord(),KEEP_ALL_PG_RECORDS);
}
} else {
// set up the output writer
setupWriter(getToolkit().getSAMFileHeader());
writerToUse = writer;
}
manager = new ConstrainedMateFixingManager(writerToUse, getToolkit().getGenomeLocParser(), MAX_ISIZE_FOR_MOVEMENT, MAX_POS_MOVE_ALLOWED, MAX_RECORDS_IN_MEMORY);

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

@ -51,8 +51,8 @@ public class MathUtils {
public static final double[] log10Cache;
public static final double[] log10FactorialCache;
private static final double[] jacobianLogTable;
private static final double JACOBIAN_LOG_TABLE_STEP = 0.001;
private static final double JACOBIAN_LOG_TABLE_INV_STEP = 1.0 / 0.001;
private static final double JACOBIAN_LOG_TABLE_STEP = 0.0001;
private static final double JACOBIAN_LOG_TABLE_INV_STEP = 1.0 / JACOBIAN_LOG_TABLE_STEP;
private static final double MAX_JACOBIAN_TOLERANCE = 8.0;
private static final int JACOBIAN_LOG_TABLE_SIZE = (int) (MAX_JACOBIAN_TOLERANCE / JACOBIAN_LOG_TABLE_STEP) + 1;
private static final int MAXN = 50000;
@ -596,7 +596,6 @@ public class MathUtils {
if (keepInLogSpace) {
for (int i = 0; i < array.length; i++) {
array[i] -= maxValue;
array[i] = Math.max(array[i], LOG10_P_OF_ZERO);
}
return array;
}
@ -613,8 +612,11 @@ public class MathUtils {
sum += normalized[i];
for (int i = 0; i < array.length; i++) {
double x = normalized[i] / sum;
if (takeLog10OfOutput)
x = Math.max(Math.log10(x), LOG10_P_OF_ZERO);
if (takeLog10OfOutput) {
x = Math.log10(x);
if ( x < LOG10_P_OF_ZERO || Double.isInfinite(x) )
x = array[i] - maxValue;
}
normalized[i] = x;
}

7
public/java/src/org/broadinstitute/sting/utils/classloader/JVMUtils.java
View File

@ -32,7 +32,6 @@ import org.reflections.util.ClasspathHelper;
import java.io.File;
import java.io.IOException;
import java.lang.annotation.Annotation;
import java.lang.reflect.*;
import java.net.URL;
import java.util.*;
@ -198,7 +197,7 @@ public class JVMUtils {
* @return the list of class path urls.
*/
public static Set<URL> getClasspathURLs() {
return ClasspathHelper.getUrlsForManifestsCurrentClasspath();
return ClasspathHelper.forManifest();
}
/**
@ -240,8 +239,8 @@ public class JVMUtils {
/**
* Returns a comma-separated list of the names of the interfaces implemented by this class
*
* @param covClass
* @return
* @param covClass class
* @return names of interfaces
*/
public static String classInterfaces(final Class covClass) {
final List<String> interfaces = new ArrayList<String>();

14
public/java/src/org/broadinstitute/sting/utils/classloader/PluginManager.java
View File

@ -25,8 +25,6 @@
package org.broadinstitute.sting.utils.classloader;
import ch.qos.logback.classic.Level;
import ch.qos.logback.classic.Logger;
import org.broadinstitute.sting.gatk.WalkerManager;
import org.broadinstitute.sting.gatk.filters.FilterManager;
import org.broadinstitute.sting.utils.exceptions.DynamicClassResolutionException;
@ -35,7 +33,6 @@ import org.broadinstitute.sting.utils.exceptions.UserException;
import org.reflections.Reflections;
import org.reflections.scanners.SubTypesScanner;
import org.reflections.util.ConfigurationBuilder;
import org.slf4j.LoggerFactory;
import java.io.File;
import java.lang.reflect.Constructor;
@ -57,9 +54,8 @@ public class PluginManager<PluginType> {
private static final Reflections defaultReflections;
static {
// turn off logging in the reflections library - they talk too much (to the wrong logger factory as well, logback)
Logger logger = (ch.qos.logback.classic.Logger) LoggerFactory.getLogger(Reflections.class);
logger.setLevel(Level.OFF);
// turn off logging in the reflections library - they talk too much
Reflections.log = null;
Set<URL> classPathUrls = new LinkedHashSet<URL>();
@ -179,9 +175,9 @@ public class PluginManager<PluginType> {
/**
* Sorts, in place, the list of plugins according to getName() on each element
*
* @param unsortedPlugins
* @param unsortedPlugins unsorted plugins
*/
private final void sortPlugins(final List<Class<? extends PluginType>> unsortedPlugins) {
private void sortPlugins(final List<Class<? extends PluginType>> unsortedPlugins) {
Collections.sort(unsortedPlugins, new ComparePluginsByName());
}
@ -235,7 +231,7 @@ public class PluginManager<PluginType> {
* @param plugin Name of the plugin for which to search.
* @return True if the plugin exists, false otherwise.
*/
public boolean exists(Class<?> plugin) {
public boolean exists(Class<? extends PluginType> plugin) {
return pluginsByName.containsValue(plugin);
}

43
public/java/src/org/broadinstitute/sting/utils/collections/LoggingNestedIntegerArray.java
View File

@ -34,7 +34,12 @@ import java.io.PrintStream;
* to the provided output stream. For testing/debugging purposes.
*
* Log entries are of the following form (fields are tab-separated):
* LABEL VALUE KEY1 KEY2 ... KEY_N
* LABEL OPERATION VALUE KEY1 KEY2 ... KEY_N
*
* A header line is written before the log entries giving the dimensions of this NestedIntegerArray.
* It has the form:
*
* # LABEL SIZE_OF_FIRST_DIMENSION SIZE_OF_SECOND_DIMENSION ... SIZE_OF_NTH_DIMENSION
*
* @author David Roazen
*/
@ -43,6 +48,9 @@ public class LoggingNestedIntegerArray<T> extends NestedIntegerArray<T> {
private PrintStream log;
private String logEntryLabel;
public static final String HEADER_LINE_PREFIX = "# ";
public enum NestedIntegerArrayOperation { GET, PUT };
/**
*
* @param log output stream to which to log update operations
@ -57,6 +65,37 @@ public class LoggingNestedIntegerArray<T> extends NestedIntegerArray<T> {
}
this.log = log;
this.logEntryLabel = logEntryLabel != null ? logEntryLabel : "";
// Write the header line recording the dimensions of this NestedIntegerArray:
StringBuilder logHeaderLine = new StringBuilder();
logHeaderLine.append(HEADER_LINE_PREFIX);
logHeaderLine.append(this.logEntryLabel);
for ( int dimension : dimensions ) {
logHeaderLine.append("\t");
logHeaderLine.append(dimension);
}
this.log.println(logHeaderLine.toString());
}
@Override
public T get( final int... keys ) {
StringBuilder logEntry = new StringBuilder();
logEntry.append(logEntryLabel);
logEntry.append("\t");
logEntry.append(NestedIntegerArrayOperation.GET);
logEntry.append("\t"); // empty field for the datum value
for ( int key : keys ) {
logEntry.append("\t");
logEntry.append(key);
}
log.println(logEntry.toString());
return super.get(keys);
}
@Override
@ -67,6 +106,8 @@ public class LoggingNestedIntegerArray<T> extends NestedIntegerArray<T> {
logEntry.append(logEntryLabel);
logEntry.append("\t");
logEntry.append(NestedIntegerArrayOperation.PUT);
logEntry.append("\t");
logEntry.append(value);
for ( int key : keys ) {
logEntry.append("\t");

18
public/java/src/org/broadinstitute/sting/utils/help/HelpFormatter.java
View File

@ -26,10 +26,7 @@
package org.broadinstitute.sting.utils.help;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.commandline.ArgumentDefinition;
import org.broadinstitute.sting.commandline.ArgumentDefinitionGroup;
import org.broadinstitute.sting.commandline.ArgumentDefinitions;
import org.broadinstitute.sting.commandline.ArgumentMatchSource;
import org.broadinstitute.sting.commandline.*;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.text.TextFormattingUtils;
@ -273,9 +270,9 @@ public class HelpFormatter {
* Generate a standard header for the logger
*
* @param applicationDetails details of the application to run.
* @param parsedArgs the command line arguments passed in
* @param parsedArgs the arguments passed in
*/
public static void generateHeaderInformation(ApplicationDetails applicationDetails, Map<ArgumentMatchSource, List<String>> parsedArgs) {
public static void generateHeaderInformation(ApplicationDetails applicationDetails, Map<ArgumentMatchSource, ParsedArgs> parsedArgs) {
DateFormat dateFormat = new SimpleDateFormat("yyyy/MM/dd HH:mm:ss");
java.util.Date date = new java.util.Date();
@ -286,19 +283,16 @@ public class HelpFormatter {
for (String headerLine : applicationDetails.applicationHeader)
logger.info(headerLine);
logger.debug("Current directory: " + System.getProperty("user.dir"));
for (Map.Entry<ArgumentMatchSource, List<String>> entry: parsedArgs.entrySet()) {
for (Map.Entry<ArgumentMatchSource, ParsedArgs> entry: parsedArgs.entrySet()) {
ArgumentMatchSource matchSource = entry.getKey();
final String sourceName;
switch (matchSource.getType()) {
case CommandLine: sourceName = "Program"; break;
case File: sourceName = matchSource.getFile().getPath(); break;
case Provider: sourceName = matchSource.getDescription(); break;
default: throw new RuntimeException("Unexpected argument match source type: " + matchSource.getType());
}
String output = sourceName + " Args:";
for (String str : entry.getValue()) {
output = output + " " + str;
}
String output = sourceName + " Args: " + entry.getValue().getDescription();
logger.info(output);
}
logger.info("Date/Time: " + dateFormat.format(date));

14
public/java/src/org/broadinstitute/sting/utils/variantcontext/Genotype.java
View File

@ -298,12 +298,16 @@ public abstract class Genotype implements Comparable<Genotype> {
* @return true if all samples PLs are equal and == 0
*/
public boolean isNonInformative() {
for ( final int PL : getPL() ) {
if ( PL != 0 )
return false;
}
if ( getPL() == null )
return true;
else {
for ( final int PL : getPL() ) {
if ( PL != 0 )
return false;
}
return true;
return true;
}
}
/**

2
public/java/test/org/broadinstitute/sting/commandline/ArgumentMatchSiteUnitTest.java
View File

@ -39,7 +39,7 @@ public class ArgumentMatchSiteUnitTest {
@Test
public void testFile() {
ArgumentMatchSource source = new ArgumentMatchSource(new File("test"));
ArgumentMatchSource source = new ArgumentMatchFileSource(new File("test"));
ArgumentMatchSite site = new ArgumentMatchSite(source, 1);
Assert.assertEquals(site.getSource(), source);
Assert.assertEquals(site.getIndex(), 1);

16
public/java/test/org/broadinstitute/sting/commandline/ArgumentMatchSourceUnitTest.java
View File

@ -35,15 +35,15 @@ public class ArgumentMatchSourceUnitTest extends BaseTest {
public void testCommandLine() {
ArgumentMatchSource source = ArgumentMatchSource.COMMAND_LINE;
Assert.assertEquals(source.getType(), ArgumentMatchSourceType.CommandLine);
Assert.assertNull(source.getFile());
Assert.assertNull(source.getDescription());
}
@Test
public void testFile() {
File f = new File("test");
ArgumentMatchSource source = new ArgumentMatchSource(f);
Assert.assertEquals(source.getType(), ArgumentMatchSourceType.File);
Assert.assertEquals(source.getFile(), f);
ArgumentMatchSource source = new ArgumentMatchFileSource(f);
Assert.assertEquals(source.getType(), ArgumentMatchSourceType.Provider);
Assert.assertEquals(source.getDescription(), "file " + f.getAbsolutePath());
}
@Test(expectedExceptions = IllegalArgumentException.class)
@ -54,8 +54,8 @@ public class ArgumentMatchSourceUnitTest extends BaseTest {
@Test
public void testEquals() {
ArgumentMatchSource cmdLine = ArgumentMatchSource.COMMAND_LINE;
ArgumentMatchSource fileA = new ArgumentMatchSource(new File("a"));
ArgumentMatchSource fileB = new ArgumentMatchSource(new File("b"));
ArgumentMatchSource fileA = new ArgumentMatchFileSource(new File("a"));
ArgumentMatchSource fileB = new ArgumentMatchFileSource(new File("b"));
Assert.assertFalse(cmdLine.equals(null));
@ -75,8 +75,8 @@ public class ArgumentMatchSourceUnitTest extends BaseTest {
@Test
public void testCompareTo() {
ArgumentMatchSource cmdLine = ArgumentMatchSource.COMMAND_LINE;
ArgumentMatchSource fileA = new ArgumentMatchSource(new File("a"));
ArgumentMatchSource fileB = new ArgumentMatchSource(new File("b"));
ArgumentMatchSource fileA = new ArgumentMatchFileSource(new File("a"));
ArgumentMatchSource fileB = new ArgumentMatchFileSource(new File("b"));
Assert.assertTrue(cmdLine.compareTo(cmdLine) == 0);
Assert.assertTrue(cmdLine.compareTo(fileA) < 0);

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

@ -28,7 +28,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
public void testMultiSamplePilot1() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
baseCommand + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10,022,000-10,025,000", 1,
Arrays.asList("fe9c0e9e4b4ee4677145748cdd2285ff"));
Arrays.asList("b3abf320f7d02d0e3b2883833419130e"));
executeTest("test MultiSample Pilot1", spec);
}
@ -60,7 +60,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
public void testMultipleSNPAlleles() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper -R " + b37KGReference + " -nosl --no_cmdline_in_header -glm BOTH --dbsnp " + b37dbSNP129 + " -I " + privateTestDir + "multiallelic.snps.bam -o %s -L " + privateTestDir + "multiallelic.snps.intervals", 1,
Arrays.asList("772e14d8c908044c04053d204bad69ef"));
Arrays.asList("26af30187316f742878c85f0ed091837"));
executeTest("test Multiple SNP alleles", spec);
}
@ -76,7 +76,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
public void testReverseTrim() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper -R " + b37KGReference + " -nosl --no_cmdline_in_header -glm INDEL -I " + validationDataLocation + "CEUTrio.HiSeq.b37.chr20.10_11mb.bam -o %s -L 20:10289124 -L 20:10090289", 1,
Arrays.asList("1fb69aa3857e921191997daa73f1b687"));
Arrays.asList("aa9cf96ab8f5aa844387e3aef1f27249"));
executeTest("test reverse trim", spec);
}
@ -249,7 +249,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
" -o %s" +
" -L 1:10,000,000-10,100,000",
1,
Arrays.asList("5d19e3077e0cabbb364f68676a09ebe0"));
Arrays.asList("04a87b87ee4323eba853c78f25551d1a"));
executeTest(String.format("test multiple technologies"), spec);
}
@ -268,7 +268,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
" -L 1:10,000,000-10,100,000" +
" -baq CALCULATE_AS_NECESSARY",
1,
Arrays.asList("8a1931095f70523ad11cb99b30df7b84"));
Arrays.asList("950fb032cc9902ae48bd21f272d2fd52"));
executeTest(String.format("test calling with BAQ"), spec);
}
@ -287,7 +287,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
" -o %s" +
" -L 1:10,000,000-10,500,000",
1,
Arrays.asList("64a491b5276fd5d1cd04260ea3e63cf7"));
Arrays.asList("b3df138254ed141b61a758df87757e0d"));
executeTest(String.format("test indel caller in SLX"), spec);
}
@ -302,7 +302,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
" -minIndelCnt 1" +
" -L 1:10,000,000-10,100,000",
1,
Arrays.asList("f63a8b8061e6c5999408d34798061895"));
Arrays.asList("63fd9488daadd4baaef0a98f02916996"));
executeTest(String.format("test indel caller in SLX with low min allele count"), spec);
}
@ -315,7 +315,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
" -o %s" +
" -L 1:10,000,000-10,500,000",
1,
Arrays.asList("c9d684ff2f2a9083480db6e962d612a9"));
Arrays.asList("52b5a432092995c92fe71e1942689ba8"));
executeTest(String.format("test indel calling, multiple technologies"), spec);
}
@ -325,7 +325,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
baseCommandIndels + " --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles " + privateTestDir + "indelAllelesForUG.vcf -I " + validationDataLocation +
"pilot2_daughters.chr20.10k-11k.bam -o %s -L 20:10,000,000-10,100,000", 1,
Arrays.asList("833fd97c6f32d7af6c9c088a78e51f68"));
Arrays.asList("7e3f67bf371112be5dbadb4fe6faa52a"));
executeTest("test MultiSample Pilot2 indels with alleles passed in", spec);
}
@ -335,7 +335,7 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
baseCommandIndels + " --output_mode EMIT_ALL_SITES --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles "
+ privateTestDir + "indelAllelesForUG.vcf -I " + validationDataLocation +
"pilot2_daughters.chr20.10k-11k.bam -o %s -L 20:10,000,000-10,100,000", 1,
Arrays.asList("95b73c24c68dc475516571d9f49dfb1e"));
Arrays.asList("bc31c4977cb7e563ddf9c8dea27f3f4f"));
executeTest("test MultiSample Pilot2 indels with alleles passed in and emitting all sites", spec);
}
@ -343,13 +343,13 @@ public class UnifiedGenotyperIntegrationTest extends WalkerTest {
public void testMultiSampleIndels1() {
WalkerTest.WalkerTestSpec spec1 = new WalkerTest.WalkerTestSpec(
baseCommandIndels + " -I " + validationDataLocation + "low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10450700-10551000", 1,
Arrays.asList("3bdbf48de30bac58f3bcbc5bf3aa63aa"));
Arrays.asList("7fc488fe16dea9f023bfcfdaa908a548"));
List<File> result = executeTest("test MultiSample Pilot1 CEU indels", spec1).getFirst();
WalkerTest.WalkerTestSpec spec2 = new WalkerTest.WalkerTestSpec(
baseCommandIndels + " --genotyping_mode GENOTYPE_GIVEN_ALLELES -alleles " + result.get(0).getAbsolutePath() + " -I " + validationDataLocation +
"low_coverage_CEU.chr1.10k-11k.bam -o %s -L 1:10450700-10551000", 1,
Arrays.asList("beee9457d7cea42006ac45400db5e873"));
Arrays.asList("f3ff7fe0f15f31eadd726c711d6bf3de"));
executeTest("test MultiSample Pilot1 CEU indels using GENOTYPE_GIVEN_ALLELES", spec2);
}

10
public/java/test/org/broadinstitute/sting/gatk/walkers/indels/IndelRealignerIntegrationTest.java
View File

@ -113,4 +113,14 @@ public class IndelRealignerIntegrationTest extends WalkerTest {
executeTest(String.format("realigner [%s]", entry.getKey()), spec);
}
}
@Test
public void testNWayOut() {
WalkerTestSpec spec1 = new WalkerTestSpec(
baseCommandPrefix + " -nWayOut .clean.bam ",
1,
Arrays.asList("d41d8cd98f00b204e9800998ecf8427e"));
executeTest("test realigner nWayOut", spec1);
}
}

153
public/java/test/org/broadinstitute/sting/utils/MathUtilsUnitTest.java
View File

@ -225,65 +225,67 @@ public class MathUtilsUnitTest extends BaseTest {
@Test
public void testApproximateLog10SumLog10() {
final double requiredPrecision = 1E-4;
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0}), 0.0, 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-5.15}), -5.15, 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {130.0}), 130.0, 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-0.145}), -0.145, 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0}), 0.0, requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-5.15}), -5.15, requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {130.0}), 130.0, requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-0.145}), -0.145, requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(0.0, 0.0), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-1.0, 0.0), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(0.0, -1.0), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-2.2, -3.5), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-1.0, -7.1), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(5.0, 6.2), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(38.1, 16.2), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-38.1, 6.2), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-19.1, -37.1), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-29.1, -27.6), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-0.12345, -0.23456), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-15.7654, -17.0101), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-0.12345, Double.NEGATIVE_INFINITY), -0.12345, 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-15.7654, Double.NEGATIVE_INFINITY), -15.7654, 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(0.0, 0.0), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-1.0, 0.0), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(0.0, -1.0), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-2.2, -3.5), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-1.0, -7.1), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(5.0, 6.2), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(38.1, 16.2), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-38.1, 6.2), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-19.1, -37.1), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-29.1, -27.6), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-0.12345, -0.23456), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-15.7654, -17.0101), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-0.12345, Double.NEGATIVE_INFINITY), -0.12345, requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-15.7654, Double.NEGATIVE_INFINITY), -15.7654, requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0, 0.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-1.0, 0.0}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0, -1.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-2.2, -3.5}), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-1.0, -7.1}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {5.0, 6.2}), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {38.1, 16.2}), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-38.1, 6.2}), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-19.1, -37.1}), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-29.1, -27.6}), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-0.12345, -0.23456}), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-15.7654, -17.0101}), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0, 0.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-1.0, 0.0}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0, -1.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-2.2, -3.5}), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-1.0, -7.1}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {5.0, 6.2}), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {38.1, 16.2}), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-38.1, 6.2}), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-19.1, -37.1}), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-29.1, -27.6}), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-0.12345, -0.23456}), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-15.7654, -17.0101}), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0, 0.0, 0.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-1.0, 0.0, 0.0}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0, -1.0, -2.5}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0) + Math.pow(10.0, -2.5)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-2.2, -3.5, -1.1}), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5) + Math.pow(10.0, -1.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-1.0, -7.1, 0.5}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1) + Math.pow(10.0, 0.5)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {5.0, 6.2, 1.3}), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2) + Math.pow(10.0, 1.3)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {38.1, 16.2, 18.1}), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2) + Math.pow(10.0, 18.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-38.1, 6.2, 26.6}), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2) + Math.pow(10.0, 26.6)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-19.1, -37.1, -45.1}), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1) + Math.pow(10.0, -45.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-29.1, -27.6, -26.2}), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6) + Math.pow(10.0, -26.2)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-0.12345, -0.23456, -0.34567}), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456) + Math.pow(10.0, -0.34567)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-15.7654, -17.0101, -17.9341}), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101) + Math.pow(10.0, -17.9341)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0, 0.0, 0.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-1.0, 0.0, 0.0}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {0.0, -1.0, -2.5}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0) + Math.pow(10.0, -2.5)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-2.2, -3.5, -1.1}), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5) + Math.pow(10.0, -1.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-1.0, -7.1, 0.5}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1) + Math.pow(10.0, 0.5)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {5.0, 6.2, 1.3}), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2) + Math.pow(10.0, 1.3)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {38.1, 16.2, 18.1}), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2) + Math.pow(10.0, 18.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-38.1, 6.2, 26.6}), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2) + Math.pow(10.0, 26.6)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-19.1, -37.1, -45.1}), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1) + Math.pow(10.0, -45.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-29.1, -27.6, -26.2}), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6) + Math.pow(10.0, -26.2)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-0.12345, -0.23456, -0.34567}), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456) + Math.pow(10.0, -0.34567)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(new double[] {-15.7654, -17.0101, -17.9341}), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101) + Math.pow(10.0, -17.9341)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(0.0, 0.0, 0.0), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-1.0, 0.0, 0.0), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(0.0, -1.0, -2.5), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0) + Math.pow(10.0, -2.5)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-2.2, -3.5, -1.1), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5) + Math.pow(10.0, -1.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-1.0, -7.1, 0.5), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1) + Math.pow(10.0, 0.5)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(5.0, 6.2, 1.3), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2) + Math.pow(10.0, 1.3)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(38.1, 16.2, 18.1), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2) + Math.pow(10.0, 18.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-38.1, 6.2, 26.6), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2) + Math.pow(10.0, 26.6)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-19.1, -37.1, -45.1), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1) + Math.pow(10.0, -45.1)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-29.1, -27.6, -26.2), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6) + Math.pow(10.0, -26.2)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-0.12345, -0.23456, -0.34567), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456) + Math.pow(10.0, -0.34567)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-15.7654, -17.0101, -17.9341), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101) + Math.pow(10.0, -17.9341)), 1e-3);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(0.0, 0.0, 0.0), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-1.0, 0.0, 0.0), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(0.0, -1.0, -2.5), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0) + Math.pow(10.0, -2.5)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-2.2, -3.5, -1.1), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5) + Math.pow(10.0, -1.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-1.0, -7.1, 0.5), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1) + Math.pow(10.0, 0.5)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(5.0, 6.2, 1.3), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2) + Math.pow(10.0, 1.3)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(38.1, 16.2, 18.1), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2) + Math.pow(10.0, 18.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-38.1, 6.2, 26.6), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2) + Math.pow(10.0, 26.6)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-19.1, -37.1, -45.1), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1) + Math.pow(10.0, -45.1)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-29.1, -27.6, -26.2), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6) + Math.pow(10.0, -26.2)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-0.12345, -0.23456, -0.34567), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456) + Math.pow(10.0, -0.34567)), requiredPrecision);
Assert.assertEquals(MathUtils.approximateLog10SumLog10(-15.7654, -17.0101, -17.9341), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101) + Math.pow(10.0, -17.9341)), requiredPrecision);
}
@Test
@ -299,14 +301,47 @@ public class MathUtilsUnitTest extends BaseTest {
@Test
public void testLog10sumLog10() {
final double requiredPrecision = 1E-14;
final double log3 = 0.477121254719662;
Assert.assertEquals(MathUtils.compareDoubles(MathUtils.log10sumLog10(new double[]{0.0, 0.0, 0.0}), log3), 0);
Assert.assertEquals(MathUtils.compareDoubles(MathUtils.log10sumLog10(new double[] {0.0, 0.0, 0.0}, 0), log3), 0);
Assert.assertEquals(MathUtils.compareDoubles(MathUtils.log10sumLog10(new double[]{0.0, 0.0, 0.0}, 0, 3), log3), 0);
Assert.assertEquals(MathUtils.log10sumLog10(new double[]{0.0, 0.0, 0.0}), log3, requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {0.0, 0.0, 0.0}, 0), log3, requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[]{0.0, 0.0, 0.0}, 0, 3), log3, requiredPrecision);
final double log2 = 0.301029995663981;
Assert.assertEquals(MathUtils.compareDoubles(MathUtils.log10sumLog10(new double[] {0.0, 0.0, 0.0}, 0, 2), log2), 0);
Assert.assertEquals(MathUtils.compareDoubles(MathUtils.log10sumLog10(new double[] {0.0, 0.0, 0.0}, 0, 1), 0.0), 0);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {0.0, 0.0, 0.0}, 0, 2), log2, requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {0.0, 0.0, 0.0}, 0, 1), 0.0, requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {0.0}), 0.0, requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-5.15}), -5.15, requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {130.0}), 130.0, requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-0.145}), -0.145, requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {0.0, 0.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-1.0, 0.0}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {0.0, -1.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-2.2, -3.5}), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-1.0, -7.1}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {5.0, 6.2}), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {38.1, 16.2}), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-38.1, 6.2}), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-19.1, -37.1}), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-29.1, -27.6}), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-0.12345, -0.23456}), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-15.7654, -17.0101}), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {0.0, 0.0, 0.0}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-1.0, 0.0, 0.0}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, 0.0) + Math.pow(10.0, 0.0)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {0.0, -1.0, -2.5}), Math.log10(Math.pow(10.0, 0.0) + Math.pow(10.0, -1.0) + Math.pow(10.0, -2.5)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-2.2, -3.5, -1.1}), Math.log10(Math.pow(10.0, -2.2) + Math.pow(10.0, -3.5) + Math.pow(10.0, -1.1)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-1.0, -7.1, 0.5}), Math.log10(Math.pow(10.0, -1.0) + Math.pow(10.0, -7.1) + Math.pow(10.0, 0.5)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {5.0, 6.2, 1.3}), Math.log10(Math.pow(10.0, 5.0) + Math.pow(10.0, 6.2) + Math.pow(10.0, 1.3)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {38.1, 16.2, 18.1}), Math.log10(Math.pow(10.0, 38.1) + Math.pow(10.0, 16.2) + Math.pow(10.0, 18.1)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-38.1, 6.2, 26.6}), Math.log10(Math.pow(10.0, -38.1) + Math.pow(10.0, 6.2) + Math.pow(10.0, 26.6)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-19.1, -37.1, -45.1}), Math.log10(Math.pow(10.0, -19.1) + Math.pow(10.0, -37.1) + Math.pow(10.0, -45.1)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-29.1, -27.6, -26.2}), Math.log10(Math.pow(10.0, -29.1) + Math.pow(10.0, -27.6) + Math.pow(10.0, -26.2)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-0.12345, -0.23456, -0.34567}), Math.log10(Math.pow(10.0, -0.12345) + Math.pow(10.0, -0.23456) + Math.pow(10.0, -0.34567)), requiredPrecision);
Assert.assertEquals(MathUtils.log10sumLog10(new double[] {-15.7654, -17.0101, -17.9341}), Math.log10(Math.pow(10.0, -15.7654) + Math.pow(10.0, -17.0101) + Math.pow(10.0, -17.9341)), requiredPrecision);
}
@Test

2
public/java/test/org/broadinstitute/sting/utils/nanoScheduler/NanoSchedulerIntegrationTest.java
View File

@ -21,7 +21,7 @@ public class NanoSchedulerIntegrationTest extends WalkerTest {
for ( final int nct : Arrays.asList(1, 2) ) {
// tests.add(new Object[]{ "SNP", "a1c7546f32a8919a3f3a70a04b2e8322", nt, nct });
//// tests.add(new Object[]{ "INDEL", "0a6d2be79f4f8a4b0eb788cc4751b31b", nt, nct });
tests.add(new Object[]{ "BOTH", "78ce72d8f9d029313f5f2ceb02bb9822", nt, nct });
tests.add(new Object[]{ "BOTH", "8cad82c3a5f5b932042933f136663c8a", nt, nct });
}
return tests.toArray(new Object[][]{});

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

@ -28,7 +28,7 @@ import function.QFunction
import java.io.File
import org.broadinstitute.sting.commandline._
import org.broadinstitute.sting.queue.util._
import org.broadinstitute.sting.queue.engine.{QGraphSettings, QGraph}
import org.broadinstitute.sting.queue.engine.{QStatusMessenger, QGraphSettings, QGraph}
import collection.JavaConversions._
import org.broadinstitute.sting.utils.classloader.PluginManager
import org.broadinstitute.sting.utils.exceptions.UserException
@ -90,29 +90,46 @@ class QCommandLine extends CommandLineProgram with Logging {
private var qScriptClasses: File = _
private var shuttingDown = false
private lazy val pluginManager = {
private lazy val qScriptPluginManager = {
qScriptClasses = IOUtils.tempDir("Q-Classes-", "", settings.qSettings.tempDirectory)
qScriptManager.loadScripts(scripts, qScriptClasses)
new PluginManager[QScript](classOf[QScript], Seq(qScriptClasses.toURI.toURL))
}
QFunction.parsingEngine = new ParsingEngine(this)
private lazy val qStatusMessengerPluginManager = {
new PluginManager[QStatusMessenger](classOf[QStatusMessenger])
}
ClassFieldCache.parsingEngine = new ParsingEngine(this)
/**
* Takes the QScripts passed in, runs their script() methods, retrieves their generated
* functions, and then builds and runs a QGraph based on the dependencies.
*/
def execute = {
val allStatusMessengers = qStatusMessengerPluginManager.createAllTypes()
if (settings.qSettings.runName == null)
settings.qSettings.runName = FilenameUtils.removeExtension(scripts.head.getName)
if (IOUtils.isDefaultTempDir(settings.qSettings.tempDirectory))
settings.qSettings.tempDirectory = IOUtils.absolute(settings.qSettings.runDirectory, ".queue/tmp")
qGraph.initializeWithSettings(settings)
val allQScripts = pluginManager.createAllTypes()
for (statusMessenger <- allStatusMessengers) {
loadArgumentsIntoObject(statusMessenger)
}
for (statusMessenger <- allStatusMessengers) {
statusMessenger.started()
}
val allQScripts = qScriptPluginManager.createAllTypes()
for (script <- allQScripts) {
logger.info("Scripting " + pluginManager.getName(script.getClass.asSubclass(classOf[QScript])))
logger.info("Scripting " + qScriptPluginManager.getName(script.getClass.asSubclass(classOf[QScript])))
loadArgumentsIntoObject(script)
// TODO: Pulling inputs can be time/io expensive! Some scripts are using the files to generate functions-- even for dry runs-- so pull it all down for now.
//if (settings.run)
script.pullInputs()
script.qSettings = settings.qSettings
try {
script.script()
@ -142,11 +159,18 @@ class QCommandLine extends CommandLineProgram with Logging {
logger.info("Writing final jobs report...")
qGraph.writeJobsReport()
if (!qGraph.success) {
if (!success) {
logger.info("Done with errors")
qGraph.logFailed()
for (statusMessenger <- allStatusMessengers)
statusMessenger.exit("Done with errors")
1
} else {
if (settings.run) {
allQScripts.foreach(_.pushOutputs())
for (statusMessenger <- allStatusMessengers)
statusMessenger.done(allQScripts.map(_.remoteOutputs))
}
0
}
}
@ -158,19 +182,30 @@ class QCommandLine extends CommandLineProgram with Logging {
override def canAddArgumentsDynamically = true
/**
* Returns the list of QScripts passed in via -S so that their
* arguments can be inspected before QScript.script is called.
* @return Array of QScripts passed in.
* Returns the list of QScripts passed in via -S and other plugins
* so that their arguments can be inspected before QScript.script is called.
* @return Array of dynamic sources
*/
override def getArgumentSources =
pluginManager.getPlugins.toIterable.toArray.asInstanceOf[Array[Class[_]]]
override def getArgumentSources = {
var plugins = Seq.empty[Class[_]]
plugins ++= qScriptPluginManager.getPlugins
plugins ++= qStatusMessengerPluginManager.getPlugins
plugins.toArray
}
/**
* Returns the name of a QScript
* @return The name of a QScript
* Returns the name of a script/plugin
* @return The name of a script/plugin
*/
override def getArgumentSourceName(source: Class[_]) =
pluginManager.getName(source.asSubclass(classOf[QScript]))
override def getArgumentSourceName(source: Class[_]) = {
if (classOf[QScript].isAssignableFrom(source))
qScriptPluginManager.getName(source.asSubclass(classOf[QScript]))
else if (classOf[QStatusMessenger].isAssignableFrom(source))
qStatusMessengerPluginManager.getName(source.asSubclass(classOf[QStatusMessenger]))
else
null
}
/**
* Returns a ScalaCompoundArgumentTypeDescriptor that can parse argument sources into scala collections.

34
public/scala/src/org/broadinstitute/sting/queue/QScript.scala
View File

@ -27,7 +27,8 @@ package org.broadinstitute.sting.queue
import engine.JobRunInfo
import org.broadinstitute.sting.queue.function.QFunction
import annotation.target.field
import util.{StringFileConversions, PrimitiveOptionConversions, Logging}
import util._
import org.broadinstitute.sting.commandline.ArgumentSource
/**
* Defines a Queue pipeline as a collection of CommandLineFunctions.
@ -106,6 +107,37 @@ trait QScript extends Logging with PrimitiveOptionConversions with StringFileCon
def addAll(functions: Seq[QFunction]) {
functions.foreach( f => add(f) )
}
def pullInputs() {
val inputs = ClassFieldCache.getFieldFiles(this, inputFields)
for (remoteFile <- filterRemoteFiles(inputs)) {
logger.info("Pulling %s from %s".format(remoteFile.getAbsolutePath, remoteFile.remoteDescription))
remoteFile.pullToLocal()
}
}
def pushOutputs() {
val outputs = ClassFieldCache.getFieldFiles(this, outputFields)
for (remoteFile <- filterRemoteFiles(outputs)) {
logger.info("Pushing %s to %s".format(remoteFile.getAbsolutePath, remoteFile.remoteDescription))
remoteFile.pushToRemote()
}
}
def remoteOutputs: Map[ArgumentSource, Seq[RemoteFile]] =
outputFields.map(field => (field -> filterRemoteFiles(ClassFieldCache.getFieldFiles(this, field)))).filter(tuple => !tuple._2.isEmpty).toMap
private def filterRemoteFiles(fields: Seq[File]): Seq[RemoteFile] =
fields.filter(field => field != null && field.isInstanceOf[RemoteFile]).map(_.asInstanceOf[RemoteFile])
/** The complete list of fields. */
def functionFields: Seq[ArgumentSource] = ClassFieldCache.classFunctionFields(this.getClass)
/** The @Input fields. */
def inputFields: Seq[ArgumentSource] = ClassFieldCache.classInputFields(this.getClass)
/** The @Output fields. */
def outputFields: Seq[ArgumentSource] = ClassFieldCache.classOutputFields(this.getClass)
/** The @Argument fields. */
def argumentFields: Seq[ArgumentSource] = ClassFieldCache.classArgumentFields(this.getClass)
}
object QScript {

7
public/scala/src/org/broadinstitute/sting/queue/QScriptManager.scala
View File

@ -11,6 +11,7 @@ import org.apache.log4j.Level
import scala.tools.nsc.util.{FakePos, NoPosition, Position}
import org.broadinstitute.sting.queue.util.TextFormatUtils._
import org.broadinstitute.sting.utils.classloader.JVMUtils
import tools.util.StringOps
/**
* Plugin manager for QScripts which loads QScripts into the current class loader.
@ -63,7 +64,7 @@ object QScriptManager extends Logging {
* Heavily based on scala/src/compiler/scala/tools/nsc/reporters/ConsoleReporter.scala
*/
private class Log4JReporter(val settings: Settings) extends AbstractReporter {
def displayPrompt { throw new UnsupportedOperationException("Unable to prompt the user. Prompting should be off.") }
def displayPrompt() { throw new UnsupportedOperationException("Unable to prompt the user. Prompting should be off.") }
/**
* Displays the message at position with severity.
@ -98,9 +99,9 @@ object QScriptManager extends Logging {
*/
def printSummary() {
if (WARNING.count > 0)
printMessage(Level.WARN, countElementsAsString(WARNING.count, "warning") + " found")
printMessage(Level.WARN, StringOps.countElementsAsString(WARNING.count, "warning") + " found")
if (ERROR.count > 0)
printMessage(Level.ERROR, countElementsAsString(ERROR.count, "error") + " found")
printMessage(Level.ERROR, StringOps.countElementsAsString(ERROR.count, "error") + " found")
}
/**

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

@ -185,7 +185,7 @@ class FunctionEdge(val function: QFunction, val inputs: QNode, val outputs: QNod
val tailLines = IOUtils.tail(errorFile, maxLines)
val nl = "%n".format()
val summary = if (tailLines.size > maxLines) "Last %d lines".format(maxLines) else "Contents"
this.function.jobErrorLines = collection.JavaConversions.asScalaIterable(tailLines).toSeq
this.function.jobErrorLines = collection.JavaConversions.collectionAsScalaIterable(tailLines).toSeq
logger.error("%s of %s:%n%s".format(summary, errorFile, StringUtils.join(tailLines, nl)))
} else {
logger.error("Unable to access log file: %s".format(errorFile))

13
public/scala/src/org/broadinstitute/sting/queue/engine/QStatusMessenger.scala 100644
View File

@ -0,0 +1,13 @@
package org.broadinstitute.sting.queue.engine
import org.broadinstitute.sting.commandline.ArgumentSource
import org.broadinstitute.sting.queue.util.RemoteFile
/**
* Plugin to sends QStatus messages
*/
trait QStatusMessenger {
def started()
def done(files: Seq[Map[ArgumentSource, Seq[RemoteFile]]])
def exit(message: String)
}

407
public/scala/src/org/broadinstitute/sting/queue/extensions/cancer/MuTect.scala 100644
View File

File diff suppressed because one or more lines are too long

7
public/scala/src/org/broadinstitute/sting/queue/extensions/gatk/BamGatherFunction.scala
View File

@ -28,6 +28,7 @@ import org.broadinstitute.sting.queue.function.scattergather.GatherFunction
import org.broadinstitute.sting.queue.extensions.picard.PicardBamFunction
import org.broadinstitute.sting.queue.function.{RetryMemoryLimit, QFunction}
import org.broadinstitute.sting.gatk.io.stubs.SAMFileWriterArgumentTypeDescriptor
import org.broadinstitute.sting.queue.util.ClassFieldCache
/**
* Merges BAM files using net.sf.picard.sam.MergeSamFiles.
@ -47,13 +48,13 @@ class BamGatherFunction extends GatherFunction with PicardBamFunction with Retry
// bam_compression and index_output_bam_on_the_fly from SAMFileWriterArgumentTypeDescriptor
// are added by the GATKExtensionsGenerator to the subclass of CommandLineGATK
val compression = QFunction.findField(originalFunction.getClass, SAMFileWriterArgumentTypeDescriptor.COMPRESSION_FULLNAME)
val compression = ClassFieldCache.findField(originalFunction.getClass, SAMFileWriterArgumentTypeDescriptor.COMPRESSION_FULLNAME)
this.compressionLevel = originalGATK.getFieldValue(compression).asInstanceOf[Option[Int]]
val disableIndex = QFunction.findField(originalFunction.getClass, SAMFileWriterArgumentTypeDescriptor.DISABLE_INDEXING_FULLNAME)
val disableIndex = ClassFieldCache.findField(originalFunction.getClass, SAMFileWriterArgumentTypeDescriptor.DISABLE_INDEXING_FULLNAME)
this.createIndex = Some(!originalGATK.getFieldValue(disableIndex).asInstanceOf[Boolean])
val enableMD5 = QFunction.findField(originalFunction.getClass, SAMFileWriterArgumentTypeDescriptor.ENABLE_MD5_FULLNAME)
val enableMD5 = ClassFieldCache.findField(originalFunction.getClass, SAMFileWriterArgumentTypeDescriptor.ENABLE_MD5_FULLNAME)
this.createMD5 = Some(originalGATK.getFieldValue(enableMD5).asInstanceOf[Boolean])
super.freezeFieldValues()

2
public/scala/src/org/broadinstitute/sting/queue/extensions/gatk/GATKIntervals.scala
View File

@ -92,6 +92,6 @@ object GATKIntervals {
}
private def createBinding(interval: String, argumentName: String, tags: Tags): IntervalBinding[Feature] = {
ArgumentTypeDescriptor.parseBinding(interval, classOf[Feature], classOf[IntervalBinding[Feature]], argumentName, tags, argumentName).asInstanceOf[IntervalBinding[Feature]]
ArgumentTypeDescriptor.parseBinding(new ArgumentMatchStringValue(interval), classOf[Feature], classOf[IntervalBinding[Feature]], argumentName, tags, argumentName).asInstanceOf[IntervalBinding[Feature]]
}
}

5
public/scala/src/org/broadinstitute/sting/queue/extensions/gatk/VcfGatherFunction.scala
View File

@ -27,6 +27,7 @@ package org.broadinstitute.sting.queue.extensions.gatk
import org.broadinstitute.sting.queue.function.scattergather.GatherFunction
import org.broadinstitute.sting.queue.function.{RetryMemoryLimit, QFunction}
import org.broadinstitute.sting.gatk.io.stubs.VCFWriterArgumentTypeDescriptor
import org.broadinstitute.sting.queue.util.ClassFieldCache
/**
* Merges a vcf text file.
@ -46,10 +47,10 @@ class VcfGatherFunction extends CombineVariants with GatherFunction with RetryMe
// NO_HEADER and sites_only from VCFWriterArgumentTypeDescriptor
// are added by the GATKExtensionsGenerator to the subclass of CommandLineGATK
val noHeader = QFunction.findField(originalFunction.getClass, VCFWriterArgumentTypeDescriptor.NO_HEADER_ARG_NAME)
val noHeader = ClassFieldCache.findField(originalFunction.getClass, VCFWriterArgumentTypeDescriptor.NO_HEADER_ARG_NAME)
this.no_cmdline_in_header = originalGATK.getFieldValue(noHeader).asInstanceOf[Boolean]
val sitesOnly = QFunction.findField(originalFunction.getClass, VCFWriterArgumentTypeDescriptor.SITES_ONLY_ARG_NAME)
val sitesOnly = ClassFieldCache.findField(originalFunction.getClass, VCFWriterArgumentTypeDescriptor.SITES_ONLY_ARG_NAME)
this.sites_only = originalGATK.getFieldValue(sitesOnly).asInstanceOf[Boolean]
// ensure that the gather function receives the same unsafe parameter as the scattered function

60
public/scala/src/org/broadinstitute/sting/queue/extensions/picard/CalculateHsMetrics.scala 100644
View File

@ -0,0 +1,60 @@
package org.broadinstitute.sting.queue.extensions.picard
import org.broadinstitute.sting.commandline.{Argument, Output, Input}
import java.io.File
/**
* Created with IntelliJ IDEA.
* User: delangel
* Date: 10/9/12
* Time: 5:59 PM
* To change this template use File | Settings | File Templates.
*/
class CalculateHsMetrics extends org.broadinstitute.sting.queue.function.JavaCommandLineFunction with PicardBamFunction {
analysisName = "CalculateHsMetrics"
javaMainClass = "net.sf.picard.sam.CalculateHsMetrics"
@Input(doc="The input SAM or BAM files to analyze. Must be coordinate sorted.", shortName = "input", fullName = "input_bam_files", required = true)
var input: Seq[File] = Nil
@Output(doc="The output file to write statistics to", shortName = "output", fullName = "output_file", required = true)
var output: File = _
@Argument(doc="Interval list with targets", shortName = "targets", fullName = "target_list", required = true)
var targets: File = _
@Argument(doc="Interval list with baits", shortName = "baits", fullName = "bait_list", required = true)
var baits: File = _
@Argument(doc="Reference file", shortName = "reference", fullName = "reference", required = true)
var reference: File = _
/*
@Argument(doc = "Maximum number of file handles to keep open when spilling read ends to disk. Set this number a little lower than the per-process maximum number of file that may be open. This number can be found by executing the 'ulimit -n' command on a Unix system.", shortName = "max_file_handles", fullName ="max_file_handles_for_read_ends_maps", required=false)
var MAX_FILE_HANDLES_FOR_READ_ENDS_MAP: Int = -1;
@Argument(doc = "This number, plus the maximum RAM available to the JVM, determine the memory footprint used by some of the sorting collections. If you are running out of memory, try reducing this number.", shortName = "sorting_ratio", fullName = "sorting_collection_size_ratio", required = false)
var SORTING_COLLECTION_SIZE_RATIO: Double = -1
*/
override def freezeFieldValues() {
super.freezeFieldValues()
// if (outputIndex == null && output != null)
// outputIndex = new File(output.getName.stripSuffix(".bam") + ".bai")
}
val level = "SAMPLE"
override def inputBams = input
override def outputBam = output
//this.sortOrder = null
//this.createIndex = Some(true)
override def commandLine = super.commandLine +
required("BAIT_INTERVALS=" + baits) +
required("TARGET_INTERVALS=" + targets) +
required("REFERENCE_SEQUENCE=" + reference) +
optional("METRIC_ACCUMULATION_LEVEL="+level)/*+
conditional(REMOVE_DUPLICATES, "REMOVE_DUPLICATES=true") +
conditional(MAX_FILE_HANDLES_FOR_READ_ENDS_MAP > 0, "MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=" + MAX_FILE_HANDLES_FOR_READ_ENDS_MAP.toString) +
conditional(SORTING_COLLECTION_SIZE_RATIO > 0, "SORTING_COLLECTION_SIZE_RATIO=" + SORTING_COLLECTION_SIZE_RATIO.toString) */
}

32
public/scala/src/org/broadinstitute/sting/queue/extensions/picard/CollectGcBiasMetrics.scala 100644
View File

@ -0,0 +1,32 @@
package org.broadinstitute.sting.queue.extensions.picard
import org.broadinstitute.sting.commandline.{Argument, Output, Input}
import java.io.File
/**
* Created with IntelliJ IDEA.
* User: delangel
* Date: 10/10/12
* Time: 10:37 AM
* To change this template use File | Settings | File Templates.
*/
class CollectGcBiasMetrics extends org.broadinstitute.sting.queue.function.JavaCommandLineFunction with PicardBamFunction {
analysisName = "CalculateGcMetrics"
javaMainClass = "net.sf.picard.sam.CalculateGcMetrics"
@Input(doc="The input SAM or BAM files to analyze. Must be coordinate sorted.", shortName = "input", fullName = "input_bam_files", required = true)
var input: Seq[File] = Nil
@Output(doc="The output file to write statistics to", shortName = "output", fullName = "output_file", required = true)
var output: File = _
@Argument(doc="Reference file", shortName = "reference", fullName = "reference", required = true)
var reference: File = _
override def inputBams = input
override def outputBam = output
override def commandLine = super.commandLine +
required("CHART_OUTPUT=" + output+".pdf") +
required("REFERENCE_SEQUENCE=" + reference) +
required("ASSUME_SORTED=true")
}

36
public/scala/src/org/broadinstitute/sting/queue/extensions/picard/CollectMultipleMetrics.scala 100644
View File

@ -0,0 +1,36 @@
package org.broadinstitute.sting.queue.extensions.picard
import org.broadinstitute.sting.commandline.{Argument, Output, Input}
import java.io.File
/**
* Created with IntelliJ IDEA.
* User: delangel
* Date: 10/10/12
* Time: 10:37 AM
* To change this template use File | Settings | File Templates.
*/
class CollectMultipleMetrics extends org.broadinstitute.sting.queue.function.JavaCommandLineFunction with PicardBamFunction{
analysisName = "CalculateMultipleMetrics"
javaMainClass = "net.sf.picard.sam.CalculateMultipleMetrics"
@Input(doc="The input SAM or BAM files to analyze. Must be coordinate sorted.", shortName = "input", fullName = "input_bam_files", required = true)
var input: Seq[File] = Nil
@Output(doc="The output file to write statistics to", shortName = "output", fullName = "output_file", required = true)
var output: File = _
@Argument(doc="Reference file", shortName = "reference", fullName = "reference", required = true)
var reference: File = _
override def inputBams = input
override def outputBam = output
override def commandLine = super.commandLine +
required("REFERENCE_SEQUENCE=" + reference) +
required("ASSUME_SORTED=true") +
required("PROGRAM=QualityScoreDistribution") +
required("PROGRAM=MeanQualityByCycle") +
required("PROGRAM=CollectAlignmentSummaryMetrics" )
}

80
public/scala/src/org/broadinstitute/sting/queue/function/QFunction.scala
View File

@ -28,7 +28,6 @@ import java.io.File
import java.lang.annotation.Annotation
import org.broadinstitute.sting.commandline._
import org.broadinstitute.sting.queue.{QException, QSettings}
import collection.JavaConversions._
import java.lang.IllegalStateException
import org.broadinstitute.sting.queue.util._
import org.broadinstitute.sting.utils.io.IOUtils
@ -194,13 +193,13 @@ trait QFunction extends Logging with QJobReport {
def failOutputs: Seq[File] = statusPrefixes.map(path => new File(path + ".fail"))
/** The complete list of fields on this CommandLineFunction. */
def functionFields = QFunction.classFields(this.functionFieldClass).functionFields
def functionFields: Seq[ArgumentSource] = ClassFieldCache.classFunctionFields(this.functionFieldClass)
/** The @Input fields on this CommandLineFunction. */
def inputFields = QFunction.classFields(this.functionFieldClass).inputFields
def inputFields: Seq[ArgumentSource] = ClassFieldCache.classInputFields(this.functionFieldClass)
/** The @Output fields on this CommandLineFunction. */
def outputFields = QFunction.classFields(this.functionFieldClass).outputFields
def outputFields: Seq[ArgumentSource] = ClassFieldCache.classOutputFields(this.functionFieldClass)
/** The @Argument fields on this CommandLineFunction. */
def argumentFields = QFunction.classFields(this.functionFieldClass).argumentFields
def argumentFields: Seq[ArgumentSource] = ClassFieldCache.classArgumentFields(this.functionFieldClass)
/**
* Returns the class that should be used for looking up fields.
@ -475,79 +474,12 @@ trait QFunction extends Logging with QJobReport {
* @param source Field to get the value for.
* @return value of the field.
*/
def getFieldValue(source: ArgumentSource) = ReflectionUtils.getValue(invokeObj(source), source.field)
def getFieldValue(source: ArgumentSource) = ClassFieldCache.getFieldValue(this, source)
/**
* Gets the value of a field.
* @param source Field to set the value for.
* @return value of the field.
*/
def setFieldValue(source: ArgumentSource, value: Any) = ReflectionUtils.setValue(invokeObj(source), source.field, value)
/**
* Walks gets the fields in this object or any collections in that object
* recursively to find the object holding the field to be retrieved or set.
* @param source Field find the invoke object for.
* @return Object to invoke the field on.
*/
private def invokeObj(source: ArgumentSource) = source.parentFields.foldLeft[AnyRef](this)(ReflectionUtils.getValue(_, _))
}
object QFunction {
var parsingEngine: ParsingEngine = _
/**
* The list of fields defined on a class
* @param clazz The class to lookup fields.
*/
private class ClassFields(clazz: Class[_]) {
/** The complete list of fields on this CommandLineFunction. */
val functionFields: Seq[ArgumentSource] = parsingEngine.extractArgumentSources(clazz).toSeq
/** The @Input fields on this CommandLineFunction. */
val inputFields = functionFields.filter(source => ReflectionUtils.hasAnnotation(source.field, classOf[Input]))
/** The @Output fields on this CommandLineFunction. */
val outputFields = functionFields.filter(source => ReflectionUtils.hasAnnotation(source.field, classOf[Output]))
/** The @Argument fields on this CommandLineFunction. */
val argumentFields = functionFields.filter(source => ReflectionUtils.hasAnnotation(source.field, classOf[Argument]))
}
/**
* The mapping from class to fields.
*/
private var classFieldsMap = Map.empty[Class[_], ClassFields]
/**
* Returns the field on clazz.
* @param clazz Class to search.
* @param name Name of the field to return.
* @return Argument source for the field.
*/
def findField(clazz: Class[_], name: String) = {
classFields(clazz).functionFields.find(_.field.getName == name) match {
case Some(source) => source
case None => throw new QException("Could not find a field on class %s with name %s".format(clazz, name))
}
}
/**
* Returns the fields for a class.
* @param clazz Class to retrieve fields for.
* @return the fields for the class.
*/
private def classFields(clazz: Class[_]) = {
classFieldsMap.get(clazz) match {
case Some(classFields) => classFields
case None =>
val classFields = new ClassFields(clazz)
classFieldsMap += clazz -> classFields
classFields
}
}
/**
* Returns the Seq of fields for a QFunction class.
* @param clazz Class to retrieve fields for.
* @return the fields of the class.
*/
def classFunctionFields(clazz: Class[_]) = classFields(clazz).functionFields
def setFieldValue(source: ArgumentSource, value: Any) = ClassFieldCache.setFieldValue(this, source, value)
}

9
public/scala/src/org/broadinstitute/sting/queue/function/scattergather/CloneFunction.scala
View File

@ -25,13 +25,14 @@
package org.broadinstitute.sting.queue.function.scattergather
import org.broadinstitute.sting.commandline.ArgumentSource
import org.broadinstitute.sting.queue.function.{QFunction, CommandLineFunction}
import org.broadinstitute.sting.queue.function.CommandLineFunction
import org.broadinstitute.sting.queue.util.ClassFieldCache
/**
* Shadow clones another command line function.
*/
object CloneFunction {
private lazy val cloneFunctionFields = QFunction.classFunctionFields(classOf[CloneFunction])
private lazy val cloneFunctionFields = ClassFieldCache.classFunctionFields(classOf[CloneFunction])
}
class CloneFunction extends CommandLineFunction {
@ -76,7 +77,7 @@ class CloneFunction extends CommandLineFunction {
def commandLine = withScatterPart(() => originalFunction.commandLine)
def getFieldValue(field: String): AnyRef = {
val source = QFunction.findField(originalFunction.getClass, field)
val source = ClassFieldCache.findField(originalFunction.getClass, field)
getFieldValue(source)
}
@ -98,7 +99,7 @@ class CloneFunction extends CommandLineFunction {
}
def setFieldValue(field: String, value: Any) {
val source = QFunction.findField(originalFunction.getClass, field)
val source = ClassFieldCache.findField(originalFunction.getClass, field)
setFieldValue(source, value)
}

183
public/scala/src/org/broadinstitute/sting/queue/util/ClassFieldCache.scala 100644
View File

@ -0,0 +1,183 @@
package org.broadinstitute.sting.queue.util
import org.broadinstitute.sting.commandline._
import scala.Some
import org.broadinstitute.sting.queue.QException
import collection.JavaConversions._
import java.io.File
/**
* Utilities and a static cache of argument fields for various classes populated by the parsingEngine.
* Because this class works with the ParsingEngine it can walk @ArgumentCollection hierarchies.
*/
object ClassFieldCache {
var parsingEngine: ParsingEngine = _
//
// Field caching
//
/**
* The list of fields defined on a class
* @param clazz The class to lookup fields.
*/
private class ClassFields(clazz: Class[_]) {
/** The complete list of fields on this CommandLineFunction. */
val functionFields: Seq[ArgumentSource] = parsingEngine.extractArgumentSources(clazz).toSeq
/** The @Input fields on this CommandLineFunction. */
val inputFields: Seq[ArgumentSource] = functionFields.filter(source => ReflectionUtils.hasAnnotation(source.field, classOf[Input]))
/** The @Output fields on this CommandLineFunction. */
val outputFields: Seq[ArgumentSource] = functionFields.filter(source => ReflectionUtils.hasAnnotation(source.field, classOf[Output]))
/** The @Argument fields on this CommandLineFunction. */
val argumentFields: Seq[ArgumentSource] = functionFields.filter(source => ReflectionUtils.hasAnnotation(source.field, classOf[Argument]))
}
/**
* The mapping from class to fields.
*/
private var classFieldsMap = Map.empty[Class[_], ClassFields]
/**
* Returns the fields for a class.
* @param clazz Class to retrieve fields for.
* @return the fields for the class.
*/
private def classFields(clazz: Class[_]): ClassFields = {
classFieldsMap.get(clazz) match {
case Some(classFields) => classFields
case None =>
val classFields = new ClassFields(clazz)
classFieldsMap += clazz -> classFields
classFields
}
}
/**
* Returns the field on clazz.
* @param clazz Class to search.
* @param name Name of the field to return.
* @return Argument source for the field.
*/
def findField(clazz: Class[_], name: String): ArgumentSource = {
classFields(clazz).functionFields.find(_.field.getName == name) match {
case Some(source) => source
case None => throw new QException("Could not find a field on class %s with name %s".format(clazz, name))
}
}
/**
* Returns the Seq of fields for a QFunction class.
* @param clazz Class to retrieve fields for.
* @return the fields of the class.
*/
def classFunctionFields(clazz: Class[_]): Seq[ArgumentSource] = classFields(clazz).functionFields
/**
* Returns the Seq of inputs for a QFunction class.
* @param clazz Class to retrieve inputs for.
* @return the inputs of the class.
*/
def classInputFields(clazz: Class[_]): Seq[ArgumentSource] = classFields(clazz).inputFields
/**
* Returns the Seq of outputs for a QFunction class.
* @param clazz Class to retrieve outputs for.
* @return the outputs of the class.
*/
def classOutputFields(clazz: Class[_]): Seq[ArgumentSource] = classFields(clazz).outputFields
/**
* Returns the Seq of arguments for a QFunction class.
* @param clazz Class to retrieve arguments for.
* @return the arguments of the class.
*/
def classArgumentFields(clazz: Class[_]): Seq[ArgumentSource] = classFields(clazz).argumentFields
//
// get/set fields as AnyRef
//
/**
* Gets the value of a field.
* @param obj Top level object storing the source info.
* @param source Field to get the value for.
* @return value of the field.
*/
def getFieldValue(obj: AnyRef, source: ArgumentSource) = ReflectionUtils.getValue(invokeObj(obj, source), source.field)
/**
* Gets the value of a field.
* @param obj Top level object storing the source info.
* @param source Field to set the value for.
* @return value of the field.
*/
def setFieldValue(obj: AnyRef, source: ArgumentSource, value: Any) = ReflectionUtils.setValue(invokeObj(obj, source), source.field, value)
/**
* Walks gets the fields in this object or any collections in that object
* recursively to find the object holding the field to be retrieved or set.
* @param obj Top level object storing the source info.
* @param source Field find the invoke object for.
* @return Object to invoke the field on.
*/
private def invokeObj(obj: AnyRef, source: ArgumentSource) = source.parentFields.foldLeft[AnyRef](obj)(ReflectionUtils.getValue(_, _))
//
// get/set fields as java.io.File
//
/**
* Gets the files from the fields. The fields must be a File, a FileExtension, or a Seq or Set of either.
* @param obj Top level object storing the source info.
* @param fields Fields to get files.
* @return for the fields.
*/
def getFieldFiles(obj: AnyRef, fields: Seq[ArgumentSource]): Seq[File] = {
var files: Seq[File] = Nil
for (field <- fields)
files ++= getFieldFiles(obj, field)
files.distinct
}
/**
* Gets the files from the field. The field must be a File, a FileExtension, or a Seq or Set of either.
* @param obj Top level object storing the source info.
* @param field Field to get files.
* @return for the field.
*/
def getFieldFiles(obj: AnyRef, field: ArgumentSource): Seq[File] = {
var files: Seq[File] = Nil
CollectionUtils.foreach(getFieldValue(obj, field), (fieldValue) => {
val file = fieldValueToFile(field, fieldValue)
if (file != null)
files :+= file
})
files.distinct
}
/**
* Gets the file from the field. The field must be a File or a FileExtension and not a Seq or Set.
* @param obj Top level object storing the source info.
* @param field Field to get the file.
* @return for the field.
*/
def getFieldFile(obj: AnyRef, field: ArgumentSource): File =
fieldValueToFile(field, getFieldValue(obj, field))
/**
* Converts the field value to a file. The field must be a File or a FileExtension.
* @param field Field to get the file.
* @param value Value of the File or FileExtension or null.
* @return Null if value is null, otherwise the File.
* @throws QException if the value is not a File or FileExtension.
*/
private def fieldValueToFile(field: ArgumentSource, value: Any): File = value match {
case file: File => file
case null => null
case unknown => throw new QException("Non-file found. Try removing the annotation, change the annotation to @Argument, or extend File with FileExtension: %s: %s".format(field.field, unknown))
}
}

3
public/scala/src/org/broadinstitute/sting/queue/util/QScriptUtils.scala
View File

@ -57,7 +57,8 @@ object QScriptUtils {
for (file <- fromFile(in).getLines())
if (!file.startsWith("#") && !file.isEmpty )
list :+= new File(file.trim())
list.sortWith(_.compareTo(_) < 0)
// list.sortWith(_.compareTo(_) < 0)
list
}
/**

14
public/scala/src/org/broadinstitute/sting/queue/util/RemoteFile.scala 100644
View File

@ -0,0 +1,14 @@
package org.broadinstitute.sting.queue.util
import java.io.File
import org.broadinstitute.sting.utils.io.FileExtension
/**
* An extension of java.io.File that can be pulled from or pushed to a remote location.
*/
trait RemoteFile extends File with FileExtension {
def pullToLocal()
def pushToRemote()
def deleteRemote()
def remoteDescription: String
}

1
settings/ivysettings.xml
View File

@ -7,7 +7,6 @@
<artifact pattern="${repository.dir}/[organisation]/[artifact]-[revision].[ext]" />
<ivy pattern="${repository.dir}/[organisation]/[module]-[revision].xml" />
</filesystem>
<ibiblio name="reflections-repo" m2compatible="true" root="http://reflections.googlecode.com/svn/repo" />
<ibiblio name="maven" root="http://repo1.maven.org/maven2" m2compatible="true" />
</chain>
</resolvers>