Merge branch 'master' of github.com:broadinstitute/gsa-unstable
Conflicts: .gitignore
This commit is contained in:
Intel Repocontact
commit
ff2a972ab5
|
|
@ -38,3 +38,5 @@ resources/
|
|||
velocity.log
|
||||
perf
|
||||
verify
|
||||
maven-metadata-local.xml
|
||||
dependency-reduced-pom.xml
|
||||
|
|
|
|||
|
|
@ -0,0 +1,173 @@
|
|||
#!/bin/sh
|
||||
|
||||
mvn_args="verify"
|
||||
mvn_properties=
|
||||
mvn_clean=
|
||||
unknown_args=
|
||||
property_regex='-D(.*)=(.*)'
|
||||
unit_test_regex='.*UnitTest'
|
||||
post_script=
|
||||
run_type="run"
|
||||
|
||||
for arg in "${@}" ; do
|
||||
if [[ "${arg}" == "dry" ]] ; then
|
||||
run_type="dry"
|
||||
|
||||
elif [[ "${arg}" == "clean" ]] ; then
|
||||
mvn_clean="clean"
|
||||
mvn_args=
|
||||
|
||||
elif [[ "${arg}" =~ ${property_regex} ]] ; then
|
||||
property_name=${BASH_REMATCH[1]}
|
||||
property_value=${BASH_REMATCH[2]}
|
||||
|
||||
if [[ "${property_name}" == "single" ]] ; then
|
||||
test_property="test"
|
||||
test_disabled="it.test"
|
||||
if [[ ! "${property_value}" =~ ${unit_test_regex} ]] ; then
|
||||
test_property="it.test"
|
||||
test_disabled="test"
|
||||
fi
|
||||
|
||||
mvn_properties="${mvn_properties} -D${test_disabled}=disabled -D${test_property}=${property_value}"
|
||||
|
||||
elif [[ "${property_name}" == "test.debug.port" ]] ; then
|
||||
mvn_properties="${mvn_properties} -Dmaven.surefire.debug=\"-Xdebug -Xrunjdwp:transport=dt_socket,server=y,suspend=y,address=${property_value}\""
|
||||
mvn_properties="${mvn_properties} -Dmaven.failsafe.debug=\"-Xdebug -Xrunjdwp:transport=dt_socket,server=y,suspend=y,address=${property_value}\""
|
||||
|
||||
elif [[ "${property_name}" == "test.default.maxmemory" ]] ; then
|
||||
mvn_properties="${mvn_properties} -Dtest.maxmemory=${property_value}"
|
||||
|
||||
else
|
||||
unknown_args="${unknown_args} \"${arg}\""
|
||||
|
||||
fi
|
||||
|
||||
else
|
||||
if [[ "${arg}" != "dist" && "${mvn_args}" != "" && "${mvn_args}" != "verify" ]] ; then
|
||||
echo "Sorry, this script does not currently support mixing targets." >&2
|
||||
exit 1
|
||||
|
||||
elif [[ "${arg}" == "dist" ]] ; then
|
||||
mvn_args="verify"
|
||||
|
||||
elif [[ "${arg}" == "gatk" ]] ; then
|
||||
mvn_args="verify '-P!queue'"
|
||||
|
||||
elif [[ "${arg}" == "test.compile" ]] ; then
|
||||
mvn_args="test-compile"
|
||||
|
||||
elif [[ "${arg}" == "gatkdocs" ]] ; then
|
||||
local_repo="sitetemprepo"
|
||||
mvn_args="install -Dmaven.repo.local=${local_repo} -Ddisable.queue && mvn site -Dmaven.repo.local=${local_repo} -Ddisable.queue"
|
||||
|
||||
elif [[ "${arg}" == "package.gatk.full" ]] ; then
|
||||
mvn_args="package '-P!private,!queue'"
|
||||
|
||||
elif [[ "${arg}" == "package.gatk.all" ]] ; then
|
||||
mvn_args="package '-P!queue'"
|
||||
|
||||
elif [[ "${arg}" == "package.queue.full" ]] ; then
|
||||
mvn_args="package '-P!private'"
|
||||
|
||||
elif [[ "${arg}" == "package.queue.all" ]] ; then
|
||||
mvn_args="package"
|
||||
|
||||
# elif [[ "${arg}" == "release.gatk.full" ]] ; then
|
||||
# mvn_args="package '-P!private,!queue'"
|
||||
# post_script=" && private/src/main/scripts/shell/copy_release.sh public/gatk-package/target/GenomeAnalysisTK-*.tar.bz2"
|
||||
|
||||
# elif [[ "${arg}" == "release.queue.full" ]] ; then
|
||||
# mvn_args="package '-P!private'"
|
||||
# post_script=" && private/src/main/scripts/shell/copy_release.sh public/queue-package/target/Queue-*.tar.bz2"
|
||||
|
||||
elif [[ "${arg}" == "build-picard-private" ]] ; then
|
||||
mvn_args="mvn install -f private/picard-maven/pom.xml"
|
||||
|
||||
# TODO: clover support
|
||||
# see ant and maven docs for clover:
|
||||
# https://confluence.atlassian.com/display/CLOVER/1.+QuickStart+Guide
|
||||
# https://confluence.atlassian.com/display/CLOVER/Clover-for-Maven+2+and+3+User%27s+Guide
|
||||
#
|
||||
#elif [[ "${arg}" == "clover.report" ]] ; then
|
||||
# mvn_args=...
|
||||
#
|
||||
#elif [[ "${arg}" == "with.clover" ]] ; then
|
||||
# mvn_args=...
|
||||
|
||||
# TODO: This runs *all* commit tests, including the few on Queue.
|
||||
elif [[ "${arg}" == "gatkfull.binary.release.tests" ]] ; then
|
||||
local_repo="sitetemprepo"
|
||||
mvn_args="install -Dmaven.repo.local=${local_repo} && mvn verify"
|
||||
mvn_args="${mvn_args} -Dmaven.repo.local=${local_repo}"
|
||||
mvn_args="${mvn_args} -Dsting.packagetests.enabled=true"
|
||||
mvn_args="${mvn_args} -Dsting.packagecommittests.skipped=false"
|
||||
|
||||
# TODO: This runs only the pipeline tests (full, non-dry run), but not the commit tests for Queue.
|
||||
elif [[ "${arg}" == "queuefull.binary.release.tests" ]] ; then
|
||||
local_repo="sitetemprepo"
|
||||
mvn_args="install -Dmaven.repo.local=${local_repo} && mvn verify"
|
||||
mvn_args="${mvn_args} -Dmaven.repo.local=${local_repo}"
|
||||
mvn_args="${mvn_args} -Dsting.packagetests.enabled=true"
|
||||
mvn_args="${mvn_args} -Dsting.packagepipelinetests.skipped=false"
|
||||
mvn_args="${mvn_args} -Dsting.pipelinetests.run=true"
|
||||
|
||||
elif [[ "${arg}" == "committests" ]] ; then
|
||||
mvn_args="verify -Dsting.committests.skipped=false"
|
||||
|
||||
elif [[ "${arg}" == "test" ]] ; then
|
||||
mvn_args="test -Dsting.unittests.skipped=false"
|
||||
|
||||
elif [[ "${arg}" == "unittest" ]] ; then
|
||||
mvn_args="test -Dsting.unittests.skipped=false"
|
||||
|
||||
elif [[ "${arg}" == "integrationtest" ]] ; then
|
||||
mvn_args="verify -Dsting.integrationtests.skipped=false"
|
||||
|
||||
elif [[ "${arg}" == "largescaletest" ]] ; then
|
||||
mvn_args="verify -Dsting.largescaletests.skipped=false"
|
||||
|
||||
elif [[ "${arg}" == "knowledgebasetest" ]] ; then
|
||||
mvn_args="verify -Dsting.knowledgebasetests.skipped=false"
|
||||
|
||||
elif [[ "${arg}" == "pipelinetest" ]] ; then
|
||||
mvn_args="verify -Dsting.pipelinetests.skipped=false"
|
||||
|
||||
elif [[ "${arg}" == "pipelinetestrun" ]] ; then
|
||||
mvn_args="verify -Dsting.pipelinetests.skipped=false -Dsting.pipelinetests.run=true"
|
||||
|
||||
elif [[ "${arg}" == "fasttest" ]] ; then
|
||||
mvn_args="verify -Dsting.committests.skipped=false -pl private/gatk-private -am -Dresource.bundle.skip=true"
|
||||
|
||||
else
|
||||
unknown_args="${unknown_args} \"${arg}\""
|
||||
|
||||
fi
|
||||
|
||||
fi
|
||||
|
||||
done
|
||||
|
||||
mvn_cmd=
|
||||
if [[ "${mvn_clean}" != "" ]] ; then
|
||||
if [[ "${mvn_args}" != "" ]] ; then
|
||||
mvn_cmd="mvn ${mvn_clean} && mvn ${mvn_args}"
|
||||
else
|
||||
mvn_cmd="mvn ${mvn_clean}"
|
||||
fi
|
||||
else
|
||||
mvn_cmd="mvn ${mvn_args}"
|
||||
fi
|
||||
|
||||
if [[ "${unknown_args}" != "" ]] ; then
|
||||
echo "Unrecognized arguments:${unknown_args}" >&2
|
||||
|
||||
else
|
||||
echo "Equivalent maven command"
|
||||
echo "${mvn_cmd}${mvn_properties}${post_script}"
|
||||
|
||||
if [[ "${run_type}" != "dry" ]] ; then
|
||||
sh -c "${mvn_cmd}${mvn_properties}${post_script}"
|
||||
fi
|
||||
|
||||
fi
|
||||
117
ivy.xml
117
ivy.xml
|
|
@ -1,117 +0,0 @@
|
|||
<!--
|
||||
~ Copyright (c) 2012, 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.
|
||||
-->
|
||||
|
||||
<ivy-module version="1.0">
|
||||
<info organisation="org.broadinstitute" module="Sting"/>
|
||||
<configurations>
|
||||
<conf name="default" description="the core dependencies for the GATK"/>
|
||||
</configurations>
|
||||
<dependencies defaultconf="default">
|
||||
<dependency org="net.sf" name="sam" rev="latest.integration"/>
|
||||
<dependency org="net.sf" name="picard" rev="latest.integration"/>
|
||||
<dependency org="edu.mit.broad" name="picard-private-parts" rev="latest.integration"/>
|
||||
|
||||
<!-- Tribble -->
|
||||
<dependency org="org.broad" name="tribble" rev="latest.integration"/>
|
||||
|
||||
<!-- Variant -->
|
||||
<dependency org="org.broadinstitute" name="variant" rev="latest.integration"/>
|
||||
|
||||
<dependency org="log4j" name="log4j" rev="1.2.15"/>
|
||||
<dependency org="javax.mail" name="mail" rev="1.4.4"/>
|
||||
<dependency org="colt" name="colt" rev="1.2.0"/>
|
||||
<dependency org="it.unimi.dsi" name="fastutil" rev="6.5.3" />
|
||||
|
||||
<!-- <dependency org="jboss" name="javassist" rev="3.7.ga"/> -->
|
||||
<dependency org="org.simpleframework" name="simple-xml" rev="2.0.4"/>
|
||||
<dependency org="org.apache.bcel" name="bcel" rev="5.2"/>
|
||||
|
||||
<!-- Dependencies for reflections mvn repository -->
|
||||
<dependency org="org.reflections" name="reflections" rev="0.9.8"/>
|
||||
<dependency org="org.slf4j" name="slf4j-log4j12" rev="1.6.1"/>
|
||||
|
||||
<!-- Matrix package from math.nist.gov -->
|
||||
<dependency org="gov.nist" name="Jama" rev="1.0.2"/>
|
||||
|
||||
<!-- Dependencies for the graph aligner -->
|
||||
<dependency org="net.sf.jgrapht" name="jgrapht" rev="0.8.3"/>
|
||||
|
||||
<!-- Dependencies for the html walker documention -->
|
||||
<dependency org="org.freemarker" name="freemarker" rev="2.3.18"/>
|
||||
|
||||
<!-- Commons Dependencies -->
|
||||
<dependency org="org.apache.commons" name="commons-email" rev="1.2"/>
|
||||
<dependency org="org.apache.commons" name="commons-jexl" rev="2.1.1"/>
|
||||
<dependency org="commons-lang" name="commons-lang" rev="2.5"/>
|
||||
<dependency org="commons-logging" name="commons-logging" rev="1.1.1"/>
|
||||
<dependency org="commons-io" name="commons-io" rev="2.1"/>
|
||||
<dependency org="commons-collections" name="commons-collections" rev="3.2.1"/>
|
||||
<dependency org="org.apache.commons" name="commons-math" rev="2.2"/>
|
||||
|
||||
<!-- Lucene core utilities -->
|
||||
<!-- <dependency org="org.apache.lucene" name="lucene-core" rev="3.0.3"/> -->
|
||||
|
||||
<!-- Dependencies for LSF, DRMAA, and other C libraries -->
|
||||
<dependency org="net.java.dev.jna" name="jna" rev="3.2.7"/>
|
||||
|
||||
<!-- Dependencies for amazon.com S3 support -->
|
||||
<dependency org="net.java.dev.jets3t" name="jets3t" rev="0.8.1"/>
|
||||
|
||||
<!-- 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.10.2"/>
|
||||
<dependency org="org.scala-lang" name="scala-library" rev="2.10.2"/>
|
||||
|
||||
<!-- testing and evaluation dependencies -->
|
||||
<dependency org="org.testng" name="testng" rev="6.8"/>
|
||||
<dependency org="org.uncommons" name="reportng" rev="1.1.2"/>
|
||||
<dependency org="com.google.caliper" name="caliper" rev="0.5-rc1"/>
|
||||
<dependency org="com.google.inject" name="guice" rev="3.0"/>
|
||||
|
||||
<!-- Contracts for Java and dependencies -->
|
||||
<dependency org="com.google.code.cofoja" name="cofoja" rev="1.0-r139"/>
|
||||
<dependency org="asm" name="asm-all" rev="3.3.1"/>
|
||||
|
||||
<!-- POI, for reading pipeline files -->
|
||||
<dependency org="org.apache.poi" name="poi" rev="3.8-beta3"/>
|
||||
<dependency org="org.apache.poi" name="poi-ooxml" rev="3.8-beta3"/>
|
||||
|
||||
<!-- snpEff annotator for pipelines -->
|
||||
<dependency org="net.sf.snpeff" name="snpeff" rev="2.0.5"/>
|
||||
|
||||
<!-- MongoDB for the GXDB project -->
|
||||
<dependency org="org.mongodb" name="mongo-java-driver" rev="2.7.3"/>
|
||||
|
||||
<!-- GSON and HTTP for talking to the REST API on Vanilla Forums -->
|
||||
<dependency org="com.google.code.gson" name="gson" rev="2.2.2"/>
|
||||
<dependency org="org.apache.httpcomponents" name="httpclient" rev="4.1.1"/>
|
||||
|
||||
<!-- Exclude dependencies on sun libraries where the downloads aren't available but included in the jvm. -->
|
||||
<exclude org="javax.servlet"/>
|
||||
<exclude org="javax.jms"/>
|
||||
<exclude org="com.sun.*"/>
|
||||
</dependencies>
|
||||
</ivy-module>
|
||||
|
|
@ -0,0 +1,862 @@
|
|||
<?xml version="1.0" encoding="UTF-8"?>
|
||||
<project xmlns="http://maven.apache.org/POM/4.0.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">
|
||||
<modelVersion>4.0.0</modelVersion>
|
||||
|
||||
<!--
|
||||
This pom is the aggregator for all sting/gatk poms
|
||||
See also:
|
||||
http://maven.apache.org/pom.html#Inheritance_v
|
||||
http://maven.apache.org/guides/introduction/introduction-to-the-pom.html#Project_Inheritance_vs_Project_Aggregation
|
||||
http://stackoverflow.com/questions/1992213/maven-parent-pom-vs-modules-pom
|
||||
-->
|
||||
|
||||
<parent>
|
||||
<groupId>org.broadinstitute.sting</groupId>
|
||||
<artifactId>sting-root</artifactId>
|
||||
<version>2.8-SNAPSHOT</version>
|
||||
<relativePath>public/sting-root</relativePath>
|
||||
</parent>
|
||||
|
||||
<artifactId>sting-aggregator</artifactId>
|
||||
<packaging>pom</packaging>
|
||||
<name>Sting Aggregator</name>
|
||||
|
||||
<modules>
|
||||
<module>public</module>
|
||||
<!-- private & protected optionally enabled as profiles -->
|
||||
</modules>
|
||||
|
||||
<properties>
|
||||
<sting.basedir>${project.basedir}</sting.basedir>
|
||||
<resource.bundle.path>StingText.properties</resource.bundle.path>
|
||||
<resource.bundle.skip>false</resource.bundle.skip>
|
||||
<!-- TODO: Need a better a way to say "don't include hidden" by default -->
|
||||
<gatkdocs.include.hidden>-build-timestamp "${maven.build.timestamp}"</gatkdocs.include.hidden>
|
||||
|
||||
<!--
|
||||
Phases of the build that may be disabled to speed up compilation.
|
||||
-->
|
||||
<sting.jar.phase>package</sting.jar.phase>
|
||||
<sting.generate-resources.phase>generate-resources</sting.generate-resources.phase>
|
||||
<sting.process-resources.phase>process-resources</sting.process-resources.phase>
|
||||
<sting.process-test-resources.phase>process-test-resources</sting.process-test-resources.phase>
|
||||
|
||||
<!--
|
||||
Package tests ensure the consistency of the packaged / shaded jars.
|
||||
It runs the tests where the monolithic jar is the only dependency on the classpath.
|
||||
-->
|
||||
<sting.packagecommittests.skipped>true</sting.packagecommittests.skipped>
|
||||
<sting.packageunittests.skipped>${sting.packagecommittests.skipped}</sting.packageunittests.skipped>
|
||||
<sting.packageintegrationtests.skipped>${sting.packagecommittests.skipped}</sting.packageintegrationtests.skipped>
|
||||
<sting.packagepipelinetests.skipped>${sting.packagecommittests.skipped}</sting.packagepipelinetests.skipped>
|
||||
<sting.packagelargescaletests.skipped>true</sting.packagelargescaletests.skipped>
|
||||
<sting.packageknowledgebasetests.skipped>true</sting.packageknowledgebasetests.skipped>
|
||||
|
||||
<!--
|
||||
Serial tests use the test jars to run tests, such that all tests are run from a single TestNG invocation.
|
||||
This is different that the invoker, that runs the test classes from the filesystem, but pointing at the packaged JAR files.
|
||||
TODO: Currently, all tests run within each package, since packages already collect dependencies for shading an uber jar.
|
||||
TODO: Should there be another level up of tests, possibly running "all tests" via this aggregator level?
|
||||
TODO: If that require the aggregator to be dependent on the child dependencies, perhaps a better approach might be another monolithic test project.
|
||||
-->
|
||||
<sting.serialcommittests.skipped>true</sting.serialcommittests.skipped>
|
||||
<sting.serialunittests.skipped>${sting.serialcommittests.skipped}</sting.serialunittests.skipped>
|
||||
<sting.serialintegrationtests.skipped>${sting.serialcommittests.skipped}</sting.serialintegrationtests.skipped>
|
||||
<sting.serialpipelinetests.skipped>${sting.serialcommittests.skipped}</sting.serialpipelinetests.skipped>
|
||||
<sting.seriallargescaletests.skipped>true</sting.seriallargescaletests.skipped>
|
||||
<sting.serialknowledgebasetests.skipped>true</sting.serialknowledgebasetests.skipped>
|
||||
</properties>
|
||||
|
||||
<dependencies>
|
||||
<dependency>
|
||||
<groupId>com.sun</groupId>
|
||||
<artifactId>tools</artifactId>
|
||||
</dependency>
|
||||
</dependencies>
|
||||
|
||||
<build>
|
||||
<!-- Plugin configuration -->
|
||||
<pluginManagement>
|
||||
<plugins>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-clean-plugin</artifactId>
|
||||
<configuration>
|
||||
<filesets>
|
||||
<!--
|
||||
TODO: Once GATKDoclet stops hard coding paths, we remove this and leave
|
||||
TODO: it to the standard maven conventions to clean up for us
|
||||
-->
|
||||
<fileset>
|
||||
<directory>gatkdocs</directory>
|
||||
</fileset>
|
||||
<fileset>
|
||||
<directory>${basedir}</directory>
|
||||
<includes>
|
||||
<include>javadoc.sh</include>
|
||||
<include>options</include>
|
||||
<include>packages</include>
|
||||
</includes>
|
||||
</fileset>
|
||||
<!--
|
||||
Shade dumps this temp file in basedir, with no good way to reconfigure.
|
||||
https://jira.codehaus.org/browse/MSHADE-145
|
||||
-->
|
||||
<fileset>
|
||||
<directory>${basedir}</directory>
|
||||
<includes>
|
||||
<include>dependency-reduced-pom.xml</include>
|
||||
</includes>
|
||||
</fileset>
|
||||
</filesets>
|
||||
</configuration>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-dependency-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>unpack-direct-dependencies</id>
|
||||
<goals>
|
||||
<goal>unpack-dependencies</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<excludeTransitive>true</excludeTransitive>
|
||||
<outputDirectory>${project.build.outputDirectory}</outputDirectory>
|
||||
<includeTypes>jar</includeTypes>
|
||||
<excludeScope>system</excludeScope>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
<!-- TODO: Change the ResourceBundleExtractorDoclet to not require log4j.properties file -->
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-resources-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>default-resources</id>
|
||||
<goals>
|
||||
<goal>resources</goal>
|
||||
</goals>
|
||||
<phase>${sting.process-resources.phase}</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>default-testResources</id>
|
||||
<goals>
|
||||
<goal>testResources</goal>
|
||||
</goals>
|
||||
<phase>${sting.process-test-resources.phase}</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>copy-resource-bundle-log4j</id>
|
||||
<goals>
|
||||
<goal>copy-resources</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<outputDirectory>${project.reporting.outputDirectory}/apidocs</outputDirectory>
|
||||
<resources>
|
||||
<resource>
|
||||
<directory>${sting.basedir}/sting-utils/src/main/config/org/broadinstitute/sting/utils/help</directory>
|
||||
</resource>
|
||||
</resources>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-javadoc-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>extract-resource-bundle</id>
|
||||
<goals>
|
||||
<goal>javadoc</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<!-- Allow skipping for "fasttest" -->
|
||||
<skip>${resource.bundle.skip}</skip>
|
||||
<doclet>org.broadinstitute.sting.utils.help.ResourceBundleExtractorDoclet</doclet>
|
||||
<!-- Required as doclet uses reflection to access classes for documentation, instead of source java-->
|
||||
<docletPath>${project.build.outputDirectory}</docletPath>
|
||||
<docletArtifact>
|
||||
<groupId>${project.groupId}</groupId>
|
||||
<!-- TODO: THIS IS SUPPOSED TO BE STING-UTILS! -->
|
||||
<artifactId>gatk-framework</artifactId>
|
||||
<version>${project.version}</version>
|
||||
</docletArtifact>
|
||||
<maxmemory>2g</maxmemory>
|
||||
<useStandardDocletOptions>false</useStandardDocletOptions>
|
||||
<quiet>true</quiet>
|
||||
<additionalparam>-build-timestamp "${maven.build.timestamp}" -absolute-version ${build.version} -out ${project.build.outputDirectory}/${resource.bundle.path}</additionalparam>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-compiler-plugin</artifactId>
|
||||
<configuration>
|
||||
<proc>none</proc>
|
||||
<annotationProcessors>
|
||||
<annotationProcessor>com.google.java.contract.core.apt.AnnotationProcessor</annotationProcessor>
|
||||
</annotationProcessors>
|
||||
</configuration>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>default-compile</id>
|
||||
<phase>none</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>default-testCompile</id>
|
||||
<phase>none</phase>
|
||||
</execution>
|
||||
<!--
|
||||
Explicit package-info creation to match existing ant output.
|
||||
-->
|
||||
<execution>
|
||||
<id>compile-package-info</id>
|
||||
<goals>
|
||||
<goal>compile</goal>
|
||||
</goals>
|
||||
<phase>compile</phase>
|
||||
<configuration>
|
||||
<compilerArgs>
|
||||
<arg>-Xpkginfo:always</arg>
|
||||
</compilerArgs>
|
||||
<includes>
|
||||
<include>**/package-info.java</include>
|
||||
</includes>
|
||||
</configuration>
|
||||
</execution>
|
||||
<!--
|
||||
TODO: Currently disabled in build.xml. Here as well.
|
||||
<execution>
|
||||
<id>compile-annotations</id>
|
||||
<phase>compile</phase>
|
||||
<goals>
|
||||
<goal>compile</goal>
|
||||
</goals>
|
||||
<configuration>
|
||||
<proc>only</proc>
|
||||
</configuration>
|
||||
</execution>
|
||||
-->
|
||||
<execution>
|
||||
<id>compile-java</id>
|
||||
<goals>
|
||||
<goal>compile</goal>
|
||||
</goals>
|
||||
<phase>compile</phase>
|
||||
<configuration>
|
||||
<!--
|
||||
Package info is supposed to be in source:
|
||||
http://maven.apache.org/plugins/maven-javadoc-plugin/examples/javadoc-resources.html
|
||||
But maven-compile-plugin doesn't auto exclude these:
|
||||
https://jira.codehaus.org/browse/MCOMPILER-205?page=com.atlassian.jira.plugin.system.issuetabpanels:comment-tabpanel&focusedCommentId=326505#comment-326505
|
||||
So, explicitly exclude them
|
||||
-->
|
||||
<excludes>
|
||||
<exclude>**/package-info.java</exclude>
|
||||
</excludes>
|
||||
</configuration>
|
||||
</execution>
|
||||
<!--
|
||||
TODO: Currently disabled in build.xml. Here as well.
|
||||
<execution>
|
||||
<id>testCompile-annotations</id>
|
||||
<phase>test-compile</phase>
|
||||
<goals>
|
||||
<goal>testCompile</goal>
|
||||
</goals>
|
||||
<configuration>
|
||||
<proc>only</proc>
|
||||
</configuration>
|
||||
</execution>
|
||||
-->
|
||||
<execution>
|
||||
<id>testCompile-java</id>
|
||||
<goals>
|
||||
<goal>testCompile</goal>
|
||||
</goals>
|
||||
<phase>test-compile</phase>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
|
||||
<plugin>
|
||||
<groupId>org.scala-tools</groupId>
|
||||
<artifactId>maven-scala-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<goals>
|
||||
<goal>compile</goal>
|
||||
<goal>testCompile</goal>
|
||||
</goals>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-jar-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>default-jar</id>
|
||||
<phase>${sting.jar.phase}</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>test-jar</id>
|
||||
<goals>
|
||||
<goal>test-jar</goal>
|
||||
</goals>
|
||||
<phase>${sting.jar.phase}</phase>
|
||||
<configuration>
|
||||
<skipIfEmpty>true</skipIfEmpty>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-shade-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>sting-executable</id>
|
||||
<goals>
|
||||
<goal>shade</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<minimizeJar>true</minimizeJar>
|
||||
<artifactSet>
|
||||
<excludes>
|
||||
<exclude>org.broadinstitute.sting:gsalib:tar.gz:*</exclude>
|
||||
<exclude>org.broadinstitute.sting:*:tar.bz2:example-resources</exclude>
|
||||
</excludes>
|
||||
</artifactSet>
|
||||
<transformers>
|
||||
<transformer implementation="org.apache.maven.plugins.shade.resource.ManifestResourceTransformer">
|
||||
<manifestEntries>
|
||||
<Main-Class>${app.main.class}</Main-Class>
|
||||
</manifestEntries>
|
||||
</transformer>
|
||||
<transformer implementation="org.apache.maven.plugins.shade.resource.AppendingTransformer">
|
||||
<resource>${resource.bundle.path}</resource>
|
||||
</transformer>
|
||||
</transformers>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-assembly-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>example-resources</id>
|
||||
<goals>
|
||||
<goal>single</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<descriptors>
|
||||
<descriptor>src/main/assembly/example-resources.xml</descriptor>
|
||||
</descriptors>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>binary-dist</id>
|
||||
<goals>
|
||||
<goal>single</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<descriptors>
|
||||
<descriptor>src/main/assembly/binary-dist.xml</descriptor>
|
||||
</descriptors>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
|
||||
<!-- TODO: Remove temporary symbolic link creation, after fixing test paths to use local resources. -->
|
||||
<plugin>
|
||||
<groupId>com.pyx4j</groupId>
|
||||
<artifactId>maven-junction-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>link-public-testdata</id>
|
||||
<goals>
|
||||
<goal>link</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${basedir}/public/testdata</dst>
|
||||
<src>${sting.basedir}/public/gatk-framework/src/test/resources</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-public-testdata</id>
|
||||
<goals>
|
||||
<goal>unlink</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${basedir}/public/testdata</dst>
|
||||
<src>${sting.basedir}/public/gatk-framework/src/test/resources</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>link-private-testdata</id>
|
||||
<goals>
|
||||
<goal>link</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${basedir}/private/testdata</dst>
|
||||
<src>${sting.basedir}/private/gatk-private/src/test/resources</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-private-testdata</id>
|
||||
<goals>
|
||||
<goal>unlink</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${basedir}/private/testdata</dst>
|
||||
<src>${sting.basedir}/private/gatk-private/src/test/resources</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>link-public-qscript</id>
|
||||
<goals>
|
||||
<goal>link</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${basedir}/public/scala/qscript</dst>
|
||||
<src>${sting.basedir}/public/queue-framework/src/main/qscripts</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-public-qscript</id>
|
||||
<goals>
|
||||
<goal>unlink</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${basedir}/public/scala/qscript</dst>
|
||||
<src>${sting.basedir}/public/queue-framework/src/main/qscripts</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>link-private-qscript</id>
|
||||
<goals>
|
||||
<goal>link</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${basedir}/private/scala/qscript</dst>
|
||||
<src>${sting.basedir}/private/queue-private/src/main/qscripts</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-private-qscript</id>
|
||||
<goals>
|
||||
<goal>unlink</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${basedir}/private/scala/qscript</dst>
|
||||
<src>${sting.basedir}/private/queue-private/src/main/qscripts</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>link-binary-jar</id>
|
||||
<goals>
|
||||
<goal>link</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${sting.basedir}/target/${sting.binary-dist.name}.${project.packaging}</dst>
|
||||
<src>${project.build.directory}/${project.build.finalName}.${project.packaging}</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>link-git-release</id>
|
||||
<goals>
|
||||
<goal>link</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<links>
|
||||
<link>
|
||||
<dst>${project.build.directory}/${sting.binary-dist.name}-${build.version}.tar.bz2</dst>
|
||||
<src>${project.build.directory}/${project.build.finalName}-binary-dist.tar.bz2</src>
|
||||
</link>
|
||||
</links>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-invoker-plugin</artifactId>
|
||||
<configuration>
|
||||
<skipInvocation>true</skipInvocation>
|
||||
<debug>false</debug>
|
||||
<pom>${sting.basedir}/public/package-tests/pom.xml</pom>
|
||||
<noLog>true</noLog>
|
||||
<streamLogs>true</streamLogs>
|
||||
<localRepositoryPath>${sting.basedir}/${maven.repo.local}</localRepositoryPath>
|
||||
<properties>
|
||||
<test>${test}</test>
|
||||
<it.test>${it.test}</it.test>
|
||||
<skip>false</skip>
|
||||
<maven.test.skip>false</maven.test.skip>
|
||||
<sting.packagetests.artifactId>${sting.packagetests.artifactId}</sting.packagetests.artifactId>
|
||||
<sting.packagetests.testClasses>${project.build.testOutputDirectory}</sting.packagetests.testClasses>
|
||||
<sting.packagetests.basedir>${project.basedir}</sting.packagetests.basedir>
|
||||
<sting.pipelinetests.run>${sting.pipelinetests.run}</sting.pipelinetests.run>
|
||||
<maven.surefire.debug>${maven.surefire.debug}</maven.surefire.debug>
|
||||
<maven.failsafe.debug>${maven.failsafe.debug}</maven.failsafe.debug>
|
||||
</properties>
|
||||
<!--
|
||||
To allow using separated integration-test and verify, each execution must use a separate reportsDirectory.
|
||||
|
||||
Otherwise, when an empty "pipeline test" later runs on an IntegrationTest class,
|
||||
it overwrites the results of the previous invocation always with a zero exit status.
|
||||
|
||||
Mixing in the test name into the reportsDirectory also avoids collisions, when different maven jobs run tests in parallel.
|
||||
|
||||
Similarly generating unique failsafe summary reports to avoid collisions.
|
||||
-->
|
||||
</configuration>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>package-unittests</id>
|
||||
<goals>
|
||||
<goal>run</goal>
|
||||
</goals>
|
||||
<configuration>
|
||||
<goals>
|
||||
<goal>test</goal>
|
||||
</goals>
|
||||
<reportsDirectory>${project.build.directory}/invoker-reports/unit/${test}</reportsDirectory>
|
||||
<skipInvocation>${sting.packageunittests.skipped}</skipInvocation>
|
||||
<properties>
|
||||
<unittests.profile.enabled>true</unittests.profile.enabled>
|
||||
<sting.packageunittests.skipped>${sting.packageunittests.skipped}</sting.packageunittests.skipped>
|
||||
</properties>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>package-integrationtests</id>
|
||||
<goals>
|
||||
<goal>integration-test</goal>
|
||||
<goal>verify</goal>
|
||||
</goals>
|
||||
<configuration>
|
||||
<goals>
|
||||
<goal>verify</goal>
|
||||
</goals>
|
||||
<reportsDirectory>${project.build.directory}/invoker-reports/integration/${it.test}</reportsDirectory>
|
||||
<skipInvocation>${sting.packageintegrationtests.skipped}</skipInvocation>
|
||||
<properties>
|
||||
<integrationtests.profile.enabled>true</integrationtests.profile.enabled>
|
||||
<sting.packageintegrationtests.skipped>${sting.packageintegrationtests.skipped}</sting.packageintegrationtests.skipped>
|
||||
<failsafe.summaryFile>${project.build.directory}/failsafe-reports/integration/failsafe-summary-${it.test}.xml</failsafe.summaryFile>
|
||||
</properties>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>package-pipelinetests</id>
|
||||
<goals>
|
||||
<goal>integration-test</goal>
|
||||
<goal>verify</goal>
|
||||
</goals>
|
||||
<configuration>
|
||||
<goals>
|
||||
<goal>verify</goal>
|
||||
</goals>
|
||||
<reportsDirectory>${project.build.directory}/invoker-reports/pipeline/${it.test}</reportsDirectory>
|
||||
<skipInvocation>${sting.packagepipelinetests.skipped}</skipInvocation>
|
||||
<properties>
|
||||
<integrationtests.profile.enabled>true</integrationtests.profile.enabled>
|
||||
<sting.packagepipelinetests.skipped>${sting.packagepipelinetests.skipped}</sting.packagepipelinetests.skipped>
|
||||
<failsafe.summaryFile>${project.build.directory}/failsafe-reports/pipeline/failsafe-summary-${it.test}.xml</failsafe.summaryFile>
|
||||
</properties>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>package-largescaletests</id>
|
||||
<goals>
|
||||
<goal>integration-test</goal>
|
||||
<goal>verify</goal>
|
||||
</goals>
|
||||
<configuration>
|
||||
<goals>
|
||||
<goal>verify</goal>
|
||||
</goals>
|
||||
<reportsDirectory>${project.build.directory}/invoker-reports/largescale/${it.test}</reportsDirectory>
|
||||
<skipInvocation>${sting.packagelargescaletests.skipped}</skipInvocation>
|
||||
<properties>
|
||||
<integrationtests.profile.enabled>true</integrationtests.profile.enabled>
|
||||
<sting.packagelargescaletests.skipped>${sting.packagelargescaletests.skipped}</sting.packagelargescaletests.skipped>
|
||||
<failsafe.summaryFile>${project.build.directory}/failsafe-reports/largescale/failsafe-summary-${it.test}.xml</failsafe.summaryFile>
|
||||
</properties>
|
||||
</configuration>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>package-knowledgebasetests</id>
|
||||
<goals>
|
||||
<goal>integration-test</goal>
|
||||
<goal>verify</goal>
|
||||
</goals>
|
||||
<configuration>
|
||||
<goals>
|
||||
<goal>verify</goal>
|
||||
</goals>
|
||||
<reportsDirectory>${project.build.directory}/invoker-reports/knowledgebase/${it.test}</reportsDirectory>
|
||||
<skipInvocation>${sting.packageknowledgebasetests.skipped}</skipInvocation>
|
||||
<properties>
|
||||
<integrationtests.profile.enabled>true</integrationtests.profile.enabled>
|
||||
<sting.packageknowledgebasetests.skipped>${sting.packageknowledgebasetests.skipped}</sting.packageknowledgebasetests.skipped>
|
||||
<failsafe.summaryFile>${project.build.directory}/failsafe-reports/knowledgebase/failsafe-summary-${it.test}.xml</failsafe.summaryFile>
|
||||
</properties>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-install-plugin</artifactId>
|
||||
<version>2.5</version>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>install-package</id>
|
||||
<goals>
|
||||
<goal>install-file</goal>
|
||||
</goals>
|
||||
<phase>none</phase>
|
||||
<configuration>
|
||||
<generatePom>true</generatePom>
|
||||
<groupId>${project.groupId}</groupId>
|
||||
<artifactId>${project.artifactId}</artifactId>
|
||||
<version>${project.version}</version>
|
||||
<packaging>${project.packaging}</packaging>
|
||||
<file>${project.build.directory}/${project.build.finalName}.${project.packaging}</file>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
</plugins>
|
||||
</pluginManagement>
|
||||
|
||||
<!-- Invoke plugins that always run -->
|
||||
<plugins>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-failsafe-plugin</artifactId>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>com.pyx4j</groupId>
|
||||
<artifactId>maven-junction-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>link-public-testdata</id>
|
||||
<phase>process-test-resources</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-public-testdata</id>
|
||||
<phase>clean</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>link-public-qscript</id>
|
||||
<phase>process-test-resources</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-public-qscript</id>
|
||||
<phase>clean</phase>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-clean-plugin</artifactId>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>com.google.code.sortpom</groupId>
|
||||
<artifactId>maven-sortpom-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>package-tests</id>
|
||||
<goals>
|
||||
<goal>sort</goal>
|
||||
</goals>
|
||||
<phase>verify</phase>
|
||||
<inherited>false</inherited>
|
||||
<configuration>
|
||||
<pomFile>public/package-tests/pom.xml</pomFile>
|
||||
</configuration>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
</plugins>
|
||||
</build>
|
||||
|
||||
<reporting>
|
||||
<plugins>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-javadoc-plugin</artifactId>
|
||||
<version>2.9.1</version>
|
||||
<reportSets>
|
||||
<!-- By not specifying an id, shut off all default javadocs from mvn site -->
|
||||
<reportSet>
|
||||
<reports />
|
||||
</reportSet>
|
||||
<reportSet>
|
||||
<id>generate-gatk-docs</id>
|
||||
<reports>
|
||||
<report>aggregate</report>
|
||||
</reports>
|
||||
<!-- Only generate the GATK Docs across the parent aggregation, not the children too. -->
|
||||
<inherited>false</inherited>
|
||||
<configuration>
|
||||
<doclet>org.broadinstitute.sting.utils.help.GATKDoclet</doclet>
|
||||
<docletArtifact>
|
||||
<groupId>${project.groupId}</groupId>
|
||||
<artifactId>gatk-package</artifactId>
|
||||
<version>${project.version}</version>
|
||||
</docletArtifact>
|
||||
<useStandardDocletOptions>false</useStandardDocletOptions>
|
||||
<quiet>true</quiet>
|
||||
<show>private</show>
|
||||
<additionalparam>-build-timestamp "${maven.build.timestamp}" -absolute-version ${build.version} ${gatkdocs.include.hidden} -settings-dir ${sting.basedir}/settings/helpTemplates -destination-dir ${project.build.directory}/gatkdocs</additionalparam>
|
||||
</configuration>
|
||||
</reportSet>
|
||||
</reportSets>
|
||||
</plugin>
|
||||
</plugins>
|
||||
</reporting>
|
||||
|
||||
<profiles>
|
||||
<!-- Optionally include protected -->
|
||||
<profile>
|
||||
<id>protected</id>
|
||||
<activation>
|
||||
<file>
|
||||
<exists>${basedir}/protected/pom.xml</exists>
|
||||
</file>
|
||||
</activation>
|
||||
<modules>
|
||||
<module>protected</module>
|
||||
</modules>
|
||||
</profile>
|
||||
|
||||
<!-- Optionally include private -->
|
||||
<profile>
|
||||
<id>private</id>
|
||||
<activation>
|
||||
<file>
|
||||
<exists>${basedir}/private/pom.xml</exists>
|
||||
</file>
|
||||
</activation>
|
||||
<modules>
|
||||
<module>private</module>
|
||||
</modules>
|
||||
<build>
|
||||
<plugins>
|
||||
<!-- TODO: All tests require access to private. For now. -->
|
||||
<plugin>
|
||||
<groupId>com.pyx4j</groupId>
|
||||
<artifactId>maven-junction-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>link-private-testdata</id>
|
||||
<phase>process-test-resources</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-private-testdata</id>
|
||||
<phase>clean</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>link-private-qscript</id>
|
||||
<phase>process-test-resources</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-private-qscript</id>
|
||||
<phase>clean</phase>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
</plugins>
|
||||
</build>
|
||||
</profile>
|
||||
|
||||
<!-- Collection of properties for use during package testing -->
|
||||
<profile>
|
||||
<id>packagetests-enabled</id>
|
||||
<activation>
|
||||
<property>
|
||||
<name>sting.packagetests.enabled</name>
|
||||
<value>true</value>
|
||||
</property>
|
||||
</activation>
|
||||
<properties>
|
||||
<maven.javadoc.skip>true</maven.javadoc.skip>
|
||||
<sting.generate-gatk-extensions.skipped>true</sting.generate-gatk-extensions.skipped>
|
||||
<sting.jar.phase>none</sting.jar.phase>
|
||||
<sting.generate-resources.phase>none</sting.generate-resources.phase>
|
||||
<sting.process-resources.phase>none</sting.process-resources.phase>
|
||||
<sting.process-test-resources.phase>none</sting.process-test-resources.phase>
|
||||
</properties>
|
||||
</profile>
|
||||
</profiles>
|
||||
|
||||
</project>
|
||||
|
|
@ -0,0 +1,139 @@
|
|||
<?xml version="1.0" encoding="UTF-8"?>
|
||||
<project xmlns="http://maven.apache.org/POM/4.0.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/maven-v4_0_0.xsd">
|
||||
<modelVersion>4.0.0</modelVersion>
|
||||
|
||||
<parent>
|
||||
<groupId>org.broadinstitute.sting</groupId>
|
||||
<artifactId>sting-aggregator</artifactId>
|
||||
<version>2.8-SNAPSHOT</version>
|
||||
<relativePath>../..</relativePath>
|
||||
</parent>
|
||||
|
||||
<artifactId>gatk-protected</artifactId>
|
||||
<packaging>jar</packaging>
|
||||
<name>GATK Protected</name>
|
||||
|
||||
<properties>
|
||||
<sting.basedir>${project.basedir}/../..</sting.basedir>
|
||||
<sting.packagetests.artifactId>gatk-package</sting.packagetests.artifactId>
|
||||
</properties>
|
||||
|
||||
<dependencies>
|
||||
<dependency>
|
||||
<groupId>${project.groupId}</groupId>
|
||||
<artifactId>gatk-framework</artifactId>
|
||||
<version>${project.version}</version>
|
||||
</dependency>
|
||||
|
||||
<dependency>
|
||||
<groupId>net.sf.jgrapht</groupId>
|
||||
<artifactId>jgrapht</artifactId>
|
||||
</dependency>
|
||||
|
||||
<dependency>
|
||||
<groupId>gov.nist.math</groupId>
|
||||
<artifactId>jama</artifactId>
|
||||
</dependency>
|
||||
|
||||
<dependency>
|
||||
<groupId>it.unimi.dsi</groupId>
|
||||
<artifactId>fastutil</artifactId>
|
||||
</dependency>
|
||||
|
||||
<dependency>
|
||||
<groupId>${project.groupId}</groupId>
|
||||
<artifactId>gatk-framework</artifactId>
|
||||
<version>${project.version}</version>
|
||||
<type>test-jar</type>
|
||||
<scope>test</scope>
|
||||
</dependency>
|
||||
|
||||
<dependency>
|
||||
<groupId>org.testng</groupId>
|
||||
<artifactId>testng</artifactId>
|
||||
<scope>test</scope>
|
||||
</dependency>
|
||||
<dependency>
|
||||
<groupId>com.google.caliper</groupId>
|
||||
<artifactId>caliper</artifactId>
|
||||
<scope>test</scope>
|
||||
</dependency>
|
||||
</dependencies>
|
||||
|
||||
<build>
|
||||
<plugins>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-resources-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>copy-resource-bundle-log4j</id>
|
||||
<phase>prepare-package</phase>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-javadoc-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>extract-resource-bundle</id>
|
||||
<phase>prepare-package</phase>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
<plugin>
|
||||
<groupId>org.apache.maven.plugins</groupId>
|
||||
<artifactId>maven-invoker-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>package-unittests</id>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>package-integrationtests</id>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>package-largescaletests</id>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>package-knowledgebasetests</id>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>package-pipelinetests</id>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
</plugins>
|
||||
</build>
|
||||
|
||||
<profiles>
|
||||
<profile>
|
||||
<id>private</id>
|
||||
<activation>
|
||||
<file>
|
||||
<exists>${basedir}/../../private/gatk-private/pom.xml</exists>
|
||||
</file>
|
||||
</activation>
|
||||
<build>
|
||||
<plugins>
|
||||
<!-- TODO: All tests require access to private. For now. -->
|
||||
<plugin>
|
||||
<groupId>com.pyx4j</groupId>
|
||||
<artifactId>maven-junction-plugin</artifactId>
|
||||
<executions>
|
||||
<execution>
|
||||
<id>link-private-testdata</id>
|
||||
<phase>process-test-resources</phase>
|
||||
</execution>
|
||||
<execution>
|
||||
<id>unlink-private-testdata</id>
|
||||
<phase>clean</phase>
|
||||
</execution>
|
||||
</executions>
|
||||
</plugin>
|
||||
</plugins>
|
||||
</build>
|
||||
</profile>
|
||||
</profiles>
|
||||
|
||||
</project>
|
||||
|
|
@ -0,0 +1,517 @@
|
|||
/*
|
||||
* By downloading the PROGRAM you agree to the following terms of use:
|
||||
*
|
||||
* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
|
||||
*
|
||||
* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
|
||||
*
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
|
||||
* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
|
||||
*
|
||||
* 1. DEFINITIONS
|
||||
* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
|
||||
*
|
||||
* 2. LICENSE
|
||||
* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
|
||||
* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
|
||||
* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
|
||||
* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
|
||||
*
|
||||
* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
|
||||
* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
|
||||
* Copyright 2012 Broad Institute, Inc.
|
||||
* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
|
||||
* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
|
||||
*
|
||||
* 4. INDEMNIFICATION
|
||||
* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
|
||||
*
|
||||
* 5. NO REPRESENTATIONS OR WARRANTIES
|
||||
* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
|
||||
* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
|
||||
*
|
||||
* 6. ASSIGNMENT
|
||||
* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
|
||||
*
|
||||
* 7. MISCELLANEOUS
|
||||
* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
|
||||
* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
|
||||
* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
|
||||
* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
|
||||
* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
|
||||
* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
|
||||
* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
|
||||
*/
|
||||
|
||||
package org.broadinstitute.sting.gatk.walkers.annotator;
|
||||
|
||||
import cern.jet.math.Arithmetic;
|
||||
import org.apache.log4j.Logger;
|
||||
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.annotator.interfaces.ActiveRegionBasedAnnotation;
|
||||
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.AnnotatorCompatible;
|
||||
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.InfoFieldAnnotation;
|
||||
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.StandardAnnotation;
|
||||
import org.broadinstitute.sting.utils.genotyper.MostLikelyAllele;
|
||||
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
|
||||
import org.broadinstitute.sting.utils.QualityUtils;
|
||||
import org.broadinstitute.variant.variantcontext.Genotype;
|
||||
import org.broadinstitute.variant.variantcontext.GenotypesContext;
|
||||
import org.broadinstitute.variant.vcf.VCFHeaderLineType;
|
||||
import org.broadinstitute.variant.vcf.VCFInfoHeaderLine;
|
||||
import org.broadinstitute.sting.utils.pileup.PileupElement;
|
||||
import org.broadinstitute.sting.utils.sam.GATKSAMRecord;
|
||||
import org.broadinstitute.sting.utils.sam.ReadUtils;
|
||||
import org.broadinstitute.variant.variantcontext.Allele;
|
||||
import org.broadinstitute.variant.variantcontext.VariantContext;
|
||||
|
||||
import java.util.*;
|
||||
|
||||
|
||||
/**
|
||||
* Phred-scaled p-value using Fisher's Exact Test to detect strand bias
|
||||
*
|
||||
* <p>Phred-scaled p-value using Fisher's Exact Test to detect strand bias (the variation
|
||||
* being seen on only the forward or only the reverse strand) in the reads. More bias is
|
||||
* indicative of false positive calls.
|
||||
* </p>
|
||||
*
|
||||
* <h3>Caveat</h3>
|
||||
* <p>The Fisher Strand test may not be calculated for certain complex indel cases or for multi-allelic sites.</p>
|
||||
*/
|
||||
public class FisherStrand extends InfoFieldAnnotation implements StandardAnnotation, ActiveRegionBasedAnnotation {
|
||||
private final static boolean ENABLE_DEBUGGING = false;
|
||||
private final static Logger logger = Logger.getLogger(FisherStrand.class);
|
||||
|
||||
private static final String FS = "FS";
|
||||
private static final double MIN_PVALUE = 1E-320;
|
||||
private static final int MIN_QUAL_FOR_FILTERED_TEST = 17;
|
||||
private static final int MIN_COUNT = 2;
|
||||
|
||||
public Map<String, Object> annotate(final RefMetaDataTracker tracker,
|
||||
final AnnotatorCompatible walker,
|
||||
final ReferenceContext ref,
|
||||
final Map<String, AlignmentContext> stratifiedContexts,
|
||||
final VariantContext vc,
|
||||
final Map<String, PerReadAlleleLikelihoodMap> stratifiedPerReadAlleleLikelihoodMap) {
|
||||
if ( !vc.isVariant() )
|
||||
return null;
|
||||
|
||||
if ( vc.hasGenotypes() ) {
|
||||
final int[][] tableFromPerSampleAnnotations = getTableFromSamples( vc.getGenotypes() );
|
||||
if ( tableFromPerSampleAnnotations != null ) {
|
||||
return pValueForBestTable(tableFromPerSampleAnnotations, null);
|
||||
}
|
||||
}
|
||||
|
||||
if (vc.isSNP() && stratifiedContexts != null) {
|
||||
final int[][] tableNoFiltering = getSNPContingencyTable(stratifiedContexts, vc.getReference(), vc.getAltAlleleWithHighestAlleleCount(), -1);
|
||||
final int[][] tableFiltering = getSNPContingencyTable(stratifiedContexts, vc.getReference(), vc.getAltAlleleWithHighestAlleleCount(), MIN_QUAL_FOR_FILTERED_TEST);
|
||||
printTable("unfiltered", tableNoFiltering);
|
||||
printTable("filtered", tableFiltering);
|
||||
return pValueForBestTable(tableFiltering, tableNoFiltering);
|
||||
}
|
||||
else if (stratifiedPerReadAlleleLikelihoodMap != null) {
|
||||
// either SNP with no alignment context, or indels: per-read likelihood map needed
|
||||
final int[][] table = getContingencyTable(stratifiedPerReadAlleleLikelihoodMap, vc);
|
||||
// logger.info("VC " + vc);
|
||||
// printTable(table, 0.0);
|
||||
return pValueForBestTable(table, null);
|
||||
}
|
||||
else
|
||||
// for non-snp variants, we need per-read likelihoods.
|
||||
// for snps, we can get same result from simple pileup
|
||||
return null;
|
||||
}
|
||||
|
||||
/**
|
||||
* Create the FisherStrand table by retrieving the per-sample strand bias annotation and adding them together
|
||||
* @param genotypes the genotypes from which to pull out the per-sample strand bias annotation
|
||||
* @return the table used for the FisherStrand p-value calculation, will be null if none of the genotypes contain the per-sample SB annotation
|
||||
*/
|
||||
private int[][] getTableFromSamples( final GenotypesContext genotypes ) {
|
||||
if( genotypes == null ) { throw new IllegalArgumentException("Genotypes cannot be null."); }
|
||||
|
||||
final int[] sbArray = {0,0,0,0}; // reference-forward-reverse -by- alternate-forward-reverse
|
||||
boolean foundData = false;
|
||||
|
||||
for( final Genotype g : genotypes ) {
|
||||
if( g.isNoCall() || ! g.hasAnyAttribute(StrandBiasBySample.STRAND_BIAS_BY_SAMPLE_KEY_NAME) )
|
||||
continue;
|
||||
|
||||
foundData = true;
|
||||
final String sbbsString = (String) g.getAnyAttribute(StrandBiasBySample.STRAND_BIAS_BY_SAMPLE_KEY_NAME);
|
||||
final int[] data = encodeSBBS(sbbsString);
|
||||
if ( passesMinimumThreshold(data) ) {
|
||||
for( int index = 0; index < sbArray.length; index++ ) {
|
||||
sbArray[index] += data[index];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return ( foundData ? decodeSBBS(sbArray) : null );
|
||||
}
|
||||
|
||||
/**
|
||||
* Does this strand data array pass the minimum threshold for inclusion?
|
||||
*
|
||||
* @param data the array
|
||||
* @return true if it passes the minimum threshold, false otherwise
|
||||
*/
|
||||
private static boolean passesMinimumThreshold(final int[] data) {
|
||||
// the ref and alt totals must each be greater than MIN_COUNT
|
||||
return data[0] + data[1] > MIN_COUNT && data[2] + data[3] > MIN_COUNT;
|
||||
}
|
||||
|
||||
/**
|
||||
* Create an annotation for the highest (i.e., least significant) p-value of table1 and table2
|
||||
*
|
||||
* @param table1 a contingency table, may be null
|
||||
* @param table2 a contingency table, may be null
|
||||
* @return annotation result for FS given tables
|
||||
*/
|
||||
private Map<String, Object> pValueForBestTable(final int[][] table1, final int[][] table2) {
|
||||
if ( table2 == null )
|
||||
return table1 == null ? null : annotationForOneTable(pValueForContingencyTable(table1));
|
||||
else if (table1 == null)
|
||||
return annotationForOneTable(pValueForContingencyTable(table2));
|
||||
else { // take the one with the best (i.e., least significant pvalue)
|
||||
double pvalue1 = pValueForContingencyTable(table1);
|
||||
double pvalue2 = pValueForContingencyTable(table2);
|
||||
return annotationForOneTable(Math.max(pvalue1, pvalue2));
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Returns an annotation result given a pValue
|
||||
*
|
||||
* @param pValue
|
||||
* @return a hash map from FS -> phred-scaled pValue
|
||||
*/
|
||||
private Map<String, Object> annotationForOneTable(final double pValue) {
|
||||
final Object value = String.format("%.3f", QualityUtils.phredScaleErrorRate(Math.max(pValue, MIN_PVALUE))); // prevent INFINITYs
|
||||
return Collections.singletonMap(FS, value);
|
||||
}
|
||||
|
||||
public List<String> getKeyNames() {
|
||||
return Collections.singletonList(FS);
|
||||
}
|
||||
|
||||
public List<VCFInfoHeaderLine> getDescriptions() {
|
||||
return Collections.singletonList(new VCFInfoHeaderLine(FS, 1, VCFHeaderLineType.Float, "Phred-scaled p-value using Fisher's exact test to detect strand bias"));
|
||||
}
|
||||
|
||||
/**
|
||||
* Helper function to turn the FisherStrand table into the SB annotation array
|
||||
* @param table the table used by the FisherStrand annotation
|
||||
* @return the array used by the per-sample Strand Bias annotation
|
||||
*/
|
||||
public static List<Integer> getContingencyArray( final int[][] table ) {
|
||||
if(table.length != 2) { throw new IllegalArgumentException("Expecting a 2x2 strand bias table."); }
|
||||
if(table[0].length != 2) { throw new IllegalArgumentException("Expecting a 2x2 strand bias table."); }
|
||||
final List<Integer> list = new ArrayList<>(4); // TODO - if we ever want to do something clever with multi-allelic sites this will need to change
|
||||
list.add(table[0][0]);
|
||||
list.add(table[0][1]);
|
||||
list.add(table[1][0]);
|
||||
list.add(table[1][1]);
|
||||
return list;
|
||||
}
|
||||
|
||||
/**
|
||||
* Helper function to parse the genotype annotation into the SB annotation array
|
||||
* @param string the string that is returned by genotype.getAnnotation("SB")
|
||||
* @return the array used by the per-sample Strand Bias annotation
|
||||
*/
|
||||
private static int[] encodeSBBS( final String string ) {
|
||||
final int[] array = new int[4];
|
||||
final StringTokenizer tokenizer = new StringTokenizer(string, ",", false);
|
||||
for( int index = 0; index < 4; index++ ) {
|
||||
array[index] = Integer.parseInt(tokenizer.nextToken());
|
||||
}
|
||||
return array;
|
||||
}
|
||||
|
||||
/**
|
||||
* Helper function to turn the SB annotation array into the FisherStrand table
|
||||
* @param array the array used by the per-sample Strand Bias annotation
|
||||
* @return the table used by the FisherStrand annotation
|
||||
*/
|
||||
private static int[][] decodeSBBS( final int[] array ) {
|
||||
if(array.length != 4) { throw new IllegalArgumentException("Expecting a length = 4 strand bias array."); }
|
||||
final int[][] table = new int[2][2];
|
||||
table[0][0] = array[0];
|
||||
table[0][1] = array[1];
|
||||
table[1][0] = array[2];
|
||||
table[1][1] = array[3];
|
||||
return table;
|
||||
}
|
||||
|
||||
private Double pValueForContingencyTable(int[][] originalTable) {
|
||||
final int[][] normalizedTable = normalizeContingencyTable(originalTable);
|
||||
|
||||
int[][] table = copyContingencyTable(normalizedTable);
|
||||
|
||||
double pCutoff = computePValue(table);
|
||||
//printTable(table, pCutoff);
|
||||
|
||||
double pValue = pCutoff;
|
||||
while (rotateTable(table)) {
|
||||
double pValuePiece = computePValue(table);
|
||||
|
||||
//printTable(table, pValuePiece);
|
||||
|
||||
if (pValuePiece <= pCutoff) {
|
||||
pValue += pValuePiece;
|
||||
}
|
||||
}
|
||||
|
||||
table = copyContingencyTable(normalizedTable);
|
||||
while (unrotateTable(table)) {
|
||||
double pValuePiece = computePValue(table);
|
||||
|
||||
//printTable(table, pValuePiece);
|
||||
|
||||
if (pValuePiece <= pCutoff) {
|
||||
pValue += pValuePiece;
|
||||
}
|
||||
}
|
||||
|
||||
//System.out.printf("P-cutoff: %f\n", pCutoff);
|
||||
//System.out.printf("P-value: %f\n\n", pValue);
|
||||
|
||||
// min is necessary as numerical precision can result in pValue being slightly greater than 1.0
|
||||
return Math.min(pValue, 1.0);
|
||||
}
|
||||
|
||||
// how large do we want the normalized table to be?
|
||||
private static final double TARGET_TABLE_SIZE = 200.0;
|
||||
|
||||
/**
|
||||
* Normalize the table so that the entries are not too large.
|
||||
* Note that this method does NOT necessarily make a copy of the table being passed in!
|
||||
*
|
||||
* @param table the original table
|
||||
* @return a normalized version of the table or the original table if it is already normalized
|
||||
*/
|
||||
private static int[][] normalizeContingencyTable(final int[][] table) {
|
||||
final int sum = table[0][0] + table[0][1] + table[1][0] + table[1][1];
|
||||
if ( sum <= TARGET_TABLE_SIZE * 2 )
|
||||
return table;
|
||||
|
||||
final double normalizationFactor = (double)sum / TARGET_TABLE_SIZE;
|
||||
|
||||
final int[][] normalized = new int[2][2];
|
||||
for ( int i = 0; i < 2; i++ ) {
|
||||
for ( int j = 0; j < 2; j++ )
|
||||
normalized[i][j] = (int)(table[i][j] / normalizationFactor);
|
||||
}
|
||||
|
||||
return normalized;
|
||||
}
|
||||
|
||||
private static int [][] copyContingencyTable(int [][] t) {
|
||||
int[][] c = new int[2][2];
|
||||
|
||||
for ( int i = 0; i < 2; i++ )
|
||||
for ( int j = 0; j < 2; j++ )
|
||||
c[i][j] = t[i][j];
|
||||
|
||||
return c;
|
||||
}
|
||||
|
||||
|
||||
private static void printTable(int[][] table, double pValue) {
|
||||
logger.info(String.format("%d %d; %d %d : %f", table[0][0], table[0][1], table[1][0], table[1][1], pValue));
|
||||
}
|
||||
|
||||
/**
|
||||
* Printing information to logger.info for debugging purposes
|
||||
*
|
||||
* @param name the name of the table
|
||||
* @param table the table itself
|
||||
*/
|
||||
private void printTable(final String name, final int[][] table) {
|
||||
if ( ENABLE_DEBUGGING ) {
|
||||
final String pValue = (String)annotationForOneTable(pValueForContingencyTable(table)).get(FS);
|
||||
logger.info(String.format("FS %s (REF+, REF-, ALT+, ALT-) = (%d, %d, %d, %d) = %s",
|
||||
name, table[0][0], table[0][1], table[1][0], table[1][1], pValue));
|
||||
}
|
||||
}
|
||||
|
||||
private static boolean rotateTable(int[][] table) {
|
||||
table[0][0] -= 1;
|
||||
table[1][0] += 1;
|
||||
|
||||
table[0][1] += 1;
|
||||
table[1][1] -= 1;
|
||||
|
||||
return (table[0][0] >= 0 && table[1][1] >= 0);
|
||||
}
|
||||
|
||||
private static boolean unrotateTable(int[][] table) {
|
||||
table[0][0] += 1;
|
||||
table[1][0] -= 1;
|
||||
|
||||
table[0][1] -= 1;
|
||||
table[1][1] += 1;
|
||||
|
||||
return (table[0][1] >= 0 && table[1][0] >= 0);
|
||||
}
|
||||
|
||||
private static double computePValue(int[][] table) {
|
||||
|
||||
int[] rowSums = { sumRow(table, 0), sumRow(table, 1) };
|
||||
int[] colSums = { sumColumn(table, 0), sumColumn(table, 1) };
|
||||
int N = rowSums[0] + rowSums[1];
|
||||
|
||||
// calculate in log space so we don't die with high numbers
|
||||
double pCutoff = Arithmetic.logFactorial(rowSums[0])
|
||||
+ Arithmetic.logFactorial(rowSums[1])
|
||||
+ Arithmetic.logFactorial(colSums[0])
|
||||
+ Arithmetic.logFactorial(colSums[1])
|
||||
- Arithmetic.logFactorial(table[0][0])
|
||||
- Arithmetic.logFactorial(table[0][1])
|
||||
- Arithmetic.logFactorial(table[1][0])
|
||||
- Arithmetic.logFactorial(table[1][1])
|
||||
- Arithmetic.logFactorial(N);
|
||||
return Math.exp(pCutoff);
|
||||
}
|
||||
|
||||
private static int sumRow(int[][] table, int column) {
|
||||
int sum = 0;
|
||||
for (int r = 0; r < table.length; r++) {
|
||||
sum += table[r][column];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
private static int sumColumn(int[][] table, int row) {
|
||||
int sum = 0;
|
||||
for (int c = 0; c < table[row].length; c++) {
|
||||
sum += table[row][c];
|
||||
}
|
||||
|
||||
return sum;
|
||||
}
|
||||
|
||||
/**
|
||||
Allocate and fill a 2x2 strand contingency table. In the end, it'll look something like this:
|
||||
* fw rc
|
||||
* allele1 # #
|
||||
* allele2 # #
|
||||
* @return a 2x2 contingency table
|
||||
*/
|
||||
public static int[][] getContingencyTable( final Map<String, PerReadAlleleLikelihoodMap> stratifiedPerReadAlleleLikelihoodMap, final VariantContext vc) {
|
||||
if( stratifiedPerReadAlleleLikelihoodMap == null ) { throw new IllegalArgumentException("stratifiedPerReadAlleleLikelihoodMap cannot be null"); }
|
||||
if( vc == null ) { throw new IllegalArgumentException("input vc cannot be null"); }
|
||||
|
||||
final Allele ref = vc.getReference();
|
||||
final Allele alt = vc.getAltAlleleWithHighestAlleleCount();
|
||||
final int[][] table = new int[2][2];
|
||||
|
||||
for (final PerReadAlleleLikelihoodMap maps : stratifiedPerReadAlleleLikelihoodMap.values() ) {
|
||||
final int[] myTable = new int[4];
|
||||
for (final Map.Entry<GATKSAMRecord,Map<Allele,Double>> el : maps.getLikelihoodReadMap().entrySet()) {
|
||||
final MostLikelyAllele mostLikelyAllele = PerReadAlleleLikelihoodMap.getMostLikelyAllele(el.getValue());
|
||||
final GATKSAMRecord read = el.getKey();
|
||||
final int representativeCount = read.isReducedRead() ? read.getReducedCount(ReadUtils.getReadCoordinateForReferenceCoordinateUpToEndOfRead(read, vc.getStart(), ReadUtils.ClippingTail.RIGHT_TAIL)) : 1;
|
||||
updateTable(myTable, mostLikelyAllele.getAlleleIfInformative(), read, ref, alt, representativeCount);
|
||||
}
|
||||
if ( passesMinimumThreshold(myTable) )
|
||||
copyToMainTable(myTable, table);
|
||||
}
|
||||
|
||||
return table;
|
||||
}
|
||||
|
||||
/**
|
||||
* Helper method to copy the per-sample table to the main table
|
||||
*
|
||||
* @param perSampleTable per-sample table (single dimension)
|
||||
* @param mainTable main table (two dimensions)
|
||||
*/
|
||||
private static void copyToMainTable(final int[] perSampleTable, final int[][] mainTable) {
|
||||
mainTable[0][0] += perSampleTable[0];
|
||||
mainTable[0][1] += perSampleTable[1];
|
||||
mainTable[1][0] += perSampleTable[2];
|
||||
mainTable[1][1] += perSampleTable[3];
|
||||
}
|
||||
|
||||
/**
|
||||
Allocate and fill a 2x2 strand contingency table. In the end, it'll look something like this:
|
||||
* fw rc
|
||||
* allele1 # #
|
||||
* allele2 # #
|
||||
* @return a 2x2 contingency table
|
||||
*/
|
||||
private static int[][] getSNPContingencyTable(final Map<String, AlignmentContext> stratifiedContexts,
|
||||
final Allele ref,
|
||||
final Allele alt,
|
||||
final int minQScoreToConsider ) {
|
||||
int[][] table = new int[2][2];
|
||||
|
||||
for ( Map.Entry<String, AlignmentContext> sample : stratifiedContexts.entrySet() ) {
|
||||
final int[] myTable = new int[4];
|
||||
for (PileupElement p : sample.getValue().getBasePileup()) {
|
||||
|
||||
if ( ! isUsableBase(p) ) // ignore deletions and bad MQ
|
||||
continue;
|
||||
|
||||
if ( p.getQual() < minQScoreToConsider || p.getMappingQual() < minQScoreToConsider )
|
||||
continue;
|
||||
|
||||
updateTable(myTable, Allele.create(p.getBase(), false), p.getRead(), ref, alt, p.getRepresentativeCount());
|
||||
}
|
||||
if ( passesMinimumThreshold(myTable) )
|
||||
copyToMainTable(myTable, table);
|
||||
}
|
||||
|
||||
return table;
|
||||
}
|
||||
|
||||
/**
|
||||
* Can the base in this pileup element be used in comparative tests?
|
||||
*
|
||||
* @param p the pileup element to consider
|
||||
*
|
||||
* @return true if this base is part of a meaningful read for comparison, false otherwise
|
||||
*/
|
||||
private static boolean isUsableBase(final PileupElement p) {
|
||||
return !( p.isDeletion() ||
|
||||
p.getMappingQual() == 0 ||
|
||||
p.getMappingQual() == QualityUtils.MAPPING_QUALITY_UNAVAILABLE ||
|
||||
((int) p.getQual()) < QualityUtils.MIN_USABLE_Q_SCORE);
|
||||
}
|
||||
|
||||
private static void updateTable(final int[] table, final Allele allele, final GATKSAMRecord read, final Allele ref, final Allele alt, final int representativeCount) {
|
||||
|
||||
final boolean matchesRef = allele.equals(ref, true);
|
||||
final boolean matchesAlt = allele.equals(alt, true);
|
||||
|
||||
if ( matchesRef || matchesAlt ) {
|
||||
final int offset = matchesRef ? 0 : 2;
|
||||
|
||||
if ( read.isStrandless() ) {
|
||||
|
||||
// ignore strandless reduced reads because they are always on the forward strand!
|
||||
if ( !read.isReducedRead() ) {
|
||||
|
||||
// a strandless read counts as observations on both strand, at 50% weight, with a minimum of 1
|
||||
// (the 1 is to ensure that a strandless read always counts as an observation on both strands, even
|
||||
// if the read is only seen once, because it's a merged read or other)
|
||||
final int toAdd = Math.max(representativeCount / 2, 1);
|
||||
table[offset] += toAdd;
|
||||
table[offset + 1] += toAdd;
|
||||
}
|
||||
} else {
|
||||
// a normal read with an actual strand
|
||||
final boolean isFW = !read.getReadNegativeStrandFlag();
|
||||
table[offset + (isFW ? 0 : 1)] += representativeCount;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,191 @@
|
|||
/*
|
||||
* By downloading the PROGRAM you agree to the following terms of use:
|
||||
*
|
||||
* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
|
||||
*
|
||||
* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
|
||||
*
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
|
||||
* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
|
||||
*
|
||||
* 1. DEFINITIONS
|
||||
* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
|
||||
*
|
||||
* 2. LICENSE
|
||||
* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
|
||||
* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
|
||||
* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
|
||||
* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
|
||||
*
|
||||
* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
|
||||
* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
|
||||
* Copyright 2012 Broad Institute, Inc.
|
||||
* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
|
||||
* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
|
||||
*
|
||||
* 4. INDEMNIFICATION
|
||||
* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
|
||||
*
|
||||
* 5. NO REPRESENTATIONS OR WARRANTIES
|
||||
* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
|
||||
* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
|
||||
*
|
||||
* 6. ASSIGNMENT
|
||||
* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
|
||||
*
|
||||
* 7. MISCELLANEOUS
|
||||
* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
|
||||
* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
|
||||
* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
|
||||
* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
|
||||
* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
|
||||
* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
|
||||
* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
|
||||
*/
|
||||
|
||||
package org.broadinstitute.sting.gatk.walkers.annotator;
|
||||
|
||||
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
|
||||
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.annotator.interfaces.ActiveRegionBasedAnnotation;
|
||||
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.AnnotatorCompatible;
|
||||
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.InfoFieldAnnotation;
|
||||
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.StandardAnnotation;
|
||||
import org.broadinstitute.sting.utils.MathUtils;
|
||||
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
|
||||
import org.broadinstitute.sting.utils.variant.GATKVariantContextUtils;
|
||||
import org.broadinstitute.variant.vcf.VCFHeaderLineType;
|
||||
import org.broadinstitute.variant.vcf.VCFInfoHeaderLine;
|
||||
import org.broadinstitute.variant.variantcontext.Genotype;
|
||||
import org.broadinstitute.variant.variantcontext.GenotypesContext;
|
||||
import org.broadinstitute.variant.variantcontext.VariantContext;
|
||||
|
||||
import java.util.*;
|
||||
|
||||
/**
|
||||
* Variant confidence (from the QUAL field) / unfiltered depth of non-reference samples. Note that the QD is also normalized by event length.
|
||||
*
|
||||
* Low scores are indicative of false positive calls and artifacts. Note that QualByDepth requires sequencing
|
||||
* reads associated with the samples with polymorphic genotypes.
|
||||
*/
|
||||
public class QualByDepth extends InfoFieldAnnotation implements StandardAnnotation, ActiveRegionBasedAnnotation {
|
||||
// private final static Logger logger = Logger.getLogger(QualByDepth.class);
|
||||
|
||||
public Map<String, Object> annotate(final RefMetaDataTracker tracker,
|
||||
final AnnotatorCompatible walker,
|
||||
final ReferenceContext ref,
|
||||
final Map<String, AlignmentContext> stratifiedContexts,
|
||||
final VariantContext vc,
|
||||
final Map<String, PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap ) {
|
||||
if ( !vc.hasLog10PError() )
|
||||
return null;
|
||||
|
||||
final GenotypesContext genotypes = vc.getGenotypes();
|
||||
if ( genotypes == null || genotypes.size() == 0 )
|
||||
return null;
|
||||
|
||||
int standardDepth = 0;
|
||||
int ADrestrictedDepth = 0;
|
||||
|
||||
for ( final Genotype genotype : genotypes ) {
|
||||
|
||||
// we care only about variant calls with likelihoods
|
||||
if ( !genotype.isHet() && !genotype.isHomVar() )
|
||||
continue;
|
||||
|
||||
// if we have the AD values for this sample, let's make sure that the variant depth is greater than 1!
|
||||
// TODO -- If we like how this is working and want to apply it to a situation other than the single sample HC pipeline,
|
||||
// TODO -- then we will need to modify the annotateContext() - and related - routines in the VariantAnnotatorEngine
|
||||
// TODO -- so that genotype-level annotations are run first (to generate AD on the samples) and then the site-level
|
||||
// TODO -- annotations must come afterwards (so that QD can use the AD).
|
||||
if ( genotype.hasAD() ) {
|
||||
final int[] AD = genotype.getAD();
|
||||
final int totalADdepth = (int)MathUtils.sum(AD);
|
||||
if ( totalADdepth - AD[0] > 1 )
|
||||
ADrestrictedDepth += totalADdepth;
|
||||
standardDepth += totalADdepth;
|
||||
continue;
|
||||
}
|
||||
|
||||
if (stratifiedContexts!= null && !stratifiedContexts.isEmpty()) {
|
||||
final AlignmentContext context = stratifiedContexts.get(genotype.getSampleName());
|
||||
if ( context == null )
|
||||
continue;
|
||||
standardDepth += context.getBasePileup().depthOfCoverage();
|
||||
|
||||
} else if (perReadAlleleLikelihoodMap != null) {
|
||||
final PerReadAlleleLikelihoodMap perReadAlleleLikelihoods = perReadAlleleLikelihoodMap.get(genotype.getSampleName());
|
||||
if (perReadAlleleLikelihoods == null || perReadAlleleLikelihoods.isEmpty())
|
||||
continue;
|
||||
|
||||
standardDepth += perReadAlleleLikelihoods.getNumberOfStoredElements();
|
||||
} else if ( genotype.hasDP() ) {
|
||||
standardDepth += genotype.getDP();
|
||||
}
|
||||
}
|
||||
|
||||
// if the AD-restricted depth is a usable value (i.e. not zero), then we should use that one going forward
|
||||
if ( ADrestrictedDepth > 0 )
|
||||
standardDepth = ADrestrictedDepth;
|
||||
|
||||
if ( standardDepth == 0 )
|
||||
return null;
|
||||
|
||||
final double altAlleleLength = GATKVariantContextUtils.getMeanAltAlleleLength(vc);
|
||||
// Hack: when refContext == null then we know we are coming from the HaplotypeCaller and do not want to do a
|
||||
// full length-based normalization (because the indel length problem is present only in the UnifiedGenotyper)
|
||||
double QD = -10.0 * vc.getLog10PError() / ((double)standardDepth * indelNormalizationFactor(altAlleleLength, ref != null));
|
||||
|
||||
// Hack: see note in the fixTooHighQD method below
|
||||
QD = fixTooHighQD(QD);
|
||||
|
||||
final Map<String, Object> map = new HashMap<>();
|
||||
map.put(getKeyNames().get(0), String.format("%.2f", QD));
|
||||
return map;
|
||||
}
|
||||
|
||||
/**
|
||||
* Generate the indel normalization factor.
|
||||
*
|
||||
* @param altAlleleLength the average alternate allele length for the call
|
||||
* @param increaseNormalizationAsLengthIncreases should we apply a normalization factor based on the allele length?
|
||||
* @return a possitive double
|
||||
*/
|
||||
private double indelNormalizationFactor(final double altAlleleLength, final boolean increaseNormalizationAsLengthIncreases) {
|
||||
return ( increaseNormalizationAsLengthIncreases ? Math.max(altAlleleLength / 3.0, 1.0) : 1.0);
|
||||
}
|
||||
|
||||
/**
|
||||
* The haplotype caller generates very high quality scores when multiple events are on the
|
||||
* same haplotype. This causes some very good variants to have unusually high QD values,
|
||||
* and VQSR will filter these out. This code looks at the QD value, and if it is above
|
||||
* threshold we map it down to the mean high QD value, with some jittering
|
||||
*
|
||||
* // TODO -- remove me when HaplotypeCaller bubble caller is live
|
||||
*
|
||||
* @param QD the raw QD score
|
||||
* @return a QD value
|
||||
*/
|
||||
private double fixTooHighQD(final double QD) {
|
||||
if ( QD < MAX_QD_BEFORE_FIXING ) {
|
||||
return QD;
|
||||
} else {
|
||||
return IDEAL_HIGH_QD + GenomeAnalysisEngine.getRandomGenerator().nextGaussian() * JITTER_SIGMA;
|
||||
}
|
||||
}
|
||||
|
||||
private final static double MAX_QD_BEFORE_FIXING = 35;
|
||||
private final static double IDEAL_HIGH_QD = 30;
|
||||
private final static double JITTER_SIGMA = 3;
|
||||
|
||||
public List<String> getKeyNames() { return Arrays.asList("QD"); }
|
||||
|
||||
public List<VCFInfoHeaderLine> getDescriptions() {
|
||||
return Arrays.asList(new VCFInfoHeaderLine(getKeyNames().get(0), 1, VCFHeaderLineType.Float, "Variant Confidence/Quality by Depth"));
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
|
@ -0,0 +1,99 @@
|
|||
/*
|
||||
* By downloading the PROGRAM you agree to the following terms of use:
|
||||
*
|
||||
* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
|
||||
*
|
||||
* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
|
||||
*
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
|
||||
* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
|
||||
*
|
||||
* 1. DEFINITIONS
|
||||
* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
|
||||
*
|
||||
* 2. LICENSE
|
||||
* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
|
||||
* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
|
||||
* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
|
||||
* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
|
||||
*
|
||||
* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
|
||||
* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
|
||||
* Copyright 2012 Broad Institute, Inc.
|
||||
* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
|
||||
* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
|
||||
*
|
||||
* 4. INDEMNIFICATION
|
||||
* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
|
||||
*
|
||||
* 5. NO REPRESENTATIONS OR WARRANTIES
|
||||
* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
|
||||
* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
|
||||
*
|
||||
* 6. ASSIGNMENT
|
||||
* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
|
||||
*
|
||||
* 7. MISCELLANEOUS
|
||||
* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
|
||||
* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
|
||||
* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
|
||||
* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
|
||||
* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
|
||||
* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
|
||||
* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
|
||||
*/
|
||||
|
||||
package org.broadinstitute.sting.gatk.walkers.annotator;
|
||||
|
||||
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.annotator.interfaces.AnnotatorCompatible;
|
||||
import org.broadinstitute.sting.gatk.walkers.annotator.interfaces.GenotypeAnnotation;
|
||||
import org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap;
|
||||
import org.broadinstitute.variant.variantcontext.Genotype;
|
||||
import org.broadinstitute.variant.variantcontext.GenotypeBuilder;
|
||||
import org.broadinstitute.variant.variantcontext.VariantContext;
|
||||
import org.broadinstitute.variant.vcf.VCFFormatHeaderLine;
|
||||
import org.broadinstitute.variant.vcf.VCFHeaderLineType;
|
||||
|
||||
import java.util.*;
|
||||
|
||||
/**
|
||||
* Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias
|
||||
* User: rpoplin
|
||||
* Date: 8/28/13
|
||||
*/
|
||||
|
||||
public class StrandBiasBySample extends GenotypeAnnotation {
|
||||
|
||||
public final static String STRAND_BIAS_BY_SAMPLE_KEY_NAME = "SB";
|
||||
|
||||
@Override
|
||||
public void annotate(final RefMetaDataTracker tracker,
|
||||
final AnnotatorCompatible walker,
|
||||
final ReferenceContext ref,
|
||||
final AlignmentContext stratifiedContext,
|
||||
final VariantContext vc,
|
||||
final Genotype g,
|
||||
final GenotypeBuilder gb,
|
||||
final PerReadAlleleLikelihoodMap alleleLikelihoodMap) {
|
||||
if ( ! isAppropriateInput(alleleLikelihoodMap, g) )
|
||||
return;
|
||||
|
||||
final int[][] table = FisherStrand.getContingencyTable(Collections.singletonMap(g.getSampleName(), alleleLikelihoodMap), vc);
|
||||
|
||||
gb.attribute(STRAND_BIAS_BY_SAMPLE_KEY_NAME, FisherStrand.getContingencyArray(table));
|
||||
}
|
||||
|
||||
@Override
|
||||
public List<String> getKeyNames() { return Collections.singletonList(STRAND_BIAS_BY_SAMPLE_KEY_NAME); }
|
||||
|
||||
@Override
|
||||
public List<VCFFormatHeaderLine> getDescriptions() { return Collections.singletonList(new VCFFormatHeaderLine(getKeyNames().get(0), 4, VCFHeaderLineType.Integer, "Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.")); }
|
||||
|
||||
private boolean isAppropriateInput(final PerReadAlleleLikelihoodMap map, final Genotype g) {
|
||||
return ! (map == null || g == null || !g.isCalled());
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,844 @@
|
|||
/*
|
||||
* By downloading the PROGRAM you agree to the following terms of use:
|
||||
*
|
||||
* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
|
||||
*
|
||||
* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
|
||||
*
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
|
||||
* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
|
||||
*
|
||||
* 1. DEFINITIONS
|
||||
* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
|
||||
*
|
||||
* 2. LICENSE
|
||||
* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
|
||||
* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
|
||||
* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
|
||||
* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
|
||||
*
|
||||
* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
|
||||
* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
|
||||
* Copyright 2012 Broad Institute, Inc.
|
||||
* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
|
||||
* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
|
||||
*
|
||||
* 4. INDEMNIFICATION
|
||||
* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
|
||||
*
|
||||
* 5. NO REPRESENTATIONS OR WARRANTIES
|
||||
* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
|
||||
* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
|
||||
*
|
||||
* 6. ASSIGNMENT
|
||||
* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
|
||||
*
|
||||
* 7. MISCELLANEOUS
|
||||
* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
|
||||
* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
|
||||
* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
|
||||
* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
|
||||
* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
|
||||
* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
|
||||
* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
|
||||
*/
|
||||
|
||||
package org.broadinstitute.sting.gatk.walkers.genotyper;
|
||||
|
||||
import com.google.java.contract.Ensures;
|
||||
import com.google.java.contract.Requires;
|
||||
import org.apache.log4j.Logger;
|
||||
import org.broadinstitute.sting.commandline.RodBinding;
|
||||
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
|
||||
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
|
||||
import org.broadinstitute.sting.gatk.contexts.AlignmentContextUtils;
|
||||
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
|
||||
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
|
||||
import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotatorEngine;
|
||||
import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.AFCalc;
|
||||
import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.AFCalcFactory;
|
||||
import org.broadinstitute.sting.gatk.walkers.genotyper.afcalc.AFCalcResult;
|
||||
import org.broadinstitute.sting.utils.*;
|
||||
import org.broadinstitute.sting.utils.baq.BAQ;
|
||||
import org.broadinstitute.sting.utils.classloader.PluginManager;
|
||||
import org.broadinstitute.sting.utils.variant.GATKVariantContextUtils;
|
||||
import org.broadinstitute.sting.utils.BaseUtils;
|
||||
import org.broadinstitute.variant.vcf.VCFConstants;
|
||||
import org.broadinstitute.sting.utils.exceptions.UserException;
|
||||
import org.broadinstitute.sting.utils.pileup.PileupElement;
|
||||
import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
|
||||
import org.broadinstitute.variant.variantcontext.*;
|
||||
|
||||
import java.io.PrintStream;
|
||||
import java.lang.reflect.Constructor;
|
||||
import java.util.*;
|
||||
|
||||
public class UnifiedGenotyperEngine {
|
||||
public static final String LOW_QUAL_FILTER_NAME = "LowQual";
|
||||
private static final String GPSTRING = "GENERALPLOIDY";
|
||||
|
||||
public static final String NUMBER_OF_DISCOVERED_ALLELES_KEY = "NDA";
|
||||
public static final String PL_FOR_ALL_SNP_ALLELES_KEY = "APL";
|
||||
|
||||
public static final double HUMAN_SNP_HETEROZYGOSITY = 1e-3;
|
||||
public static final double HUMAN_INDEL_HETEROZYGOSITY = 1e-4;
|
||||
|
||||
private static final int SNP_MODEL = 0;
|
||||
private static final int INDEL_MODEL = 1;
|
||||
|
||||
public enum OUTPUT_MODE {
|
||||
/** produces calls only at variant sites */
|
||||
EMIT_VARIANTS_ONLY,
|
||||
/** produces calls at variant sites and confident reference sites */
|
||||
EMIT_ALL_CONFIDENT_SITES,
|
||||
/** produces calls at any callable site regardless of confidence; this argument is intended only for point
|
||||
* mutations (SNPs) in DISCOVERY mode or generally when running in GENOTYPE_GIVEN_ALLELES mode; it will by
|
||||
* no means produce a comprehensive set of indels in DISCOVERY mode */
|
||||
EMIT_ALL_SITES
|
||||
}
|
||||
|
||||
// the unified argument collection
|
||||
private final UnifiedArgumentCollection UAC;
|
||||
public UnifiedArgumentCollection getUAC() { return UAC; }
|
||||
|
||||
// the annotation engine
|
||||
private final VariantAnnotatorEngine annotationEngine;
|
||||
|
||||
// the model used for calculating genotypes
|
||||
private ThreadLocal<Map<String, GenotypeLikelihoodsCalculationModel>> glcm = new ThreadLocal<Map<String, GenotypeLikelihoodsCalculationModel>>();
|
||||
private final List<GenotypeLikelihoodsCalculationModel.Model> modelsToUse = new ArrayList<GenotypeLikelihoodsCalculationModel.Model>(2);
|
||||
|
||||
// the model used for calculating p(non-ref)
|
||||
private ThreadLocal<AFCalc> afcm = new ThreadLocal<AFCalc>();
|
||||
|
||||
// because the allele frequency priors are constant for a given i, we cache the results to avoid having to recompute everything
|
||||
private final double[] log10AlleleFrequencyPriorsSNPs;
|
||||
private final double[] log10AlleleFrequencyPriorsIndels;
|
||||
|
||||
// samples in input
|
||||
private final Set<String> samples;
|
||||
|
||||
// the various loggers and writers
|
||||
private final Logger logger;
|
||||
private final PrintStream verboseWriter;
|
||||
|
||||
// number of chromosomes (ploidy * samples) in input
|
||||
private final int ploidy;
|
||||
private final int N;
|
||||
|
||||
// the standard filter to use for calls below the confidence threshold but above the emit threshold
|
||||
private static final Set<String> filter = new HashSet<String>(1);
|
||||
|
||||
private final GenomeLocParser genomeLocParser;
|
||||
private final boolean BAQEnabledOnCMDLine;
|
||||
|
||||
// ---------------------------------------------------------------------------------------------------------
|
||||
//
|
||||
// Public interface functions
|
||||
//
|
||||
// ---------------------------------------------------------------------------------------------------------
|
||||
@Requires({"toolkit != null", "UAC != null"})
|
||||
public UnifiedGenotyperEngine(GenomeAnalysisEngine toolkit, UnifiedArgumentCollection UAC) {
|
||||
this(toolkit, UAC, Logger.getLogger(UnifiedGenotyperEngine.class), null, null, SampleUtils.getSAMFileSamples(toolkit.getSAMFileHeader()), GATKVariantContextUtils.DEFAULT_PLOIDY);
|
||||
}
|
||||
|
||||
protected UnifiedGenotyperEngine(GenomeAnalysisEngine toolkit, Set<String> samples, UnifiedArgumentCollection UAC) {
|
||||
this(toolkit, UAC, Logger.getLogger(UnifiedGenotyperEngine.class), null, null, samples, GATKVariantContextUtils.DEFAULT_PLOIDY);
|
||||
}
|
||||
|
||||
@Requires({"toolkit != null", "UAC != null", "logger != null", "samples != null && samples.size() > 0","ploidy>0"})
|
||||
public UnifiedGenotyperEngine(GenomeAnalysisEngine toolkit, UnifiedArgumentCollection UAC, Logger logger, PrintStream verboseWriter, VariantAnnotatorEngine engine, Set<String> samples, int ploidy) {
|
||||
this.BAQEnabledOnCMDLine = toolkit.getArguments().BAQMode != BAQ.CalculationMode.OFF;
|
||||
genomeLocParser = toolkit.getGenomeLocParser();
|
||||
this.samples = new TreeSet<String>(samples);
|
||||
// note that, because we cap the base quality by the mapping quality, minMQ cannot be less than minBQ
|
||||
this.UAC = UAC;
|
||||
|
||||
this.logger = logger;
|
||||
this.verboseWriter = verboseWriter;
|
||||
this.annotationEngine = engine;
|
||||
|
||||
this.ploidy = ploidy;
|
||||
this.N = samples.size() * ploidy;
|
||||
log10AlleleFrequencyPriorsSNPs = new double[N+1];
|
||||
log10AlleleFrequencyPriorsIndels = new double[N+1];
|
||||
computeAlleleFrequencyPriors(N, log10AlleleFrequencyPriorsSNPs, UAC.heterozygosity,UAC.inputPrior);
|
||||
computeAlleleFrequencyPriors(N, log10AlleleFrequencyPriorsIndels, UAC.INDEL_HETEROZYGOSITY, UAC.inputPrior);
|
||||
|
||||
filter.add(LOW_QUAL_FILTER_NAME);
|
||||
|
||||
determineGLModelsToUse();
|
||||
|
||||
// do argument checking
|
||||
if (UAC.annotateAllSitesWithPLs) {
|
||||
if (!modelsToUse.contains(GenotypeLikelihoodsCalculationModel.Model.SNP))
|
||||
throw new IllegalArgumentException("Invalid genotype likelihood model specification: Only diploid SNP model can be used in conjunction with option allSitePLs");
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* @see #calculateLikelihoodsAndGenotypes(org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker, org.broadinstitute.sting.gatk.contexts.ReferenceContext, org.broadinstitute.sting.gatk.contexts.AlignmentContext, java.util.Set)
|
||||
*
|
||||
* same as the full call but with allSamples == null
|
||||
*
|
||||
* @param tracker the meta data tracker
|
||||
* @param refContext the reference base
|
||||
* @param rawContext contextual information around the locus
|
||||
* @return the VariantCallContext object
|
||||
*/
|
||||
public List<VariantCallContext> calculateLikelihoodsAndGenotypes(final RefMetaDataTracker tracker,
|
||||
final ReferenceContext refContext,
|
||||
final AlignmentContext rawContext) {
|
||||
return calculateLikelihoodsAndGenotypes(tracker, refContext, rawContext, null);
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Compute full calls at a given locus. Entry point for engine calls from the UnifiedGenotyper.
|
||||
*
|
||||
* If allSamples != null, then the output variantCallContext is guarenteed to contain a genotype
|
||||
* for every sample in allSamples. If it's null there's no such guarentee. Providing this
|
||||
* argument is critical when the resulting calls will be written to a VCF file.
|
||||
*
|
||||
* @param tracker the meta data tracker
|
||||
* @param refContext the reference base
|
||||
* @param rawContext contextual information around the locus
|
||||
* @param allSamples set of all sample names that we might call (i.e., those in the VCF header)
|
||||
* @return the VariantCallContext object
|
||||
*/
|
||||
public List<VariantCallContext> calculateLikelihoodsAndGenotypes(final RefMetaDataTracker tracker,
|
||||
final ReferenceContext refContext,
|
||||
final AlignmentContext rawContext,
|
||||
final Set<String> allSamples) {
|
||||
final List<VariantCallContext> results = new ArrayList<VariantCallContext>(2);
|
||||
|
||||
final List<GenotypeLikelihoodsCalculationModel.Model> models = getGLModelsToUse(tracker, refContext, rawContext);
|
||||
|
||||
final Map<String, org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap = new HashMap<String, org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap>();
|
||||
|
||||
if ( models.isEmpty() ) {
|
||||
results.add(UAC.OutputMode == OUTPUT_MODE.EMIT_ALL_SITES && UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ? generateEmptyContext(tracker, refContext, null, rawContext) : null);
|
||||
}
|
||||
else {
|
||||
for ( final GenotypeLikelihoodsCalculationModel.Model model : models ) {
|
||||
perReadAlleleLikelihoodMap.clear();
|
||||
final Map<String, AlignmentContext> stratifiedContexts = getFilteredAndStratifiedContexts(UAC, refContext, rawContext, model);
|
||||
if ( stratifiedContexts == null ) {
|
||||
results.add(UAC.OutputMode == OUTPUT_MODE.EMIT_ALL_SITES && UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ? generateEmptyContext(tracker, refContext, null, rawContext) : null);
|
||||
}
|
||||
else {
|
||||
final VariantContext vc = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.COMPLETE, null, true, model, perReadAlleleLikelihoodMap);
|
||||
if ( vc != null )
|
||||
results.add(calculateGenotypes(tracker, refContext, rawContext, stratifiedContexts, vc, model, true, perReadAlleleLikelihoodMap));
|
||||
// todo - uncomment if we want to also emit a null ref call (with no QUAL) if there's no evidence for REF and if EMIT_ALL_SITES is set
|
||||
// else if (UAC.OutputMode == OUTPUT_MODE.EMIT_ALL_SITES)
|
||||
// results.add(generateEmptyContext(tracker, refContext, null, rawContext));
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return results;
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute GLs at a given locus. Entry point for engine calls from UGCalcLikelihoods.
|
||||
*
|
||||
* @param tracker the meta data tracker
|
||||
* @param refContext the reference base
|
||||
* @param rawContext contextual information around the locus
|
||||
* @param perReadAlleleLikelihoodMap Map to store per-sample, per-read, per-allele likelihoods (only used for indels)
|
||||
* @return the VariantContext object
|
||||
*/
|
||||
public VariantContext calculateLikelihoods(final RefMetaDataTracker tracker,
|
||||
final ReferenceContext refContext,
|
||||
final AlignmentContext rawContext,
|
||||
final Map<String, org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap) {
|
||||
final List<GenotypeLikelihoodsCalculationModel.Model> models = getGLModelsToUse(tracker, refContext, rawContext);
|
||||
if ( models.isEmpty() ) {
|
||||
return null;
|
||||
}
|
||||
|
||||
for ( final GenotypeLikelihoodsCalculationModel.Model model : models ) {
|
||||
final Map<String, AlignmentContext> stratifiedContexts = getFilteredAndStratifiedContexts(UAC, refContext, rawContext, model);
|
||||
// return the first valid one we encounter
|
||||
if ( stratifiedContexts != null )
|
||||
return calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.COMPLETE, null, true, model, perReadAlleleLikelihoodMap);
|
||||
|
||||
}
|
||||
|
||||
return null;
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute genotypes at a given locus. Entry point for engine calls from UGCallVariants.
|
||||
*
|
||||
* @param tracker the meta data tracker
|
||||
* @param refContext the reference base
|
||||
* @param rawContext contextual information around the locus
|
||||
* @param vc the GL-annotated variant context
|
||||
* @return the VariantCallContext object
|
||||
*/
|
||||
public VariantCallContext calculateGenotypes(final RefMetaDataTracker tracker,
|
||||
final ReferenceContext refContext,
|
||||
final AlignmentContext rawContext,
|
||||
final VariantContext vc) {
|
||||
final List<GenotypeLikelihoodsCalculationModel.Model> models = getGLModelsToUse(tracker, refContext, rawContext);
|
||||
if ( models.isEmpty() ) {
|
||||
return null;
|
||||
}
|
||||
|
||||
// return the first one
|
||||
final GenotypeLikelihoodsCalculationModel.Model model = models.get(0);
|
||||
final Map<String, AlignmentContext> stratifiedContexts = getFilteredAndStratifiedContexts(UAC, refContext, rawContext, model);
|
||||
return calculateGenotypes(tracker, refContext, rawContext, stratifiedContexts, vc, model, null);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute genotypes at a given locus.
|
||||
*
|
||||
* @param vc the GL-annotated variant context
|
||||
* @return the VariantCallContext object
|
||||
*/
|
||||
public VariantCallContext calculateGenotypes(VariantContext vc) {
|
||||
return calculateGenotypes(null, null, null, null, vc, GenotypeLikelihoodsCalculationModel.Model.valueOf("SNP"), null);
|
||||
}
|
||||
|
||||
|
||||
// ---------------------------------------------------------------------------------------------------------
|
||||
//
|
||||
// Private implementation helpers
|
||||
//
|
||||
// ---------------------------------------------------------------------------------------------------------
|
||||
|
||||
// private method called by both UnifiedGenotyper and UGCalcLikelihoods entry points into the engine
|
||||
private VariantContext calculateLikelihoods(final RefMetaDataTracker tracker,
|
||||
final ReferenceContext refContext,
|
||||
final Map<String, AlignmentContext> stratifiedContexts,
|
||||
final AlignmentContextUtils.ReadOrientation type,
|
||||
final List<Allele> alternateAllelesToUse,
|
||||
final boolean useBAQedPileup,
|
||||
final GenotypeLikelihoodsCalculationModel.Model model,
|
||||
final Map<String, org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap) {
|
||||
|
||||
// initialize the data for this thread if that hasn't been done yet
|
||||
if ( glcm.get() == null ) {
|
||||
glcm.set(getGenotypeLikelihoodsCalculationObject(logger, UAC));
|
||||
}
|
||||
|
||||
return glcm.get().get(model.name()).getLikelihoods(tracker, refContext, stratifiedContexts, type, alternateAllelesToUse, useBAQedPileup && BAQEnabledOnCMDLine, genomeLocParser, perReadAlleleLikelihoodMap);
|
||||
}
|
||||
|
||||
private VariantCallContext generateEmptyContext(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, AlignmentContext> stratifiedContexts, AlignmentContext rawContext) {
|
||||
VariantContext vc;
|
||||
if ( UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ) {
|
||||
VariantContext vcInput = UnifiedGenotyperEngine.getVCFromAllelesRod(tracker, ref, rawContext.getLocation(), false, logger, UAC.alleles);
|
||||
if ( vcInput == null )
|
||||
return null;
|
||||
vc = new VariantContextBuilder("UG_call", ref.getLocus().getContig(), vcInput.getStart(), vcInput.getEnd(), vcInput.getAlleles()).make();
|
||||
} else {
|
||||
// deal with bad/non-standard reference bases
|
||||
if ( !Allele.acceptableAlleleBases(new byte[]{ref.getBase()}) )
|
||||
return null;
|
||||
|
||||
Set<Allele> alleles = new HashSet<Allele>();
|
||||
alleles.add(Allele.create(ref.getBase(), true));
|
||||
vc = new VariantContextBuilder("UG_call", ref.getLocus().getContig(), ref.getLocus().getStart(), ref.getLocus().getStart(), alleles).make();
|
||||
}
|
||||
|
||||
if ( annotationEngine != null ) {
|
||||
// Note: we want to use the *unfiltered* and *unBAQed* context for the annotations
|
||||
final ReadBackedPileup pileup = rawContext.getBasePileup();
|
||||
stratifiedContexts = AlignmentContextUtils.splitContextBySampleName(pileup);
|
||||
|
||||
vc = annotationEngine.annotateContext(tracker, ref, stratifiedContexts, vc);
|
||||
}
|
||||
|
||||
return new VariantCallContext(vc, false);
|
||||
}
|
||||
|
||||
public VariantCallContext calculateGenotypes(final VariantContext vc, final GenotypeLikelihoodsCalculationModel.Model model, final Map<String, org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap) {
|
||||
return calculateGenotypes(null, null, null, null, vc, model, perReadAlleleLikelihoodMap);
|
||||
}
|
||||
|
||||
public VariantCallContext calculateGenotypes(final VariantContext vc, final GenotypeLikelihoodsCalculationModel.Model model) {
|
||||
return calculateGenotypes(null, null, null, null, vc, model, null);
|
||||
}
|
||||
|
||||
public VariantCallContext calculateGenotypes(final RefMetaDataTracker tracker,
|
||||
final ReferenceContext refContext,
|
||||
final AlignmentContext rawContext,
|
||||
final Map<String, AlignmentContext> stratifiedContexts,
|
||||
final VariantContext vc,
|
||||
final GenotypeLikelihoodsCalculationModel.Model model,
|
||||
final Map<String, org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap) {
|
||||
return calculateGenotypes(tracker, refContext, rawContext, stratifiedContexts, vc, model, false, perReadAlleleLikelihoodMap);
|
||||
}
|
||||
|
||||
/**
|
||||
* Main entry function to calculate genotypes of a given VC with corresponding GL's
|
||||
* @param tracker Tracker
|
||||
* @param refContext Reference context
|
||||
* @param rawContext Raw context
|
||||
* @param stratifiedContexts Stratified alignment contexts
|
||||
* @param vc Input VC
|
||||
* @param model GL calculation model
|
||||
* @param inheritAttributesFromInputVC Output VC will contain attributes inherited from input vc
|
||||
* @return VC with assigned genotypes
|
||||
*/
|
||||
public VariantCallContext calculateGenotypes(final RefMetaDataTracker tracker, final ReferenceContext refContext,
|
||||
final AlignmentContext rawContext, Map<String, AlignmentContext> stratifiedContexts,
|
||||
final VariantContext vc, final GenotypeLikelihoodsCalculationModel.Model model,
|
||||
final boolean inheritAttributesFromInputVC,
|
||||
final Map<String, org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap) {
|
||||
|
||||
boolean limitedContext = tracker == null || refContext == null || rawContext == null || stratifiedContexts == null;
|
||||
|
||||
// TODO TODO TODO TODO
|
||||
// REFACTOR THIS FUNCTION, TOO UNWIELDY!!
|
||||
|
||||
// initialize the data for this thread if that hasn't been done yet
|
||||
if ( afcm.get() == null ) {
|
||||
afcm.set(AFCalcFactory.createAFCalc(UAC, N, logger));
|
||||
}
|
||||
|
||||
// if input VC can't be genotyped, exit with either null VCC or, in case where we need to emit all sites, an empty call
|
||||
if (!canVCbeGenotyped(vc)) {
|
||||
if (UAC.OutputMode == OUTPUT_MODE.EMIT_ALL_SITES && !limitedContext)
|
||||
return generateEmptyContext(tracker, refContext, stratifiedContexts, rawContext);
|
||||
else
|
||||
return null;
|
||||
|
||||
}
|
||||
|
||||
// estimate our confidence in a reference call and return
|
||||
if ( vc.getNSamples() == 0 ) {
|
||||
if ( limitedContext )
|
||||
return null;
|
||||
return (UAC.OutputMode != OUTPUT_MODE.EMIT_ALL_SITES ?
|
||||
estimateReferenceConfidence(vc, stratifiedContexts, getTheta(model), false, 1.0) :
|
||||
generateEmptyContext(tracker, refContext, stratifiedContexts, rawContext));
|
||||
}
|
||||
|
||||
final AFCalcResult AFresult = afcm.get().getLog10PNonRef(vc, getAlleleFrequencyPriors(model));
|
||||
|
||||
// is the most likely frequency conformation AC=0 for all alternate alleles?
|
||||
boolean bestGuessIsRef = true;
|
||||
|
||||
// determine which alternate alleles have AF>0
|
||||
final List<Allele> myAlleles = new ArrayList<>(vc.getAlleles().size());
|
||||
final List<Integer> alleleCountsofMLE = new ArrayList<>(vc.getAlleles().size());
|
||||
myAlleles.add(vc.getReference());
|
||||
for ( int i = 0; i < AFresult.getAllelesUsedInGenotyping().size(); i++ ) {
|
||||
final Allele alternateAllele = AFresult.getAllelesUsedInGenotyping().get(i);
|
||||
if ( alternateAllele.isReference() )
|
||||
continue;
|
||||
|
||||
// Compute if the site is considered polymorphic with sufficient confidence relative to our
|
||||
// phred-scaled emission QUAL
|
||||
final boolean isNonRef = AFresult.isPolymorphicPhredScaledQual(alternateAllele, UAC.STANDARD_CONFIDENCE_FOR_EMITTING);
|
||||
final boolean toInclude = isNonRef || alternateAllele == GATKVariantContextUtils.NON_REF_SYMBOLIC_ALLELE ||
|
||||
UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ||
|
||||
UAC.annotateAllSitesWithPLs;
|
||||
|
||||
bestGuessIsRef &= !isNonRef;
|
||||
|
||||
if (toInclude) {
|
||||
myAlleles.add(alternateAllele);
|
||||
alleleCountsofMLE.add(AFresult.getAlleleCountAtMLE(alternateAllele));
|
||||
}
|
||||
}
|
||||
|
||||
final double PoFGT0 = Math.pow(10, AFresult.getLog10PosteriorOfAFGT0());
|
||||
|
||||
// note the math.abs is necessary because -10 * 0.0 => -0.0 which isn't nice
|
||||
final double phredScaledConfidence =
|
||||
Math.abs(! bestGuessIsRef || UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES || UAC.annotateAllSitesWithPLs
|
||||
? -10 * AFresult.getLog10PosteriorOfAFEq0()
|
||||
: -10 * AFresult.getLog10PosteriorOfAFGT0());
|
||||
|
||||
// return a null call if we don't pass the confidence cutoff or the most likely allele frequency is zero
|
||||
if ( UAC.OutputMode != OUTPUT_MODE.EMIT_ALL_SITES && !passesEmitThreshold(phredScaledConfidence, bestGuessIsRef) ) {
|
||||
// technically, at this point our confidence in a reference call isn't accurately estimated
|
||||
// because it didn't take into account samples with no data, so let's get a better estimate
|
||||
return limitedContext ? null : estimateReferenceConfidence(vc, stratifiedContexts, getTheta(model), true, PoFGT0);
|
||||
}
|
||||
|
||||
// start constructing the resulting VC
|
||||
final GenomeLoc loc = genomeLocParser.createGenomeLoc(vc);
|
||||
final VariantContextBuilder builder = new VariantContextBuilder("UG_call", loc.getContig(), loc.getStart(), loc.getStop(), myAlleles);
|
||||
builder.log10PError(phredScaledConfidence/-10.0);
|
||||
if ( ! passesCallThreshold(phredScaledConfidence) )
|
||||
builder.filters(filter);
|
||||
|
||||
// create the genotypes
|
||||
final GenotypesContext genotypes = afcm.get().subsetAlleles(vc, myAlleles, true,ploidy);
|
||||
builder.genotypes(genotypes);
|
||||
|
||||
// print out stats if we have a writer
|
||||
if ( verboseWriter != null && !limitedContext )
|
||||
printVerboseData(refContext.getLocus().toString(), vc, PoFGT0, phredScaledConfidence, model);
|
||||
|
||||
// *** note that calculating strand bias involves overwriting data structures, so we do that last
|
||||
final HashMap<String, Object> attributes = new HashMap<String, Object>();
|
||||
|
||||
// inherit attributed from input vc if requested
|
||||
if (inheritAttributesFromInputVC)
|
||||
attributes.putAll(vc.getAttributes());
|
||||
// if the site was downsampled, record that fact
|
||||
if ( !limitedContext && rawContext.hasPileupBeenDownsampled() )
|
||||
attributes.put(VCFConstants.DOWNSAMPLED_KEY, true);
|
||||
|
||||
if ( UAC.ANNOTATE_NUMBER_OF_ALLELES_DISCOVERED )
|
||||
attributes.put(NUMBER_OF_DISCOVERED_ALLELES_KEY, vc.getAlternateAlleles().size());
|
||||
|
||||
// add the MLE AC and AF annotations
|
||||
if ( alleleCountsofMLE.size() > 0 ) {
|
||||
attributes.put(VCFConstants.MLE_ALLELE_COUNT_KEY, alleleCountsofMLE);
|
||||
final int AN = builder.make().getCalledChrCount();
|
||||
final ArrayList<Double> MLEfrequencies = new ArrayList<Double>(alleleCountsofMLE.size());
|
||||
// the MLEAC is allowed to be larger than the AN (e.g. in the case of all PLs being 0, the GT is ./. but the exact model may arbitrarily choose an AC>1)
|
||||
for ( int AC : alleleCountsofMLE )
|
||||
MLEfrequencies.add(Math.min(1.0, (double)AC / (double)AN));
|
||||
attributes.put(VCFConstants.MLE_ALLELE_FREQUENCY_KEY, MLEfrequencies);
|
||||
}
|
||||
|
||||
if ( UAC.COMPUTE_SLOD && !limitedContext && !bestGuessIsRef ) {
|
||||
//final boolean DEBUG_SLOD = false;
|
||||
|
||||
// the overall lod
|
||||
//double overallLog10PofNull = AFresult.log10AlleleFrequencyPosteriors[0];
|
||||
final double overallLog10PofF = AFresult.getLog10LikelihoodOfAFGT0();
|
||||
//if ( DEBUG_SLOD ) System.out.println("overallLog10PofF=" + overallLog10PofF);
|
||||
|
||||
final List<Allele> allAllelesToUse = builder.make().getAlleles();
|
||||
|
||||
// the forward lod
|
||||
final VariantContext vcForward = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.FORWARD, allAllelesToUse, false, model, perReadAlleleLikelihoodMap);
|
||||
final AFCalcResult forwardAFresult = afcm.get().getLog10PNonRef(vcForward, getAlleleFrequencyPriors(model));
|
||||
//double[] normalizedLog10Posteriors = MathUtils.normalizeFromLog10(AFresult.log10AlleleFrequencyPosteriors, true);
|
||||
final double forwardLog10PofNull = forwardAFresult.getLog10LikelihoodOfAFEq0();
|
||||
final double forwardLog10PofF = forwardAFresult.getLog10LikelihoodOfAFGT0();
|
||||
//if ( DEBUG_SLOD ) System.out.println("forwardLog10PofNull=" + forwardLog10PofNull + ", forwardLog10PofF=" + forwardLog10PofF);
|
||||
|
||||
// the reverse lod
|
||||
final VariantContext vcReverse = calculateLikelihoods(tracker, refContext, stratifiedContexts, AlignmentContextUtils.ReadOrientation.REVERSE, allAllelesToUse, false, model, perReadAlleleLikelihoodMap);
|
||||
final AFCalcResult reverseAFresult = afcm.get().getLog10PNonRef(vcReverse, getAlleleFrequencyPriors(model));
|
||||
//normalizedLog10Posteriors = MathUtils.normalizeFromLog10(AFresult.log10AlleleFrequencyPosteriors, true);
|
||||
final double reverseLog10PofNull = reverseAFresult.getLog10LikelihoodOfAFEq0();
|
||||
final double reverseLog10PofF = reverseAFresult.getLog10LikelihoodOfAFGT0();
|
||||
//if ( DEBUG_SLOD ) System.out.println("reverseLog10PofNull=" + reverseLog10PofNull + ", reverseLog10PofF=" + reverseLog10PofF);
|
||||
|
||||
final double forwardLod = forwardLog10PofF + reverseLog10PofNull - overallLog10PofF;
|
||||
final double reverseLod = reverseLog10PofF + forwardLog10PofNull - overallLog10PofF;
|
||||
//if ( DEBUG_SLOD ) System.out.println("forward lod=" + forwardLod + ", reverse lod=" + reverseLod);
|
||||
|
||||
// strand score is max bias between forward and reverse strands
|
||||
double strandScore = Math.max(forwardLod, reverseLod);
|
||||
// rescale by a factor of 10
|
||||
strandScore *= 10.0;
|
||||
//logger.debug(String.format("SLOD=%f", strandScore));
|
||||
|
||||
if ( !Double.isNaN(strandScore) )
|
||||
attributes.put("SB", strandScore);
|
||||
}
|
||||
|
||||
// finish constructing the resulting VC
|
||||
builder.attributes(attributes);
|
||||
VariantContext vcCall = builder.make();
|
||||
|
||||
if ( annotationEngine != null && !limitedContext ) { // limitedContext callers need to handle annotations on their own by calling their own annotationEngine
|
||||
// Note: we want to use the *unfiltered* and *unBAQed* context for the annotations
|
||||
final ReadBackedPileup pileup = rawContext.getBasePileup();
|
||||
stratifiedContexts = AlignmentContextUtils.splitContextBySampleName(pileup);
|
||||
|
||||
vcCall = annotationEngine.annotateContext(tracker, refContext, stratifiedContexts, vcCall, perReadAlleleLikelihoodMap);
|
||||
}
|
||||
|
||||
// if we are subsetting alleles (either because there were too many or because some were not polymorphic)
|
||||
// then we may need to trim the alleles (because the original VariantContext may have had to pad at the end).
|
||||
if ( myAlleles.size() != vc.getAlleles().size() && !limitedContext ) // limitedContext callers need to handle allele trimming on their own to keep their perReadAlleleLikelihoodMap alleles in sync
|
||||
vcCall = GATKVariantContextUtils.reverseTrimAlleles(vcCall);
|
||||
|
||||
return new VariantCallContext(vcCall, confidentlyCalled(phredScaledConfidence, PoFGT0));
|
||||
}
|
||||
|
||||
/**
|
||||
* Determine whether input VC to calculateGenotypes() can be genotyped and AF can be computed.
|
||||
* @param vc Input VC
|
||||
* @return Status check
|
||||
*/
|
||||
@Requires("vc != null")
|
||||
protected boolean canVCbeGenotyped(final VariantContext vc) {
|
||||
// protect against too many alternate alleles that we can't even run AF on:
|
||||
if (vc.getNAlleles()> GenotypeLikelihoods.MAX_ALT_ALLELES_THAT_CAN_BE_GENOTYPED) {
|
||||
logger.warn("Attempting to genotype more than "+GenotypeLikelihoods.MAX_ALT_ALLELES_THAT_CAN_BE_GENOTYPED +
|
||||
" alleles. Site will be skipped at location "+vc.getChr()+":"+vc.getStart());
|
||||
return false;
|
||||
}
|
||||
else return true;
|
||||
|
||||
}
|
||||
|
||||
private Map<String, AlignmentContext> getFilteredAndStratifiedContexts(UnifiedArgumentCollection UAC, ReferenceContext refContext, AlignmentContext rawContext, final GenotypeLikelihoodsCalculationModel.Model model) {
|
||||
|
||||
if ( !BaseUtils.isRegularBase(refContext.getBase()) )
|
||||
return null;
|
||||
|
||||
Map<String, AlignmentContext> stratifiedContexts = null;
|
||||
|
||||
if ( model.name().contains("INDEL") ) {
|
||||
|
||||
final ReadBackedPileup pileup = rawContext.getBasePileup().getMappingFilteredPileup(UAC.MIN_BASE_QUALTY_SCORE);
|
||||
// don't call when there is no coverage
|
||||
if ( pileup.getNumberOfElements() == 0 && UAC.OutputMode != OUTPUT_MODE.EMIT_ALL_SITES )
|
||||
return null;
|
||||
|
||||
// stratify the AlignmentContext and cut by sample
|
||||
stratifiedContexts = AlignmentContextUtils.splitContextBySampleName(pileup);
|
||||
|
||||
} else if ( model.name().contains("SNP") ) {
|
||||
|
||||
// stratify the AlignmentContext and cut by sample
|
||||
stratifiedContexts = AlignmentContextUtils.splitContextBySampleName(rawContext.getBasePileup());
|
||||
|
||||
if ( !(UAC.OutputMode == OUTPUT_MODE.EMIT_ALL_SITES && UAC.GenotypingMode != GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES) ) {
|
||||
int numDeletions = 0;
|
||||
for ( final PileupElement p : rawContext.getBasePileup() ) {
|
||||
if ( p.isDeletion() )
|
||||
numDeletions += p.getRepresentativeCount();
|
||||
}
|
||||
if ( ((double) numDeletions) / ((double) rawContext.getBasePileup().depthOfCoverage()) > UAC.MAX_DELETION_FRACTION ) {
|
||||
return null;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return stratifiedContexts;
|
||||
}
|
||||
|
||||
private final double getRefBinomialProbLog10(final int depth) {
|
||||
return MathUtils.log10BinomialProbability(depth, 0);
|
||||
}
|
||||
|
||||
private VariantCallContext estimateReferenceConfidence(VariantContext vc, Map<String, AlignmentContext> contexts, double theta, boolean ignoreCoveredSamples, double initialPofRef) {
|
||||
if ( contexts == null )
|
||||
return null;
|
||||
|
||||
double log10POfRef = Math.log10(initialPofRef);
|
||||
|
||||
// for each sample that we haven't examined yet
|
||||
for ( String sample : samples ) {
|
||||
final AlignmentContext context = contexts.get(sample);
|
||||
if ( ignoreCoveredSamples && context != null )
|
||||
continue;
|
||||
final int depth = context == null ? 0 : context.getBasePileup().depthOfCoverage();
|
||||
log10POfRef += estimateLog10ReferenceConfidenceForOneSample(depth, theta);
|
||||
}
|
||||
|
||||
return new VariantCallContext(vc, QualityUtils.phredScaleLog10CorrectRate(log10POfRef) >= UAC.STANDARD_CONFIDENCE_FOR_CALLING, false);
|
||||
}
|
||||
|
||||
/**
|
||||
* Compute the log10 probability of a sample with sequencing depth and no alt allele is actually truly homozygous reference
|
||||
*
|
||||
* Assumes the sample is diploid
|
||||
*
|
||||
* @param depth the depth of the sample
|
||||
* @param theta the heterozygosity of this species (between 0 and 1)
|
||||
* @return a valid log10 probability of the sample being hom-ref
|
||||
*/
|
||||
@Requires({"depth >= 0", "theta >= 0.0 && theta <= 1.0"})
|
||||
@Ensures("MathUtils.goodLog10Probability(result)")
|
||||
protected double estimateLog10ReferenceConfidenceForOneSample(final int depth, final double theta) {
|
||||
final double log10PofNonRef = Math.log10(theta / 2.0) + getRefBinomialProbLog10(depth);
|
||||
return MathUtils.log10OneMinusX(Math.pow(10.0, log10PofNonRef));
|
||||
}
|
||||
|
||||
protected void printVerboseData(String pos, VariantContext vc, double PofF, double phredScaledConfidence, final GenotypeLikelihoodsCalculationModel.Model model) {
|
||||
Allele refAllele = null, altAllele = null;
|
||||
for ( Allele allele : vc.getAlleles() ) {
|
||||
if ( allele.isReference() )
|
||||
refAllele = allele;
|
||||
else
|
||||
altAllele = allele;
|
||||
}
|
||||
|
||||
for (int i = 0; i <= N; i++) {
|
||||
StringBuilder AFline = new StringBuilder("AFINFO\t");
|
||||
AFline.append(pos);
|
||||
AFline.append("\t");
|
||||
AFline.append(refAllele);
|
||||
AFline.append("\t");
|
||||
if ( altAllele != null )
|
||||
AFline.append(altAllele);
|
||||
else
|
||||
AFline.append("N/A");
|
||||
AFline.append("\t");
|
||||
AFline.append(i + "/" + N + "\t");
|
||||
AFline.append(String.format("%.2f\t", ((float)i)/N));
|
||||
AFline.append(String.format("%.8f\t", getAlleleFrequencyPriors(model)[i]));
|
||||
verboseWriter.println(AFline.toString());
|
||||
}
|
||||
|
||||
verboseWriter.println("P(f>0) = " + PofF);
|
||||
verboseWriter.println("Qscore = " + phredScaledConfidence);
|
||||
verboseWriter.println();
|
||||
}
|
||||
|
||||
protected boolean passesEmitThreshold(double conf, boolean bestGuessIsRef) {
|
||||
return (UAC.OutputMode == OUTPUT_MODE.EMIT_ALL_CONFIDENT_SITES || !bestGuessIsRef) && conf >= Math.min(UAC.STANDARD_CONFIDENCE_FOR_CALLING, UAC.STANDARD_CONFIDENCE_FOR_EMITTING);
|
||||
}
|
||||
|
||||
protected boolean passesCallThreshold(double conf) {
|
||||
return conf >= UAC.STANDARD_CONFIDENCE_FOR_CALLING;
|
||||
}
|
||||
|
||||
protected boolean confidentlyCalled(double conf, double PofF) {
|
||||
return conf >= UAC.STANDARD_CONFIDENCE_FOR_CALLING ||
|
||||
(UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES && QualityUtils.phredScaleErrorRate(PofF) >= UAC.STANDARD_CONFIDENCE_FOR_CALLING);
|
||||
}
|
||||
|
||||
private void determineGLModelsToUse() {
|
||||
String modelPrefix = "";
|
||||
if ( !UAC.GLmodel.name().contains(GPSTRING) && UAC.samplePloidy != GATKVariantContextUtils.DEFAULT_PLOIDY )
|
||||
modelPrefix = GPSTRING;
|
||||
|
||||
// GGA mode => must initialize both the SNP and indel models
|
||||
if ( UAC.GenotypingMode == GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES ||
|
||||
UAC.GLmodel.name().toUpperCase().contains("BOTH") ) {
|
||||
modelsToUse.add(GenotypeLikelihoodsCalculationModel.Model.valueOf(modelPrefix+"SNP"));
|
||||
modelsToUse.add(GenotypeLikelihoodsCalculationModel.Model.valueOf(modelPrefix+"INDEL"));
|
||||
}
|
||||
else {
|
||||
modelsToUse.add(GenotypeLikelihoodsCalculationModel.Model.valueOf(modelPrefix+UAC.GLmodel.name().toUpperCase()));
|
||||
}
|
||||
}
|
||||
|
||||
// decide whether we are currently processing SNPs, indels, neither, or both
|
||||
private List<GenotypeLikelihoodsCalculationModel.Model> getGLModelsToUse(final RefMetaDataTracker tracker,
|
||||
final ReferenceContext refContext,
|
||||
final AlignmentContext rawContext) {
|
||||
if ( UAC.GenotypingMode != GenotypeLikelihoodsCalculationModel.GENOTYPING_MODE.GENOTYPE_GIVEN_ALLELES )
|
||||
return modelsToUse;
|
||||
|
||||
if ( modelsToUse.size() != 2 )
|
||||
throw new IllegalStateException("GGA mode assumes that we have initialized both the SNP and indel models but found " + modelsToUse);
|
||||
|
||||
// if we're genotyping given alleles then we need to choose the model corresponding to the variant type requested
|
||||
final VariantContext vcInput = getVCFromAllelesRod(tracker, refContext, rawContext.getLocation(), false, logger, UAC.alleles);
|
||||
|
||||
if ( vcInput == null ) {
|
||||
return Collections.emptyList(); // no work to be done
|
||||
} else if ( vcInput.isSNP() ) {
|
||||
return Collections.singletonList(modelsToUse.get(SNP_MODEL));
|
||||
} else if ( vcInput.isIndel() || vcInput.isMixed() ) {
|
||||
return Collections.singletonList(modelsToUse.get(INDEL_MODEL));
|
||||
} else {
|
||||
return Collections.emptyList(); // No support for other types yet
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Function that fills vector with allele frequency priors. By default, infinite-sites, neutral variation prior is used,
|
||||
* where Pr(AC=i) = theta/i where theta is heterozygosity
|
||||
* @param N Number of chromosomes
|
||||
* @param priors (output) array to be filled with priors
|
||||
* @param heterozygosity default heterozygosity to use, if inputPriors is empty
|
||||
* @param inputPriors Input priors to use (in which case heterozygosity is ignored)
|
||||
*/
|
||||
public static void computeAlleleFrequencyPriors(final int N, final double[] priors, final double heterozygosity, final List<Double> inputPriors) {
|
||||
|
||||
|
||||
double sum = 0.0;
|
||||
|
||||
if (!inputPriors.isEmpty()) {
|
||||
// user-specified priors
|
||||
if (inputPriors.size() != N)
|
||||
throw new UserException.BadArgumentValue("inputPrior","Invalid length of inputPrior vector: vector length must be equal to # samples +1 ");
|
||||
|
||||
int idx = 1;
|
||||
for (final double prior: inputPriors) {
|
||||
if (prior < 0.0)
|
||||
throw new UserException.BadArgumentValue("Bad argument: negative values not allowed","inputPrior");
|
||||
priors[idx++] = Math.log10(prior);
|
||||
sum += prior;
|
||||
}
|
||||
}
|
||||
else {
|
||||
// for each i
|
||||
for (int i = 1; i <= N; i++) {
|
||||
final double value = heterozygosity / (double)i;
|
||||
priors[i] = Math.log10(value);
|
||||
sum += value;
|
||||
}
|
||||
}
|
||||
|
||||
// protection against the case of heterozygosity too high or an excessive number of samples (which break population genetics assumptions)
|
||||
if (sum > 1.0) {
|
||||
throw new UserException.BadArgumentValue("heterozygosity","The heterozygosity value is set too high relative to the number of samples to be processed, or invalid values specified if input priors were provided - try reducing heterozygosity value or correct input priors.");
|
||||
}
|
||||
// null frequency for AF=0 is (1 - sum(all other frequencies))
|
||||
priors[0] = Math.log10(1.0 - sum);
|
||||
}
|
||||
|
||||
protected double[] getAlleleFrequencyPriors( final GenotypeLikelihoodsCalculationModel.Model model ) {
|
||||
if (model.name().toUpperCase().contains("SNP"))
|
||||
return log10AlleleFrequencyPriorsSNPs;
|
||||
else if (model.name().toUpperCase().contains("INDEL"))
|
||||
return log10AlleleFrequencyPriorsIndels;
|
||||
else
|
||||
throw new IllegalArgumentException("Unexpected GenotypeCalculationModel " + model);
|
||||
|
||||
}
|
||||
|
||||
protected double getTheta( final GenotypeLikelihoodsCalculationModel.Model model ) {
|
||||
if( model.name().contains("SNP") )
|
||||
return HUMAN_SNP_HETEROZYGOSITY;
|
||||
if( model.name().contains("INDEL") )
|
||||
return HUMAN_INDEL_HETEROZYGOSITY;
|
||||
else throw new IllegalArgumentException("Unexpected GenotypeCalculationModel " + model);
|
||||
}
|
||||
|
||||
private static Map<String, GenotypeLikelihoodsCalculationModel> getGenotypeLikelihoodsCalculationObject(Logger logger, UnifiedArgumentCollection UAC) {
|
||||
|
||||
final Map<String, GenotypeLikelihoodsCalculationModel> glcm = new HashMap<String, GenotypeLikelihoodsCalculationModel>();
|
||||
final List<Class<? extends GenotypeLikelihoodsCalculationModel>> glmClasses = new PluginManager<GenotypeLikelihoodsCalculationModel>(GenotypeLikelihoodsCalculationModel.class).getPlugins();
|
||||
|
||||
for (int i = 0; i < glmClasses.size(); i++) {
|
||||
final Class<? extends GenotypeLikelihoodsCalculationModel> glmClass = glmClasses.get(i);
|
||||
final String key = glmClass.getSimpleName().replaceAll("GenotypeLikelihoodsCalculationModel","").toUpperCase();
|
||||
try {
|
||||
final Object args[] = new Object[]{UAC,logger};
|
||||
final Constructor c = glmClass.getDeclaredConstructor(UnifiedArgumentCollection.class, Logger.class);
|
||||
glcm.put(key, (GenotypeLikelihoodsCalculationModel)c.newInstance(args));
|
||||
}
|
||||
catch (Exception e) {
|
||||
throw new UserException("The likelihoods model provided for the -glm argument (" + UAC.GLmodel + ") is not a valid option: " + e.getMessage());
|
||||
}
|
||||
}
|
||||
|
||||
return glcm;
|
||||
}
|
||||
|
||||
public static VariantContext getVCFromAllelesRod(RefMetaDataTracker tracker, ReferenceContext ref, GenomeLoc loc, boolean requireSNP, Logger logger, final RodBinding<VariantContext> allelesBinding) {
|
||||
if ( tracker == null || ref == null || logger == null )
|
||||
return null;
|
||||
VariantContext vc = null;
|
||||
|
||||
// search for usable record
|
||||
for ( final VariantContext vc_input : tracker.getValues(allelesBinding, loc) ) {
|
||||
if ( vc_input != null && ! vc_input.isFiltered() && (! requireSNP || vc_input.isSNP() )) {
|
||||
if ( vc == null ) {
|
||||
vc = vc_input;
|
||||
} else {
|
||||
logger.warn("Multiple valid VCF records detected in the alleles input file at site " + ref.getLocus() + ", only considering the first record");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return vc;
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,333 @@
|
|||
/*
|
||||
* By downloading the PROGRAM you agree to the following terms of use:
|
||||
*
|
||||
* BROAD INSTITUTE - SOFTWARE LICENSE AGREEMENT - FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
|
||||
*
|
||||
* This Agreement is made between the Broad Institute, Inc. with a principal address at 7 Cambridge Center, Cambridge, MA 02142 (BROAD) and the LICENSEE and is effective at the date the downloading is completed (EFFECTIVE DATE).
|
||||
*
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
|
||||
* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
|
||||
* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
|
||||
*
|
||||
* 1. DEFINITIONS
|
||||
* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK2 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute/GATK on the EFFECTIVE DATE.
|
||||
*
|
||||
* 2. LICENSE
|
||||
* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM.
|
||||
* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
|
||||
* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
|
||||
* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
|
||||
*
|
||||
* 3. OWNERSHIP OF INTELLECTUAL PROPERTY
|
||||
* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
|
||||
* Copyright 2012 Broad Institute, Inc.
|
||||
* Notice of attribution: The GATK2 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
|
||||
* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
|
||||
*
|
||||
* 4. INDEMNIFICATION
|
||||
* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
|
||||
*
|
||||
* 5. NO REPRESENTATIONS OR WARRANTIES
|
||||
* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
|
||||
* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
|
||||
*
|
||||
* 6. ASSIGNMENT
|
||||
* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
|
||||
*
|
||||
* 7. MISCELLANEOUS
|
||||
* 7.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
|
||||
* 7.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
|
||||
* 7.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
|
||||
* 7.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
|
||||
* 7.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
|
||||
* 7.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
|
||||
* 7.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
|
||||
*/
|
||||
|
||||
package org.broadinstitute.sting.gatk.walkers.genotyper.afcalc;
|
||||
|
||||
import org.broadinstitute.sting.utils.MathUtils;
|
||||
import org.broadinstitute.sting.utils.variant.GATKVariantContextUtils;
|
||||
import org.broadinstitute.variant.variantcontext.Allele;
|
||||
import org.broadinstitute.variant.variantcontext.GenotypeLikelihoods;
|
||||
import org.broadinstitute.variant.variantcontext.GenotypesContext;
|
||||
import org.broadinstitute.variant.variantcontext.VariantContext;
|
||||
|
||||
import java.util.*;
|
||||
|
||||
public abstract class DiploidExactAFCalc extends ExactAFCalc {
|
||||
public DiploidExactAFCalc(final int nSamples, final int maxAltAlleles, final int ploidy) {
|
||||
super(nSamples, maxAltAlleles, ploidy);
|
||||
if ( ploidy != 2 ) throw new IllegalArgumentException("ploidy must be two for DiploidExactAFCalc and subclasses but saw " + ploidy);
|
||||
}
|
||||
|
||||
@Override
|
||||
protected AFCalcResult computeLog10PNonRef(final VariantContext vc,
|
||||
final double[] log10AlleleFrequencyPriors) {
|
||||
final int numAlternateAlleles = vc.getNAlleles() - 1;
|
||||
final ArrayList<double[]> genotypeLikelihoods = getGLs(vc.getGenotypes(), true);
|
||||
final int numSamples = genotypeLikelihoods.size()-1;
|
||||
final int numChr = 2*numSamples;
|
||||
|
||||
// queue of AC conformations to process
|
||||
final LinkedList<ExactACset> ACqueue = new LinkedList<>();
|
||||
|
||||
// mapping of ExactACset indexes to the objects
|
||||
final HashMap<ExactACcounts, ExactACset> indexesToACset = new HashMap<>(numChr+1);
|
||||
|
||||
// add AC=0 to the queue
|
||||
final int[] zeroCounts = new int[numAlternateAlleles];
|
||||
ExactACset zeroSet = new ExactACset(numSamples+1, new ExactACcounts(zeroCounts));
|
||||
ACqueue.add(zeroSet);
|
||||
indexesToACset.put(zeroSet.getACcounts(), zeroSet);
|
||||
|
||||
while ( !ACqueue.isEmpty() ) {
|
||||
getStateTracker().incNEvaluations(); // keep track of the number of evaluations
|
||||
|
||||
// compute log10Likelihoods
|
||||
final ExactACset set = ACqueue.remove();
|
||||
|
||||
calculateAlleleCountConformation(set, genotypeLikelihoods, numChr, ACqueue, indexesToACset, log10AlleleFrequencyPriors);
|
||||
|
||||
// clean up memory
|
||||
indexesToACset.remove(set.getACcounts());
|
||||
//if ( DEBUG )
|
||||
// System.out.printf(" *** removing used set=%s%n", set.ACcounts);
|
||||
}
|
||||
|
||||
return getResultFromFinalState(vc, log10AlleleFrequencyPriors);
|
||||
}
|
||||
|
||||
|
||||
@Override
|
||||
protected GenotypesContext reduceScopeGenotypes(final VariantContext vc, final List<Allele> allelesToUse) {
|
||||
return GATKVariantContextUtils.subsetDiploidAlleles(vc, allelesToUse, GATKVariantContextUtils.GenotypeAssignmentMethod.SET_TO_NO_CALL);
|
||||
}
|
||||
|
||||
@Override
|
||||
protected void reduceScopeCalculateLikelihoodSums(final VariantContext vc, final LikelihoodSum[] likelihoodSums) {
|
||||
final ArrayList<double[]> GLs = getGLs(vc.getGenotypes(), true);
|
||||
for ( final double[] likelihoods : GLs ) {
|
||||
final int PLindexOfBestGL = MathUtils.maxElementIndex(likelihoods);
|
||||
if ( PLindexOfBestGL != PL_INDEX_OF_HOM_REF ) {
|
||||
final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair alleles = GenotypeLikelihoods.getAllelePair(PLindexOfBestGL);
|
||||
final int alleleLikelihoodIndex1 = alleles.alleleIndex1 - 1;
|
||||
final int alleleLikelihoodIndex2 = alleles.alleleIndex2 - 1;
|
||||
if ( alleles.alleleIndex1 != 0 )
|
||||
likelihoodSums[alleleLikelihoodIndex1].sum += likelihoods[PLindexOfBestGL] - likelihoods[PL_INDEX_OF_HOM_REF];
|
||||
// don't double-count it
|
||||
if ( alleles.alleleIndex2 != 0 && alleles.alleleIndex2 != alleles.alleleIndex1 )
|
||||
likelihoodSums[alleleLikelihoodIndex2].sum += likelihoods[PLindexOfBestGL] - likelihoods[PL_INDEX_OF_HOM_REF];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private static final class DependentSet {
|
||||
public final int[] ACcounts;
|
||||
public final int PLindex;
|
||||
|
||||
public DependentSet(final int[] ACcounts, final int PLindex) {
|
||||
this.ACcounts = ACcounts;
|
||||
this.PLindex = PLindex;
|
||||
}
|
||||
}
|
||||
|
||||
private double calculateAlleleCountConformation(final ExactACset set,
|
||||
final ArrayList<double[]> genotypeLikelihoods,
|
||||
final int numChr,
|
||||
final LinkedList<ExactACset> ACqueue,
|
||||
final HashMap<ExactACcounts, ExactACset> indexesToACset,
|
||||
final double[] log10AlleleFrequencyPriors) {
|
||||
|
||||
//if ( DEBUG )
|
||||
// System.out.printf(" *** computing LofK for set=%s%n", set.ACcounts);
|
||||
|
||||
// compute the log10Likelihoods
|
||||
computeLofK(set, genotypeLikelihoods, log10AlleleFrequencyPriors);
|
||||
|
||||
final double log10LofK = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1];
|
||||
|
||||
// can we abort early because the log10Likelihoods are so small?
|
||||
if ( getStateTracker().abort(log10LofK, set.getACcounts(), true) ) {
|
||||
//if ( DEBUG )
|
||||
// System.out.printf(" *** breaking early set=%s log10L=%.2f maxLog10L=%.2f%n", set.ACcounts, log10LofK, maxLog10L);
|
||||
return log10LofK;
|
||||
}
|
||||
|
||||
// iterate over higher frequencies if possible
|
||||
final int ACwiggle = numChr - set.getACsum();
|
||||
if ( ACwiggle == 0 ) // all alternate alleles already sum to 2N so we cannot possibly go to higher frequencies
|
||||
return log10LofK;
|
||||
|
||||
final int numAltAlleles = set.getACcounts().getCounts().length;
|
||||
|
||||
// add conformations for the k+1 case
|
||||
for ( int allele = 0; allele < numAltAlleles; allele++ ) {
|
||||
final int[] ACcountsClone = set.getACcounts().getCounts().clone();
|
||||
ACcountsClone[allele]++;
|
||||
// to get to this conformation, a sample would need to be AB (remember that ref=0)
|
||||
final int PLindex = GenotypeLikelihoods.calculatePLindex(0, allele+1);
|
||||
updateACset(ACcountsClone, numChr, set, PLindex, ACqueue, indexesToACset, genotypeLikelihoods);
|
||||
}
|
||||
|
||||
// add conformations for the k+2 case if it makes sense; note that the 2 new alleles may be the same or different
|
||||
if ( ACwiggle > 1 ) {
|
||||
final ArrayList<DependentSet> differentAlleles = new ArrayList<>(numAltAlleles * numAltAlleles);
|
||||
final ArrayList<DependentSet> sameAlleles = new ArrayList<>(numAltAlleles);
|
||||
|
||||
for ( int allele_i = 0; allele_i < numAltAlleles; allele_i++ ) {
|
||||
for ( int allele_j = allele_i; allele_j < numAltAlleles; allele_j++ ) {
|
||||
final int[] ACcountsClone = set.getACcounts().getCounts().clone();
|
||||
ACcountsClone[allele_i]++;
|
||||
ACcountsClone[allele_j]++;
|
||||
|
||||
// to get to this conformation, a sample would need to be BB or BC (remember that ref=0, so add one to the index)
|
||||
final int PLindex = GenotypeLikelihoods.calculatePLindex(allele_i+1, allele_j+1);
|
||||
if ( allele_i == allele_j )
|
||||
sameAlleles.add(new DependentSet(ACcountsClone, PLindex));
|
||||
else
|
||||
differentAlleles.add(new DependentSet(ACcountsClone, PLindex));
|
||||
}
|
||||
}
|
||||
|
||||
// IMPORTANT: we must first add the cases where the 2 new alleles are different so that the queue maintains its ordering
|
||||
for ( DependentSet dependent : differentAlleles )
|
||||
updateACset(dependent.ACcounts, numChr, set, dependent.PLindex, ACqueue, indexesToACset, genotypeLikelihoods);
|
||||
for ( DependentSet dependent : sameAlleles )
|
||||
updateACset(dependent.ACcounts, numChr, set, dependent.PLindex, ACqueue, indexesToACset, genotypeLikelihoods);
|
||||
}
|
||||
|
||||
return log10LofK;
|
||||
}
|
||||
|
||||
// adds the ExactACset represented by the ACcounts to the ACqueue if not already there (creating it if needed) and
|
||||
// also pushes its value to the given callingSetIndex.
|
||||
private void updateACset(final int[] newSetCounts,
|
||||
final int numChr,
|
||||
final ExactACset dependentSet,
|
||||
final int PLsetIndex,
|
||||
final Queue<ExactACset> ACqueue,
|
||||
final HashMap<ExactACcounts, ExactACset> indexesToACset,
|
||||
final ArrayList<double[]> genotypeLikelihoods) {
|
||||
final ExactACcounts index = new ExactACcounts(newSetCounts);
|
||||
if ( !indexesToACset.containsKey(index) ) {
|
||||
ExactACset set = new ExactACset(numChr/2 +1, index);
|
||||
indexesToACset.put(index, set);
|
||||
ACqueue.add(set);
|
||||
}
|
||||
|
||||
// push data from the dependency to the new set
|
||||
//if ( DEBUG )
|
||||
// System.out.println(" *** pushing data from " + index + " to " + dependencySet.ACcounts);
|
||||
pushData(indexesToACset.get(index), dependentSet, PLsetIndex, genotypeLikelihoods);
|
||||
}
|
||||
|
||||
private void computeLofK(final ExactACset set,
|
||||
final ArrayList<double[]> genotypeLikelihoods,
|
||||
final double[] log10AlleleFrequencyPriors) {
|
||||
|
||||
set.getLog10Likelihoods()[0] = 0.0; // the zero case
|
||||
final int totalK = set.getACsum();
|
||||
|
||||
// special case for k = 0 over all k
|
||||
if ( totalK == 0 ) {
|
||||
for ( int j = 1; j < set.getLog10Likelihoods().length; j++ )
|
||||
set.getLog10Likelihoods()[j] = set.getLog10Likelihoods()[j-1] + genotypeLikelihoods.get(j)[HOM_REF_INDEX];
|
||||
|
||||
final double log10Lof0 = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1];
|
||||
getStateTracker().setLog10LikelihoodOfAFzero(log10Lof0);
|
||||
getStateTracker().setLog10PosteriorOfAFzero(log10Lof0 + log10AlleleFrequencyPriors[0]);
|
||||
return;
|
||||
}
|
||||
|
||||
// if we got here, then k > 0 for at least one k.
|
||||
// the non-AA possible conformations were already dealt with by pushes from dependent sets;
|
||||
// now deal with the AA case (which depends on previous cells in this column) and then update the L(j,k) value
|
||||
for ( int j = 1; j < set.getLog10Likelihoods().length; j++ ) {
|
||||
|
||||
if ( totalK < 2*j-1 ) {
|
||||
final double[] gl = genotypeLikelihoods.get(j);
|
||||
final double conformationValue = MathUtils.log10Cache[2*j-totalK] + MathUtils.log10Cache[2*j-totalK-1] + set.getLog10Likelihoods()[j-1] + gl[HOM_REF_INDEX];
|
||||
set.getLog10Likelihoods()[j] = MathUtils.approximateLog10SumLog10(set.getLog10Likelihoods()[j], conformationValue);
|
||||
}
|
||||
|
||||
final double logDenominator = MathUtils.log10Cache[2*j] + MathUtils.log10Cache[2*j-1];
|
||||
set.getLog10Likelihoods()[j] = set.getLog10Likelihoods()[j] - logDenominator;
|
||||
}
|
||||
|
||||
double log10LofK = set.getLog10Likelihoods()[set.getLog10Likelihoods().length-1];
|
||||
|
||||
// update the MLE if necessary
|
||||
getStateTracker().updateMLEifNeeded(log10LofK, set.getACcounts().getCounts());
|
||||
|
||||
// apply the priors over each alternate allele
|
||||
for ( final int ACcount : set.getACcounts().getCounts() ) {
|
||||
if ( ACcount > 0 )
|
||||
log10LofK += log10AlleleFrequencyPriors[ACcount];
|
||||
}
|
||||
|
||||
getStateTracker().updateMAPifNeeded(log10LofK, set.getACcounts().getCounts());
|
||||
}
|
||||
|
||||
private void pushData(final ExactACset targetSet,
|
||||
final ExactACset dependentSet,
|
||||
final int PLsetIndex,
|
||||
final ArrayList<double[]> genotypeLikelihoods) {
|
||||
final int totalK = targetSet.getACsum();
|
||||
|
||||
for ( int j = 1; j < targetSet.getLog10Likelihoods().length; j++ ) {
|
||||
|
||||
if ( totalK <= 2*j ) { // skip impossible conformations
|
||||
final double[] gl = genotypeLikelihoods.get(j);
|
||||
final double conformationValue =
|
||||
determineCoefficient(PLsetIndex, j, targetSet.getACcounts().getCounts(), totalK) + dependentSet.getLog10Likelihoods()[j-1] + gl[PLsetIndex];
|
||||
targetSet.getLog10Likelihoods()[j] = MathUtils.approximateLog10SumLog10(targetSet.getLog10Likelihoods()[j], conformationValue);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
private double determineCoefficient(int PLindex, final int j, final int[] ACcounts, final int totalK) {
|
||||
// the closed form representation generalized for multiple alleles is as follows:
|
||||
// AA: (2j - totalK) * (2j - totalK - 1)
|
||||
// AB: 2k_b * (2j - totalK)
|
||||
// AC: 2k_c * (2j - totalK)
|
||||
// BB: k_b * (k_b - 1)
|
||||
// BC: 2 * k_b * k_c
|
||||
// CC: k_c * (k_c - 1)
|
||||
|
||||
// find the 2 alleles that are represented by this PL index
|
||||
GenotypeLikelihoods.GenotypeLikelihoodsAllelePair alleles = GenotypeLikelihoods.getAllelePair(PLindex);
|
||||
|
||||
// *** note that throughout this method we subtract one from the alleleIndex because ACcounts ***
|
||||
// *** doesn't consider the reference allele whereas the GenotypeLikelihoods PL cache does. ***
|
||||
|
||||
// the AX het case
|
||||
if ( alleles.alleleIndex1 == 0 )
|
||||
return MathUtils.log10Cache[2*ACcounts[alleles.alleleIndex2-1]] + MathUtils.log10Cache[2*j-totalK];
|
||||
|
||||
final int k_i = ACcounts[alleles.alleleIndex1-1];
|
||||
|
||||
// the hom var case (e.g. BB, CC, DD)
|
||||
final double coeff;
|
||||
if ( alleles.alleleIndex1 == alleles.alleleIndex2 ) {
|
||||
coeff = MathUtils.log10Cache[k_i] + MathUtils.log10Cache[k_i - 1];
|
||||
}
|
||||
// the het non-ref case (e.g. BC, BD, CD)
|
||||
else {
|
||||
final int k_j = ACcounts[alleles.alleleIndex2-1];
|
||||
coeff = MathUtils.log10Cache[2] + MathUtils.log10Cache[k_i] + MathUtils.log10Cache[k_j];
|
||||
}
|
||||
|
||||
return coeff;
|
||||
}
|
||||
|
||||
public GenotypesContext subsetAlleles(final VariantContext vc,
|
||||
final List<Allele> allelesToUse,
|
||||
final boolean assignGenotypes,
|
||||
final int ploidy) {
|
||||
return allelesToUse.size() == 1
|
||||
? GATKVariantContextUtils.subsetToRefOnly(vc, ploidy)
|
||||
: GATKVariantContextUtils.subsetDiploidAlleles(vc, allelesToUse,
|
||||
assignGenotypes ? GATKVariantContextUtils.GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN : GATKVariantContextUtils.GenotypeAssignmentMethod.SET_TO_NO_CALL);
|
||||
}
|
||||
}
|
||||
Some files were not shown because too many files have changed in this diff Show More
Loading…
Reference in New Issue