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

This commit is contained in:
Guillermo del Angel 2011-07-12 12:37:43 -04:00
commit a0d074a71c
62 changed files with 2130 additions and 304 deletions

101
build.xml
View File

@ -28,6 +28,7 @@
<property name="build.dir" value="build" />
<property name="dist.dir" value="dist" />
<property name="lib.dir" value="lib" />
<property name="external.dir" value="external" />
<property name="public.dir" value="public" />
<property name="private.dir" value="private" />
@ -44,11 +45,11 @@
<property name="queue-extensions.source.dir" value="${build.dir}/queue-extensions/src" />
<!-- Contracts for Java -->
<!-- uncomment out to enable building contracts -->
<property name="use.contracts" value="true" />
<!-- To disable, run with -Duse.contracts=false -->
<property name="use.contracts" value="true" />
<property name="java.contracts" value="${build.dir}/java/contracts" />
<property name="cofojaDir" value="settings/repository/com.google/"/>
<property name="cofoja.jar" value="${cofojaDir}/cofoja-1.0-20110609.jar"/>
<property name="contracts.version" value="1.0-20110609" />
<property name="cofoja.jar" value="${lib.dir}/cofoja-${contracts.version}.jar"/>
<!-- where to find the tribble distro -->
<property name="tribble.dir" value="tribble" />
@ -69,8 +70,6 @@
<property environment="env"/>
<property name="drmaa.jar" value="${env.SGE_ROOT}/lib/drmaa.jar" />
<patternset id="java.source.pattern">
<include name="${java.public.source.dir}/**/*.java" />
<include name="${java.private.source.dir}/**/*.java" if="include.private" />
@ -103,7 +102,7 @@
</patternset>
<path id="external.dependencies">
<fileset dir="lib">
<fileset dir="${lib.dir}">
<patternset refid="dependency.mask" />
</fileset>
</path>
@ -128,14 +127,14 @@
<!-- ivy properties -->
<property name="ivy.install.version" value="2.2.0"/>
<property name="ivy.home" value="${user.home}/.ant"/>
<property name="ivy.jar.dir" value="${ivy.home}/lib"/>
<property name="ivy.jar.dir" value="${ivy.home}/${lib.dir}"/>
<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"/>
<mkdir dir="lib"/>
<mkdir dir="${lib.dir}"/>
<mkdir dir="${ivy.jar.dir}"/>
<copy file="${cofoja.jar}" toFile="lib/cofoja.jar"/>
<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"/>
@ -146,11 +145,7 @@
<property name="init.resolve.done" value="true"/>
</target>
<target name="init.gridengine" depends="init" if="include.gridengine">
<copy todir="lib" file="${drmaa.jar}"/>
</target>
<target name="resolve" depends="init.resolve,init,init.gridengine"
<target name="resolve" depends="init.resolve,init"
description="locate and download library dependencies">
<property name="ivy.conf" value="default"/>
<ivy:retrieve file="ivy.xml" conf="${ivy.conf}" />
@ -178,13 +173,23 @@
<property name="build.version" value="${git.describe.output}" />
</target>
<target name="untagged.build.version" depends="git.describe" unless="git.describe.succeeded">
<exec executable="git" outputproperty="build.version" failonerror="true">
<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" />
</exec>
<condition property="git.rev-parse.succeeded">
<equals arg1="${git.rev-parse.exit.value}" arg2="0" />
</condition>
</target>
<target name="generate.build.version" depends="tagged.build.version, untagged.build.version" />
<target name="untagged.build.version" depends="git.rev-parse" if="git.rev-parse.succeeded">
<property name="build.version" value="${git.rev-parse.output}" />
</target>
<target name="generate.build.version" depends="tagged.build.version, untagged.build.version">
<!-- Set build.version to exported if no other value has been set -->
<property name="build.version" value="exported" />
</target>
<!-- define some key locations that might change based on how the build is run -->
<target name="init" depends="generate.build.version">
@ -214,12 +219,6 @@
</or>
</condition>
<!-- Include Grid Engine in the compile if SGE_ROOT is available. -->
<!-- Based off of http://wikis.sun.com/display/GridEngine/Automating+Grid+Engine+Functions+Through+DRMAA -->
<condition property="include.gridengine">
<available file="${drmaa.jar}"/>
</condition>
<echo message="GATK build : ${gatk.target}"/>
<echo message="Scala build : ${scala.target}"/>
<echo message="source revision : ${build.version}"/>
@ -229,6 +228,10 @@
<equals arg1="${gatk.target}" arg2="private" casesensitive="false" />
</condition>
<condition property="include.contracts">
<equals arg1="${use.contracts}" arg2="true" />
</condition>
<!-- Create the build directory structure used by compile -->
<mkdir dir="${build.dir}"/>
<mkdir dir="${java.classes}"/>
@ -252,7 +255,7 @@
<target name="init.scala.compile" depends="resolve"
description="Initializes the scala ant tasks from scala-compiler.jar">
<path id="scala.classpath">
<fileset dir="lib">
<fileset dir="${lib.dir}">
<include name="scala-compiler-*.jar"/>
<include name="scala-library-*.jar"/>
</fileset>
@ -287,7 +290,7 @@
depends="gatk.compile.public.source,gatk.compile.private.source,gatk.compile.external.source"
description="compile the GATK source" />
<target name="gatk.contracts.public" depends="gatk.compile.source">
<target name="gatk.contracts.public" depends="gatk.compile.source" if="include.contracts">
<javac fork="true" memoryMaximumSize="512m" srcdir="${java.public.source.dir}" destdir="${java.contracts}" debug="true" debuglevel="lines,vars,source" tempdir="${java.io.tmpdir}" >
<classpath>
<path refid="external.dependencies" />
@ -299,7 +302,16 @@
</javac>
</target>
<target name="gatk.contracts.private" depends="gatk.compile.source" if="include.private">
<target name="check.contracts.private" depends="gatk.contracts.public">
<condition property="include.contracts.private">
<and>
<isset property="include.contracts" />
<isset property="include.private" />
</and>
</condition>
</target>
<target name="gatk.contracts.private" depends="check.contracts.private" if="include.contracts.private">
<javac fork="true" memoryMaximumSize="512m" srcdir="${java.private.source.dir}" destdir="${java.contracts}" debug="true" debuglevel="lines,vars,source" tempdir="${java.io.tmpdir}" >
<classpath>
<path refid="external.dependencies" />
@ -312,7 +324,7 @@
</target>
<target name="gatk.contracts" depends="gatk.contracts.public,gatk.contracts.private"
description="create GATK contracts" if="use.contracts" />
description="create GATK contracts" if="include.contracts" />
<target name="gatk.compile" depends="tribble,init,resolve,gatk.compile.source,gatk.contracts" />
@ -357,7 +369,6 @@
<src path="${scala.public.source.dir}" />
<src path="${queue-extensions.source.dir}" />
<include name="**/*.scala"/>
<exclude name="**/gridengine/**" unless="include.gridengine" />
</scalac>
</target>
@ -374,7 +385,6 @@
<scalac fork="true" jvmargs="-Xmx512m" destdir="${scala.classes}" classpathref="scala.dependencies" deprecation="yes" unchecked="yes">
<src path="${scala.private.source.dir}" />
<include name="**/*.scala"/>
<exclude name="**/gridengine/**" unless="include.gridengine" />
</scalac>
</target>
@ -452,7 +462,7 @@
<target name="init.jar" depends="sting.compile,extracthelp">
<mkdir dir="${dist.dir}"/>
<copy todir="${dist.dir}">
<fileset dir="lib" includes="*.jar"/>
<fileset dir="${lib.dir}" includes="*.jar"/>
</copy>
</target>
@ -464,7 +474,7 @@
<exclude name="**/utils/variantcontext/**/*.class"/>
</fileset>
<fileset dir="${java.classes}" includes="**/commandline/**/*.class"/>
<fileset dir="${java.classes}" includes="**/sting/datasources/**/*.class"/>
<fileset dir="${java.classes}" includes="**/sting/pipeline/**/*.class"/>
<fileset dir="${java.classes}" includes="**/sting/jna/**/*.class"/>
<fileset dir="${java.classes}" includes="net/sf/picard/**/*.class"/>
<fileset dir="${java.classes}" includes="net/sf/samtools/**/*.class"/>
@ -663,7 +673,7 @@
<path refid="external.dependencies" />
<pathelement location="${java.classes}"/>
<pathelement location="${java.contracts}"/>
<pathelement location="lib/testng-5.14.1.jar"/>
<pathelement location="${lib.dir}/testng-5.14.1.jar"/>
</classpath>
<compilerarg value="-proc:none"/>
<!--
@ -682,12 +692,11 @@
<src path="${scala.public.test.sources}" />
<src path="${scala.private.test.sources}" />
<include name="**/*.scala"/>
<exclude name="**/gridengine/**" unless="include.gridengine" />
<classpath>
<path refid="scala.dependencies"/>
<pathelement location="${scala.test.classes}"/>
<pathelement location="${java.test.classes}"/>
<pathelement location="lib/testng-5.14.1.jar"/>
<pathelement location="${lib.dir}/testng-5.14.1.jar"/>
</classpath>
</scalac>
</target>
@ -727,9 +736,13 @@
<equals arg1="${pipeline.run}" arg2="$${pipeline.run}" />
</condition>
<condition property="cofoja.jvm.args" value="-javaagent:${cofoja.jar} -Dcom.google.java.contract.log.contract=false" else="">
<isset property="include.contracts" />
</condition>
<mkdir dir="${report}/@{testtype}"/>
<echo message="Sting: Running @{testtype} test cases!"/>
<taskdef resource="testngtasks" classpath="lib/testng-5.14.1.jar"/>
<taskdef resource="testngtasks" classpath="${lib.dir}/testng-5.14.1.jar"/>
<testng outputDir="${report}/@{testtype}"
haltOnFailure="false" failureProperty="test.failure"
verbose="2"
@ -740,9 +753,7 @@
<jvmarg value="-Djava.awt.headless=true" />
<jvmarg value="-Dpipeline.run=${pipeline.run}" />
<jvmarg value="-Djava.io.tmpdir=${java.io.tmpdir}" />
<!-- needs to be if -->
<!--jvmarg value="-javaagent:lib/cofoja.jar"/-->
<!--jvmarg value="-Dcom.google.java.contract.log.contract=false"/-->
<jvmarg line="${cofoja.jvm.args}"/>
<!-- <jvmarg value="-Xdebug"/> -->
<!-- <jvmarg value="-Xrunjdwp:transport=dt_socket,server=y,suspend=y,address=5005"/> -->
<classpath>
@ -820,7 +831,7 @@
<!-- copy the compiled library -->
<target name="tribble.compile.copy" description="Copies the compiled tribble library" depends="tribble.compile" if="tribble.compile.exists">
<copy todir="lib">
<copy todir="${lib.dir}">
<fileset dir="${tribble.dir}/dist" includes="*.jar"/>
</copy>
</target>
@ -828,7 +839,7 @@
<!-- copy the precompiled library -->
<target name="tribble.library.copy" description="Copies the precompiled tribble library" depends="tribble.init" unless="tribble.compile.exists">
<echo message="Copying the Tribble Library..."/>
<copy todir="lib">
<copy todir="${lib.dir}">
<fileset dir="settings/repository/org.broad" includes="tribble*.jar"/>
</copy>
</target>
@ -921,8 +932,8 @@
<pathconvert property="required.picard.jars" pathsep=":">
<fileset dir="${basedir}">
<include name="staging" />
<include name="lib/picard-*.*.*.jar" />
<include name="lib/sam-*.jar" />
<include name="${lib.dir}/picard-*.*.*.jar" />
<include name="${lib.dir}/sam-*.jar" />
</fileset>
</pathconvert>
<echo message="required.picard.jars=${required.picard.jars}" />
@ -944,7 +955,7 @@
<!-- Find bug tasks -->
<!-- ******************************************************************************** -->
<path id="findbugs.classpath">
<fileset dir="lib" erroronmissingdir="true" includes="*.jar"/>
<fileset dir="${lib.dir}" erroronmissingdir="true" includes="*.jar"/>
</path>
<target name="findbugs" depends="dist">
<antcall target ="resolve">
@ -969,7 +980,7 @@
<target name="clean" description="clean up" depends="tribble.clean,clean.javadoc">
<delete dir="out"/>
<delete dir="${build.dir}"/>
<delete dir="lib"/>
<delete dir="${lib.dir}"/>
<delete dir="staging"/>
<delete dir="${dist.dir}"/>
<delete dir="pipelinetests"/>

View File

@ -48,6 +48,9 @@
<!-- Dependencies for amazon.com S3 support -->
<dependency org="net.java.dev.jets3t" name="jets3t" rev="0.8.0"/>
<!-- Dependencies for GridEngine -->
<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"/>
@ -60,6 +63,10 @@
<dependency org="net.sourceforge.findbugs" name="jsr305" rev="1.3.2" conf="test"/>
<dependency org="com.google.code.caliper" name="caliper" rev="1.0-SNAPSHOT" conf="test" />
<!-- Contracts for Java and dependencies -->
<dependency org="com.google.code.cofoja" name="cofoja" rev="1.0-20110609" />
<dependency org="asm" name="asm-all" rev="3.3.1" />
<!-- POI, for reading pipeline files -->
<dependency org="org.apache.poi" name="poi" rev="3.8-beta3" />
<dependency org="org.apache.poi" name="poi-ooxml" rev="3.8-beta3" />

View File

@ -61,7 +61,7 @@ public class AnalyzeCovariates extends CommandLineProgram {
@Argument(fullName = "path_to_Rscript", shortName = "Rscript", doc = "The path to your implementation of Rscript. For Broad users this is maybe /broad/tools/apps/R-2.6.0/bin/Rscript", required = false)
private String PATH_TO_RSCRIPT = "Rscript";
@Argument(fullName = "path_to_resources", shortName = "resources", doc = "Path to resources folder holding the Sting R scripts.", required = false)
private String PATH_TO_RESOURCES = "R/";
private String PATH_TO_RESOURCES = "public/R/";
@Argument(fullName = "ignoreQ", shortName = "ignoreQ", doc = "Ignore bases with reported quality less than this number.", required = false)
private int IGNORE_QSCORES_LESS_THAN = 5;
@Argument(fullName = "numRG", shortName = "numRG", doc = "Only process N read groups. Default value: -1 (process all read groups)", required = false)

View File

@ -55,10 +55,14 @@ public class TableFeature implements Feature {
}
public List<String> getAllValues() {
return getValuesTo(values.size()-1);
return getValuesTo(values.size());
}
public List<String> getValuesTo(int columnPosition) {
return values.subList(0,columnPosition);
}
public List<String> getHeader() {
return keys;
}
}

View File

@ -25,6 +25,7 @@
package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeaderLineCount;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeaderLineType;
import org.broadinstitute.sting.utils.codecs.vcf.VCFInfoHeaderLine;
@ -41,8 +42,8 @@ import java.util.*;
public class ChromosomeCounts implements InfoFieldAnnotation, StandardAnnotation {
private String[] keyNames = { VCFConstants.ALLELE_NUMBER_KEY, VCFConstants.ALLELE_COUNT_KEY, VCFConstants.ALLELE_FREQUENCY_KEY };
private VCFInfoHeaderLine[] descriptions = { new VCFInfoHeaderLine(VCFConstants.ALLELE_FREQUENCY_KEY, -1, VCFHeaderLineType.Float, "Allele Frequency, for each ALT allele, in the same order as listed"),
new VCFInfoHeaderLine(VCFConstants.ALLELE_COUNT_KEY, -1, VCFHeaderLineType.Integer, "Allele count in genotypes, for each ALT allele, in the same order as listed"),
private VCFInfoHeaderLine[] descriptions = { new VCFInfoHeaderLine(VCFConstants.ALLELE_FREQUENCY_KEY, VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.Float, "Allele Frequency, for each ALT allele, in the same order as listed"),
new VCFInfoHeaderLine(VCFConstants.ALLELE_COUNT_KEY, VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.Integer, "Allele count in genotypes, for each ALT allele, in the same order as listed"),
new VCFInfoHeaderLine(VCFConstants.ALLELE_NUMBER_KEY, 1, VCFHeaderLineType.Integer, "Total number of alleles in called genotypes") };
public Map<String, Object> annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, AlignmentContext> stratifiedContexts, VariantContext vc) {

View File

@ -1,5 +1,6 @@
package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeaderLineCount;
import org.broadinstitute.sting.utils.variantcontext.Allele;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
@ -142,5 +143,5 @@ public class DepthPerAlleleBySample implements GenotypeAnnotation, StandardAnnot
// public String getIndelBases()
public List<String> getKeyNames() { return Arrays.asList("AD"); }
public List<VCFFormatHeaderLine> getDescriptions() { return Arrays.asList(new VCFFormatHeaderLine(getKeyNames().get(0), VCFCompoundHeaderLine.UNBOUNDED, VCFHeaderLineType.Integer, "Allelic depths for the ref and alt alleles in the order listed")); }
public List<VCFFormatHeaderLine> getDescriptions() { return Arrays.asList(new VCFFormatHeaderLine(getKeyNames().get(0), VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.Integer, "Allelic depths for the ref and alt alleles in the order listed")); }
}

View File

@ -29,6 +29,7 @@ import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.GenotypeAnnotation;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeaderLineCount;
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
import org.broadinstitute.sting.utils.pileup.PileupElement;
import org.broadinstitute.sting.utils.pileup.ReadBackedExtendedEventPileup;
@ -200,8 +201,8 @@ public class ReadDepthAndAllelicFractionBySample implements GenotypeAnnotation {
1,
VCFHeaderLineType.Integer,
"Total read depth per sample, including MQ0"),
new VCFFormatHeaderLine(getKeyNames().get(1),
VCFCompoundHeaderLine.UNBOUNDED,
new VCFFormatHeaderLine(getKeyNames().get(1),
VCFHeaderLineCount.UNBOUNDED,
VCFHeaderLineType.Float,
"Fractions of reads (excluding MQ0 from both ref and alt) supporting each reported alternative allele, per sample"));
}

View File

@ -25,6 +25,7 @@
package org.broadinstitute.sting.gatk.walkers.annotator;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeaderLineCount;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeaderLineType;
@ -65,5 +66,5 @@ public class SampleList implements InfoFieldAnnotation {
public List<String> getKeyNames() { return Arrays.asList("Samples"); }
public List<VCFInfoHeaderLine> getDescriptions() { return Arrays.asList(new VCFInfoHeaderLine("Samples", VCFInfoHeaderLine.UNBOUNDED, VCFHeaderLineType.String, "List of polymorphic samples")); }
public List<VCFInfoHeaderLine> getDescriptions() { return Arrays.asList(new VCFInfoHeaderLine("Samples", VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.String, "List of polymorphic samples")); }
}

View File

@ -0,0 +1,122 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
import net.sf.samtools.*;
import net.sf.samtools.util.BlockCompressedInputStream;
import org.broad.tribble.readers.AsciiLineReader;
import org.broad.tribble.readers.LineReader;
import org.broadinstitute.sting.utils.codecs.vcf.VCFCodec;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeader;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import java.io.DataInputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.util.Arrays;
import java.util.Map;
import java.util.zip.GZIPInputStream;
/**
* Created by IntelliJ IDEA.
* User: depristo
* Date: 7/4/11
* Time: 1:09 PM
*
* Class implementing diffnode reader for VCF
*/
public class BAMDiffableReader implements DiffableReader {
private final static int MAX_RECORDS_TO_READ = 1000;
@Override
public String getName() { return "BAM"; }
@Override
public DiffElement readFromFile(File file) {
final SAMFileReader reader = new SAMFileReader(file, null); // null because we don't want it to look for the index
reader.setValidationStringency(SAMFileReader.ValidationStringency.SILENT);
DiffNode root = DiffNode.rooted(file.getName());
SAMRecordIterator iterator = reader.iterator();
int count = 0;
while ( iterator.hasNext() ) {
if ( count++ > MAX_RECORDS_TO_READ )
break;
final SAMRecord record = iterator.next();
// name is the read name + first of pair
String name = record.getReadName().replace('.', '_');
if ( record.getReadPairedFlag() ) {
name += record.getFirstOfPairFlag() ? "_1" : "_2";
}
DiffNode readRoot = DiffNode.empty(name, root);
// add fields
readRoot.add("NAME", record.getReadName());
readRoot.add("FLAGS", record.getFlags());
readRoot.add("RNAME", record.getReferenceName());
readRoot.add("POS", record.getAlignmentStart());
readRoot.add("MAPQ", record.getMappingQuality());
readRoot.add("CIGAR", record.getCigarString());
readRoot.add("RNEXT", record.getMateReferenceName());
readRoot.add("PNEXT", record.getMateAlignmentStart());
readRoot.add("TLEN", record.getInferredInsertSize());
readRoot.add("SEQ", record.getReadString());
readRoot.add("QUAL", record.getBaseQualityString());
for ( SAMRecord.SAMTagAndValue xt : record.getAttributes() ) {
readRoot.add(xt.tag, xt.value);
}
// add record to root
if ( ! root.hasElement(name) )
// protect ourselves from malformed files
root.add(readRoot);
}
reader.close();
return root.getBinding();
}
@Override
public boolean canRead(File file) {
final byte[] BAM_MAGIC = "BAM\1".getBytes();
final byte[] buffer = new byte[BAM_MAGIC.length];
try {
FileInputStream fstream = new FileInputStream(file);
new BlockCompressedInputStream(fstream).read(buffer,0,BAM_MAGIC.length);
return Arrays.equals(buffer, BAM_MAGIC);
} catch ( IOException e ) {
return false;
} catch ( net.sf.samtools.FileTruncatedException e ) {
return false;
}
}
}

View File

@ -0,0 +1,118 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
import com.google.java.contract.*;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
/**
* Created by IntelliJ IDEA.
* User: depristo
* Date: 7/4/11
* Time: 12:55 PM
*
* An interface that must be implemented to allow us to calculate differences
* between structured objects
*/
@Invariant({
"name != null",
"value != null",
"parent != null || name.equals(\"ROOT\")",
"value == null || value.getBinding() == this"})
public class DiffElement {
public final static DiffElement ROOT = new DiffElement();
final private String name;
final private DiffElement parent;
final private DiffValue value;
/**
* For ROOT only
*/
private DiffElement() {
this.name = "ROOT";
this.parent = null;
this.value = new DiffValue(this, "ROOT");
}
@Requires({"name != null", "parent != null", "value != null"})
public DiffElement(String name, DiffElement parent, DiffValue value) {
if ( name.equals("ROOT") ) throw new IllegalArgumentException("Cannot use reserved name ROOT");
this.name = name;
this.parent = parent;
this.value = value;
this.value.setBinding(this);
}
@Ensures({"result != null"})
public String getName() {
return name;
}
public DiffElement getParent() {
return parent;
}
@Ensures({"result != null"})
public DiffValue getValue() {
return value;
}
public boolean isRoot() { return this == ROOT; }
@Ensures({"result != null"})
@Override
public String toString() {
return getName() + "=" + getValue().toString();
}
public String toString(int offset) {
return (offset > 0 ? Utils.dupString(' ', offset) : 0) + getName() + "=" + getValue().toString(offset);
}
@Ensures({"result != null"})
public final String fullyQualifiedName() {
if ( isRoot() )
return "";
else if ( parent.isRoot() )
return name;
else
return parent.fullyQualifiedName() + "." + name;
}
@Ensures({"result != null"})
public String toOneLineString() {
return getName() + "=" + getValue().toOneLineString();
}
@Ensures({"result != null"})
public DiffNode getValueAsNode() {
if ( getValue().isCompound() )
return (DiffNode)getValue();
else
throw new ReviewedStingException("Illegal request conversion of a DiffValue into a DiffNode: " + this);
}
}

View File

@ -0,0 +1,423 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
import com.google.java.contract.Requires;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.gatk.report.GATKReport;
import org.broadinstitute.sting.gatk.report.GATKReportTable;
import org.broadinstitute.sting.gatk.walkers.varianteval.stratifications.VariantStratifier;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.classloader.PluginManager;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.exceptions.UserException;
import java.io.File;
import java.io.PrintStream;
import java.util.*;
/**
* Created by IntelliJ IDEA.
* User: depristo
* Date: 7/4/11
* Time: 12:51 PM
* A generic engine for comparing tree-structured objects
*/
public class DiffEngine {
final protected static Logger logger = Logger.getLogger(DiffEngine.class);
private final Map<String, DiffableReader> readers = new HashMap<String, DiffableReader>();
public DiffEngine() {
loadDiffableReaders();
}
// --------------------------------------------------------------------------------
//
// difference calculation
//
// --------------------------------------------------------------------------------
public List<Difference> diff(DiffElement master, DiffElement test) {
DiffValue masterValue = master.getValue();
DiffValue testValue = test.getValue();
if ( masterValue.isCompound() && masterValue.isCompound() ) {
return diff(master.getValueAsNode(), test.getValueAsNode());
} else if ( masterValue.isAtomic() && testValue.isAtomic() ) {
return diff(masterValue, testValue);
} else {
// structural difference in types. one is node, other is leaf
return Arrays.asList(new Difference(master, test));
}
}
public List<Difference> diff(DiffNode master, DiffNode test) {
Set<String> allNames = new HashSet<String>(master.getElementNames());
allNames.addAll(test.getElementNames());
List<Difference> diffs = new ArrayList<Difference>();
for ( String name : allNames ) {
DiffElement masterElt = master.getElement(name);
DiffElement testElt = test.getElement(name);
if ( masterElt == null && testElt == null ) {
throw new ReviewedStingException("BUG: unexceptedly got two null elements for field: " + name);
} else if ( masterElt == null || testElt == null ) { // if either is null, we are missing a value
// todo -- should one of these be a special MISSING item?
diffs.add(new Difference(masterElt, testElt));
} else {
diffs.addAll(diff(masterElt, testElt));
}
}
return diffs;
}
public List<Difference> diff(DiffValue master, DiffValue test) {
if ( master.getValue().equals(test.getValue()) ) {
return Collections.emptyList();
} else {
return Arrays.asList(new Difference(master.getBinding(), test.getBinding()));
}
}
// --------------------------------------------------------------------------------
//
// Summarizing differences
//
// --------------------------------------------------------------------------------
/**
* Emits a summary of the diffs to out. Suppose you have the following three differences:
*
* A.X.Z:1!=2
* A.Y.Z:3!=4
* B.X.Z:5!=6
*
* The above is the itemized list of the differences. The summary looks for common differences
* in the name hierarchy, counts those shared elements, and emits the differences that occur
* in order of decreasing counts.
*
* So, in the above example, what are the shared elements?
*
* A.X.Z and B.X.Z share X.Z, so there's a *.X.Z with count 2
* A.X.Z, A.Y.Z, and B.X.Z all share *.*.Z, with count 3
* Each of A.X.Z, A.Y.Z, and B.X.Z are individually unique, with count 1
*
* So we would emit the following summary:
*
* *.*.Z: 3
* *.X.Z: 2
* A.X.Z: 1 [specific difference: 1!=2]
* A.Y.Z: 1 [specific difference: 3!=4]
* B.X.Z: 1 [specific difference: 5!=6]
*
* The algorithm to accomplish this calculation is relatively simple. Start with all of the
* concrete differences. For each pair of differences A1.A2....AN and B1.B2....BN:
*
* find the longest common subsequence Si.Si+1...SN where Ai = Bi = Si
* If i == 0, then there's no shared substructure
* If i > 0, then generate the summarized value X = *.*...Si.Si+1...SN
* if X is a known summary, increment it's count, otherwise set its count to 1
*
* Not that only pairs of the same length are considered as potentially equivalent
*
* @param params determines how we display the items
* @param diffs
*/
public void reportSummarizedDifferences(List<Difference> diffs, SummaryReportParams params ) {
printSummaryReport(summarizeDifferences(diffs), params );
}
public List<SummarizedDifference> summarizeDifferences(List<Difference> diffs) {
List<String[]> diffPaths = new ArrayList<String[]>(diffs.size());
for ( Difference diff1 : diffs ) {
diffPaths.add(diffNameToPath(diff1.getFullyQualifiedName()));
}
return summarizedDifferencesOfPaths(diffPaths);
}
final protected static String[] diffNameToPath(String diffName) {
return diffName.split("\\.");
}
protected List<SummarizedDifference> summarizedDifferencesOfPaths(List<String[]> diffPaths) {
Map<String, SummarizedDifference> summaries = new HashMap<String, SummarizedDifference>();
// create the initial set of differences
for ( int i = 0; i < diffPaths.size(); i++ ) {
for ( int j = 0; j <= i; j++ ) {
String[] diffPath1 = diffPaths.get(i);
String[] diffPath2 = diffPaths.get(j);
if ( diffPath1.length == diffPath2.length ) {
int lcp = longestCommonPostfix(diffPath1, diffPath2);
String path = lcp > 0 ? summarizedPath(diffPath2, lcp) : Utils.join(".", diffPath2);
addSummary(summaries, path, true);
}
}
}
// count differences
for ( String[] diffPath : diffPaths ) {
for ( SummarizedDifference sumDiff : summaries.values() ) {
if ( sumDiff.matches(diffPath) )
addSummary(summaries, sumDiff.getPath(), false);
}
}
List<SummarizedDifference> sortedSummaries = new ArrayList<SummarizedDifference>(summaries.values());
Collections.sort(sortedSummaries);
return sortedSummaries;
}
private static void addSummary(Map<String, SummarizedDifference> summaries, String path, boolean onlyCatalog) {
if ( summaries.containsKey(path) ) {
if ( ! onlyCatalog )
summaries.get(path).incCount();
} else {
SummarizedDifference sumDiff = new SummarizedDifference(path);
summaries.put(sumDiff.getPath(), sumDiff);
}
}
protected void printSummaryReport(List<SummarizedDifference> sortedSummaries, SummaryReportParams params ) {
GATKReport report = new GATKReport();
final String tableName = "diffences";
report.addTable(tableName, "Summarized differences between the master and test files.\nSee http://www.broadinstitute.org/gsa/wiki/index.php/DiffObjectsWalker_and_SummarizedDifferences for more information");
GATKReportTable table = report.getTable(tableName);
table.addPrimaryKey("Difference", true);
table.addColumn("NumberOfOccurrences", 0);
int count = 0, count1 = 0;
for ( SummarizedDifference diff : sortedSummaries ) {
if ( diff.getCount() < params.minSumDiffToShow )
// in order, so break as soon as the count is too low
break;
if ( params.maxItemsToDisplay != 0 && count++ > params.maxItemsToDisplay )
break;
if ( diff.getCount() == 1 ) {
count1++;
if ( params.maxCountOneItems != 0 && count1 > params.maxCountOneItems )
break;
}
table.set(diff.getPath(), "NumberOfOccurrences", diff.getCount());
}
table.write(params.out);
}
protected static int longestCommonPostfix(String[] diffPath1, String[] diffPath2) {
int i = 0;
for ( ; i < diffPath1.length; i++ ) {
int j = diffPath1.length - i - 1;
if ( ! diffPath1[j].equals(diffPath2[j]) )
break;
}
return i;
}
/**
* parts is [A B C D]
* commonPostfixLength: how many parts are shared at the end, suppose its 2
* We want to create a string *.*.C.D
*
* @param parts
* @param commonPostfixLength
* @return
*/
protected static String summarizedPath(String[] parts, int commonPostfixLength) {
int stop = parts.length - commonPostfixLength;
if ( stop > 0 ) parts = parts.clone();
for ( int i = 0; i < stop; i++ ) {
parts[i] = "*";
}
return Utils.join(".", parts);
}
/**
* TODO -- all of the algorithms above should use SummarizedDifference instead
* TODO -- of some SummarizedDifferences and some low-level String[]
*/
public static class SummarizedDifference implements Comparable<SummarizedDifference> {
final String path; // X.Y.Z
final String[] parts;
int count = 0;
public SummarizedDifference(String path) {
this.path = path;
this.parts = diffNameToPath(path);
}
public void incCount() { count++; }
public int getCount() {
return count;
}
/**
* The fully qualified path object A.B.C etc
* @return
*/
public String getPath() {
return path;
}
/**
* @return the length of the parts of this summary
*/
public int length() {
return this.parts.length;
}
/**
* Returns true if the string parts matches this summary. Matches are
* must be equal() everywhere where this summary isn't *.
* @param otherParts
* @return
*/
public boolean matches(String[] otherParts) {
if ( otherParts.length != length() )
return false;
// TODO optimization: can start at right most non-star element
for ( int i = 0; i < length(); i++ ) {
String part = parts[i];
if ( ! part.equals("*") && ! part.equals(otherParts[i]) )
return false;
}
return true;
}
@Override
public String toString() {
return String.format("%s:%d", getPath(), getCount());
}
@Override
public int compareTo(SummarizedDifference other) {
// sort first highest to lowest count, then by lowest to highest path
int countCmp = Integer.valueOf(count).compareTo(other.count);
return countCmp != 0 ? -1 * countCmp : path.compareTo(other.path);
}
}
// --------------------------------------------------------------------------------
//
// plugin manager
//
// --------------------------------------------------------------------------------
public void loadDiffableReaders() {
List<Class<? extends DiffableReader>> drClasses = new PluginManager<DiffableReader>( DiffableReader.class ).getPlugins();
logger.info("Loading diffable modules:");
for (Class<? extends DiffableReader> drClass : drClasses ) {
logger.info("\t" + drClass.getSimpleName());
try {
DiffableReader dr = drClass.newInstance();
readers.put(dr.getName(), dr);
} catch (InstantiationException e) {
throw new ReviewedStingException("Unable to instantiate module '" + drClass.getSimpleName() + "'");
} catch (IllegalAccessException e) {
throw new ReviewedStingException("Illegal access error when trying to instantiate '" + drClass.getSimpleName() + "'");
}
}
}
protected Map<String, DiffableReader> getReaders() {
return readers;
}
protected DiffableReader getReader(String name) {
return readers.get(name);
}
/**
* Returns a reader appropriate for this file, or null if no such reader exists
* @param file
* @return
*/
public DiffableReader findReaderForFile(File file) {
for ( DiffableReader reader : readers.values() )
if (reader.canRead(file) )
return reader;
return null;
}
/**
* Returns true if reader appropriate for this file, or false if no such reader exists
* @param file
* @return
*/
public boolean canRead(File file) {
return findReaderForFile(file) != null;
}
public DiffElement createDiffableFromFile(File file) {
DiffableReader reader = findReaderForFile(file);
if ( reader == null )
throw new UserException("Unsupported file type: " + file);
else
return reader.readFromFile(file);
}
public static boolean simpleDiffFiles(File masterFile, File testFile, DiffEngine.SummaryReportParams params) {
DiffEngine diffEngine = new DiffEngine();
if ( diffEngine.canRead(masterFile) && diffEngine.canRead(testFile) ) {
DiffElement master = diffEngine.createDiffableFromFile(masterFile);
DiffElement test = diffEngine.createDiffableFromFile(testFile);
List<Difference> diffs = diffEngine.diff(master, test);
diffEngine.reportSummarizedDifferences(diffs, params);
return true;
} else {
return false;
}
}
public static class SummaryReportParams {
PrintStream out = System.out;
int maxItemsToDisplay = 0;
int maxCountOneItems = 0;
int minSumDiffToShow = 0;
public SummaryReportParams(PrintStream out, int maxItemsToDisplay, int maxCountOneItems, int minSumDiffToShow) {
this.out = out;
this.maxItemsToDisplay = maxItemsToDisplay;
this.maxCountOneItems = maxCountOneItems;
this.minSumDiffToShow = minSumDiffToShow;
}
}
}

View File

@ -0,0 +1,239 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
import com.google.java.contract.Requires;
import org.broadinstitute.sting.utils.Utils;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import java.util.*;
/**
* Created by IntelliJ IDEA.
* User: depristo
* Date: 7/4/11
* Time: 12:55 PM
*
* An interface that must be implemented to allow us to calculate differences
* between structured objects
*/
public class DiffNode extends DiffValue {
private Map<String, DiffElement> getElementMap() {
return (Map<String, DiffElement>)super.getValue();
}
private static Map<String, DiffElement> emptyElements() { return new HashMap<String, DiffElement>(); }
private DiffNode(Map<String, DiffElement> elements) {
super(elements);
}
private DiffNode(DiffElement binding, Map<String, DiffElement> elements) {
super(binding, elements);
}
// ---------------------------------------------------------------------------
//
// constructors
//
// ---------------------------------------------------------------------------
public static DiffNode rooted(String name) {
return empty(name, DiffElement.ROOT);
}
public static DiffNode empty(String name, DiffElement parent) {
DiffNode df = new DiffNode(emptyElements());
DiffElement elt = new DiffElement(name, parent, df);
df.setBinding(elt);
return df;
}
public static DiffNode empty(String name, DiffValue parent) {
return empty(name, parent.getBinding());
}
// ---------------------------------------------------------------------------
//
// accessors
//
// ---------------------------------------------------------------------------
@Override
public boolean isAtomic() { return false; }
public Collection<String> getElementNames() {
return getElementMap().keySet();
}
public Collection<DiffElement> getElements() {
return getElementMap().values();
}
private Collection<DiffElement> getElements(boolean atomicOnly) {
List<DiffElement> elts = new ArrayList<DiffElement>();
for ( DiffElement elt : getElements() )
if ( (atomicOnly && elt.getValue().isAtomic()) || (! atomicOnly && elt.getValue().isCompound()))
elts.add(elt);
return elts;
}
public Collection<DiffElement> getAtomicElements() {
return getElements(true);
}
public Collection<DiffElement> getCompoundElements() {
return getElements(false);
}
public DiffElement getElement(String name) {
for ( DiffElement elt : getElements() )
if ( elt.getName().equals(name) )
return elt;
return null;
}
/**
* Returns true if name is bound in this node
* @param name
* @return
*/
public boolean hasElement(String name) {
return getElement(name) != null;
}
// ---------------------------------------------------------------------------
//
// add
//
// ---------------------------------------------------------------------------
@Requires("elt != null")
public void add(DiffElement elt) {
if ( getElementMap().containsKey(elt.getName()) )
throw new IllegalArgumentException("Attempting to rebind already existing binding: " + elt + " node=" + this);
getElementMap().put(elt.getName(), elt);
}
@Requires("elt != null")
public void add(DiffValue elt) {
add(elt.getBinding());
}
@Requires("elts != null")
public void add(Collection<DiffElement> elts) {
for ( DiffElement e : elts )
add(e);
}
public void add(String name, Object value) {
add(new DiffElement(name, this.getBinding(), new DiffValue(value)));
}
// ---------------------------------------------------------------------------
//
// toString
//
// ---------------------------------------------------------------------------
@Override
public String toString() {
return toString(0);
}
@Override
public String toString(int offset) {
String off = offset > 0 ? Utils.dupString(' ', offset) : "";
StringBuilder b = new StringBuilder();
b.append("(").append("\n");
Collection<DiffElement> atomicElts = getAtomicElements();
for ( DiffElement elt : atomicElts ) {
b.append(elt.toString(offset + 2)).append('\n');
}
for ( DiffElement elt : getCompoundElements() ) {
b.append(elt.toString(offset + 4)).append('\n');
}
b.append(off).append(")").append("\n");
return b.toString();
}
@Override
public String toOneLineString() {
StringBuilder b = new StringBuilder();
b.append('(');
List<String> parts = new ArrayList<String>();
for ( DiffElement elt : getElements() )
parts.add(elt.toOneLineString());
b.append(Utils.join(" ", parts));
b.append(')');
return b.toString();
}
// --------------------------------------------------------------------------------
//
// fromString and toOneLineString
//
// --------------------------------------------------------------------------------
public static DiffElement fromString(String tree) {
return fromString(tree, DiffElement.ROOT);
}
/**
* Doesn't support full tree structure parsing
* @param tree
* @param parent
* @return
*/
private static DiffElement fromString(String tree, DiffElement parent) {
// X=(A=A B=B C=(D=D))
String[] parts = tree.split("=", 2);
if ( parts.length != 2 )
throw new ReviewedStingException("Unexpected tree structure: " + tree + " parts=" + parts);
String name = parts[0];
String value = parts[1];
if ( value.length() == 0 )
throw new ReviewedStingException("Illegal tree structure: " + value + " at " + tree);
if ( value.charAt(0) == '(' ) {
if ( ! value.endsWith(")") )
throw new ReviewedStingException("Illegal tree structure. Missing ): " + value + " at " + tree);
String subtree = value.substring(1, value.length()-1);
DiffNode rec = DiffNode.empty(name, parent);
String[] subParts = subtree.split(" ");
for ( String subPart : subParts ) {
rec.add(fromString(subPart, rec.getBinding()));
}
return rec.getBinding();
} else {
return new DiffValue(name, parent, value).getBinding();
}
}
}

View File

@ -0,0 +1,113 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
import org.apache.xmlbeans.impl.tool.Diff;
import org.broadinstitute.sting.commandline.Argument;
import org.broadinstitute.sting.commandline.Output;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.Requires;
import org.broadinstitute.sting.gatk.walkers.RodWalker;
import java.io.File;
import java.io.PrintStream;
import java.util.List;
/**
* Compares two record-oriented files, itemizing specific difference between equivalent
* records in the two files. Reports both itemized and summarized differences.
* @author Mark DePristo
* @version 0.1
*/
@Requires(value={})
public class DiffObjectsWalker extends RodWalker<Integer, Integer> {
@Output(doc="File to which results should be written",required=true)
protected PrintStream out;
@Argument(fullName="maxRecords", shortName="M", doc="Max. number of records to process", required=false)
int MAX_RECORDS = 0;
@Argument(fullName="maxCount1Records", shortName="M1", doc="Max. number of records occuring exactly once in the file to process", required=false)
int MAX_COUNT1_RECORDS = 0;
@Argument(fullName="minCountForDiff", shortName="MCFD", doc="Min number of observations for a records to display", required=false)
int minCountForDiff = 1;
@Argument(fullName="showItemizedDifferences", shortName="SID", doc="Should we enumerate all differences between the files?", required=false)
boolean showItemizedDifferences = false;
@Argument(fullName="master", shortName="m", doc="Master file: expected results", required=true)
File masterFile;
@Argument(fullName="test", shortName="t", doc="Test file: new results to compare to the master file", required=true)
File testFile;
final DiffEngine diffEngine = new DiffEngine();
@Override
public void initialize() {
}
@Override
public Integer map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
return 0;
}
@Override
public Integer reduceInit() {
return 0;
}
@Override
public Integer reduce(Integer counter, Integer sum) {
return counter + sum;
}
@Override
public void onTraversalDone(Integer sum) {
out.printf("Reading master file %s%n", masterFile);
DiffElement master = diffEngine.createDiffableFromFile(masterFile);
out.printf("Reading test file %s%n", testFile);
DiffElement test = diffEngine.createDiffableFromFile(testFile);
// out.printf("Master diff objects%n");
// out.println(master.toString());
// out.printf("Test diff objects%n");
// out.println(test.toString());
List<Difference> diffs = diffEngine.diff(master, test);
if ( showItemizedDifferences ) {
out.printf("Itemized results%n");
for ( Difference diff : diffs )
out.printf("DIFF: %s%n", diff.toString());
}
DiffEngine.SummaryReportParams params = new DiffEngine.SummaryReportParams(out, MAX_RECORDS, MAX_COUNT1_RECORDS, minCountForDiff);
diffEngine.reportSummarizedDifferences(diffs, params);
}
}

View File

@ -0,0 +1,90 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
import org.broadinstitute.sting.utils.Utils;
/**
* Created by IntelliJ IDEA.
* User: depristo
* Date: 7/4/11
* Time: 12:55 PM
*
* An interface that must be implemented to allow us to calculate differences
* between structured objects
*/
public class DiffValue {
private DiffElement binding = null;
final private Object value;
public DiffValue(Object value) {
this.value = value;
}
public DiffValue(DiffElement binding, Object value) {
this.binding = binding;
this.value = value;
}
public DiffValue(DiffValue parent, Object value) {
this(parent.getBinding(), value);
}
public DiffValue(String name, DiffElement parent, Object value) {
this.binding = new DiffElement(name, parent, this);
this.value = value;
}
public DiffValue(String name, DiffValue parent, Object value) {
this(name, parent.getBinding(), value);
}
public DiffElement getBinding() {
return binding;
}
protected void setBinding(DiffElement binding) {
this.binding = binding;
}
public Object getValue() {
return value;
}
public String toString() {
return getValue().toString();
}
public String toString(int offset) {
return toString();
}
public String toOneLineString() {
return getValue().toString();
}
public boolean isAtomic() { return true; }
public boolean isCompound() { return ! isAtomic(); }
}

View File

@ -0,0 +1,50 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
import com.google.java.contract.Ensures;
import com.google.java.contract.Requires;
import java.io.File;
/**
* Created by IntelliJ IDEA.
* User: depristo
* Date: 7/4/11
* Time: 1:09 PM
*
* Interface for readers creating diffable objects from a file
*/
public interface DiffableReader {
@Ensures("result != null")
public String getName();
@Ensures("result != null")
@Requires("file != null")
public DiffElement readFromFile(File file);
@Requires("file != null")
public boolean canRead(File file);
}

View File

@ -0,0 +1,58 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
/**
* Created by IntelliJ IDEA.
* User: depristo
* Date: 7/4/11
* Time: 12:53 PM
*
* Represents a specific difference between two specific DiffElements
*/
public class Difference {
DiffElement master, test;
public Difference(DiffElement master, DiffElement test) {
if ( master == null && test == null ) throw new IllegalArgumentException("Master and test both cannot be null");
this.master = master;
this.test = test;
}
public String toString() {
return String.format("%s:%s!=%s",
getFullyQualifiedName(),
getOneLineString(master),
getOneLineString(test));
}
public String getFullyQualifiedName() {
return (master == null ? test : master).fullyQualifiedName();
}
private static String getOneLineString(DiffElement elt) {
return elt == null ? "MISSING" : elt.getValue().toOneLineString();
}
}

View File

@ -0,0 +1,119 @@
/*
* Copyright (c) 2011, The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.diffengine;
import org.broad.tribble.readers.AsciiLineReader;
import org.broad.tribble.readers.LineReader;
import org.broadinstitute.sting.utils.codecs.vcf.VCFCodec;
import org.broadinstitute.sting.utils.codecs.vcf.VCFConstants;
import org.broadinstitute.sting.utils.codecs.vcf.VCFHeader;
import org.broadinstitute.sting.utils.variantcontext.Genotype;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import java.io.*;
import java.util.Arrays;
import java.util.Map;
import java.util.zip.GZIPInputStream;
/**
* Created by IntelliJ IDEA.
* User: depristo
* Date: 7/4/11
* Time: 1:09 PM
*
* Class implementing diffnode reader for VCF
*/
public class VCFDiffableReader implements DiffableReader {
@Override
public String getName() { return "VCF"; }
@Override
public DiffElement readFromFile(File file) {
DiffNode root = DiffNode.rooted(file.getName());
try {
LineReader lineReader = new AsciiLineReader(new FileInputStream(file));
VCFCodec vcfCodec = new VCFCodec();
VCFHeader header = (VCFHeader)vcfCodec.readHeader(lineReader);
String line = lineReader.readLine();
while ( line != null ) {
VariantContext vc = (VariantContext)vcfCodec.decode(line);
String name = vc.getChr() + ":" + vc.getStart();
DiffNode vcRoot = DiffNode.empty(name, root);
// add fields
vcRoot.add("CHROM", vc.getChr());
vcRoot.add("POS", vc.getStart());
vcRoot.add("ID", vc.hasID() ? vc.getID() : VCFConstants.MISSING_VALUE_v4);
vcRoot.add("REF", vc.getReference());
vcRoot.add("ALT", vc.getAlternateAlleles());
vcRoot.add("QUAL", vc.hasNegLog10PError() ? vc.getNegLog10PError() * 10 : VCFConstants.MISSING_VALUE_v4);
vcRoot.add("FILTER", vc.getFilters());
// add info fields
for (Map.Entry<String, Object> attribute : vc.getAttributes().entrySet()) {
if ( ! attribute.getKey().startsWith("_") && ! attribute.getKey().equals(VariantContext.ID_KEY))
vcRoot.add(attribute.getKey(), attribute.getValue());
}
for (Genotype g : vc.getGenotypes().values() ) {
DiffNode gRoot = DiffNode.empty(g.getSampleName(), vcRoot);
gRoot.add("GT", g.getGenotypeString());
gRoot.add("GQ", g.hasNegLog10PError() ? g.getNegLog10PError() * 10 : VCFConstants.MISSING_VALUE_v4 );
for (Map.Entry<String, Object> attribute : g.getAttributes().entrySet()) {
if ( ! attribute.getKey().startsWith("_") )
gRoot.add(attribute.getKey(), attribute.getValue());
}
vcRoot.add(gRoot);
}
root.add(vcRoot);
line = lineReader.readLine();
}
lineReader.close();
} catch ( IOException e ) {
return null;
}
return root.getBinding();
}
@Override
public boolean canRead(File file) {
try {
final String VCF4_HEADER = "##fileformat=VCFv4";
char[] buff = new char[VCF4_HEADER.length()];
new FileReader(file).read(buff, 0, VCF4_HEADER.length());
String firstLine = new String(buff);
return firstLine.startsWith(VCF4_HEADER);
} catch ( IOException e ) {
return false;
}
}
}

View File

@ -92,25 +92,31 @@ public class UnifiedArgumentCollection {
@Argument(fullName = "indel_heterozygosity", shortName = "indelHeterozygosity", doc = "Heterozygosity for indel calling", required = false)
public double INDEL_HETEROZYGOSITY = 1.0/8000;
@Hidden
@Argument(fullName = "indelGapContinuationPenalty", shortName = "indelGCP", doc = "Indel gap continuation penalty", required = false)
public double INDEL_GAP_CONTINUATION_PENALTY = 10.0;
@Hidden
@Argument(fullName = "indelGapOpenPenalty", shortName = "indelGOP", doc = "Indel gap open penalty", required = false)
public double INDEL_GAP_OPEN_PENALTY = 45.0;
@Hidden
@Argument(fullName = "indelHaplotypeSize", shortName = "indelHSize", doc = "Indel haplotype size", required = false)
public int INDEL_HAPLOTYPE_SIZE = 80;
@Hidden
@Argument(fullName = "doContextDependentGapPenalties", shortName = "doCDP", doc = "Vary gap penalties by context", required = false)
public boolean DO_CONTEXT_DEPENDENT_PENALTIES = true;
//gdebug+
@Hidden
// experimental arguments, NOT TO BE USED BY ANYONE WHOSE INITIALS AREN'T GDA!!!
@Hidden
@Argument(fullName = "getGapPenaltiesFromData", shortName = "dataGP", doc = "Vary gap penalties by context - EXPERIMENTAL, DO NO USE", required = false)
public boolean GET_GAP_PENALTIES_FROM_DATA = false;
@Hidden
@Argument(fullName="indel_recal_file", shortName="recalFile", required=false, doc="Filename for the input covariates table recalibration .csv file - EXPERIMENTAL, DO NO USE")
public File INDEL_RECAL_FILE = new File("indel.recal_data.csv");
@Hidden
@Argument(fullName = "indelDebug", shortName = "indelDebug", doc = "Output indel debug info", required = false)
public boolean OUTPUT_DEBUG_INDEL_INFO = false;
@Hidden

View File

@ -37,7 +37,6 @@ import org.broadinstitute.sting.gatk.datasources.rmd.ReferenceOrderedDataSource;
import org.broadinstitute.sting.utils.*;
import org.broadinstitute.sting.utils.baq.BAQ;
import org.broadinstitute.sting.commandline.*;
import org.broadinstitute.sting.utils.codecs.vcf.VCFUtils;
import java.util.*;
import java.io.PrintStream;
@ -158,7 +157,7 @@ public class UnifiedGenotyper extends LocusWalker<VariantCallContext, UnifiedGen
}
// FORMAT and INFO fields
headerInfo.addAll(VCFUtils.getSupportedHeaderStrings());
headerInfo.addAll(getSupportedHeaderStrings());
// FILTER fields
if ( UAC.STANDARD_CONFIDENCE_FOR_EMITTING < UAC.STANDARD_CONFIDENCE_FOR_CALLING )
@ -167,6 +166,20 @@ public class UnifiedGenotyper extends LocusWalker<VariantCallContext, UnifiedGen
return headerInfo;
}
/**
* return a set of supported format lines; what we currently support for output in the genotype fields of a VCF
* @return a set of VCF format lines
*/
private static Set<VCFFormatHeaderLine> getSupportedHeaderStrings() {
Set<VCFFormatHeaderLine> result = new HashSet<VCFFormatHeaderLine>();
result.add(new VCFFormatHeaderLine(VCFConstants.GENOTYPE_KEY, 1, VCFHeaderLineType.String, "Genotype"));
result.add(new VCFFormatHeaderLine(VCFConstants.GENOTYPE_QUALITY_KEY, 1, VCFHeaderLineType.Float, "Genotype Quality"));
result.add(new VCFFormatHeaderLine(VCFConstants.DEPTH_KEY, 1, VCFHeaderLineType.Integer, "Read Depth (only filtered reads used for calling)"));
result.add(new VCFFormatHeaderLine(VCFConstants.PHRED_GENOTYPE_LIKELIHOODS_KEY, VCFHeaderLineCount.UNBOUNDED, VCFHeaderLineType.Float, "Normalized, Phred-scaled likelihoods for AA,AB,BB genotypes where A=ref and B=alt; if site is not biallelic, number of likelihoods if n*(n+1)/2"));
return result;
}
/**
* Compute at a given locus.
*

View File

@ -5,6 +5,7 @@ import net.sf.samtools.*;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.broadinstitute.sting.utils.exceptions.UserException;
import java.util.*;
@ -113,9 +114,10 @@ public class ConstrainedMateFixingManager {
HashMap<String, SAMRecordHashObject> forMateMatching = new HashMap<String, SAMRecordHashObject>();
TreeSet<SAMRecord> waitingReads = new TreeSet<SAMRecord>(comparer);
private <T> T remove(TreeSet<T> treeSet) {
final T first = treeSet.first();
treeSet.remove(first);
private SAMRecord remove(TreeSet<SAMRecord> treeSet) {
final SAMRecord first = treeSet.first();
if ( !treeSet.remove(first) )
throw new UserException("Error caching SAM record " + first.getReadName() + ", which is usually caused by malformed SAM/BAM files in which multiple identical copies of a read are present.");
return first;
}

View File

@ -43,9 +43,9 @@ public class AlleleCount extends VariantStratifier {
if (eval != null) {
int AC = -1;
if ( eval.hasAttribute("AC") )
if ( eval.hasAttribute("AC") && eval.getAttribute("AC") instanceof Integer ) {
AC = eval.getAttributeAsInt("AC");
else if ( eval.isVariant() ) {
} else if ( eval.isVariant() ) {
for (Allele allele : eval.getAlternateAlleles())
AC = Math.max(AC, eval.getChromosomeCount(allele));
} else

View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 2010 The Broad Institute
* Copyright (c) 2011 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation

View File

@ -1,3 +1,28 @@
/*
* Copyright (c) 2011 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.variantrecalibration;
import Jama.Matrix;
@ -72,7 +97,7 @@ public class GaussianMixtureModel {
int ttt = 0;
while( ttt++ < numIterations ) {
// Estep: assign each variant to the nearest cluster
// E step: assign each variant to the nearest cluster
for( final VariantDatum datum : data ) {
double minDistance = Double.MAX_VALUE;
MultivariateGaussian minGaussian = null;
@ -87,7 +112,7 @@ public class GaussianMixtureModel {
datum.assignment = minGaussian;
}
// Mstep: update gaussian means based on assigned variants
// M step: update gaussian means based on assigned variants
for( final MultivariateGaussian gaussian : gaussians ) {
gaussian.zeroOutMu();
int numAssigned = 0;
@ -204,26 +229,29 @@ public class GaussianMixtureModel {
}
public double evaluateDatumMarginalized( final VariantDatum datum ) {
int numVals = 0;
int numSamples = 0;
double sumPVarInGaussian = 0.0;
int numIter = 10;
final int numIterPerMissingAnnotation = 10; // Trade off here between speed of computation and accuracy of the marginalization
final double[] pVarInGaussianLog10 = new double[gaussians.size()];
// for each dimension
for( int iii = 0; iii < datum.annotations.length; iii++ ) {
// marginalize over the missing dimension by drawing X random values for the missing annotation and averaging the lod
// if it is missing marginalize over the missing dimension by drawing X random values for the missing annotation and averaging the lod
if( datum.isNull[iii] ) {
for( int ttt = 0; ttt < numIter; ttt++ ) {
datum.annotations[iii] = Normal.staticNextDouble(0.0, 1.0);
for( int ttt = 0; ttt < numIterPerMissingAnnotation; ttt++ ) {
datum.annotations[iii] = GenomeAnalysisEngine.getRandomGenerator().nextGaussian(); // draw a random sample from the standard normal distribution
// evaluate this random data point
int gaussianIndex = 0;
for( final MultivariateGaussian gaussian : gaussians ) {
pVarInGaussianLog10[gaussianIndex++] = gaussian.pMixtureLog10 + gaussian.evaluateDatumLog10( datum );
}
sumPVarInGaussian += Math.pow(10.0, MathUtils.log10sumLog10(pVarInGaussianLog10));
numVals++;
// add this sample's probability to the pile in order to take an average in the end
sumPVarInGaussian += Math.pow(10.0, MathUtils.log10sumLog10(pVarInGaussianLog10)); // p = 10 ^ Sum(pi_k * p(v|n,k))
numSamples++;
}
}
}
return Math.log10( sumPVarInGaussian / ((double) numVals) );
return Math.log10( sumPVarInGaussian / ((double) numSamples) );
}
}

View File

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

View File

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

View File

@ -1,5 +1,5 @@
/*
* Copyright (c) 2010 The Broad Institute
* Copyright (c) 2011 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation

View File

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

View File

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

View File

@ -1,6 +1,30 @@
/*
* Copyright (c) 2011 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.variantrecalibration;
import cern.jet.random.Normal;
import org.apache.log4j.Logger;
import org.broadinstitute.sting.utils.variantcontext.VariantContext;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
@ -58,19 +82,11 @@ public class VariantDataManager {
}
foundZeroVarianceAnnotation = foundZeroVarianceAnnotation || (theSTD < 1E-6);
if( annotationKeys.get(iii).toLowerCase().contains("ranksum") ) { // BUGBUG: to clean up
for( final VariantDatum datum : data ) {
if( datum.annotations[iii] > 0.0 ) { datum.annotations[iii] /= 3.0; }
}
}
meanVector[iii] = theMean;
varianceVector[iii] = theSTD;
for( final VariantDatum datum : data ) {
datum.annotations[iii] = ( datum.isNull[iii] ? Normal.staticNextDouble(0.0, 1.0) : ( datum.annotations[iii] - theMean ) / theSTD );
// Each data point is now [ (x - mean) / standard deviation ]
if( annotationKeys.get(iii).toLowerCase().contains("ranksum") && datum.isNull[iii] && datum.annotations[iii] > 0.0 ) {
datum.annotations[iii] /= 3.0;
}
// Transform each data point via: (x - mean) / standard deviation
datum.annotations[iii] = ( datum.isNull[iii] ? GenomeAnalysisEngine.getRandomGenerator().nextGaussian() : ( datum.annotations[iii] - theMean ) / theSTD );
}
}
if( foundZeroVarianceAnnotation ) {
@ -139,7 +155,7 @@ public class VariantDataManager {
final int numBadSitesAdded = trainingData.size();
logger.info( "Found " + numBadSitesAdded + " variants overlapping bad sites training tracks." );
// Next, sort the variants by the LOD coming from the positive model and add to the list the bottom X percent of variants
// Next sort the variants by the LOD coming from the positive model and add to the list the bottom X percent of variants
Collections.sort( data );
final int numToAdd = Math.max( minimumNumber - trainingData.size(), Math.round((float)bottomPercentage * data.size()) );
if( numToAdd > data.size() ) {
@ -217,23 +233,15 @@ public class VariantDataManager {
double value;
try {
if( annotationKey.equalsIgnoreCase("QUAL") ) {
value = vc.getPhredScaledQual();
} else if( annotationKey.equalsIgnoreCase("DP") ) {
value = Double.parseDouble( (String)vc.getAttribute( "DP" ) ) / Double.parseDouble( (String)vc.getAttribute( "AN" ) );
} else {
value = Double.parseDouble( (String)vc.getAttribute( annotationKey ) );
if( Double.isInfinite(value) ) { value = Double.NaN; }
if( annotationKey.equalsIgnoreCase("InbreedingCoeff") && value > 0.05 ) { value = Double.NaN; }
if( jitter && annotationKey.equalsIgnoreCase("HRUN") ) { // Integer valued annotations must be jittered a bit to work in this GMM
value += -0.25 + 0.5 * GenomeAnalysisEngine.getRandomGenerator().nextDouble();
}
if( annotationKey.equalsIgnoreCase("HaplotypeScore") && MathUtils.compareDoubles(value, 0.0, 0.0001) == 0 ) { value = -0.2 + 0.4*GenomeAnalysisEngine.getRandomGenerator().nextDouble(); }
if( annotationKey.equalsIgnoreCase("FS") && MathUtils.compareDoubles(value, 0.0, 0.01) == 0 ) { value = -0.2 + 0.4*GenomeAnalysisEngine.getRandomGenerator().nextDouble(); }
value = Double.parseDouble( (String)vc.getAttribute( annotationKey ) );
if( Double.isInfinite(value) ) { value = Double.NaN; }
if( jitter && annotationKey.equalsIgnoreCase("HRUN") ) { // Integer valued annotations must be jittered a bit to work in this GMM
value += -0.25 + 0.5 * GenomeAnalysisEngine.getRandomGenerator().nextDouble();
}
if( jitter && annotationKey.equalsIgnoreCase("HaplotypeScore") && MathUtils.compareDoubles(value, 0.0, 0.0001) == 0 ) { value = -0.2 + 0.4*GenomeAnalysisEngine.getRandomGenerator().nextDouble(); }
if( jitter && annotationKey.equalsIgnoreCase("FS") && MathUtils.compareDoubles(value, 0.0, 0.001) == 0 ) { value = -0.2 + 0.4*GenomeAnalysisEngine.getRandomGenerator().nextDouble(); }
} catch( Exception e ) {
value = Double.NaN; // The VQSR works with missing data now by marginalizing over the missing dimension when evaluating Gaussians
value = Double.NaN; // The VQSR works with missing data by marginalizing over the missing dimension when evaluating the Gaussian mixture model
}
return value;

View File

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

View File

@ -88,7 +88,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
@Argument(fullName="rscript_file", shortName="rscriptFile", doc="The output rscript file generated by the VQSR to aid in visualization of the input data and learned model", required=false)
private String RSCRIPT_FILE = null;
@Argument(fullName = "path_to_resources", shortName = "resources", doc = "Path to resources folder holding the Sting R scripts.", required=false)
private String PATH_TO_RESOURCES = "R/";
private String PATH_TO_RESOURCES = "public/R/";
@Argument(fullName="ts_filter_level", shortName="ts_filter_level", doc="The truth sensitivity level at which to start filtering, used here to indicate filtered variants in plots", required=false)
private double TS_FILTER_LEVEL = 99.0;
@ -118,6 +118,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
//---------------------------------------------------------------------------------------------------------------
public void initialize() {
if( !PATH_TO_RESOURCES.endsWith("/") ) { PATH_TO_RESOURCES = PATH_TO_RESOURCES + "/"; }
dataManager = new VariantDataManager( new ArrayList<String>(Arrays.asList(USE_ANNOTATIONS)), VRAC );
if( IGNORE_INPUT_FILTERS != null ) {
@ -228,19 +229,23 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
public void onTraversalDone( final ExpandingArrayList<VariantDatum> reduceSum ) {
dataManager.setData( reduceSum );
dataManager.normalizeData(); // Each data point is now (x - mean) / standard deviation
// Generate the positive model using the training data and evaluate each variant
final GaussianMixtureModel goodModel = engine.generateModel( dataManager.getTrainingData() );
engine.evaluateData( dataManager.getData(), goodModel, false );
// Generate the negative model using the worst performing data and evaluate each variant contrastively
final GaussianMixtureModel badModel = engine.generateModel( dataManager.selectWorstVariants( VRAC.PERCENT_BAD_VARIANTS, VRAC.MIN_NUM_BAD_VARIANTS ) );
engine.evaluateData( dataManager.getData(), badModel, true );
engine.calculateWorstPerformingAnnotation( dataManager.getData(), goodModel, badModel );
final ExpandingArrayList<VariantDatum> randomData = dataManager.getRandomDataForPlotting( 6000 );
// Find the VQSLOD cutoff values which correspond to the various tranches of calls requested by the user
final int nCallsAtTruth = TrancheManager.countCallsAtTruth( dataManager.getData(), Double.NEGATIVE_INFINITY );
final TrancheManager.SelectionMetric metric = new TrancheManager.TruthSensitivityMetric( nCallsAtTruth );
final List<Tranche> tranches = TrancheManager.findTranches( dataManager.getData(), TS_TRANCHES, metric );
tranchesStream.print(Tranche.tranchesString( tranches ));
// Find the filtering lodCutoff for display on the model PDFs. Red variants are those which were below the cutoff and filtered out of the final callset.
double lodCutoff = 0.0;
for( final Tranche tranche : tranches ) {
if( MathUtils.compareDoubles(tranche.ts, TS_FILTER_LEVEL, 0.0001)==0 ) {
@ -252,7 +257,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
dataManager.writeOutRecalibrationTable( RECAL_FILE );
if( RSCRIPT_FILE != null ) {
logger.info( "Writing out visualization Rscript file...");
createVisualizationScript( randomData, goodModel, badModel, lodCutoff );
createVisualizationScript( dataManager.getRandomDataForPlotting( 6000 ), goodModel, badModel, lodCutoff );
}
// Execute Rscript command to create the tranche plot
@ -278,6 +283,8 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
} catch( FileNotFoundException e ) {
throw new UserException.CouldNotCreateOutputFile(RSCRIPT_FILE, "", e);
}
// We make extensive use of the ggplot2 R library: http://had.co.nz/ggplot2/
stream.println("library(ggplot2)");
createArrangeFunction( stream );
@ -371,6 +378,7 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
}
}
// The Arrange function is how we place the 4 model plots on one page
// from http://gettinggeneticsdone.blogspot.com/2010/03/arrange-multiple-ggplot2-plots-in-same.html
private void createArrangeFunction( final PrintStream stream ) {
stream.println("vp.layout <- function(x, y) viewport(layout.pos.row=x, layout.pos.col=y)");
@ -395,5 +403,4 @@ public class VariantRecalibrator extends RodWalker<ExpandingArrayList<VariantDat
stream.println("}");
stream.println("}");
}
}
}

View File

@ -1,3 +1,28 @@
/*
* Copyright (c) 2011 The Broad Institute
*
* Permission is hereby granted, free of charge, to any person
* obtaining a copy of this software and associated documentation
* files (the "Software"), to deal in the Software without
* restriction, including without limitation the rights to use,
* copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following
* conditions:
*
* The above copyright notice and this permission notice shall be
* included in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR
* THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
package org.broadinstitute.sting.gatk.walkers.variantrecalibration;
import org.broadinstitute.sting.commandline.Argument;
@ -35,7 +60,7 @@ public class VariantRecalibratorArgumentCollection {
@Argument(fullName="priorCounts", shortName="priorCounts", doc="The number of prior counts to use in variational Bayes algorithm.", required=false)
public double PRIOR_COUNTS = 20.0;
@Argument(fullName="percentBadVariants", shortName="percentBad", doc="What percentage of the worst scoring variants to use when building the Gaussian mixture model of bad variants. 0.07 means bottom 7 percent.", required=false)
public double PERCENT_BAD_VARIANTS = 0.015;
public double PERCENT_BAD_VARIANTS = 0.03;
@Argument(fullName="minNumBadVariants", shortName="minNumBad", doc="The minimum amount of worst scoring variants to use when building the Gaussian mixture model of bad variants. Will override -percentBad arugment if necessary.", required=false)
public int MIN_NUM_BAD_VARIANTS = 2000;
}

View File

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

View File

@ -225,7 +225,7 @@ public abstract class AbstractVCFCodec implements FeatureCodec, NameAwareCodec,
loc = pos + alleles.get(0).length() - 1;
} else if ( !isSingleNucleotideEvent(alleles) ) {
ArrayList<Allele> newAlleles = new ArrayList<Allele>();
loc = clipAlleles(pos, ref, alleles, newAlleles);
loc = clipAlleles(pos, ref, alleles, newAlleles, lineNo);
alleles = newAlleles;
}
@ -504,7 +504,7 @@ public abstract class AbstractVCFCodec implements FeatureCodec, NameAwareCodec,
* @param clippedAlleles output list of clipped alleles
* @return a list of alleles, clipped to the reference
*/
protected static long clipAlleles(long position, String ref, List<Allele> unclippedAlleles, List<Allele> clippedAlleles) {
protected static long clipAlleles(long position, String ref, List<Allele> unclippedAlleles, List<Allele> clippedAlleles, int lineNo) {
// Note that the computation of forward clipping here is meant only to see whether there is a common
// base to all alleles, and to correctly compute reverse clipping,
@ -522,6 +522,8 @@ public abstract class AbstractVCFCodec implements FeatureCodec, NameAwareCodec,
}
if (a.length() - reverseClipped <= forwardClipping || a.length() - forwardClipping == 0)
clipping = false;
else if (ref.length() == reverseClipped)
generateException("bad alleles encountered", lineNo);
else if (a.getBases()[a.length()-reverseClipped-1] != ref.getBytes()[ref.length()-reverseClipped-1])
clipping = false;
}

View File

@ -358,16 +358,8 @@ public class StandardVCFWriter implements VCFWriter {
mWriter.write(key);
if ( !entry.getValue().equals("") ) {
int numVals = 1;
VCFInfoHeaderLine metaData = mHeader.getInfoHeaderLine(key);
if ( metaData != null )
numVals = metaData.getCount();
// take care of unbounded encoding
if ( numVals == VCFInfoHeaderLine.UNBOUNDED )
numVals = 1;
if ( numVals > 0 ) {
if ( metaData == null || metaData.getCountType() != VCFHeaderLineCount.INTEGER || metaData.getCount() != 0 ) {
mWriter.write("=");
mWriter.write(entry.getValue());
}
@ -423,7 +415,7 @@ public class StandardVCFWriter implements VCFWriter {
VCFFormatHeaderLine metaData = mHeader.getFormatHeaderLine(key);
if ( metaData != null ) {
int numInFormatField = metaData.getCount();
int numInFormatField = metaData.getCount(vc.getAlternateAlleles().size());
if ( numInFormatField > 1 && val.equals(VCFConstants.MISSING_VALUE_v4) ) {
// If we have a missing field but multiple values are expected, we need to construct a new string with all fields.
// For example, if Number=2, the string has to be ".,."

View File

@ -0,0 +1,28 @@
package org.broadinstitute.sting.utils.codecs.vcf;
/**
* @author ebanks
* A class representing a key=value entry for ALT fields in the VCF header
*/
public class VCFAltHeaderLine extends VCFSimpleHeaderLine {
/**
* create a VCF filter header line
*
* @param name the name for this header line
* @param description the description for this header line
*/
public VCFAltHeaderLine(String name, String description) {
super(name, description, SupportedHeaderLineType.ALT);
}
/**
* create a VCF info header line
*
* @param line the header line
* @param version the vcf header version
*/
protected VCFAltHeaderLine(String line, VCFHeaderVersion version) {
super(line, version, SupportedHeaderLineType.ALT);
}
}

View File

@ -24,6 +24,8 @@
package org.broadinstitute.sting.utils.codecs.vcf;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import java.util.Arrays;
import java.util.LinkedHashMap;
import java.util.Map;
@ -43,26 +45,43 @@ public abstract class VCFCompoundHeaderLine extends VCFHeaderLine implements VCF
// the field types
private String name;
private int count;
private int count = -1;
private VCFHeaderLineCount countType;
private String description;
private VCFHeaderLineType type;
// access methods
public String getName() { return name; }
public int getCount() { return count; }
public String getDescription() { return description; }
public VCFHeaderLineType getType() { return type; }
public VCFHeaderLineCount getCountType() { return countType; }
public int getCount() {
if ( countType != VCFHeaderLineCount.INTEGER )
throw new ReviewedStingException("Asking for header line count when type is not an integer");
return count;
}
//
public void setNumberToUnbounded() { this.count = UNBOUNDED; }
// utility method
public int getCount(int numAltAlleles) {
int myCount;
switch ( countType ) {
case INTEGER: myCount = count; break;
case UNBOUNDED: myCount = -1; break;
case A: myCount = numAltAlleles; break;
case G: myCount = ((numAltAlleles + 1) * (numAltAlleles + 2) / 2); break;
default: throw new ReviewedStingException("Unknown count type: " + countType);
}
return myCount;
}
public void setNumberToUnbounded() {
countType = VCFHeaderLineCount.UNBOUNDED;
count = -1;
}
// our type of line, i.e. format, info, etc
private final SupportedHeaderLineType lineType;
// line numerical values are allowed to be unbounded (or unknown), which is
// marked with a dot (.)
public static final int UNBOUNDED = -1; // the value we store internally for unbounded types
/**
* create a VCF format header line
*
@ -70,10 +89,12 @@ public abstract class VCFCompoundHeaderLine extends VCFHeaderLine implements VCF
* @param count the count for this header line
* @param type the type for this header line
* @param description the description for this header line
* @param lineType the header line type
*/
protected VCFCompoundHeaderLine(String name, int count, VCFHeaderLineType type, String description, SupportedHeaderLineType lineType) {
super(lineType.toString(), "");
this.name = name;
this.countType = VCFHeaderLineCount.INTEGER;
this.count = count;
this.type = type;
this.description = description;
@ -81,20 +102,53 @@ public abstract class VCFCompoundHeaderLine extends VCFHeaderLine implements VCF
validate();
}
/**
* create a VCF format header line
*
* @param name the name for this header line
* @param count the count type for this header line
* @param type the type for this header line
* @param description the description for this header line
* @param lineType the header line type
*/
protected VCFCompoundHeaderLine(String name, VCFHeaderLineCount count, VCFHeaderLineType type, String description, SupportedHeaderLineType lineType) {
super(lineType.toString(), "");
this.name = name;
this.countType = count;
this.type = type;
this.description = description;
this.lineType = lineType;
validate();
}
/**
* create a VCF format header line
*
* @param line the header line
* @param version the VCF header version
* @param lineType the header line type
*
*/
protected VCFCompoundHeaderLine(String line, VCFHeaderVersion version, SupportedHeaderLineType lineType) {
super(lineType.toString(), "");
Map<String,String> mapping = VCFHeaderLineTranslator.parseLine(version,line, Arrays.asList("ID","Number","Type","Description"));
name = mapping.get("ID");
count = (version == VCFHeaderVersion.VCF4_0 || version == VCFHeaderVersion.VCF4_1) ?
mapping.get("Number").equals(VCFConstants.UNBOUNDED_ENCODING_v4) ? UNBOUNDED : Integer.valueOf(mapping.get("Number")) :
mapping.get("Number").equals(VCFConstants.UNBOUNDED_ENCODING_v3) ? UNBOUNDED : Integer.valueOf(mapping.get("Number"));
count = -1;
final String numberStr = mapping.get("Number");
if ( numberStr.equals(VCFConstants.PER_ALLELE_COUNT) ) {
countType = VCFHeaderLineCount.A;
} else if ( numberStr.equals(VCFConstants.PER_GENOTYPE_COUNT) ) {
countType = VCFHeaderLineCount.G;
} else if ( ((version == VCFHeaderVersion.VCF4_0 || version == VCFHeaderVersion.VCF4_1) &&
numberStr.equals(VCFConstants.UNBOUNDED_ENCODING_v4)) ||
((version == VCFHeaderVersion.VCF3_2 || version == VCFHeaderVersion.VCF3_3) &&
numberStr.equals(VCFConstants.UNBOUNDED_ENCODING_v3)) ) {
countType = VCFHeaderLineCount.UNBOUNDED;
} else {
countType = VCFHeaderLineCount.INTEGER;
count = Integer.valueOf(numberStr);
}
type = VCFHeaderLineType.valueOf(mapping.get("Type"));
if (type == VCFHeaderLineType.Flag && !allowFlagValues())
throw new IllegalArgumentException("Flag is an unsupported type for this kind of field");
@ -121,7 +175,15 @@ public abstract class VCFCompoundHeaderLine extends VCFHeaderLine implements VCF
protected String toStringEncoding() {
Map<String,Object> map = new LinkedHashMap<String,Object>();
map.put("ID", name);
map.put("Number", count == UNBOUNDED ? VCFConstants.UNBOUNDED_ENCODING_v4 : count);
Object number;
switch ( countType ) {
case A: number = VCFConstants.PER_ALLELE_COUNT; break;
case G: number = VCFConstants.PER_GENOTYPE_COUNT; break;
case UNBOUNDED: number = VCFConstants.UNBOUNDED_ENCODING_v4; break;
case INTEGER:
default: number = count;
}
map.put("Number", number);
map.put("Type", type);
map.put("Description", description);
return lineType.toString() + "=" + VCFHeaderLine.toStringEncoding(map);
@ -136,15 +198,13 @@ public abstract class VCFCompoundHeaderLine extends VCFHeaderLine implements VCF
if ( !(o instanceof VCFCompoundHeaderLine) )
return false;
VCFCompoundHeaderLine other = (VCFCompoundHeaderLine)o;
return name.equals(other.name) &&
count == other.count &&
description.equals(other.description) &&
type == other.type &&
lineType == other.lineType;
return equalsExcludingDescription(other) &&
description.equals(other.description);
}
public boolean equalsExcludingDescription(VCFCompoundHeaderLine other) {
return count == other.count &&
countType == other.countType &&
type == other.type &&
lineType == other.lineType &&
name.equals(other.name);

View File

@ -99,6 +99,8 @@ public final class VCFConstants {
public static final String MISSING_DEPTH_v3 = "-1";
public static final String UNBOUNDED_ENCODING_v4 = ".";
public static final String UNBOUNDED_ENCODING_v3 = "-1";
public static final String PER_ALLELE_COUNT = "A";
public static final String PER_GENOTYPE_COUNT = "G";
public static final String EMPTY_ALLELE = ".";
public static final String EMPTY_GENOTYPE = "./.";
public static final double MAX_GENOTYPE_QUAL = 99.0;

View File

@ -1,19 +1,10 @@
package org.broadinstitute.sting.utils.codecs.vcf;
import java.util.Arrays;
import java.util.LinkedHashMap;
import java.util.Map;
/**
* @author ebanks
* A class representing a key=value entry for FILTER fields in the VCF header
*/
public class VCFFilterHeaderLine extends VCFHeaderLine implements VCFNamedHeaderLine {
private String name;
private String description;
public class VCFFilterHeaderLine extends VCFSimpleHeaderLine {
/**
* create a VCF filter header line
@ -22,12 +13,7 @@ public class VCFFilterHeaderLine extends VCFHeaderLine implements VCFNamedHeader
* @param description the description for this header line
*/
public VCFFilterHeaderLine(String name, String description) {
super("FILTER", "");
this.name = name;
this.description = description;
if ( name == null || description == null )
throw new IllegalArgumentException(String.format("Invalid VCFCompoundHeaderLine: key=%s name=%s desc=%s", super.getKey(), name, description ));
super(name, description, SupportedHeaderLineType.FILTER);
}
/**
@ -37,34 +23,6 @@ public class VCFFilterHeaderLine extends VCFHeaderLine implements VCFNamedHeader
* @param version the vcf header version
*/
protected VCFFilterHeaderLine(String line, VCFHeaderVersion version) {
super("FILTER", "");
Map<String,String> mapping = VCFHeaderLineTranslator.parseLine(version,line, Arrays.asList("ID","Description"));
name = mapping.get("ID");
description = mapping.get("Description");
if ( description == null && ALLOW_UNBOUND_DESCRIPTIONS ) // handle the case where there's no description provided
description = UNBOUND_DESCRIPTION;
}
protected String toStringEncoding() {
Map<String,Object> map = new LinkedHashMap<String,Object>();
map.put("ID", name);
map.put("Description", description);
return "FILTER=" + VCFHeaderLine.toStringEncoding(map);
}
public boolean equals(Object o) {
if ( !(o instanceof VCFFilterHeaderLine) )
return false;
VCFFilterHeaderLine other = (VCFFilterHeaderLine)o;
return name.equals(other.name) &&
description.equals(other.description);
}
public String getName() {
return name;
}
public String getDescription() {
return description;
super(line, version, SupportedHeaderLineType.FILTER);
}
}

View File

@ -16,6 +16,10 @@ public class VCFFormatHeaderLine extends VCFCompoundHeaderLine {
throw new IllegalArgumentException("Flag is an unsupported type for format fields");
}
public VCFFormatHeaderLine(String name, VCFHeaderLineCount count, VCFHeaderLineType type, String description) {
super(name, count, type, description, SupportedHeaderLineType.FORMAT);
}
protected VCFFormatHeaderLine(String line, VCFHeaderVersion version) {
super(line, version, SupportedHeaderLineType.FORMAT);
}

View File

@ -0,0 +1,8 @@
package org.broadinstitute.sting.utils.codecs.vcf;
/**
* the count encodings we use for fields in VCF header lines
*/
public enum VCFHeaderLineCount {
INTEGER, A, G, UNBOUNDED;
}

View File

@ -13,6 +13,10 @@ public class VCFInfoHeaderLine extends VCFCompoundHeaderLine {
super(name, count, type, description, SupportedHeaderLineType.INFO);
}
public VCFInfoHeaderLine(String name, VCFHeaderLineCount count, VCFHeaderLineType type, String description) {
super(name, count, type, description, SupportedHeaderLineType.INFO);
}
protected VCFInfoHeaderLine(String line, VCFHeaderVersion version) {
super(line, version, SupportedHeaderLineType.INFO);
}

View File

@ -0,0 +1,81 @@
package org.broadinstitute.sting.utils.codecs.vcf;
import java.util.Arrays;
import java.util.LinkedHashMap;
import java.util.Map;
/**
* @author ebanks
* A class representing a key=value entry for simple VCF header types
*/
public abstract class VCFSimpleHeaderLine extends VCFHeaderLine implements VCFNamedHeaderLine {
public enum SupportedHeaderLineType {
FILTER, ALT;
}
private String name;
private String description;
// our type of line, i.e. filter, alt, etc
private final SupportedHeaderLineType lineType;
/**
* create a VCF filter header line
*
* @param name the name for this header line
* @param description the description for this header line
* @param lineType the header line type
*/
public VCFSimpleHeaderLine(String name, String description, SupportedHeaderLineType lineType) {
super(lineType.toString(), "");
this.lineType = lineType;
this.name = name;
this.description = description;
if ( name == null || description == null )
throw new IllegalArgumentException(String.format("Invalid VCFSimpleHeaderLine: key=%s name=%s desc=%s", super.getKey(), name, description ));
}
/**
* create a VCF info header line
*
* @param line the header line
* @param version the vcf header version
* @param lineType the header line type
*/
protected VCFSimpleHeaderLine(String line, VCFHeaderVersion version, SupportedHeaderLineType lineType) {
super(lineType.toString(), "");
this.lineType = lineType;
Map<String,String> mapping = VCFHeaderLineTranslator.parseLine(version,line, Arrays.asList("ID","Description"));
name = mapping.get("ID");
description = mapping.get("Description");
if ( description == null && ALLOW_UNBOUND_DESCRIPTIONS ) // handle the case where there's no description provided
description = UNBOUND_DESCRIPTION;
}
protected String toStringEncoding() {
Map<String,Object> map = new LinkedHashMap<String,Object>();
map.put("ID", name);
map.put("Description", description);
return lineType.toString() + "=" + VCFHeaderLine.toStringEncoding(map);
}
public boolean equals(Object o) {
if ( !(o instanceof VCFSimpleHeaderLine) )
return false;
VCFSimpleHeaderLine other = (VCFSimpleHeaderLine)o;
return name.equals(other.name) &&
description.equals(other.description);
}
public String getName() {
return name;
}
public String getDescription() {
return description;
}
}

View File

@ -180,19 +180,4 @@ public class VCFUtils {
return new HashSet<VCFHeaderLine>(map.values());
}
/**
* return a set of supported format lines; what we currently support for output in the genotype fields of a VCF
* @return a set of VCF format lines
*/
public static Set<VCFFormatHeaderLine> getSupportedHeaderStrings() {
Set<VCFFormatHeaderLine> result = new HashSet<VCFFormatHeaderLine>();
result.add(new VCFFormatHeaderLine(VCFConstants.GENOTYPE_KEY, 1, VCFHeaderLineType.String, "Genotype"));
result.add(new VCFFormatHeaderLine(VCFConstants.GENOTYPE_QUALITY_KEY, 1, VCFHeaderLineType.Float, "Genotype Quality"));
result.add(new VCFFormatHeaderLine(VCFConstants.DEPTH_KEY, 1, VCFHeaderLineType.Integer, "Read Depth (only filtered reads used for calling)"));
result.add(new VCFFormatHeaderLine(VCFConstants.PHRED_GENOTYPE_LIKELIHOODS_KEY, -1, VCFHeaderLineType.Float, "Normalized, Phred-scaled likelihoods for AA,AB,BB genotypes where A=ref and B=alt; if site is not biallelic, number of likelihoods if n*(n+1)/2"));
return result;
}
}

View File

@ -133,8 +133,12 @@ public class Haplotype {
byte[] basesBeforeVariant = Arrays.copyOfRange(refBases,startIdxInReference,startIdxInReference+numPrefBases);
int startAfter = startIdxInReference+numPrefBases+ refAllele.getBases().length;
// protect against long events that overrun available reference context
if (startAfter > refBases.length)
startAfter = refBases.length;
byte[] basesAfterVariant = Arrays.copyOfRange(refBases,
startIdxInReference+numPrefBases+ refAllele.getBases().length, refBases.length);
startAfter, refBases.length);
// Create location for all haplotypes

View File

@ -108,7 +108,7 @@ public class Allele implements Comparable<Allele> {
this.bases = bases;
if ( ! acceptableAlleleBases(bases) )
throw new IllegalArgumentException("Unexpected base in allele bases " + new String(bases));
throw new IllegalArgumentException("Unexpected base in allele bases \'" + new String(bases)+"\'");
}
private Allele(String bases, boolean isRef) {

View File

@ -1343,6 +1343,15 @@ public class VariantContext implements Feature { // to enable tribble intergrati
return (int)stop;
}
private boolean hasSymbolicAlleles() {
for (Allele a: getAlleles()) {
if (a.isSymbolic()) {
return true;
}
}
return false;
}
public static VariantContext createVariantContextWithPaddedAlleles(VariantContext inputVC, byte inputRefBase, boolean refBaseShouldBeAppliedToEndOfAlleles) {
Allele refAllele = inputVC.getReference();
@ -1352,7 +1361,9 @@ public class VariantContext implements Feature { // to enable tribble intergrati
// We need to pad a VC with a common base if the length of the reference allele is less than the length of the VariantContext.
// This happens because the position of e.g. an indel is always one before the actual event (as per VCF convention).
long locLength = (inputVC.getEnd() - inputVC.getStart()) + 1;
if (refAllele.length() == locLength)
if (inputVC.hasSymbolicAlleles())
padVC = true;
else if (refAllele.length() == locLength)
padVC = false;
else if (refAllele.length() == locLength-1)
padVC = true;

View File

@ -4,6 +4,7 @@ import org.apache.commons.io.FileUtils;
import org.apache.log4j.*;
import org.apache.log4j.spi.LoggingEvent;
import org.broadinstitute.sting.commandline.CommandLineUtils;
import org.broadinstitute.sting.gatk.walkers.diffengine.DiffEngine;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import org.testng.Assert;
@ -12,6 +13,10 @@ import java.io.*;
import java.math.BigInteger;
import java.security.MessageDigest;
import java.security.NoSuchAlgorithmException;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
/**
*
@ -107,6 +112,57 @@ public abstract class BaseTest {
}
}
/**
* Simple generic utility class to creating TestNG data providers:
*
* 1: inherit this class, as in
*
* private class SummarizeDifferenceTest extends TestDataProvider {
* public SummarizeDifferenceTest() {
* super(SummarizeDifferenceTest.class);
* }
* ...
* }
*
* Provide a reference to your class to the TestDataProvider constructor.
*
* 2: Create instances of your subclass. Return from it the call to getTests, providing
* the class type of your test
*
* @DataProvider(name = "summaries")
* public Object[][] createSummaries() {
* new SummarizeDifferenceTest().addDiff("A", "A").addSummary("A:2");
* new SummarizeDifferenceTest().addDiff("A", "B").addSummary("A:1", "B:1");
* return SummarizeDifferenceTest.getTests(SummarizeDifferenceTest.class);
* }
*
* This class magically tracks created objects of this
*/
public static class TestDataProvider {
private static final Map<Class, List<Object>> tests = new HashMap<Class, List<Object>>();
/**
* Create a new TestDataProvider instance bound to the class variable C
* @param c
*/
public TestDataProvider(Class c) {
if ( ! tests.containsKey(c) )
tests.put(c, new ArrayList<Object>());
tests.get(c).add(this);
}
/**
* Return all of the data providers in the form expected by TestNG of type class C
* @param c
* @return
*/
public static Object[][] getTests(Class c) {
List<Object[]> params2 = new ArrayList<Object[]>();
for ( Object x : tests.get(c) ) params2.add(new Object[]{x});
return params2.toArray(new Object[][]{});
}
}
/**
* test if the file exists
*
@ -279,11 +335,14 @@ public abstract class BaseTest {
if (parameterize || expectedMD5.equals("")) {
// Don't assert
} else {
Assert.assertEquals(filemd5sum, expectedMD5, name + " Mismatching MD5s");
} else if ( filemd5sum.equals(expectedMD5) ) {
System.out.println(String.format(" => %s PASSED", name));
} else {
Assert.fail(String.format("%s has mismatching MD5s: expected=%s observed=%s", name, expectedMD5, filemd5sum));
}
return filemd5sum;
}
@ -326,7 +385,12 @@ public abstract class BaseTest {
System.out.printf("##### Path to calculated file (MD5=%s): %s%n", filemd5sum, pathToFileMD5File);
System.out.printf("##### Diff command: diff %s %s%n", pathToExpectedMD5File, pathToFileMD5File);
// todo -- add support for simple inline display of the first N differences for text file
// inline differences
DiffEngine.SummaryReportParams params = new DiffEngine.SummaryReportParams(System.out, 20, 10, 0);
boolean success = DiffEngine.simpleDiffFiles(new File(pathToExpectedMD5File), new File(pathToFileMD5File), params);
if ( success )
System.out.printf("Note that the above list is not comprehensive. At most 20 lines of output, and 10 specific differences will be listed. Please use -T DiffObjects -R public/testdata/exampleFASTA.fasta -m %s -t %s to explore the differences more freely%n",
pathToExpectedMD5File, pathToFileMD5File);
}
}

View File

@ -34,76 +34,76 @@ import java.util.Arrays;
* Tests CombineVariants
*/
public class CombineVariantsIntegrationTest extends WalkerTest {
// public static String baseTestString(String args) {
// return "-T CombineVariants -NO_HEADER -L 1:1-50,000,000 -o %s -R " + b36KGReference + args;
// }
//
// public void test1InOut(String file, String md5, boolean vcf3) {
// test1InOut(file, md5, "", vcf3);
// }
//
// public void test1InOut(String file, String md5, String args, boolean vcf3) {
// WalkerTestSpec spec = new WalkerTestSpec(
// baseTestString(" -priority v1 -B:v1,VCF" + (vcf3 ? "3 " : " ") + validationDataLocation + file + args),
// 1,
// Arrays.asList(md5));
// executeTest("testInOut1--" + file, spec);
// }
//
// public void combine2(String file1, String file2, String args, String md5, boolean vcf3) {
// WalkerTestSpec spec = new WalkerTestSpec(
// baseTestString(" -priority v1,v2 -B:v1,VCF" + (vcf3 ? "3 " : " ") + validationDataLocation + file1 + " -B:v2,VCF" + (vcf3 ? "3 " : " ") + validationDataLocation + file2 + args),
// 1,
// Arrays.asList(md5));
// executeTest("combine2 1:" + new File(file1).getName() + " 2:" + new File(file2).getName(), spec);
// }
//
// public void combineSites(String args, String md5) {
// String file1 = "1000G_omni2.5.b37.sites.vcf";
// String file2 = "hapmap_3.3.b37.sites.vcf";
// WalkerTestSpec spec = new WalkerTestSpec(
// "-T CombineVariants -NO_HEADER -o %s -R " + b37KGReference
// + " -L 1:1-10,000,000 -B:omni,VCF " + validationDataLocation + file1
// + " -B:hm3,VCF " + validationDataLocation + file2 + args,
// 1,
// Arrays.asList(md5));
// executeTest("combineSites 1:" + new File(file1).getName() + " 2:" + new File(file2).getName() + " args = " + args, spec);
// }
//
//
// @Test public void test1SNP() { test1InOut("pilot2.snps.vcf4.genotypes.vcf", "2117fff6e0d182cd20be508e9661829c", true); }
// @Test public void test2SNP() { test1InOut("pilot2.snps.vcf4.genotypes.vcf", "2cfaf7af3dd119df08b8a9c1f72e2f93", " -setKey foo", true); }
// @Test public void test3SNP() { test1InOut("pilot2.snps.vcf4.genotypes.vcf", "1474ac0fde2ce42a3c24f1c97eab333e", " -setKey null", true); }
// @Test public void testOfficialCEUPilotCalls() { test1InOut("CEU.trio.2010_03.genotypes.vcf.gz", "7fc66df048a0ab08cf507906e1d4a308", false); } // official project VCF files in tabix format
//
// @Test public void test1Indel1() { test1InOut("CEU.dindel.vcf4.trio.2010_06.indel.genotypes.vcf", "ec9715f53dbf4531570557c212822f12", false); }
// @Test public void test1Indel2() { test1InOut("CEU.dindel.vcf4.low_coverage.2010_06.indel.genotypes.vcf", "f1072be5f5c6ee810276d9ca6537224d", false); }
//
// @Test public void combineTrioCalls() { combine2("CEU.trio.2010_03.genotypes.vcf.gz", "YRI.trio.2010_03.genotypes.vcf.gz", "", "b77a1eec725201d9d8e74ee0c45638d3", false); } // official project VCF files in tabix format
// @Test public void combineTrioCallsMin() { combine2("CEU.trio.2010_03.genotypes.vcf.gz", "YRI.trio.2010_03.genotypes.vcf.gz", " -minimalVCF", "802977fdfd2f4905b501bb06800f60af", false); } // official project VCF files in tabix format
// @Test public void combine2Indels() { combine2("CEU.dindel.vcf4.trio.2010_06.indel.genotypes.vcf", "CEU.dindel.vcf4.low_coverage.2010_06.indel.genotypes.vcf", "", "a67157287dd2b24b5cdf7ebf8fcbbe9a", false); }
//
// @Test public void combineSNPsAndIndels() { combine2("CEU.trio.2010_03.genotypes.vcf.gz", "CEU.dindel.vcf4.low_coverage.2010_06.indel.genotypes.vcf", "", "e1f4718a179f1196538a33863da04f53", false); }
//
// @Test public void uniqueSNPs() { combine2("pilot2.snps.vcf4.genotypes.vcf", "yri.trio.gatk_glftrio.intersection.annotated.filtered.chr1.vcf", "", "b3783384b7c8e877b971033e90beba48", true); }
//
// @Test public void omniHM3Union() { combineSites(" -filteredRecordsMergeType KEEP_IF_ANY_UNFILTERED", "902e541c87caa72134db6293fc46f0ad"); }
// @Test public void omniHM3Intersect() { combineSites(" -filteredRecordsMergeType KEEP_IF_ALL_UNFILTERED", "f339ad4bb5863b58b9c919ce7d040bb9"); }
//
// @Test public void threeWayWithRefs() {
// WalkerTestSpec spec = new WalkerTestSpec(
// baseTestString(" -B:NA19240_BGI,VCF "+validationDataLocation+"NA19240.BGI.RG.vcf" +
// " -B:NA19240_ILLUMINA,VCF "+validationDataLocation+"NA19240.ILLUMINA.RG.vcf" +
// " -B:NA19240_WUGSC,VCF "+validationDataLocation+"NA19240.WUGSC.RG.vcf" +
// " -B:denovoInfo,VCF "+validationDataLocation+"yri_merged_validation_data_240610.annotated.b36.vcf" +
// " -setKey centerSet" +
// " -filteredRecordsMergeType KEEP_IF_ANY_UNFILTERED" +
// " -priority NA19240_BGI,NA19240_ILLUMINA,NA19240_WUGSC,denovoInfo" +
// " -genotypeMergeOptions UNIQUIFY -L 1"),
// 1,
// Arrays.asList("a07995587b855f3214fb71940bf23c0f"));
// executeTest("threeWayWithRefs", spec);
// }
public static String baseTestString(String args) {
return "-T CombineVariants -NO_HEADER -L 1:1-50,000,000 -o %s -R " + b36KGReference + args;
}
public void test1InOut(String file, String md5, boolean vcf3) {
test1InOut(file, md5, "", vcf3);
}
public void test1InOut(String file, String md5, String args, boolean vcf3) {
WalkerTestSpec spec = new WalkerTestSpec(
baseTestString(" -priority v1 -B:v1,VCF" + (vcf3 ? "3 " : " ") + validationDataLocation + file + args),
1,
Arrays.asList(md5));
executeTest("testInOut1--" + file, spec);
}
public void combine2(String file1, String file2, String args, String md5, boolean vcf3) {
WalkerTestSpec spec = new WalkerTestSpec(
baseTestString(" -priority v1,v2 -B:v1,VCF" + (vcf3 ? "3 " : " ") + validationDataLocation + file1 + " -B:v2,VCF" + (vcf3 ? "3 " : " ") + validationDataLocation + file2 + args),
1,
Arrays.asList(md5));
executeTest("combine2 1:" + new File(file1).getName() + " 2:" + new File(file2).getName(), spec);
}
public void combineSites(String args, String md5) {
String file1 = "1000G_omni2.5.b37.sites.vcf";
String file2 = "hapmap_3.3.b37.sites.vcf";
WalkerTestSpec spec = new WalkerTestSpec(
"-T CombineVariants -NO_HEADER -o %s -R " + b37KGReference
+ " -L 1:1-10,000,000 -B:omni,VCF " + validationDataLocation + file1
+ " -B:hm3,VCF " + validationDataLocation + file2 + args,
1,
Arrays.asList(md5));
executeTest("combineSites 1:" + new File(file1).getName() + " 2:" + new File(file2).getName() + " args = " + args, spec);
}
@Test public void test1SNP() { test1InOut("pilot2.snps.vcf4.genotypes.vcf", "2117fff6e0d182cd20be508e9661829c", true); }
@Test public void test2SNP() { test1InOut("pilot2.snps.vcf4.genotypes.vcf", "2cfaf7af3dd119df08b8a9c1f72e2f93", " -setKey foo", true); }
@Test public void test3SNP() { test1InOut("pilot2.snps.vcf4.genotypes.vcf", "1474ac0fde2ce42a3c24f1c97eab333e", " -setKey null", true); }
@Test public void testOfficialCEUPilotCalls() { test1InOut("CEU.trio.2010_03.genotypes.vcf.gz", "7fc66df048a0ab08cf507906e1d4a308", false); } // official project VCF files in tabix format
@Test public void test1Indel1() { test1InOut("CEU.dindel.vcf4.trio.2010_06.indel.genotypes.vcf", "ec9715f53dbf4531570557c212822f12", false); }
@Test public void test1Indel2() { test1InOut("CEU.dindel.vcf4.low_coverage.2010_06.indel.genotypes.vcf", "f1072be5f5c6ee810276d9ca6537224d", false); }
@Test public void combineTrioCalls() { combine2("CEU.trio.2010_03.genotypes.vcf.gz", "YRI.trio.2010_03.genotypes.vcf.gz", "", "b77a1eec725201d9d8e74ee0c45638d3", false); } // official project VCF files in tabix format
@Test public void combineTrioCallsMin() { combine2("CEU.trio.2010_03.genotypes.vcf.gz", "YRI.trio.2010_03.genotypes.vcf.gz", " -minimalVCF", "802977fdfd2f4905b501bb06800f60af", false); } // official project VCF files in tabix format
@Test public void combine2Indels() { combine2("CEU.dindel.vcf4.trio.2010_06.indel.genotypes.vcf", "CEU.dindel.vcf4.low_coverage.2010_06.indel.genotypes.vcf", "", "a67157287dd2b24b5cdf7ebf8fcbbe9a", false); }
@Test public void combineSNPsAndIndels() { combine2("CEU.trio.2010_03.genotypes.vcf.gz", "CEU.dindel.vcf4.low_coverage.2010_06.indel.genotypes.vcf", "", "e1f4718a179f1196538a33863da04f53", false); }
@Test public void uniqueSNPs() { combine2("pilot2.snps.vcf4.genotypes.vcf", "yri.trio.gatk_glftrio.intersection.annotated.filtered.chr1.vcf", "", "b3783384b7c8e877b971033e90beba48", true); }
@Test public void omniHM3Union() { combineSites(" -filteredRecordsMergeType KEEP_IF_ANY_UNFILTERED", "902e541c87caa72134db6293fc46f0ad"); }
@Test public void omniHM3Intersect() { combineSites(" -filteredRecordsMergeType KEEP_IF_ALL_UNFILTERED", "f339ad4bb5863b58b9c919ce7d040bb9"); }
@Test public void threeWayWithRefs() {
WalkerTestSpec spec = new WalkerTestSpec(
baseTestString(" -B:NA19240_BGI,VCF "+validationDataLocation+"NA19240.BGI.RG.vcf" +
" -B:NA19240_ILLUMINA,VCF "+validationDataLocation+"NA19240.ILLUMINA.RG.vcf" +
" -B:NA19240_WUGSC,VCF "+validationDataLocation+"NA19240.WUGSC.RG.vcf" +
" -B:denovoInfo,VCF "+validationDataLocation+"yri_merged_validation_data_240610.annotated.b36.vcf" +
" -setKey centerSet" +
" -filteredRecordsMergeType KEEP_IF_ANY_UNFILTERED" +
" -priority NA19240_BGI,NA19240_ILLUMINA,NA19240_WUGSC,denovoInfo" +
" -genotypeMergeOptions UNIQUIFY -L 1"),
1,
Arrays.asList("a07995587b855f3214fb71940bf23c0f"));
executeTest("threeWayWithRefs", spec);
}
// complex examples with filtering, indels, and multiple alleles
@ -119,7 +119,7 @@ public class CombineVariantsIntegrationTest extends WalkerTest {
executeTest("combineComplexSites 1:" + new File(file1).getName() + " 2:" + new File(file2).getName() + " args = " + args, spec);
}
@Test public void complexTestFull() { combineComplexSites("", "64b991fd3850f83614518f7d71f0532f"); }
// @Test public void complexTestFull() { combineComplexSites("", "64b991fd3850f83614518f7d71f0532f"); }
@Test public void complexTestMinimal() { combineComplexSites(" -minimalVCF", "0db9ef50fe54b60426474273d7c7fa99"); }
@Test public void complexTestSitesOnly() { combineComplexSites(" -sites_only", "d20acb3d53ba0a02ce92d540ebeda2a9"); }
@Test public void complexTestSitesOnlyMinimal() { combineComplexSites(" -sites_only -minimalVCF", "8d1b3d120515f8b56b5a0d10bc5da713"); }

View File

@ -10,7 +10,6 @@
<class name="org.broadinstitute.sting.gatk.walkers.recalibration.QualityScoreCovariate" />
<class name="org.broadinstitute.sting.gatk.walkers.recalibration.CycleCovariate" />
<class name="org.broadinstitute.sting.gatk.walkers.recalibration.ReadGroupCovariate" />
<class name="org.broadinstitute.sting.gatk.walkers.recalibration.TileCovariate" />
</dependencies>
</executable>
<resources>

View File

@ -19,10 +19,7 @@ import org.broadinstitute.sting.gatk.phonehome.GATKRunReport
class MethodsDevelopmentCallingPipeline extends QScript {
qscript =>
@Argument(shortName="gatk", doc="gatk jar file", required=true)
var gatkJarFile: File = _
@Argument(shortName="outputDir", doc="output directory", required=true)
@Argument(shortName="outputDir", doc="output directory", required=false)
var outputDir: String = "./"
@Argument(shortName="skipCalling", doc="skip the calling part of the pipeline and only run VQSR on preset, gold standard VCF files", required=false)
@ -185,7 +182,6 @@ class MethodsDevelopmentCallingPipeline extends QScript {
trait UNIVERSAL_GATK_ARGS extends CommandLineGATK {
logging_level = "INFO";
jarFile = gatkJarFile;
memoryLimit = 4;
phone_home = if ( LOCAL_ET ) GATKRunReport.PhoneHomeOption.STANDARD else GATKRunReport.PhoneHomeOption.AWS_S3
}

View File

@ -138,30 +138,32 @@ class QGraph extends Logging {
validate()
if (running && numMissingValues == 0) {
logger.info("Generating scatter gather jobs.")
val scatterGathers = jobGraph.edgeSet.filter(edge => scatterGatherable(edge))
if (!scatterGathers.isEmpty) {
logger.info("Generating scatter gather jobs.")
var addedFunctions = List.empty[QFunction]
for (scatterGather <- scatterGathers) {
val functions = scatterGather.asInstanceOf[FunctionEdge]
.function.asInstanceOf[ScatterGatherableFunction]
.generateFunctions()
addedFunctions ++= functions
var addedFunctions = List.empty[QFunction]
for (scatterGather <- scatterGathers) {
val functions = scatterGather.asInstanceOf[FunctionEdge]
.function.asInstanceOf[ScatterGatherableFunction]
.generateFunctions()
addedFunctions ++= functions
}
logger.info("Removing original jobs.")
this.jobGraph.removeAllEdges(scatterGathers)
prune()
logger.info("Adding scatter gather jobs.")
addedFunctions.foreach(function => if (running) this.add(function))
logger.info("Regenerating graph.")
fill
val scatterGatherDotFile = if (settings.expandedDotFile != null) settings.expandedDotFile else settings.dotFile
if (scatterGatherDotFile != null)
renderToDot(scatterGatherDotFile)
validate()
}
logger.info("Removing original jobs.")
this.jobGraph.removeAllEdges(scatterGathers)
prune()
logger.info("Adding scatter gather jobs.")
addedFunctions.foreach(function => if (running) this.add(function))
logger.info("Regenerating graph.")
fill
val scatterGatherDotFile = if (settings.expandedDotFile != null) settings.expandedDotFile else settings.dotFile
if (scatterGatherDotFile != null)
renderToDot(scatterGatherDotFile)
validate()
}
}

View File

@ -286,11 +286,11 @@ object Lsf706JobRunner extends Logging {
// LSB_SHAREDIR/cluster_name/logdir/lsb.acct (man bacct)
// LSB_SHAREDIR/cluster_name/logdir/lsb.events (man bhist)
logger.debug("Job Id %s status / exitStatus / exitInfo: ??? / ??? / ???".format(runner.jobId))
val unknownStatusSeconds = (System.currentTimeMillis - runner.lastStatusUpdate)
if (unknownStatusSeconds > (unknownStatusMaxSeconds * 1000L)) {
val unknownStatusMillis = (System.currentTimeMillis - runner.lastStatusUpdate)
if (unknownStatusMillis > (unknownStatusMaxSeconds * 1000L)) {
// Unknown status has been returned for a while now.
runner.updateStatus(RunnerStatus.FAILED)
logger.error("Unable to read LSF status for %d minutes: job id %d: %s".format(unknownStatusSeconds/60, runner.jobId, runner.function.description))
logger.error("Unable to read LSF status for %0.2f minutes: job id %d: %s".format(unknownStatusMillis/(60 * 1000D), runner.jobId, runner.function.description))
}
}

View File

@ -0,0 +1,11 @@
##fileformat=VCFv4.0
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878
chr1 2646 rs62635284 G A 0.15 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:53,75:3:-12.40,-0.90,-0.00:9.03
chr1 2979 rs62635286 T G 83.67 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:31,32:9:-33.61,-2.71,-0.00:27.09
chr1 2981 rs62028691 A G 14.69 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:31,33:9:-32.12,-2.71,-0.00:27.08
chr1 4536 rs11582131 G C 0.18 PASS AC=1;AF=0.50;AN=2 GT:AD:DP:GL:GQ 0/1:42,33:16:-41.67,-4.82,-26.29:99
chr1 4562 rs11490464 C G 0.14 PASS AC=1;AF=0.50;AN=2 GT:AD:DP:GL:GQ 0/1:26,30:9:-19.64,-2.72,-14.87:99
chr1 4770 rs6682375 A G 0.32 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:9,111:84:-306.27,-28.58,-3.46:99
chr1 4793 rs6682385 A G 0.15 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:4,115:109:-350.74,-32.88,-0.10:99
chr1 5074 rs11586607 T G 0.01 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:29,97:39:-130.41,-11.75,-3.82:79.31
chr1 5137 rs62636497 A T 140.49 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:0,74:39:-148.99,-11.75,-0.01:99

View File

@ -0,0 +1,11 @@
##fileformat=VCFv4.0
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA12878
chr1 2646 rs62635284 G A 0.15 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:53,75:3:-12.40,-0.90,-0.00:9.03
chr1 2979 rs62635286 T G 83.67 CHANGED_FILTER AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:31,32:9:-33.61,-2.71,-0.00:27.09
chr1 2981 rs62028691 A G 14.69 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:31,33:9:-32.12,-2.71,-0.00:27.08
chr1 4536 rs11582131 G C 0.18 PASS AC=2;AF=0.50;AN=2 GT:AD:DP:GL:GQ 0/1:42,33:16:-41.67,-4.82,-26.29:99
chr1 4562 rs11490464 C G 0.14 PASS AC=1;AF=0.50;AN=2 GT:AD:DP:GL:GQ 1/1:26,30:9:-19.64,-2.72,-14.87:99
chr1 4770 rs6682375 A G 0.32 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 0/1:9,111:84:-306.27,-28.58,-3.46:99
chr1 4793 rs6682385 A G 0.15 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:4,114:109:-350.74,-32.88,-0.10:99
chr1 5074 rs11586607 T G 0.01 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:29,97:39:-130.41,-11.74,-3.82:79.31
chr1 5137 rs62636497 A T 140.49 PASS AC=2;AF=1.00;AN=2 GT:AD:DP:GL:GQ 1/1:0,74:39:-148.99,-11.75,-0.01:9

View File

@ -25,5 +25,7 @@
<module organisation="javax.activation" resolver="java.net" />
<module organisation="net.java.dev.jna" resolver="maven2-repository.dev.java.net" />
<module organisation="com.google.code.caliper" resolver="projects" />
<module organisation="net.sf.gridscheduler" resolver="projects" />
<module organisation="com.google.code.cofoja" resolver="projects" />
</modules>
</ivysettings>

View File

@ -0,0 +1,3 @@
<ivy-module version="1.0">
<info organisation="com.google.code.cofoja" module="cofoja" revision="1.0-20110609" status="integration" publication="20110609114800" />
</ivy-module>

View File

@ -0,0 +1,3 @@
<ivy-module version="1.0">
<info organisation="net.sf.gridscheduler" module="drmaa" revision="6.2u5p2" status="release" />
</ivy-module>