* To download and extract the data sources, you can invoke {@link FuncotatorDataSourceDownloader} in the following ways: *
{@code ./gatk FuncotatorDataSourceDownloader --somatic --validate-integrity --extract-after-download}{@code ./gatk FuncotatorDataSourceDownloader --germline --validate-integrity --extract-after-download}{@code ./gatk FuncotatorDataSourceDownloader --somatic --validate-integrity --hg38 --extract-after-download}{@code ./gatk FuncotatorDataSourceDownloader --germline --validate-integrity --hg19 --extract-after-download}
* sample1 sample1.vcf.gz
* sample2 sample2.vcf.gz sample2.vcf.gz.tbi
@@ -218,12 +216,14 @@ public final class GenomicsDBImport extends GATKTool {
public static final String MAX_NUM_INTERVALS_TO_IMPORT_IN_PARALLEL = "max-num-intervals-to-import-in-parallel";
public static final String MERGE_CONTIGS_INTO_NUM_PARTITIONS = "merge-contigs-into-num-partitions";
public static final String BYPASS_FEATURE_READER = "bypass-feature-reader";
+ public static final String VCF_HEADER_OVERRIDE = "header";
public static final int INTERVAL_LIST_SIZE_WARNING_THRESHOLD = 100;
public static final int ARRAY_COLUMN_BOUNDS_START = 0;
public static final int ARRAY_COLUMN_BOUNDS_END = 1;
public static final String SHARED_POSIXFS_OPTIMIZATIONS = GenomicsDBArgumentCollection.SHARED_POSIXFS_OPTIMIZATIONS;
public static final String USE_GCS_HDFS_CONNECTOR = GenomicsDBArgumentCollection.USE_GCS_HDFS_CONNECTOR;
+ public static final String AVOID_NIO = "avoid-nio";
@Argument(fullName = WORKSPACE_ARG_LONG_NAME,
doc = "Workspace for GenomicsDB. Can be a POSIX file system absolute or relative path or a HDFS/GCS URL. " +
@@ -239,7 +239,7 @@ public final class GenomicsDBImport extends GATKTool {
"when using the "+INTERVAL_LIST_LONG_NAME+" option. " +
"Either this or "+WORKSPACE_ARG_LONG_NAME+" must be specified. " +
"Must point to an existing workspace.",
- mutex = {WORKSPACE_ARG_LONG_NAME})
+ mutex = {WORKSPACE_ARG_LONG_NAME, VCF_HEADER_OVERRIDE})
private String incrementalImportWorkspace;
@Argument(fullName = SEGMENT_SIZE_ARG_LONG_NAME,
@@ -254,7 +254,7 @@ public final class GenomicsDBImport extends GATKTool {
" data for only a single sample. Either this or " + SAMPLE_NAME_MAP_LONG_NAME +
" must be specified.",
optional = true,
- mutex = {SAMPLE_NAME_MAP_LONG_NAME})
+ mutex = {SAMPLE_NAME_MAP_LONG_NAME, AVOID_NIO})
private List variantPaths;
@Argument(fullName = VCF_BUFFER_SIZE_ARG_NAME,
@@ -364,6 +364,13 @@ public final class GenomicsDBImport extends GATKTool {
optional = true)
private boolean sharedPosixFSOptimizations = false;
+ @Argument(fullName = VCF_HEADER_OVERRIDE,
+ doc = "Specify a vcf file to use instead of reading and combining headers from the input vcfs",
+ optional = true,
+ mutex ={INCREMENTAL_WORKSPACE_ARG_LONG_NAME}
+ )
+ private FeatureInput headerOverride = null;
+
@Argument(fullName = BYPASS_FEATURE_READER,
doc = "Use htslib to read input VCFs instead of GATK's FeatureReader. This will reduce memory usage and potentially speed up " +
"the import. Lower memory requirements may also enable parallelism through " + MAX_NUM_INTERVALS_TO_IMPORT_IN_PARALLEL +
@@ -371,6 +378,15 @@ public final class GenomicsDBImport extends GATKTool {
optional = true)
private boolean bypassFeatureReader = false;
+ @Argument(fullName = AVOID_NIO,
+ doc = "Do not attempt to open the input vcf file paths in java. This can only be used with " + BYPASS_FEATURE_READER
+ + ". It allows operating on file systems which GenomicsDB understands how to open but GATK does not. This will disable "
+ + "many of the sanity checks.",
+ mutex = {StandardArgumentDefinitions.VARIANT_LONG_NAME}
+ )
+ @Advanced
+ private boolean avoidNio = false;
+
@Argument(fullName = USE_GCS_HDFS_CONNECTOR,
doc = "Use the GCS HDFS Connector instead of the native GCS SDK client with gs:// URLs.",
optional = true)
@@ -440,10 +456,6 @@ public int getDefaultCloudIndexPrefetchBufferSize() {
// Path to combined VCF header file to be written by GenomicsDBImporter
private String vcfHeaderFile;
- // GenomicsDB callset map protobuf structure containing all callset names
- // used to write the callset json file on traversal success
- private GenomicsDBCallsetsMapProto.CallsetMappingPB callsetMappingPB;
-
//in-progress batchCount
private int batchCount = 1;
@@ -463,11 +475,14 @@ public void onStartup() {
initializeWorkspaceAndToolMode();
assertVariantPathsOrSampleNameFileWasSpecified();
assertOverwriteWorkspaceAndIncrementalImportMutuallyExclusive();
+ assertAvoidNioConditionsAreValid();
initializeHeaderAndSampleMappings();
initializeIntervals();
super.onStartup();
}
+
+
private void initializeWorkspaceAndToolMode() {
if (incrementalImportWorkspace != null && !incrementalImportWorkspace.isEmpty()) {
doIncrementalImport = true;
@@ -495,6 +510,24 @@ private void assertVariantPathsOrSampleNameFileWasSpecified(){
}
}
+ private void assertAvoidNioConditionsAreValid() {
+ if (avoidNio && (!bypassFeatureReader || headerOverride == null) ){
+ final List missing = new ArrayList<>();
+ if(!bypassFeatureReader){
+ missing.add(BYPASS_FEATURE_READER);
+ }
+ if(headerOverride == null){
+ missing.add(VCF_HEADER_OVERRIDE);
+ }
+ final String missingArgs = String.join(" and ", missing);
+
+ // this potentially produces and exception with bad grammar but that's probably ok
+ throw new CommandLineException.MissingArgument(missingArgs, "If --" +AVOID_NIO + " is set then --" + BYPASS_FEATURE_READER
+ + " and --" + VCF_HEADER_OVERRIDE + " must also be specified.");
+
+ }
+ }
+
private static void assertIntervalsCoverEntireContigs(GenomicsDBImporter importer,
List intervals) {
GenomicsDBVidMapProto.VidMappingPB vidMapPB = importer.getProtobufVidMapping();
@@ -523,32 +556,37 @@ private static void assertIntervalsCoverEntireContigs(GenomicsDBImporter importe
*/
private void initializeHeaderAndSampleMappings() {
// Only one of -V and --sampleNameMapFile may be specified
- if (variantPaths != null && variantPaths.size() > 0) {
+ if (variantPaths != null && !variantPaths.isEmpty()) {
// -V was specified
final List headers = new ArrayList<>(variantPaths.size());
sampleNameMap = new SampleNameMap();
- for (final String variantPathString : variantPaths) {
- final Path variantPath = IOUtils.getPath(variantPathString);
- if (bypassFeatureReader) {
- GATKGenomicsDBUtils.assertVariantFileIsCompressedAndIndexed(variantPath);
- }
- final VCFHeader header = getHeaderFromPath(variantPath, null);
- Utils.validate(header != null, "Null header was found in " + variantPath + ".");
- assertGVCFHasOnlyOneSample(variantPathString, header);
- headers.add(header);
- final String sampleName = header.getGenotypeSamples().get(0);
- try {
- sampleNameMap.addSample(sampleName, new URI(variantPathString));
- }
- catch(final URISyntaxException e) {
- throw new UserException("Malformed URI "+e.toString(), e);
+ if(headerOverride == null) {
+ for (final String variantPathString : variantPaths) {
+ final Path variantPath = IOUtils.getPath(variantPathString);
+ if (bypassFeatureReader) { // avoid-nio can't be set here because it requires headerOverride
+ GATKGenomicsDBUtils.assertVariantFileIsCompressedAndIndexed(variantPath);
+ }
+ final VCFHeader header = getHeaderFromPath(variantPath);
+ Utils.validate(header != null, "Null header was found in " + variantPath + ".");
+ assertGVCFHasOnlyOneSample(variantPathString, header);
+ headers.add(header);
+
+ final String sampleName = header.getGenotypeSamples().get(0);
+ try {
+ sampleNameMap.addSample(sampleName, new URI(variantPathString));
+ } catch (final URISyntaxException e) {
+ throw new UserException("Malformed URI " + e.getMessage(), e);
+ }
}
+ mergedHeaderLines = VCFUtils.smartMergeHeaders(headers, true);
+ mergedHeaderSequenceDictionary = new VCFHeader(mergedHeaderLines).getSequenceDictionary();
+ } else {
+ final VCFHeader header = getHeaderFromPath(headerOverride.toPath());
+ mergedHeaderLines = new LinkedHashSet<>(header.getMetaDataInInputOrder());
+ mergedHeaderSequenceDictionary = header.getSequenceDictionary();
}
- mergedHeaderLines = VCFUtils.smartMergeHeaders(headers, true);
- mergedHeaderSequenceDictionary = new VCFHeader(mergedHeaderLines).getSequenceDictionary();
mergedHeaderLines.addAll(getDefaultToolVCFHeaderLines());
-
} else if (sampleNameMapFile != null) {
// --sampleNameMap was specified
@@ -556,31 +594,34 @@ private void initializeHeaderAndSampleMappings() {
//the resulting database will have incorrect sample names
//see https://github.com/broadinstitute/gatk/issues/3682 for more information
// The SampleNameMap class guarantees that the samples will be sorted correctly.
- sampleNameMap = new SampleNameMap(IOUtils.getPath(sampleNameMapFile), bypassFeatureReader);
+ sampleNameMap = new SampleNameMap(IOUtils.getPath(sampleNameMapFile),
+ bypassFeatureReader && !avoidNio);
final String firstSample = sampleNameMap.getSampleNameToVcfPath().entrySet().iterator().next().getKey();
- final Path firstVCFPath = sampleNameMap.getVCFForSampleAsPath(firstSample);
- final Path firstVCFIndexPath = sampleNameMap.getVCFIndexForSampleAsPath(firstSample);
- final VCFHeader header = getHeaderFromPath(firstVCFPath, firstVCFIndexPath);
+ final VCFHeader header;
+ if(headerOverride == null){
+ final Path firstVCFPath = sampleNameMap.getVCFForSampleAsPath(firstSample);
+ header = getHeaderFromPath(firstVCFPath);
+ } else {
+ header = getHeaderFromPath(headerOverride.toPath());
+ }
//getMetaDataInInputOrder() returns an ImmutableSet - LinkedHashSet is mutable and preserves ordering
- mergedHeaderLines = new LinkedHashSet(header.getMetaDataInInputOrder());
+ mergedHeaderLines = new LinkedHashSet<>(header.getMetaDataInInputOrder());
mergedHeaderSequenceDictionary = header.getSequenceDictionary();
mergedHeaderLines.addAll(getDefaultToolVCFHeaderLines());
- }
- else if (getIntervalsFromExistingWorkspace){
+ } else if (getIntervalsFromExistingWorkspace){
final String vcfHeader = IOUtils.appendPathToDir(workspace, GenomicsDBConstants.DEFAULT_VCFHEADER_FILE_NAME);
IOUtils.assertPathsAreReadable(vcfHeader);
final String header = GenomicsDBUtils.readEntireFile(vcfHeader);
try {
File tempHeader = IOUtils.createTempFile("tempheader", ".vcf");
- Files.write(tempHeader.toPath(), header.getBytes(StandardCharsets.UTF_8));
+ Files.writeString(tempHeader.toPath(), header);
mergedHeaderSequenceDictionary = VCFFileReader.getSequenceDictionary(tempHeader);
} catch (final IOException e) {
- throw new UserException("Unable to create temporary header file to get sequence dictionary");
+ throw new UserException("Unable to create temporary header file to get sequence dictionary", e);
}
- }
- else {
+ } else {
throw new UserException(StandardArgumentDefinitions.VARIANT_LONG_NAME+" or "+
SAMPLE_NAME_MAP_LONG_NAME+" must be specified unless "+
INTERVAL_LIST_LONG_NAME+" is specified");
@@ -599,8 +640,12 @@ else if (getIntervalsFromExistingWorkspace){
}
}
- private VCFHeader getHeaderFromPath(final Path variantPath, final Path variantIndexPath) {
- try(final FeatureReader reader = getReaderFromPath(variantPath, variantIndexPath)) {
+ private VCFHeader getHeaderFromPath(final Path variantPath) {
+ //TODO make this mangling unecessary
+ final String variantURI = variantPath.toAbsolutePath().toUri().toString();
+ try(final FeatureReader reader = AbstractFeatureReader.getFeatureReader(variantURI, null, new VCFCodec(), false,
+ BucketUtils.getPrefetchingWrapper(cloudPrefetchBuffer),
+ BucketUtils.getPrefetchingWrapper(cloudIndexPrefetchBuffer))) {
return (VCFHeader) reader.getHeader();
} catch (final IOException e) {
throw new UserException("Error while reading vcf header from " + variantPath.toUri(), e);
@@ -633,9 +678,9 @@ private void writeIntervalListToFile() {
@Override
public void onTraversalStart() {
String workspaceDir = overwriteCreateOrCheckWorkspace();
- vidMapJSONFile = genomicsDBApppendPaths(workspaceDir, GenomicsDBConstants.DEFAULT_VIDMAP_FILE_NAME);
- callsetMapJSONFile = genomicsDBApppendPaths(workspaceDir, GenomicsDBConstants.DEFAULT_CALLSETMAP_FILE_NAME);
- vcfHeaderFile = genomicsDBApppendPaths(workspaceDir, GenomicsDBConstants.DEFAULT_VCFHEADER_FILE_NAME);
+ vidMapJSONFile = GATKGenomicsDBUtils.genomicsDBApppendPaths(workspaceDir, GenomicsDBConstants.DEFAULT_VIDMAP_FILE_NAME);
+ callsetMapJSONFile = GATKGenomicsDBUtils.genomicsDBApppendPaths(workspaceDir, GenomicsDBConstants.DEFAULT_CALLSETMAP_FILE_NAME);
+ vcfHeaderFile = GATKGenomicsDBUtils.genomicsDBApppendPaths(workspaceDir, GenomicsDBConstants.DEFAULT_VCFHEADER_FILE_NAME);
if (getIntervalsFromExistingWorkspace) {
// intervals may be null if merge-contigs-into-num-partitions was used to create the workspace
// if so, we need to wait for vid to be generated before writing out the interval list
@@ -775,7 +820,7 @@ private List generateIntervalListFromWorkspace() {
final int start = Integer.parseInt(partitionInfo[1]);
final int end = Integer.parseInt(partitionInfo[2]);
return new SimpleInterval(contig, start, end);
- }).filter(o -> o != null).collect(Collectors.toList());
+ }).filter(Objects::nonNull).collect(Collectors.toList());
}
private ImportConfig createImportConfig(final int batchSize) {
@@ -785,7 +830,7 @@ private ImportConfig createImportConfig(final int batchSize) {
GenomicsDBImportConfiguration.ImportConfiguration.newBuilder();
importConfigurationBuilder.addAllColumnPartitions(partitions);
importConfigurationBuilder.setSizePerColumnPartition(vcfBufferSizePerSample);
- importConfigurationBuilder.setFailIfUpdating(true && !doIncrementalImport);
+ importConfigurationBuilder.setFailIfUpdating(!doIncrementalImport);
importConfigurationBuilder.setSegmentSize(segmentSize);
importConfigurationBuilder.setConsolidateTiledbArrayAfterLoad(doConsolidation);
importConfigurationBuilder.setEnableSharedPosixfsOptimizations(sharedPosixFSOptimizations);
@@ -936,7 +981,7 @@ private FeatureReader getReaderFromPath(final Path variantPath,
/* Anonymous FeatureReader subclass that wraps returned iterators to ensure that the GVCFs do not
* contain MNPs.
*/
- return new FeatureReader() {
+ return new FeatureReader<>() {
/** Iterator that asserts that variants are not MNPs. */
class NoMnpIterator implements CloseableTribbleIterator {
private final CloseableTribbleIterator inner;
@@ -971,7 +1016,8 @@ public VariantContext next() {
return new NoMnpIterator(reader.query(chr, start, end));
}
- @Override public CloseableTribbleIterator iterator() throws IOException {
+ @Override
+ public CloseableTribbleIterator iterator() throws IOException {
return new NoMnpIterator(reader.iterator());
}
};
@@ -990,7 +1036,7 @@ public VariantContext next() {
* @return The workspace directory
*/
private String overwriteCreateOrCheckWorkspace() {
- String workspaceDir = genomicsDBGetAbsolutePath(workspace);
+ String workspaceDir = GATKGenomicsDBUtils.genomicsDBGetAbsolutePath(workspace);
// From JavaDoc for GATKGenomicsDBUtils.createTileDBWorkspacevid
// returnCode = 0 : OK. If overwriteExistingWorkspace is true and the workspace exists, it is deleted first.
// returnCode = -1 : path was not a directory
@@ -1016,7 +1062,7 @@ private String overwriteCreateOrCheckWorkspace() {
}
static class UnableToCreateGenomicsDBWorkspace extends UserException {
- private static final long serialVersionUID = 1L;
+ @Serial private static final long serialVersionUID = 1L;
UnableToCreateGenomicsDBWorkspace(final String message){
super(message);
@@ -1028,7 +1074,7 @@ static class UnableToCreateGenomicsDBWorkspace extends UserException {
* dictionary (as returned by {@link #getBestAvailableSequenceDictionary})
* to parse/verify them. Does nothing if no intervals were specified.
*/
- protected void initializeIntervals() {
+ void initializeIntervals() {
if (intervalArgumentCollection.intervalsSpecified()) {
if (getIntervalsFromExistingWorkspace || doIncrementalImport) {
logger.warn(INCREMENTAL_WORKSPACE_ARG_LONG_NAME+" was set, so ignoring specified intervals." +
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/gvs/common/CommonCode.java b/src/main/java/org/broadinstitute/hellbender/tools/gvs/common/CommonCode.java
index 234fee68b3b..83f1048f535 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/gvs/common/CommonCode.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/gvs/common/CommonCode.java
@@ -3,7 +3,7 @@
import htsjdk.variant.variantcontext.Allele;
import htsjdk.variant.variantcontext.VariantContext;
import htsjdk.variant.vcf.VCFConstants;
-import org.apache.commons.lang.StringUtils;
+import org.apache.commons.lang3.StringUtils;
import org.broadinstitute.hellbender.utils.genotyper.IndexedAlleleList;
import java.util.ArrayList;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortEngine.java
index 935a52a159f..6cd7246fe92 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortEngine.java
@@ -5,8 +5,8 @@
import htsjdk.variant.vcf.VCFConstants;
import htsjdk.variant.vcf.VCFHeader;
import org.apache.avro.generic.GenericRecord;
-import org.apache.commons.lang.StringUtils;
-import org.apache.commons.lang.math.LongRange;
+import org.apache.commons.lang3.StringUtils;
+import org.apache.commons.lang3.LongRange;
import org.apache.logging.log4j.LogManager;
import org.apache.logging.log4j.Logger;
import org.broadinstitute.hellbender.engine.FeatureContext;
@@ -39,10 +39,10 @@ public class ExtractCohortEngine {
// a better, ref-dependent different manner if/when we get around to supporting other references
// See: https://en.wikipedia.org/wiki/Pseudoautosomal_region
private final List PARRegions = List.of(
- new LongRange(23000000010001L, 23000002781479L), // PAR1, X
- new LongRange(24000000010001L, 24000002781479L), // PAR1, Y
- new LongRange(23000155701383L, 23000156030895L), // PAR2, X
- new LongRange(24000056887903L, 24000057217415L)); // PAR2, Y
+ LongRange.of(23000000010001L, 23000002781479L), // PAR1, X
+ LongRange.of(24000000010001L, 24000002781479L), // PAR1, Y
+ LongRange.of(23000155701383L, 23000156030895L), // PAR2, X
+ LongRange.of(24000056887903L, 24000057217415L)); // PAR2, Y
private final boolean printDebugInformation;
private final int localSortMaxRecordsInRam;
@@ -679,7 +679,7 @@ private String makePloidyLookupKey(final String chromosome, final String sampleI
private boolean isInPAR(final long location) {
for (LongRange region : PARRegions) {
- if (region.containsLong(location)) {
+ if (region.contains(location)) {
return true;
}
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortFilterRecord.java b/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortFilterRecord.java
index a81a57d8be3..8ab3b6b9d6f 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortFilterRecord.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortFilterRecord.java
@@ -2,7 +2,7 @@
import htsjdk.samtools.util.Locatable;
import org.apache.avro.generic.GenericRecord;
-import org.apache.commons.lang.builder.ReflectionToStringBuilder;
+import org.apache.commons.lang3.builder.ReflectionToStringBuilder;
import org.broadinstitute.hellbender.tools.gvs.common.SchemaUtils;
public class ExtractCohortFilterRecord implements Locatable {
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortRecord.java b/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortRecord.java
index ffa1e44e77e..623c0678902 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortRecord.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ExtractCohortRecord.java
@@ -2,7 +2,7 @@
import htsjdk.samtools.util.Locatable;
import org.apache.avro.generic.GenericRecord;
-import org.apache.commons.lang.builder.ReflectionToStringBuilder;
+import org.apache.commons.lang3.builder.ReflectionToStringBuilder;
import org.broadinstitute.hellbender.tools.gvs.common.SchemaUtils;
import java.util.Objects;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ReferenceRecord.java b/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ReferenceRecord.java
index 8f929262b7a..e5cda0ef47e 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ReferenceRecord.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/gvs/extract/ReferenceRecord.java
@@ -3,7 +3,7 @@
import htsjdk.samtools.util.Locatable;
import org.apache.avro.generic.GenericRecord;
import org.apache.avro.util.Utf8;
-import org.apache.commons.lang.builder.ReflectionToStringBuilder;
+import org.apache.commons.lang3.builder.ReflectionToStringBuilder;
import org.broadinstitute.hellbender.tools.gvs.common.SchemaUtils;
public class ReferenceRecord implements Locatable, Comparable {
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/gvs/ingest/VetFieldEnum.java b/src/main/java/org/broadinstitute/hellbender/tools/gvs/ingest/VetFieldEnum.java
index a82ad989d26..6df437f01e2 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/gvs/ingest/VetFieldEnum.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/gvs/ingest/VetFieldEnum.java
@@ -4,7 +4,7 @@
import htsjdk.variant.variantcontext.Genotype;
import htsjdk.variant.variantcontext.VariantContext;
import htsjdk.variant.vcf.VCFConstants;
-import org.apache.commons.lang.StringUtils;
+import org.apache.commons.lang3.StringUtils;
import org.apache.logging.log4j.LogManager;
import org.apache.logging.log4j.Logger;
import org.broadinstitute.hellbender.exceptions.UserException;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/spark/CreateHadoopBamSplittingIndex.java b/src/main/java/org/broadinstitute/hellbender/tools/spark/CreateHadoopBamSplittingIndex.java
index de6a914e42c..f1f71a6eb7f 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/spark/CreateHadoopBamSplittingIndex.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/spark/CreateHadoopBamSplittingIndex.java
@@ -1,5 +1,6 @@
package org.broadinstitute.hellbender.tools.spark;
+import com.google.common.io.Files;
import htsjdk.samtools.*;
import htsjdk.samtools.BAMSBIIndexer;
import htsjdk.samtools.seekablestream.SeekableFileStream;
@@ -18,7 +19,6 @@
import org.broadinstitute.hellbender.exceptions.UserException;
import org.broadinstitute.hellbender.utils.io.IOUtils;
import org.broadinstitute.hellbender.utils.read.ReadConstants;
-import org.codehaus.plexus.util.FileUtils;
import picard.cmdline.programgroups.OtherProgramGroup;
import java.io.*;
@@ -166,7 +166,7 @@ private static void assertBamIsCoordinateSorted(final SAMFileHeader header) {
private static void assertIsBam(final File inputBam) {
if(!BamFileIoUtils.isBamFile(inputBam)) {
throw new UserException.BadInput("A splitting index is only relevant for a bam file, but a "
- + "file with extension "+ FileUtils.getExtension(inputBam.getName()) + " was specified.");
+ + "file with extension "+ Files.getFileExtension(inputBam.getName()) + " was specified.");
}
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java b/src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java
index 95bb083f03b..de91be5e5b6 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/GATKSVVCFConstants.java
@@ -89,7 +89,10 @@ public enum ComplexVariantSubtype {
piDUP_RF,
dDUP,
dDUP_iDEL,
- INS_iDEL
+ INS_iDEL,
+ CTX_PP_QQ,
+ CTX_PQ_QP,
+ CTX_INV
}
// not defined in output vcf header but used in internal id that is currently output in the ID column
@@ -163,6 +166,7 @@ public enum ComplexVariantSubtype {
public static final String NONCODING_BREAKPOINT = "PREDICTED_NONCODING_BREAKPOINT";
public static final String NEAREST_TSS = "PREDICTED_NEAREST_TSS";
public static final String TSS_DUP = "PREDICTED_TSS_DUP";
+ public static final String PARTIAL_DISPERSED_DUP = "PREDICTED_PARTIAL_DISPERSED_DUP";
// SVTYPE classes
public enum StructuralVariantAnnotationType {
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/SVFastqUtils.java b/src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/SVFastqUtils.java
index 6c016f84436..2f51a7df114 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/SVFastqUtils.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/spark/sv/utils/SVFastqUtils.java
@@ -4,7 +4,7 @@
import com.esotericsoftware.kryo.Kryo;
import com.esotericsoftware.kryo.io.Input;
import com.esotericsoftware.kryo.io.Output;
-import com.netflix.servo.util.VisibleForTesting;
+import com.google.common.annotations.VisibleForTesting;
import htsjdk.samtools.*;
import htsjdk.samtools.util.Locatable;
import htsjdk.samtools.util.SequenceUtil;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java b/src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java
index 8edcd595881..0f95878b389 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/sv/SVCallRecord.java
@@ -132,8 +132,12 @@ private void validateCoordinates(final SAMSequenceDictionary dictionary) {
Utils.nonNull(dictionary);
validatePosition(contigA, positionA, dictionary);
validatePosition(contigB, positionB, dictionary);
- Utils.validateArg(IntervalUtils.compareLocatables(getPositionAInterval(), getPositionBInterval(), dictionary) <= 0,
- "End coordinate cannot precede start");
+ // CPX types may have position B precede A, such as dispersed duplications where A is the insertion point and
+ // B references the source sequence.
+ if (type != GATKSVVCFConstants.StructuralVariantAnnotationType.CPX) {
+ Utils.validateArg(IntervalUtils.compareLocatables(getPositionAInterval(), getPositionBInterval(), dictionary) <= 0,
+ "End coordinate cannot precede start");
+ }
}
private static void validatePosition(final String contig, final int position, final SAMSequenceDictionary dictionary) {
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinder.java b/src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinder.java
index da4e66f3fb6..9ca5056a679 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinder.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/sv/concordance/ClosestSVFinder.java
@@ -31,6 +31,7 @@
*/
public class ClosestSVFinder {
+ protected final boolean sortOutput;
protected final Map truthIdToItemMap;
protected final Map idToClusterMap;
private final SVConcordanceLinkage linkage;
@@ -49,7 +50,9 @@ public class ClosestSVFinder {
*/
public ClosestSVFinder(final SVConcordanceLinkage linkage,
final Function collapser,
+ final boolean sortOutput,
final SAMSequenceDictionary dictionary) {
+ this.sortOutput = sortOutput;
this.linkage = Utils.nonNull(linkage);
this.collapser = Utils.nonNull(collapser);
outputBuffer = new PriorityQueue<>(SVCallRecordUtils.getCallComparator(dictionary));
@@ -60,6 +63,15 @@ public ClosestSVFinder(final SVConcordanceLinkage linkage,
lastItemContig = null;
}
+ /**
+ * Sorts output by default
+ */
+ public ClosestSVFinder(final SVConcordanceLinkage linkage,
+ final Function collapser,
+ final SAMSequenceDictionary dictionary) {
+ this(linkage, collapser, true, dictionary);
+ }
+
/**
* Should be run frequently to reduce memory usage. Forced flushes must be run when a new contig is encountered.
* @param force flushes all variants in the output buffer regardless of position
@@ -70,8 +82,12 @@ public List flush(final boolean force) {
.map(c -> new ClosestPair(c.getItem(), c.getClosest()))
.map(collapser)
.collect(Collectors.toList());
- outputBuffer.addAll(collapsedRecords);
- return flushBuffer(force);
+ if (sortOutput) {
+ outputBuffer.addAll(collapsedRecords);
+ return flushBuffer(force);
+ } else {
+ return collapsedRecords;
+ }
}
/**
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/CombineGVCFs.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/CombineGVCFs.java
index 50c98a43f3d..23942237008 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/CombineGVCFs.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/CombineGVCFs.java
@@ -461,7 +461,7 @@ private VariantContext referenceBlockMerge(final List vcs, final
final GenotypeBuilder gBuilder = new GenotypeBuilder(g);
GATKVariantContextUtils.makeGenotypeCall(g.getPloidy(),
gBuilder, assignmentMethod,
- g.hasLikelihoods() ? g.getLikelihoods().getAsVector() : null, allelesToUse, g.getAlleles(), null);
+ g.hasLikelihoods() ? g.getLikelihoods().getAsVector() : null, allelesToUse, g, null);
genotypes.add(gBuilder.make());
}
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsEngine.java
index 312a90872e7..9036ecdd44d 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsEngine.java
@@ -150,8 +150,9 @@ private VariantContext regenotypeVC(final VariantContext originalVC, final Refer
final VariantContext result;
+ // only re-genotype polymorphic sites
if ( originalVC.isVariant() && originalVC.getAttributeAsInt(VCFConstants.DEPTH_KEY,0) > 0 ) {
- // only re-genotype polymorphic sites
+ // note that the calculateGenotypes method also calculates the QUAL score
final VariantContext regenotypedVC = calculateGenotypes(originalVC, includeNonVariants);
if (regenotypedVC == null) {
return null;
@@ -186,7 +187,7 @@ private VariantContext regenotypeVC(final VariantContext originalVC, final Refer
//don't count sites with no depth and no confidence towards things like AN and InbreedingCoeff
vcBuilder.genotypes(assignNoCallsAnnotationExcludedGenotypes(result.getGenotypes()));
VariantContext annotated = annotationEngine.annotateContext(vcBuilder.make(), features, ref, null, a -> true);
- return new VariantContextBuilder(annotated).genotypes(cleanupGenotypeAnnotations(result, false, keepSB)).make();
+ return new VariantContextBuilder(annotated).genotypes(cleanupGenotypeAnnotations(annotated, false, keepSB)).make();
} else if (includeNonVariants) {
// For monomorphic sites we need to make sure e.g. the hom ref genotypes are created and only then are passed to the annotation engine.
VariantContext preannotated = new VariantContextBuilder(result).genotypes(cleanupGenotypeAnnotations(result, true, false)).make();
@@ -461,24 +462,29 @@ static List cleanupGenotypeAnnotations(final VariantContext vc, final
attrs.put(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY, GenotypeGVCFs.PHASED_HOM_VAR_STRING);
}
- // create AD if it's not there
- if ( !oldGT.hasAD() && vc.isVariant() ) {
+ // create AD if it's not there, but only if there's data
+ if ( !oldGT.hasAD() && vc.isVariant() && depth > 0) {
final int[] AD = new int[vc.getNAlleles()];
AD[0] = depth;
builder.AD(AD);
}
if ( createRefGTs ) {
- // move the GQ to RGQ
- if (oldGT.hasGQ()) {
+ //keep 0 depth samples and 0 GQ samples as no-call
+ if (depth > 0 && oldGT.hasGQ()) {
+ if (oldGT.getGQ() > 0) {
+ final List refAlleles = Collections.nCopies(oldGT.getPloidy(), vc.getReference());
+ builder.alleles(refAlleles);
+ } else {
+ builder.alleles(Collections.nCopies(oldGT.getPloidy(),Allele.NO_CALL));
+ }
+
+ // move the GQ to RGQ
builder.noGQ();
attrs.put(GATKVCFConstants.REFERENCE_GENOTYPE_QUALITY, oldGT.getGQ());
- }
-
- //keep 0 depth samples and 0 GQ samples as no-call
- if (depth > 0 && oldGT.hasGQ() && oldGT.getGQ() > 0) {
- final List refAlleles = Collections.nCopies(oldGT.getPloidy(), vc.getReference());
- builder.alleles(refAlleles);
+ } else {
+ builder.alleles(Collections.nCopies(oldGT.getPloidy(),Allele.NO_CALL));
+ builder.noGQ().noDP();
}
// also, the PLs are technically no longer usable
@@ -494,8 +500,8 @@ static List cleanupGenotypeAnnotations(final VariantContext vc, final
* Does this genotype look like it has no reads and should be excluded from annotations?
*/
static boolean excludeFromAnnotations(Genotype oldGT) {
- return oldGT.isHomRef() && !oldGT.hasPL()
- && ((oldGT.hasDP() && oldGT.getDP() == 0) || !oldGT.hasDP())
+ return (oldGT.isHomRef() || oldGT.isNoCall())
+ && (!oldGT.hasDP() || oldGT.getDP() == 0)
&& oldGT.hasGQ() && oldGT.getGQ() == 0;
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/ReadAnonymizer.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/ReadAnonymizer.java
index 6c36592b6ae..16b5142fb84 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/ReadAnonymizer.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/ReadAnonymizer.java
@@ -4,7 +4,7 @@
import htsjdk.samtools.CigarElement;
import htsjdk.samtools.CigarOperator;
import org.aeonbits.owner.util.Collections;
-import org.apache.commons.lang.ArrayUtils;
+import org.apache.commons.lang3.ArrayUtils;
import org.broadinstitute.barclay.argparser.Argument;
import org.broadinstitute.barclay.argparser.CommandLineProgramProperties;
import org.broadinstitute.barclay.argparser.ExperimentalFeature;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/ReferenceConfidenceVariantContextMerger.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/ReferenceConfidenceVariantContextMerger.java
index f57d50a7cad..9484da25a9a 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/ReferenceConfidenceVariantContextMerger.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/ReferenceConfidenceVariantContextMerger.java
@@ -439,7 +439,8 @@ protected static void removeStaleAttributesAfterMerge(final Map
attributes.remove(GATKVCFConstants.MLE_ALLELE_COUNT_KEY);
attributes.remove(GATKVCFConstants.MLE_ALLELE_FREQUENCY_KEY);
attributes.remove(VCFConstants.END_KEY);
- attributes.remove(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY); //median doesn't make sense here so drop it; used for ClusteredEventFilter, which doesn't apply to MT
+ attributes.remove(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY);
+ attributes.remove(GATKVCFConstants.EVENT_COUNT_IN_REGION_KEY); //median doesn't make sense here so drop it; used for ClusteredEventFilter, which doesn't apply to MT
}
/**
@@ -579,6 +580,7 @@ private GenotypesContext mergeRefConfidenceGenotypes(final VariantContext vc,
final int ploidy = g.getPloidy();
final GenotypeBuilder genotypeBuilder = new GenotypeBuilder(g);
if (!doSomaticMerge) {
+ //do attribute subsetting
if (g.hasPL() || g.hasAD()) {
int[] perSampleIndexesOfRelevantAlleles = AlleleSubsettingUtils.getIndexesOfRelevantAllelesForGVCF(remappedAlleles, targetAlleles, vc.getStart(), g, false);
if (g.hasPL()) {
@@ -590,8 +592,10 @@ private GenotypesContext mergeRefConfidenceGenotypes(final VariantContext vc,
if (g.hasAD()) {
genotypeBuilder.AD(AlleleSubsettingUtils.generateAD(g.getAD(), perSampleIndexesOfRelevantAlleles));
}
+
+ }
//clean up low confidence hom refs for better annotations later
- } else if (GenotypeGVCFsEngine.excludeFromAnnotations(g)) {
+ if (GenotypeGVCFsEngine.excludeFromAnnotations(g)) {
genotypeBuilder.alleles(Collections.nCopies(ploidy, Allele.NO_CALL));
}
}
@@ -658,7 +662,7 @@ private GenotypesContext mergeRefConfidenceGenotypes(final VariantContext vc,
GATKVariantContextUtils.makeGenotypeCall(g.getPloidy(),
genotypeBuilder, assignmentMethod,
g.hasLikelihoods() ? g.getLikelihoods().getAsVector() : null,
- targetAlleles, originalGTAlleles, null);
+ targetAlleles, new GenotypeBuilder(g).alleles(originalGTAlleles).make(), null);
mergedGenotypes.add(genotypeBuilder.make());
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AnnotationUtils.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AnnotationUtils.java
index b0ff4985cec..7ba0a0baa24 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AnnotationUtils.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AnnotationUtils.java
@@ -3,14 +3,19 @@
import htsjdk.variant.variantcontext.Allele;
import htsjdk.variant.variantcontext.Genotype;
import htsjdk.variant.variantcontext.VariantContext;
-import org.apache.commons.lang.StringUtils;
+import org.apache.commons.lang3.StringUtils;
+import htsjdk.variant.vcf.VCFHeaderLineCount;
+import htsjdk.variant.vcf.VCFInfoHeaderLine;
import org.broadinstitute.hellbender.tools.walkers.annotator.allelespecific.*;
import org.broadinstitute.hellbender.utils.genotyper.AlleleLikelihoods;
import org.broadinstitute.hellbender.utils.read.GATKRead;
+import org.broadinstitute.hellbender.utils.variant.GATKVCFHeaderLines;
import java.util.*;
public final class AnnotationUtils {
+
+ public static final String ALLELE_SPECIFIC_ANNOTATION_KEY_PREFIX = "AS_";
public static final String ALLELE_SPECIFIC_RAW_DELIM = "|";
public static final String ALLELE_SPECIFIC_REDUCED_DELIM = ",";
public static final String ALLELE_SPECIFIC_SPLIT_REGEX = "\\|"; //String.split takes a regex, so we need to escape the pipe
@@ -74,16 +79,23 @@ public static List decodeAnyASList( final String somethingList) {
* @param annotation the annotation to be tested
* @return true if the annotation is expected to have values per-allele
*/
- public static boolean isAlleleSpecific(final InfoFieldAnnotation annotation) {
+ public static boolean isAlleleSpecific(final VariantAnnotation annotation) {
return annotation instanceof AlleleSpecificAnnotation;
}
+ public static boolean isAlleleSpecificGatkKey(final String annotationKey) {
+ final VCFInfoHeaderLine header = GATKVCFHeaderLines.getInfoLine(annotationKey);
+ return header.getCountType().equals(VCFHeaderLineCount.A) ||
+ header.getCountType().equals(VCFHeaderLineCount.R) ||
+ annotationKey.startsWith(ALLELE_SPECIFIC_ANNOTATION_KEY_PREFIX);
+ }
+
/**
- * Handles all the Java and htsjdk parsing shenanigans
- * @param rawDataString should not have surrounding brackets
+ * Handles all the Java and htsjdk parsing shenanigans from getAttributeAsString
+ * @param rawDataString may have surrounding brackets, with raw delimiter
* @return
*/
- public static List getAlleleLengthListOfString(String rawDataString) {
+ public static List getAlleleLengthListOfStringFromRawData(String rawDataString) {
if (rawDataString == null) {
return Collections.emptyList();
}
@@ -93,6 +105,21 @@ public static List getAlleleLengthListOfString(String rawDataString) {
return Arrays.asList(rawDataString.split(ALLELE_SPECIFIC_SPLIT_REGEX, -1)); //-1 to keep empty data
}
+ /**
+ * Handles all the Java and htsjdk parsing shenanigans from getAttributeAsString
+ * @param dataString may have surrounding brackets, with reduced delimieter
+ * @return
+ */
+ public static List getAlleleLengthListOfString(String dataString) {
+ if (dataString == null) {
+ return Collections.emptyList();
+ }
+ if (dataString.startsWith("[")) {
+ dataString = dataString.substring(1, dataString.length() - 1).replaceAll("\\s", "");
+ }
+ return Arrays.asList(dataString.split(ALLELE_SPECIFIC_REDUCED_DELIM, -1)); //-1 to keep empty data
+ }
+
static public String generateMissingDataWarning(final VariantContext vc, final Genotype g, final AlleleLikelihoods likelihoods) {
final StringBuilder outString = new StringBuilder("Annotation will not be calculated at position " + vc.getContig() + ":" + vc.getStart() +
" and possibly subsequent");
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AssemblyComplexity.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AssemblyComplexity.java
index 931a6bb5256..8ea5771b78e 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AssemblyComplexity.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AssemblyComplexity.java
@@ -4,12 +4,13 @@
import htsjdk.samtools.util.Locatable;
import htsjdk.variant.variantcontext.Allele;
import htsjdk.variant.variantcontext.VariantContext;
-import org.apache.commons.lang.mutable.MutableInt;
+import org.apache.commons.lang3.mutable.MutableInt;
import org.broadinstitute.barclay.argparser.Argument;
import org.apache.commons.lang3.tuple.Triple;
import org.broadinstitute.barclay.help.DocumentedFeature;
import org.broadinstitute.hellbender.engine.FeatureContext;
import org.broadinstitute.hellbender.engine.ReferenceContext;
+import org.broadinstitute.hellbender.tools.walkers.annotator.allelespecific.AlleleSpecificAnnotation;
import org.broadinstitute.hellbender.utils.MathUtils;
import org.broadinstitute.hellbender.utils.genotyper.AlleleLikelihoods;
import org.broadinstitute.hellbender.utils.haplotype.Event;
@@ -27,7 +28,7 @@
*/
@DocumentedFeature(groupName= HelpConstants.DOC_CAT_ANNOTATORS, groupSummary=HelpConstants.DOC_CAT_ANNOTATORS_SUMMARY,
summary="Describe the complexity of an assembly region")
-public class AssemblyComplexity implements JumboInfoAnnotation {
+public class AssemblyComplexity implements JumboInfoAnnotation, AlleleSpecificAnnotation {
@Argument(fullName = "assembly-complexity-reference-mode",
doc="If enabled will treat the reference as the basis for assembly complexity as opposed to estimated germline haplotypes",
@@ -189,5 +190,4 @@ private static int uniqueVariants(final Haplotype hap1, final Haplotype hap2, fi
private static int editDistance(final Haplotype hap1, final Haplotype hap2, final int excludedPosition) {
return uniqueVariants(hap1, hap2, excludedPosition) + uniqueVariants(hap2, hap1, excludedPosition);
}
-
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/FeaturizedReadSets.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/FeaturizedReadSets.java
deleted file mode 100644
index c870d9693ee..00000000000
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/FeaturizedReadSets.java
+++ /dev/null
@@ -1,175 +0,0 @@
-package org.broadinstitute.hellbender.tools.walkers.annotator;
-
-
-import htsjdk.variant.variantcontext.Allele;
-import htsjdk.variant.variantcontext.Genotype;
-import htsjdk.variant.variantcontext.GenotypeBuilder;
-import htsjdk.variant.variantcontext.VariantContext;
-import org.apache.commons.lang.StringUtils;
-import org.broadinstitute.barclay.help.DocumentedFeature;
-import org.broadinstitute.gatk.nativebindings.smithwaterman.SWOverhangStrategy;
-import org.broadinstitute.hellbender.engine.FeatureContext;
-import org.broadinstitute.hellbender.engine.ReferenceContext;
-import org.broadinstitute.hellbender.utils.Utils;
-import org.broadinstitute.hellbender.utils.genotyper.AlleleLikelihoods;
-import org.broadinstitute.hellbender.utils.haplotype.Haplotype;
-import org.broadinstitute.hellbender.utils.help.HelpConstants;
-import org.broadinstitute.hellbender.utils.read.AlignmentUtils;
-import org.broadinstitute.hellbender.utils.read.Fragment;
-import org.broadinstitute.hellbender.utils.read.GATKRead;
-import org.broadinstitute.hellbender.utils.smithwaterman.SmithWatermanAligner;
-import org.broadinstitute.hellbender.utils.smithwaterman.SmithWatermanAlignment;
-import org.broadinstitute.hellbender.utils.smithwaterman.SmithWatermanAlignmentConstants;
-import org.broadinstitute.hellbender.utils.variant.GATKVCFConstants;
-
-import java.util.*;
-import java.util.stream.Collectors;
-
-/**
- * For each sample and for each allele a list feature vectors of supporting reads
- * In order to reduce the number of delimiter characters, we flatten featurized reads. For example, suppose allele 1 has featurized reads
- * [1,2] and [3,4] and allele 2 has featurized reads [5,6] and [7,8], the annotation is
- * 1,2,3,4|5,6,7,8
- */
-@DocumentedFeature(groupName= HelpConstants.DOC_CAT_ANNOTATORS, groupSummary=HelpConstants.DOC_CAT_ANNOTATORS_SUMMARY,
- summary="Featurized read sets for Mutect3 training data")
-public class FeaturizedReadSets implements JumboGenotypeAnnotation {
- public static final int DEFAULT_BASE_QUALITY = 25;
-
- private static final int DEFAULT_MAX_REF_COUNT = Integer.MAX_VALUE;
-
- public static final int FEATURES_PER_READ = 11;
-
- private static final SmithWatermanAligner aligner = SmithWatermanAligner.getAligner(SmithWatermanAligner.Implementation.JAVA);
-
- // downsample ref reads to this count if needed
- private final int maxRefCount;
-
- public FeaturizedReadSets(final int maxRefCount) {
- this.maxRefCount = maxRefCount;
- }
-
- public FeaturizedReadSets() {
- this(DEFAULT_MAX_REF_COUNT);
- }
-
- @Override
- public void annotate(final ReferenceContext ref,
- final FeatureContext features,
- final VariantContext vc,
- final Genotype g,
- final GenotypeBuilder gb,
- final AlleleLikelihoods likelihoods,
- final AlleleLikelihoods fragmentLikelihoods,
- final AlleleLikelihoods haplotypeLikelihoods) {
- Utils.nonNull(vc);
- Utils.nonNull(g);
- Utils.nonNull(gb);
-
- if ( likelihoods == null) {
- return;
- }
-
- final List>> readVectorsByAllele = getReadVectors(vc, Collections.singletonList(g.getSampleName()),
- likelihoods, haplotypeLikelihoods, maxRefCount, Integer.MAX_VALUE);
-
- // flatten twice: over all reads supporting each allele and over all alleles
- // we can partition by allele with the countsInAlleleOrder annotation
- // and by read using the constant feature vector size
- final int[] flattenedTensorInAlleleOrder = readVectorsByAllele.stream()
- .flatMap(listOfLists -> listOfLists.stream().flatMap(List::stream))
- .mapToInt(n -> n)
- .toArray();
-
- final int[] countsInAlleleOrder = readVectorsByAllele.stream().mapToInt(List::size).toArray();
-
- gb.attribute(GATKVCFConstants.FEATURIZED_READ_SETS_KEY, flattenedTensorInAlleleOrder);
- gb.attribute(GATKVCFConstants.FEATURIZED_READ_SETS_COUNTS_KEY, countsInAlleleOrder);
- }
-
- public static List>> getReadVectors(final VariantContext vc,
- final Collection samples,
- final AlleleLikelihoods likelihoods,
- final AlleleLikelihoods haplotypeLikelihoods,
- final int refDownsample,
- final int altDownsample) {
- return getReadVectors(vc, samples, likelihoods, haplotypeLikelihoods, refDownsample, altDownsample, Collections.emptyMap());
- }
-
- // returns Lists (in allele order) of lists of read vectors supporting each allele
- public static List>> getReadVectors(final VariantContext vc,
- final Collection samples,
- final AlleleLikelihoods likelihoods,
- final AlleleLikelihoods haplotypeLikelihoods,
- final int refDownsample,
- final int altDownsample,
- final Map altDownsampleMap) {
- final Map> readsByAllele = likelihoods.alleles().stream()
- .collect(Collectors.toMap(a -> a, a -> new ArrayList<>()));
-
- samples.stream().flatMap(s -> likelihoods.bestAllelesBreakingTies(s).stream())
- .filter(ba -> ba.isInformative())
- .forEach(ba -> readsByAllele.get(ba.allele).add(ba.evidence));
-
- // downsample if necessary
- final Allele refAllele = likelihoods.alleles().stream().filter(Allele::isReference).findFirst().get();
- for (final Allele allele : likelihoods.alleles()) {
- final int downsample = allele.isReference() ? refDownsample : altDownsampleMap.getOrDefault(allele, altDownsample);
- if (readsByAllele.get(allele).size() > downsample) {
- Collections.shuffle(readsByAllele.get(allele));
- readsByAllele.put(allele, readsByAllele.get(allele).subList(0, downsample));
- }
- }
-
- final Map bestHaplotypes = new HashMap<>();
- samples.stream().flatMap(s -> haplotypeLikelihoods.bestAllelesBreakingTies(s).stream())
- .forEach(ba -> ba.evidence.getReads().forEach(read -> bestHaplotypes.put(read, ba.allele)));
-
- return vc.getAlleles().stream()
- .map(allele -> readsByAllele.get(allele).stream().map(read -> featurize(read, vc, bestHaplotypes)).collect(Collectors.toList()))
- .collect(Collectors.toList());
- }
-
-
- @Override
- public List getKeyNames() {
- return Arrays.asList(GATKVCFConstants.FEATURIZED_READ_SETS_KEY, GATKVCFConstants.FEATURIZED_READ_SETS_COUNTS_KEY);
- }
-
- private static List featurize(final GATKRead read, final VariantContext vc, final Map bestHaplotypes) {
- final List result = new ArrayList<>();
- result.add(read.getMappingQuality());
- result.add(BaseQuality.getBaseQuality(read, vc).orElse(DEFAULT_BASE_QUALITY));
- result.add(read.isFirstOfPair() ? 1 : 0);
- result.add(read.isReverseStrand() ? 1 : 0);
-
- // distances from ends of read
- final int readPosition = ReadPosition.getPosition(read, vc).orElse(0);
- result.add(readPosition);
- result.add(read.getLength() - readPosition);
-
-
- result.add(Math.abs(read.getFragmentLength()));
-
- // distances from ends of fragment
- final int fragmentStart = Math.min(read.getMateStart(), read.getUnclippedStart());
- final int fragmentEnd = fragmentStart + Math.abs(read.getFragmentLength());
- result.add(vc.getStart() - fragmentStart);
- result.add(fragmentEnd - vc.getEnd());
-
- // mismatches versus best haplotype
- final Haplotype haplotype = bestHaplotypes.get(read);
- final SmithWatermanAlignment readToHaplotypeAlignment = aligner.align(haplotype.getBases(), read.getBases(), SmithWatermanAlignmentConstants.ALIGNMENT_TO_BEST_HAPLOTYPE_SW_PARAMETERS, SWOverhangStrategy.SOFTCLIP);
- final GATKRead copy = read.copy();
- copy.setCigar(readToHaplotypeAlignment.getCigar());
- final int mismatchCount = AlignmentUtils.getMismatchCount(copy, haplotype.getBases(), readToHaplotypeAlignment.getAlignmentOffset()).numMismatches;
- result.add(mismatchCount);
-
- final long indelsVsBestHaplotype = readToHaplotypeAlignment.getCigar().getCigarElements().stream().filter(el -> el.getOperator().isIndel()).count();
- result.add((int) indelsVsBestHaplotype);
- Utils.validate(result.size() == FEATURES_PER_READ, "Wrong number of features");
-
- return result;
- }
-
-}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/OrientationBiasReadCounts.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/OrientationBiasReadCounts.java
index addbd29dfd7..96a749de50c 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/OrientationBiasReadCounts.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/OrientationBiasReadCounts.java
@@ -4,7 +4,7 @@
import htsjdk.variant.variantcontext.Genotype;
import htsjdk.variant.variantcontext.GenotypeBuilder;
import htsjdk.variant.variantcontext.VariantContext;
-import org.apache.commons.lang.mutable.MutableInt;
+import org.apache.commons.lang3.mutable.MutableInt;
import org.broadinstitute.barclay.help.DocumentedFeature;
import org.broadinstitute.hellbender.engine.FeatureContext;
import org.broadinstitute.hellbender.engine.ReferenceContext;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/ReferenceBases.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/ReferenceBases.java
index 9c2d99b97b7..3ef5ade7148 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/ReferenceBases.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/ReferenceBases.java
@@ -4,7 +4,7 @@
import htsjdk.variant.variantcontext.VariantContext;
import htsjdk.variant.vcf.VCFHeaderLineType;
import htsjdk.variant.vcf.VCFInfoHeaderLine;
-import org.apache.commons.lang.StringUtils;
+import org.apache.commons.lang3.StringUtils;
import org.broadinstitute.barclay.help.DocumentedFeature;
import org.broadinstitute.hellbender.engine.ReferenceContext;
import org.broadinstitute.hellbender.utils.Utils;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/VariantAnnotatorEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/VariantAnnotatorEngine.java
index e83be35f50f..930a6a51033 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/VariantAnnotatorEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/VariantAnnotatorEngine.java
@@ -145,6 +145,10 @@ public List getInfoAnnotations() {
return Collections.unmodifiableList(infoAnnotations);
}
+ public List getJumboInfoAnnotations() {
+ return Collections.unmodifiableList(jumboInfoAnnotations);
+ }
+
/**
*
* @param infoAnnotationClassName the name of the Java class, NOT the annotation VCF key
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_QualByDepth.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_QualByDepth.java
index 9d78efd6d1b..f8bb3b3cea9 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_QualByDepth.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_QualByDepth.java
@@ -7,7 +7,7 @@
import htsjdk.variant.vcf.VCFCompoundHeaderLine;
import htsjdk.variant.vcf.VCFHeaderLine;
import htsjdk.variant.vcf.VCFInfoHeaderLine;
-import org.apache.commons.lang.StringUtils;
+import org.apache.commons.lang3.StringUtils;
import org.broadinstitute.barclay.help.DocumentedFeature;
import org.broadinstitute.hellbender.engine.ReferenceContext;
import org.broadinstitute.hellbender.tools.walkers.annotator.AnnotationUtils;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_RankSumTest.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_RankSumTest.java
index 5b4af7ec0ff..ef01d6fbf62 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_RankSumTest.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_RankSumTest.java
@@ -5,7 +5,7 @@
import htsjdk.variant.variantcontext.GenotypesContext;
import htsjdk.variant.variantcontext.VariantContext;
import htsjdk.variant.vcf.VCFConstants;
-import org.apache.commons.lang.ArrayUtils;
+import org.apache.commons.lang3.ArrayUtils;
import org.apache.logging.log4j.LogManager;
import org.apache.logging.log4j.Logger;
import org.broadinstitute.hellbender.engine.ReferenceContext;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_StrandBiasTest.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_StrandBiasTest.java
index d8466d0d997..b521f157d28 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_StrandBiasTest.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_StrandBiasTest.java
@@ -153,7 +153,7 @@ public Map finalizeRawData(final VariantContext vc, final Varia
}
protected void parseRawDataString(ReducibleAnnotationData> myData) {
- List values = AnnotationUtils.getAlleleLengthListOfString(myData.getRawData());
+ List values = AnnotationUtils.getAlleleLengthListOfStringFromRawData(myData.getRawData());
if (values.size() != myData.getAlleles().size()) {
throw new IllegalStateException("Number of alleles and number of allele-specific entries do not match. " +
"Allele-specific annotations should have an entry for each allele including the reference.");
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/flow/FlowAnnotatorBase.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/flow/FlowAnnotatorBase.java
index 0c4a07f015a..11890925065 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/flow/FlowAnnotatorBase.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/flow/FlowAnnotatorBase.java
@@ -4,6 +4,7 @@
import htsjdk.variant.variantcontext.VariantContext;
import org.apache.logging.log4j.LogManager;
import org.apache.logging.log4j.Logger;
+import org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions;
import org.broadinstitute.hellbender.engine.ReferenceContext;
import org.broadinstitute.hellbender.exceptions.GATKException;
import org.broadinstitute.hellbender.tools.walkers.annotator.InfoFieldAnnotation;
@@ -33,7 +34,6 @@
* are accumulated as well.
*/
public abstract class FlowAnnotatorBase implements InfoFieldAnnotation {
- private final static Logger logger = LogManager.getLogger(FlowAnnotatorBase.class);
protected final OneShotLogger flowMissingOneShotLogger = new OneShotLogger(FlowAnnotatorBase.class);
@@ -149,7 +149,7 @@ private String establishReadGroupFlowOrder(final LocalContext localContext, fina
// if here, no flow order was found. may we use a default?
if ( isActualFlowOrderRequired() ) {
localContext.generateAnnotation = false;
- flowMissingOneShotLogger.warn("this.getClass().getSimpleName() + \" annotation will not be calculated, no '\" + StandardArgumentDefinitions.FLOW_ORDER_FOR_ANNOTATIONS + \"' argument provided\"");
+ flowMissingOneShotLogger.warn(this.getClass().getSimpleName() + " annotation will not be calculated, no '" + StandardArgumentDefinitions.FLOW_ORDER_FOR_ANNOTATIONS + "' argument provided");
}
return FlowBasedRead.DEFAULT_FLOW_ORDER;
@@ -201,7 +201,7 @@ protected void variantType(final VariantContext vc, final LocalContext localCont
if (isSpecial(alleles.get(i))){
continue;
}
- if ((localContext.hmerIndelLength.get(i)==null) || (localContext.hmerIndelLength.get(i)==0)){
+ if ((localContext.hmerIndelLength==null) || (localContext.hmerIndelLength.get(i)==null) || (localContext.hmerIndelLength.get(i)==0)){
isHmer=false;
}
}
@@ -249,7 +249,12 @@ protected void isHmerIndel(final VariantContext vc, final LocalContext localCont
// get byte before and after
final byte before = getReferenceNucleotide(localContext, vc.getStart() - 1);
final byte[] after = getReferenceHmerPlus(localContext, vc.getEnd() + 1, MOTIF_SIZE);
-
+ if (after.length==0){
+ flowMissingOneShotLogger.warn("Failed to get hmer from reference, probably because the " +
+ "variant is very close to the end of the chromosome, isHmerIndel and RightMotif annotations will not be calculated. " +
+ "Variant position: " + vc.getContig() + ":" + vc.getEnd() + 1);
+ return;
+ }
// build two haplotypes. add byte before and after
final byte[] refHap = buildHaplotype(before, ref.getBases(), after);
final byte[] altHap = buildHaplotype(before, alt.getBases(), after);
@@ -337,6 +342,12 @@ protected void getLeftMotif(final VariantContext vc, final LocalContext localCon
}
String motif = getRefMotif(localContext, vc.getStart() - MOTIF_SIZE, MOTIF_SIZE);
+ if (motif.length() != MOTIF_SIZE){
+ flowMissingOneShotLogger.warn("Failed to get motif from reference, probably because the variant is very close to the " +
+ "start of the chromosome. LeftMotif annotation will not be calculated. " +
+ "Variant position: " + vc.getContig() + ":" + vc.getStart());
+ return;
+ }
if ( a.length() != refLength ) {
motif = motif.substring(1) + vc.getReference().getBaseString().substring(0, 1);
}
@@ -349,8 +360,13 @@ protected void getLeftMotif(final VariantContext vc, final LocalContext localCon
protected void getRightMotif(final VariantContext vc, final LocalContext localContext) {
final int refLength = vc.getReference().length();
- final String motif = getRefMotif(localContext, vc.getStart() + refLength, MOTIF_SIZE);
-
+ String motif = getRefMotif(localContext, vc.getStart() + refLength, MOTIF_SIZE);
+ if (motif.length() != MOTIF_SIZE){
+ flowMissingOneShotLogger.warn("Failed to get motif from reference, probably because " +
+ "the variant is close to the end of the chromosome. RightMotif annotation will not be calculated. " +
+ "Variant position: " + vc.getContig() + ":" + vc.getStart());
+ return;
+ }
// fill empty entries (non indel alelles)
for ( int i = 0 ; i < localContext.rightMotif.size() ; i++ ) {
if ( localContext.rightMotif.get(i) == null ) {
@@ -365,6 +381,11 @@ protected void gcContent(final VariantContext vc, final LocalContext localContex
final int begin = vc.getStart() - (GC_CONTENT_SIZE / 2);
final String seq = getRefMotif(localContext, begin + 1, GC_CONTENT_SIZE);
+ if ( seq.length() != GC_CONTENT_SIZE ) {
+ flowMissingOneShotLogger.warn("gcContent will not be calculated at position " + vc.getContig() + ":" + vc.getStart() +
+ " because the variant is too close to the edge of the chromosome");
+ return;
+ }
int gcCount = 0;
for ( byte b : seq.getBytes() ) {
if ( b == 'G' || b == 'C' ) {
@@ -423,11 +444,11 @@ protected void cycleSkip(final VariantContext vc, final LocalContext localContex
localContext.attributes.put(GATKVCFConstants.FLOW_CYCLESKIP_STATUS, css);
}
- // get a single nucleoid from reference
+ // get a single nucleotid from reference
private byte getReferenceNucleotide(final LocalContext localContext, final int start) {
final int index = start - localContext.ref.getWindow().getStart();
final byte[] bases = localContext.ref.getBases();
- Utils.validIndex(index, bases.length);
+ Utils.validIndex(index, bases.length); // do not catch, if here the location is outside of the reference, there is a problem!
return bases[index];
}
@@ -435,7 +456,13 @@ private byte getReferenceNucleotide(final LocalContext localContext, final int s
private byte[] getReferenceHmerPlus(final LocalContext localContext, final int start, final int additional) {
int index = start - localContext.ref.getWindow().getStart();
final byte[] bases = localContext.ref.getBases();
- Utils.validIndex(index, bases.length);
+ try {
+ Utils.validIndex(index, bases.length);
+ } catch (IllegalArgumentException e) {
+ flowMissingOneShotLogger.warn("Failed to get hmer from reference, too close to the edge. " +
+ "Position: " + localContext.ref.getContig() + ":" + index);
+ return new byte[0];
+ }
// get hmer
final StringBuilder sb = new StringBuilder();
@@ -457,8 +484,12 @@ private String getRefMotif(final LocalContext localContext, final int start, fin
final byte[] bases = localContext.ref.getBases();
final int startIndex = start - localContext.ref.getWindow().getStart();
final int endIndex = startIndex + length;
- Utils.validIndex(startIndex, bases.length);
- Utils.validIndex(endIndex-1, bases.length);
+ try {
+ Utils.validIndex(startIndex, bases.length);
+ Utils.validIndex(endIndex-1, bases.length);
+ } catch (IllegalArgumentException e) {
+ return "";
+ }
return new String(Arrays.copyOfRange(bases, startIndex, endIndex));
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/ContaminationModel.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/ContaminationModel.java
index 44b07f6fcfe..93420133913 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/ContaminationModel.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/ContaminationModel.java
@@ -1,7 +1,7 @@
package org.broadinstitute.hellbender.tools.walkers.contamination;
import htsjdk.samtools.util.OverlapDetector;
-import org.apache.commons.lang.mutable.MutableDouble;
+import org.apache.commons.lang3.mutable.MutableDouble;
import org.apache.commons.lang3.tuple.ImmutablePair;
import org.apache.commons.lang3.tuple.Pair;
import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/ContaminationSegmenter.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/ContaminationSegmenter.java
index 1a023488270..607f6691c58 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/ContaminationSegmenter.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/ContaminationSegmenter.java
@@ -13,7 +13,7 @@
import java.util.stream.IntStream;
public class ContaminationSegmenter {
- public static final Range ALT_FRACTIONS_FOR_SEGMENTATION = Range.between(0.1, 0.9);
+ public static final Range ALT_FRACTIONS_FOR_SEGMENTATION = Range.of(0.1, 0.9);
public static final double KERNEL_SEGMENTER_LINEAR_COST = 1.0;
public static final double KERNEL_SEGMENTER_LOG_LINEAR_COST = 1.0;
public static final int KERNEL_SEGMENTER_DIMENSION = 100;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/GetPileupSummaries.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/GetPileupSummaries.java
index 9c4f09d1bcd..9d6c40a691c 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/GetPileupSummaries.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/contamination/GetPileupSummaries.java
@@ -9,6 +9,7 @@
import org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions;
import org.broadinstitute.hellbender.cmdline.programgroups.CoverageAnalysisProgramGroup;
import org.broadinstitute.hellbender.engine.*;
+import org.broadinstitute.hellbender.engine.filters.MappingQualityReadFilter;
import org.broadinstitute.hellbender.engine.filters.ReadFilter;
import org.broadinstitute.hellbender.engine.filters.ReadFilterLibrary;
import org.broadinstitute.hellbender.engine.filters.WellformedReadFilter;
@@ -100,6 +101,13 @@
* file that have AF of 0.01 or more.
*
*
+ *
+ * Finally for those using mappers other than bwa mem or dragen-os, {@code --minimum-mapping-quality} threshold is
+ * set to 50, which limits the usable reads that tool considers for generating pileups. Certain mappers are known to
+ * assign scores less that this threshold even for the unique mappings. If you observe all empty results in your
+ * summary file please adjust the {@code --minimum-mapping-quality} parameter according to your input files.
+ *
+ *
*/
@CommandLineProgramProperties(
summary = "Tabulates pileup metrics for inferring contamination",
@@ -112,8 +120,7 @@ public class GetPileupSummaries extends LocusWalker {
public static final String MIN_SITE_AF_LONG_NAME = "minimum-population-allele-frequency";
public static final String MAX_SITE_AF_SHORT_NAME = "max-af";
public static final String MIN_SITE_AF_SHORT_NAME = "min-af";
- public static final String MIN_MAPPING_QUALITY_LONG_NAME = "min-mapping-quality";
- public static final String MIN_MAPPING_QUALITY_SHORT_NAME = "mmq";
+
private static final double DEFAULT_MIN_POPULATION_AF = 0.01;
private static final double DEFAULT_MAX_POPULATION_AF = 0.2;
@@ -137,9 +144,6 @@ public class GetPileupSummaries extends LocusWalker {
doc = "Maximum population allele frequency of sites to consider.", optional = true)
private double maxPopulationAlleleFrequency = DEFAULT_MAX_POPULATION_AF;
- @Argument(fullName = MIN_MAPPING_QUALITY_LONG_NAME, shortName = MIN_MAPPING_QUALITY_SHORT_NAME, doc = "Minimum read mapping quality", optional = true)
- private int minMappingQuality = DEFAULT_MINIMUM_MAPPING_QUALITY;
-
private boolean sawVariantsWithoutAlleleFrequency = false;
private boolean sawVariantsWithAlleleFrequency = false;
@@ -168,6 +172,7 @@ public boolean requiresFeatures() {
@Override
public List getDefaultReadFilters() {
final List filters = new ArrayList<>();
+ filters.add(new MappingQualityReadFilter(DEFAULT_MINIMUM_MAPPING_QUALITY));
filters.add(ReadFilterLibrary.MAPPING_QUALITY_AVAILABLE);
filters.add(ReadFilterLibrary.MAPPING_QUALITY_NOT_ZERO);
filters.add(ReadFilterLibrary.MAPPED);
@@ -208,8 +213,7 @@ public void apply(AlignmentContext alignmentContext, ReferenceContext referenceC
final VariantContext vc = vcs.get(0);
if ( vc.isBiallelic() && vc.isSNP() && alleleFrequencyInRange(vc) ) {
- final ReadPileup pileup = alignmentContext.getBasePileup()
- .makeFilteredPileup(pe -> pe.getRead().getMappingQuality() >= minMappingQuality);
+ final ReadPileup pileup = alignmentContext.getBasePileup();
try {
writer.writeRecord(new PileupSummary(vc, pileup));
} catch (final IOException ex) {
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQuality.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQuality.java
new file mode 100644
index 00000000000..83ac040011c
--- /dev/null
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQuality.java
@@ -0,0 +1,626 @@
+package org.broadinstitute.hellbender.tools.walkers.featuremapping;
+
+import com.google.common.annotations.VisibleForTesting;
+import htsjdk.samtools.SAMFileHeader;
+import htsjdk.samtools.SAMUtils;
+import org.apache.commons.lang3.StringUtils;
+import org.broadinstitute.barclay.argparser.*;
+import org.broadinstitute.barclay.help.DocumentedFeature;
+import org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions;
+import org.broadinstitute.hellbender.cmdline.programgroups.FlowBasedProgramGroup;
+import org.broadinstitute.hellbender.engine.FeatureContext;
+import org.broadinstitute.hellbender.engine.GATKPath;
+import org.broadinstitute.hellbender.engine.ReadWalker;
+import org.broadinstitute.hellbender.engine.ReferenceContext;
+import org.broadinstitute.hellbender.exceptions.GATKException;
+import org.broadinstitute.hellbender.tools.FlowBasedArgumentCollection;
+import org.broadinstitute.hellbender.tools.walkers.groundtruth.SeriesStats;
+import org.broadinstitute.hellbender.utils.read.*;
+
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.util.*;
+
+@CommandLineProgramProperties(
+ summary = "This program converts the flow qualities that Ultima Genomics CRAM reports to more conventional base qualities. " +
+ "Specifically, the generated quality will report the probability that a specific base is a sequencing error mismatch, " +
+ "while auxilary tags qa, qt, qg, qc report specific probability that a specific base X is a A->X error. " +
+ "Since mismatch error in flow-based chemistries can only occur as a result of several indel errors, " +
+ "we implemented several strategies to estimate the probability of a mismatch which can be specified" +
+ "using the svnq-mode parameter: " +
+ "Legacy - the quality value from flow matrix is used. " +
+ "Optimistic - assuming that the probability of the indel errors are p1 and p2, then snvq=p1*p2 - assuming they always coincide. " +
+ "Pessimistic - snvq=(1-p1)*(1-p2) - assuming they never coincide. " +
+ "Geometric - snvq=sqrt(Optimistic*Pessimistic) - i.e. the geometric mean of the optimistic and Pessimistic modes. " +
+ "The Geometric is set as the default mode",
+ oneLineSummary = "Add SNV Quality to the flow-based CRAM",
+ programGroup = FlowBasedProgramGroup.class
+)
+
+@DocumentedFeature
+@ExperimentalFeature
+public final class AddFlowSNVQuality extends ReadWalker {
+
+ @Argument(fullName = StandardArgumentDefinitions.OUTPUT_LONG_NAME,
+ shortName = StandardArgumentDefinitions.OUTPUT_SHORT_NAME,
+ doc = "File to which reads should be written")
+ @WorkflowOutput(optionalCompanions={StandardArgumentDefinitions.OUTPUT_INDEX_COMPANION})
+ public GATKPath output = null;
+ private SAMFileGATKReadWriter outputWriter;
+
+ @ArgumentCollection
+ public FlowBasedArgumentCollection fbargs = new FlowBasedArgumentCollection();
+
+ @ArgumentCollection
+ public AddFlowSNVQualityArgumentCollection aqArgs = new AddFlowSNVQualityArgumentCollection();
+
+ public static final char PHRED_ASCII_BASE = '!';
+
+ public static final int ERROR_PROB_BAND_1LESS = 0;
+ public static final int ERROR_PROB_BAND_KEY = 1;
+ public static final int ERROR_PROB_BAND_1MORE = 2;
+ public static final int ERROR_PROB_BANDS = 3;
+
+ public double minLikelihoodProbRate = 1e-6;
+ public int maxQualityScore = 60;
+
+ // locals
+ private SeriesStats inputQualStats = new SeriesStats();
+ private SeriesStats outputBQStats = new SeriesStats();
+ private SeriesStats outputQAltStats = new SeriesStats();
+ private SeriesStats outputQCalledStats = new SeriesStats();
+ private SeriesStats outputSumPStats = new SeriesStats();
+
+ // private class to hold the base probabilities and SNVQ probabilties for a read
+ class ReadProbs {
+ double[] baseProbs;
+ double[][] snvqProbs; // length of first dimension is flow order length
+ }
+
+ @Override
+ public void onTraversalStart() {
+ super.onTraversalStart();
+ outputWriter = createSAMWriter(output, true);
+ }
+
+ @Override
+ public void closeTool() {
+ super.closeTool();
+ if ( outputWriter != null ) {
+ outputWriter.close();
+ }
+
+ try {
+ if ( aqArgs.debugCollectStatsInto != null )
+ printStats(aqArgs.debugCollectStatsInto);
+ } catch (IOException e) {
+ throw new GATKException("", e);
+ }
+ }
+
+ @Override
+ public void apply(final GATKRead read, final ReferenceContext referenceContext, final FeatureContext featureContext) {
+
+ // include supplementary alignments?
+ if ( read.isSupplementaryAlignment() && !aqArgs.keepSupplementaryAlignments ) {
+ return;
+ }
+
+ // include qc-failed reads?
+ if ( read.failsVendorQualityCheck() && !aqArgs.includeQcFailedReads ) {
+ return;
+ }
+
+ // collect input stats
+ if ( aqArgs.debugCollectStatsInto != null ) {
+ collectInputStats(read);
+ }
+
+ // add SNVQ attributes
+ addBaseQuality(read, getHeaderForReads(), aqArgs.maxPhredScore, fbargs);
+
+ // collect output stats
+ if ( aqArgs.debugCollectStatsInto != null ) {
+ collectOutputStats(read);
+ if ( aqArgs.debugReadName.size() != 0 && aqArgs.debugReadName.contains(read.getName()) ) {
+ dumpOutputRead(read);
+ }
+ }
+
+ // write read to output
+ outputWriter.addRead(read);
+ }
+
+ private void collectInputStats(GATKRead read) {
+ for ( byte q : read.getBaseQualitiesNoCopy() ) {
+ inputQualStats.add(q);
+ }
+ }
+
+ private void collectOutputStats(GATKRead read) {
+ if ( aqArgs.outputQualityAttribute != null ) {
+ if (read.hasAttribute(aqArgs.outputQualityAttribute)) {
+ for (byte q : read.getAttributeAsString(aqArgs.outputQualityAttribute).getBytes()) {
+ outputBQStats.add(SAMUtils.fastqToPhred((char) q));
+ }
+ }
+ } else {
+ for (byte q : read.getBaseQualitiesNoCopy()) {
+ outputBQStats.add(q);
+ }
+ }
+ final FlowBasedReadUtils.ReadGroupInfo rgInfo = FlowBasedReadUtils.getReadGroupInfo(getHeaderForReads(), read);
+ final byte[] bases = read.getBasesNoCopy();
+ final double[] sumP = new double[bases.length];
+ for ( int i = 0 ; i < 4 ; i++ ) {
+ byte altBase = rgInfo.flowOrder.getBytes()[i];
+ String attrValue = read.getAttributeAsString(attrNameForNonCalledBase(altBase));
+ int ofs = 0;
+ for ( byte q : attrValue.getBytes() ) {
+ if ( bases[ofs] != altBase ) {
+ outputQAltStats.add(SAMUtils.fastqToPhred((char)q));
+ } else {
+ outputQCalledStats.add(SAMUtils.fastqToPhred((char)q));
+ }
+ sumP[ofs] += Math.pow(10.0, SAMUtils.fastqToPhred((char)q) / -10.0);
+ ofs++;
+
+ }
+ }
+ for ( double p : sumP ) {
+ outputSumPStats.add(p);
+ }
+ }
+
+ // dump read as a csv for easier analysis
+ private void dumpOutputRead(GATKRead read) {
+
+ try {
+ // open file
+ final String fname = aqArgs.debugCollectStatsInto + "." + read.getName() + ".csv";
+ logger.info("dumping read into: " + fname);
+ final PrintWriter pw = new PrintWriter(fname);
+
+ // write header
+ final StringBuilder hdr = new StringBuilder();
+ hdr.append("pos,base,qual,tp,t0,bq");
+ final FlowBasedReadUtils.ReadGroupInfo rgInfo = FlowBasedReadUtils.getReadGroupInfo(getHeaderForReads(), read);
+ for (int i = 0; i < 4; i++) {
+ hdr.append(",");
+ hdr.append(attrNameForNonCalledBase(rgInfo.flowOrder.charAt(i)));
+ }
+ hdr.append(",qCalled");
+ pw.println(hdr);
+
+ // access data
+ final byte[] bases = read.getBasesNoCopy();
+ final byte[] quals = read.getBaseQualitiesNoCopy();
+ final byte[] tp = read.getAttributeAsByteArray(FlowBasedRead.FLOW_MATRIX_TAG_NAME);
+ final byte[] t0 = read.getAttributeAsByteArray(FlowBasedRead.FLOW_MATRIX_T0_TAG_NAME);
+ final byte[] bq = (aqArgs.outputQualityAttribute != null)
+ ? read.getAttributeAsString(aqArgs.outputQualityAttribute).getBytes()
+ : null;
+ final byte[][] qX = new byte[4][];
+ for (int i = 0; i < 4; i++) {
+ qX[i] = read.getAttributeAsString(attrNameForNonCalledBase(rgInfo.flowOrder.charAt(i))).getBytes();
+ }
+
+ // write lines
+ List line = new LinkedList<>();
+ for (int pos = 0; pos < bases.length; pos++) {
+ line.clear();
+
+ // position
+ line.add(Integer.toString(pos));
+
+ // base, qual
+ line.add(Character.toString(bases[pos]));
+ line.add(Integer.toString(quals[pos]));
+
+ // tp,t0,bq
+ line.add(Integer.toString(tp[pos]));
+ line.add(Integer.toString(SAMUtils.fastqToPhred((char)t0[pos])));
+ if ( bq != null ) {
+ line.add(Integer.toString(SAMUtils.fastqToPhred((char) bq[pos])));
+ }
+
+ // qX
+ int calledIndex = -1;
+ for (int i = 0; i < 4; i++) {
+ line.add(Integer.toString(SAMUtils.fastqToPhred((char)qX[i][pos])));
+ if ( bases[pos] == rgInfo.flowOrder.charAt(i) ) {
+ calledIndex = i;
+ }
+ }
+
+ // qCalled
+ if ( calledIndex >= 0 ) {
+ line.add(Integer.toString(SAMUtils.fastqToPhred((char)qX[calledIndex][pos])));
+ } else {
+ line.add("-1");
+ }
+
+ // write the line
+ pw.println(StringUtils.join(line, ','));
+ }
+
+ // close file
+ pw.close();
+ } catch (IOException e) {
+ throw new GATKException("", e);
+ }
+ }
+
+ private void printStats(final String fname) throws IOException {
+
+ inputQualStats.csvWrite(fname + ".inputQual.csv");
+ outputBQStats.csvWrite(fname + ".outputBQ.csv");
+ outputQAltStats.csvWrite(fname + ".outputQAlt.csv");
+ outputQCalledStats.csvWrite(fname + ".outputQCalled.csv");
+ outputSumPStats.csvWrite(fname + ".outputSumP.csv");
+ }
+
+ static public String attrNameForNonCalledBase(byte nonCalledBase) {
+ return attrNameForNonCalledBase((char)nonCalledBase);
+ }
+
+ static public String attrNameForNonCalledBase(char nonCalledBase) {
+ return "q" + Character.toLowerCase(nonCalledBase);
+ }
+
+ public void addBaseQuality(final GATKRead read, final SAMFileHeader hdr, double maxPhredScore, FlowBasedArgumentCollection fbargs) {
+
+ // take in phred score limit
+ if ( !Double.isNaN(maxPhredScore) ) {
+ maxQualityScore = (int)maxPhredScore;
+ minLikelihoodProbRate = Math.pow(10, -maxPhredScore / 10.0);
+ }
+
+ // convert to a flow base read
+ final FlowBasedReadUtils.ReadGroupInfo rgInfo = FlowBasedReadUtils.getReadGroupInfo(hdr, read);
+ final FlowBasedRead fbRead = new FlowBasedRead(read, rgInfo.flowOrder, rgInfo.maxClass, fbargs);
+ final int flowOrderLength = FlowBasedReadUtils.calcFlowOrderLength(rgInfo.flowOrder);
+
+ // generate base and snvq probabilities for the read
+ final ReadProbs readProbs = generateFlowReadBaseAndSNVQErrorProbabilities(fbRead, flowOrderLength, rgInfo.flowOrder.getBytes());
+
+ // install in read
+ if ( aqArgs.outputQualityAttribute != null ) {
+ read.setAttribute(aqArgs.outputQualityAttribute, new String(convertErrorProbToFastq(readProbs.baseProbs)));
+ } else {
+ read.setBaseQualities(convertErrorProbToPhred(readProbs.baseProbs));
+ }
+ for ( int i = 0 ; i < flowOrderLength ; i++ ) {
+ final String name = AddFlowSNVQuality.attrNameForNonCalledBase(rgInfo.flowOrder.charAt(i));
+ read.setAttribute(name, new String(convertErrorProbToFastq(readProbs.snvqProbs[i])));
+ }
+ }
+
+ // Not using SamUtils function since normally an error probability can not be zero.
+ // still, this method is called to convert base quality as well as snvq, which is computed.
+ // the following check is a safety, in case snvq produces a zero.
+ private char[] convertErrorProbToFastq(double[] errorProb) {
+
+ byte[] phred = convertErrorProbToPhred(errorProb);
+ return SAMUtils.phredToFastq(phred).toCharArray();
+ }
+
+ // Not using SamUtils function since normally an error probability can not be zero.
+ // still, this method is called to convert base quality as well as snvq, which is computed.
+ // the following check is a safety, in case snvq produces a zero.
+ private byte[] convertErrorProbToPhred(double[] errorProb) {
+
+ final byte[] phred = new byte[errorProb.length];
+ for ( int i = 0 ; i < errorProb.length ; i++ ) {
+
+ if ( errorProb[i] == 0 ) {
+ phred[i] = (byte)maxQualityScore;
+ } else {
+ final double p = errorProb[i];
+ phred[i] = (byte)Math.round(-10 * Math.log10(p));
+ }
+ }
+ return phred;
+ }
+
+ /**
+ * generate base and snvq probabilties for a read.
+ *
+ * @param fbRead a flow based read
+ * @param flowOrderLength number of bases in flow order (essentially number of valid base values)
+ * @param flowOrder the flow order itself (which can be the size of flowOrderLength or a repeat of it
+ *
+ * @return an instance of a private class containing the base probabilities as well as the snvq probabilities
+ */
+ private ReadProbs generateFlowReadBaseAndSNVQErrorProbabilities(final FlowBasedRead fbRead, final int flowOrderLength, byte[] flowOrder) {
+
+ /**
+ * access key and error probabilities
+ * for a description of the flow probabilities see {@link FlowBasedRead#flowMatrix}
+ */
+ final int[] key = fbRead.getKey();
+ final double[][] errorProbBands = extractErrorProbBands(fbRead, minLikelihoodProbRate);
+
+ // allocate returned prob arrays
+ final double[] baseProbs = new double[fbRead.getBasesNoCopy().length];
+ final double[][] snvqProbs = new double[flowOrderLength][];
+ for ( int i = 0 ; i < snvqProbs.length ; i++ ) {
+ snvqProbs[i] = new double[baseProbs.length];
+ }
+
+ // loop over hmers via flow key
+ int base = 0;
+ Map allBaseProb0 = new LinkedHashMap<>();
+ Map allBaseProb1 = new LinkedHashMap<>();
+
+ for ( int flow = 0 ; flow < key.length ; flow++ ) {
+ if ( key[flow] != 0 ) {
+
+ // establish initial stat
+ allBaseProb0.clear();
+ allBaseProb1.clear();
+ int flow_i = (flow % flowOrderLength);
+
+ // establish hmer quality
+ final int hmerLength = key[flow];
+ final double[] hmerBaseErrorProbs = generateHmerBaseErrorProbabilities(key, errorProbBands, flow, flowOrderLength, flowOrder, allBaseProb0, allBaseProb1);
+
+ // install value in first byte of the hmer
+ baseProbs[base++] = hmerBaseErrorProbs[0]; // first base, or only base in case of a single base hmer
+ for ( int i = 0 ; i < flowOrderLength ; i++ ) {
+ if ( allBaseProb0.containsKey(flowOrder[i]) ) {
+ snvqProbs[i][base - 1] = allBaseProb0.get(flowOrder[i]);
+ } else if ( i != flow_i ) {
+ snvqProbs[i][base - 1] = minLikelihoodProbRate;
+ }
+ }
+
+ // for hmers longer than 1
+ if ( hmerLength > 1 ) {
+
+ // skip all but last (leave with zero error probability)
+ base += (hmerLength - 2);
+
+ // fill last base from computed error probability
+ baseProbs[base++] = hmerBaseErrorProbs[1]; // last base, if hmer is longer than 1
+
+ for ( int i = 0 ; i < flowOrderLength ; i++ ) {
+ if ( allBaseProb1.containsKey(flowOrder[i]) ) {
+ final double p = allBaseProb1.get(flowOrder[i]);
+ for ( int j = 0 ; j < hmerLength - 1 ; j++ ) {
+ snvqProbs[i][base - 1 - j] = (j == 0) ? p : minLikelihoodProbRate; // all but last get the min prob
+ }
+ } else if ( i != flow_i ) {
+ for ( int j = 0 ; j < hmerLength - 1 ; j++ ) {
+ snvqProbs[i][base - 1 - j] = minLikelihoodProbRate;
+ }
+ }
+ }
+ }
+
+ // override result for the last base with the original hmer error probability
+ if ( base == baseProbs.length ) {
+ baseProbs[base - 1] = errorProbBands[ERROR_PROB_BAND_KEY][flow];
+ }
+ }
+ }
+
+ // adjust probability of called bases so that sum will be 1, also enforce min prob
+ final byte[] bases = fbRead.getBasesNoCopy();
+ for ( int ofs = 0 ; ofs < bases.length ; ofs++ ) {
+
+ // go through alt bases and accumulate p, find out index of called bases (in flow order)
+ final byte calledBase = bases[ofs];
+ double altP = 0;
+ int calledIndex = -1;
+ for (int i = 0; i < flowOrderLength; i++) {
+ if ( calledBase != flowOrder[i] ) {
+ snvqProbs[i][ofs] = Math.max(minLikelihoodProbRate, snvqProbs[i][ofs]);
+ altP += snvqProbs[i][ofs];
+ } else {
+ calledIndex = i;
+ }
+ }
+ if ( calledBase < 0 ) {
+ throw new GATKException(String.format("failed to locate called base %c in flow order %s", (char)calledBase, flowOrder));
+ }
+
+ // install probability in called base slot
+ snvqProbs[calledIndex][ofs] = Math.max(0, 1 - altP);
+
+ // at this point, bq becomes trivial (?)
+ baseProbs[ofs] = 1 - snvqProbs[calledIndex][ofs];
+ }
+
+ // build return value
+ ReadProbs readProbs = new ReadProbs();
+ readProbs.baseProbs = baseProbs;
+ readProbs.snvqProbs = snvqProbs;
+ return readProbs;
+ }
+
+ // extract error probability bands. middle (1) band is the key prob.
+ // lower (0) and high (2) are corresponding to -1 and +1 in hmer lengths
+ private static double[][] extractErrorProbBands(final FlowBasedRead flowRead, final double minValue) {
+
+ // access key
+ final int[] key = flowRead.getKey();
+
+ // allocate result
+ double[][] result = new double[ERROR_PROB_BANDS][];
+ for ( int i = 0 ; i < result.length ; i++ ) {
+ result[i] = new double[key.length];
+ }
+
+ for ( int i = 0 ; i < key.length ; i++ ) {
+
+ // extract key probability
+ result[ERROR_PROB_BAND_KEY][i] = Math.max(flowRead.getProb(i, key[i]), minValue);
+
+ // -1
+ if ( key[i] > 0 ) {
+ result[ERROR_PROB_BAND_1LESS][i] = Math.max(flowRead.getProb(i, key[i] - 1), minValue);
+ } else {
+ result[ERROR_PROB_BAND_1LESS][i] = minValue;
+ }
+
+ // +1
+ if ( key[i] < flowRead.getMaxHmer() ) {
+ result[ERROR_PROB_BAND_1MORE][i] = Math.max(flowRead.getProb(i, key[i] + 1), minValue);
+ } else {
+ result[ERROR_PROB_BAND_1MORE][i] = minValue;
+ }
+ }
+
+ return result;
+ }
+
+ @VisibleForTesting
+ protected double[] generateHmerBaseErrorProbabilities(final int[] key, final double[][] errorProbBands, final int flow,
+ final int flowOrderLength, byte[] flowOrder,
+ Map allBaseProb0, Map allBaseProb1) {
+
+ // result is left/right error probabilities
+ final double[] errorProbs = new double[2];
+ final int hmerLength = key[flow];
+
+ errorProbs[0] = generateSidedHmerBaseErrorProbability(key, errorProbBands, flow, -1, flowOrderLength, flowOrder, allBaseProb0);
+ if ( hmerLength != 1 ) {
+ errorProbs[1] = generateSidedHmerBaseErrorProbability(key, errorProbBands, flow, 1, flowOrderLength, flowOrder, allBaseProb1);
+ }
+
+ return errorProbs;
+ }
+
+ private double generateSidedHmerBaseErrorProbability(final int[] key, final double[][] errorProbBands, final int flow, final int sideIncr,
+ final int flowOrderLength, final byte[] flowOrder, final Map allBaseProb) {
+
+ // create a key slice of the area around the flow/hmer.
+ final int minIndex = Math.max(flow - (flowOrderLength - 1), 0);
+ final int maxIndex = Math.min(flow + (flowOrderLength - 1), key.length - 1);
+ final int[] slice = Arrays.copyOfRange(key, minIndex, maxIndex + 1);
+ final int hmerLength = key[flow];
+
+ // walk the flows towards the side until (and including) the first non-zero key
+ // on hmers of length 1 we walk both sides
+ final class SliceInfo {
+ int[] slice;
+ byte altByte;
+ int sideFlow;
+ }
+ final List slices = new LinkedList<>();
+ final int[] incrs = (hmerLength != 1)
+ ? new int[] { sideIncr }
+ : new int[] { sideIncr, -sideIncr};
+ for (int incr : incrs) {
+ for (int sideFlow = flow + incr; sideFlow >= 0 && sideFlow < key.length; sideFlow += incr) {
+
+ // side flow can no overflow the slice
+ if ( sideFlow < minIndex || sideFlow > maxIndex ) {
+ break;
+ }
+
+ // create a alternative key slice by incrementing sideFlow and decrementing flow
+ final int[] altSlice = Arrays.copyOf(slice, slice.length);
+ altSlice[sideFlow - minIndex] += 1;
+ altSlice[flow - minIndex] -= 1;
+ if ( sliceIsValidForConsideration(altSlice, flowOrderLength) ) {
+ SliceInfo si = new SliceInfo();
+ si.slice = altSlice;
+ si.altByte = flowOrder[sideFlow % flowOrderLength];
+ si.sideFlow = sideFlow;
+ slices.add(si);
+ }
+
+ // is the sideFlow (the first encountered) non-zero? if so, break
+ if (key[sideFlow] != 0) {
+ break;
+ }
+ }
+ }
+
+ // at this point, we have a list of valid slices. figure out the error probability for each of them
+ // and compute the base quality
+ final double keyP = sliceProbs(slice, minIndex, key, errorProbBands, flow, flow)[0];
+ double sumP = keyP;
+ for ( final SliceInfo si : slices ) {
+ final double[] sliceP = sliceProbs(si.slice, minIndex, key, errorProbBands, flow, si.sideFlow);
+ if ( allBaseProb != null ) {
+ allBaseProb.put(si.altByte, getSnvq(sliceP[0], sliceP[1], sliceP[2], aqArgs.snvMode));
+ }
+ sumP += sliceP[0];
+ }
+ final double ep = 1 - (keyP / sumP);
+
+ return ep;
+ }
+
+ static double getSnvq(final double sliceP, final double p1, final double p2, AddFlowSNVQualityArgumentCollection.SnvqModeEnum snvMode) {
+ if ( snvMode == AddFlowSNVQualityArgumentCollection.SnvqModeEnum.Legacy ) {
+ return sliceP;
+ } else if ( snvMode == AddFlowSNVQualityArgumentCollection.SnvqModeEnum.Optimistic ) {
+ return (p1 * p2);
+ } else if ( snvMode == AddFlowSNVQualityArgumentCollection.SnvqModeEnum.Pessimistic ) {
+ return (1 - (1 - p1) * (1 - p2));
+ } else if ( snvMode == AddFlowSNVQualityArgumentCollection.SnvqModeEnum.Geometric ) {
+ return Math.sqrt((p1 * p2) * (1 - (1 - p1) * (1 - p2)));
+ } else {
+ throw new GATKException("unknown snvqMode: " + snvMode);
+ }
+ }
+
+ // compute probability for a slice
+ private static double[] sliceProbs(final int[] slice, final int minIndex, final int[] key, final double[][] errorProbBands,
+ final int flow, final int sideFlow) {
+
+ double accumulatedP = 1.0;
+ double p1 = 0.0;
+ double p2 = 0.0;
+ int key_i = minIndex;
+ for ( int i = 0 ; i < slice.length ; i++, key_i++ ) {
+ final int hmer = key[key_i];
+ final int band;
+ if ( slice[i] == (hmer - 1) ) {
+ band = ERROR_PROB_BAND_1LESS;
+ } else if ( slice[i] == (hmer + 1) ) {
+ band = ERROR_PROB_BAND_1MORE;
+ } else if ( slice[i] == hmer ){
+ band = ERROR_PROB_BAND_KEY;
+ } else {
+ throw new GATKException("slice[i] and hmer are too far apart: " + slice[i] + " " + hmer);
+ }
+ final double p = errorProbBands[band][key_i];
+ accumulatedP *= p;
+
+ // collect p1/p2 (flow and sideFlow probs)
+ if ( key_i == flow ) {
+ p1 = p;
+ }
+ if ( key_i == sideFlow ) {
+ p2 = p;
+ }
+ }
+
+ return new double[] {accumulatedP, p1, p2};
+ }
+
+ static boolean sliceIsValidForConsideration(final int[] slice, final int flowOrderLength) {
+
+ // look for strings of consecutive zeros in length of flowOrderLength - 1
+ int consecutiveZeros = 0;
+ for ( int key : slice ) {
+ if ( key != 0 ) {
+ consecutiveZeros = 0;
+ } else {
+ consecutiveZeros++;
+ if ( consecutiveZeros >= (flowOrderLength - 1) ) {
+ return false;
+ }
+ }
+ }
+
+ // if here, not found -> valid
+ return true;
+ }
+}
+
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQualityArgumentCollection.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQualityArgumentCollection.java
new file mode 100644
index 00000000000..4c20b14e59c
--- /dev/null
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQualityArgumentCollection.java
@@ -0,0 +1,70 @@
+package org.broadinstitute.hellbender.tools.walkers.featuremapping;
+
+import org.broadinstitute.barclay.argparser.Advanced;
+import org.broadinstitute.barclay.argparser.Argument;
+import org.broadinstitute.barclay.argparser.Hidden;
+
+import java.io.Serializable;
+import java.util.List;
+
+/**
+ * Set of arguments for the {@link AddFlowSNVQuality}
+ */
+public class AddFlowSNVQualityArgumentCollection implements Serializable{
+ private static final long serialVersionUID = 1L;
+ public static final String MAX_PHRED_SCORE_FULL_NAME = "max-phred-score";
+ public static final String KEEP_SUPPLEMENTARY_ALIGNMENTS_FULL_NAME = "keep-supplementary-alignments";
+ public static final String INCLUDE_QC_FAILED_READ_FULL_NAME = "include-qc-failed-read";
+ public static final String SNVQ_MODE_FULL_NAME = "snvq-mode";
+ public static final String OUTPUT_QUALITY_ATTRIBUTE_FULL_NAME = "output-quality-attribute";
+ public static final String DEBUG_READ_NAME_FULL_NAME = "debug-read-name";
+ public static final String DEBUG_COLLECT_STATS_INTO_FULL_NAME = "debug-collect-stats-into";
+
+ public enum SnvqModeEnum {
+ Legacy,
+ Optimistic,
+ Pessimistic,
+ Geometric
+ };
+
+ /**
+ * maximum value for
+ * delta in score
+ **/
+ @Argument(fullName = MAX_PHRED_SCORE_FULL_NAME, doc = "Limit value for phred scores", optional = true)
+ public double maxPhredScore = Double.NaN;
+
+ /**
+ * keep supplementary alignments?
+ **/
+ @Argument(fullName = KEEP_SUPPLEMENTARY_ALIGNMENTS_FULL_NAME, doc = "keep supplementary alignments ?", optional = true)
+ public boolean keepSupplementaryAlignments = true;
+
+ @Advanced
+ @Argument(fullName= INCLUDE_QC_FAILED_READ_FULL_NAME, doc = "include reads with QC failed flag", optional = true)
+ public boolean includeQcFailedReads = true;
+
+ /**
+ * snvq computation mode
+ */
+ @Argument(fullName = SNVQ_MODE_FULL_NAME, doc = "snvq calculation mode.", optional = true)
+ public SnvqModeEnum snvMode = SnvqModeEnum.Geometric;
+
+ /**
+ * By default this tool overwrites the QUAL field with the new qualities. Setting this argument saves the original qualities in the specified SAM tag.
+ */
+ @Argument(fullName = OUTPUT_QUALITY_ATTRIBUTE_FULL_NAME, doc = "alternate SAM tag to put original quality scores instead of overwriting the QUAL field. If not used, QUAL will be overwritten.", optional = true)
+ public String outputQualityAttribute = null;
+
+ /**
+ * debug read names?
+ **/
+ @Hidden
+ @Argument(fullName = DEBUG_READ_NAME_FULL_NAME, doc = "Read names of reads to output details of as part of debugging. ", optional = true)
+ public List debugReadName = null;
+
+ @Advanced
+ @Hidden
+ @Argument(fullName= DEBUG_COLLECT_STATS_INTO_FULL_NAME, doc = "Statistics about the reads will be output to given filename.", optional = true)
+ public String debugCollectStatsInto = null;
+}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/FlowFeatureMapper.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/FlowFeatureMapper.java
index b3ddc3e21b6..2014ce3030e 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/FlowFeatureMapper.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/FlowFeatureMapper.java
@@ -2,12 +2,14 @@
import htsjdk.samtools.*;
import htsjdk.samtools.util.Locatable;
+import htsjdk.samtools.util.SequenceUtil;
import htsjdk.variant.variantcontext.Allele;
import htsjdk.variant.variantcontext.VariantContext;
import htsjdk.variant.variantcontext.VariantContextBuilder;
import htsjdk.variant.variantcontext.writer.Options;
import htsjdk.variant.variantcontext.writer.VariantContextWriter;
import htsjdk.variant.vcf.*;
+import org.apache.commons.lang3.StringUtils;
import org.apache.commons.math3.util.Precision;
import org.apache.logging.log4j.LogManager;
import org.apache.logging.log4j.Logger;
@@ -81,6 +83,19 @@
@ExperimentalFeature
public final class FlowFeatureMapper extends ReadWalker {
+ static class CopyAttrInfo {
+ public final String name;
+ public final VCFHeaderLineType type;
+ public final String desc;
+
+ public CopyAttrInfo(final String spec) {
+ final String[] toks = spec.split(",");
+ name = toks[0];
+ type = toks.length > 1 ? VCFHeaderLineType.valueOf(toks[1]) : VCFHeaderLineType.String;
+ desc = toks.length > 2 ? StringUtils.join(Arrays.copyOfRange(toks, 2, toks.length), ",") : ("copy-attr: " + name);
+ }
+ }
+
private static final Logger logger = LogManager.getLogger(FlowFeatureMapper.class);
private static final String VCB_SOURCE = "fm";
@@ -101,8 +116,18 @@ public final class FlowFeatureMapper extends ReadWalker {
private static final String VCF_SMQ_RIGHT = "X_SMQ_RIGHT";
private static final String VCF_SMQ_LEFT_MEAN = "X_SMQ_LEFT_MEAN";
private static final String VCF_SMQ_RIGHT_MEAN = "X_SMQ_RIGHT_MEAN";
+ private static final String VCF_ADJACENT_REF_DIFF = "X_ADJACENT_REF_DIFF";
+
private static final Double LOWEST_PROB = 0.0001;
+ private static final int VENDOR_QUALITY_CHECK_FLAG = 0x200;
+
+ private static final String INCLUDE_QC_FAILED_READS_FULL_NAME = "include-qc-failed-reads";
+
+ final private List copyAttrInfo = new LinkedList<>();
+
+ // order here is according to SequenceUtil.VALID_BASES_UPPER
+ final private String scoreForBaseNames[] = new String[SequenceUtil.VALID_BASES_UPPER.length];
@Argument(fullName = StandardArgumentDefinitions.OUTPUT_LONG_NAME,
shortName = StandardArgumentDefinitions.OUTPUT_SHORT_NAME,
@@ -131,6 +156,10 @@ public final class FlowFeatureMapper extends ReadWalker {
@Argument(fullName=HaplotypeCallerArgumentCollection.OUTPUT_BLOCK_LOWER_BOUNDS, doc = "Output the band lower bound for each GQ block regardless of the data it represents", optional = true)
public boolean floorBlocks = false;
+ @Advanced
+ @Argument(fullName=INCLUDE_QC_FAILED_READS_FULL_NAME, doc = "include reads with QC failed flag", optional = true)
+ public boolean includeQcFailedReads = false;
+
@ArgumentCollection
public FlowBasedArgumentCollection fbargs = new FlowBasedArgumentCollection();
@@ -175,6 +204,9 @@ protected static class MappedFeature implements Comparable {
int smqLeftMean;
int smqRightMean;
+ double scoreForBase[];
+ boolean adjacentRefDiff;
+
public MappedFeature(GATKRead read, FlowFeatureMapperArgumentCollection.MappingFeatureEnum type, byte[] readBases,
byte[] refBases, int readBasesOffset, int start, int offsetDelta) {
this.read = read;
@@ -315,8 +347,20 @@ public VCFHeader makeVCFHeader(final SAMSequenceDictionary sequenceDictionary, f
headerInfo.add(new VCFInfoHeaderLine(VCF_SMQ_RIGHT, 1, VCFHeaderLineType.Integer, "Ordinal Median quality of N bases to the right of the feature"));
headerInfo.add(new VCFInfoHeaderLine(VCF_SMQ_LEFT_MEAN, 1, VCFHeaderLineType.Integer, "Ordinal Mean quality of N bases to the left of the feature"));
headerInfo.add(new VCFInfoHeaderLine(VCF_SMQ_RIGHT_MEAN, 1, VCFHeaderLineType.Integer, "Ordinal Mean quality of N bases to the right of the feature"));
- for ( String name : fmArgs.copyAttr ) {
- headerInfo.add(new VCFInfoHeaderLine(fmArgs.copyAttrPrefix + name, 1, VCFHeaderLineType.String, "copy-attr: " + name));
+ for ( String spec : fmArgs.copyAttr ) {
+ final CopyAttrInfo info = new CopyAttrInfo(spec);
+ headerInfo.add(new VCFInfoHeaderLine(fmArgs.copyAttrPrefix + info.name, 1, info.type, info.desc));
+ copyAttrInfo.add(info);
+ }
+
+ // validation mode?
+ if ( fmArgs.reportAllAlts ) {
+ for ( int baseIndex = 0 ; baseIndex < SequenceUtil.VALID_BASES_UPPER.length ; baseIndex++ ) {
+ headerInfo.add(new VCFInfoHeaderLine(scoreNameForBase(baseIndex), 1, VCFHeaderLineType.Float, "Base specific mapping score"));
+ }
+ }
+ if ( fmArgs.reportAllAlts || fmArgs.tagBasesWithAdjacentRefDiff ) {
+ headerInfo.add(new VCFInfoHeaderLine(VCF_ADJACENT_REF_DIFF, 1, VCFHeaderLineType.Flag, "Adjacent base filter indication: indel in the adjacent 5 bases to the considered base on the read"));
}
final VCFHeader vcfHeader = new VCFHeader(headerInfo);
@@ -324,6 +368,13 @@ public VCFHeader makeVCFHeader(final SAMSequenceDictionary sequenceDictionary, f
return vcfHeader;
}
+ private String scoreNameForBase(int baseIndex) {
+ if ( scoreForBaseNames[baseIndex] == null ) {
+ scoreForBaseNames[baseIndex] = VCF_SCORE + "_" + new String(new byte[]{SequenceUtil.VALID_BASES_UPPER[baseIndex]});
+ }
+ return scoreForBaseNames[baseIndex];
+ }
+
@Override
public void apply(final GATKRead read, final ReferenceContext referenceContext, final FeatureContext featureContext) {
@@ -337,6 +388,11 @@ public void apply(final GATKRead read, final ReferenceContext referenceContext,
return;
}
+ // include qc-failed reads?
+ if ( ((read.getFlags() & VENDOR_QUALITY_CHECK_FLAG) != 0) && !includeQcFailedReads ) {
+ return;
+ }
+
// flush qeues up to this read
flushQueue(read, referenceContext);
@@ -349,6 +405,20 @@ public void apply(final GATKRead read, final ReferenceContext referenceContext,
// score the feature
fr.score = scoreFeature(fr);
+ // score for validation mode?
+ if ( fmArgs.reportAllAlts) {
+ fr.scoreForBase = new double[SequenceUtil.VALID_BASES_UPPER.length];
+ for ( int baseIndex = 0 ; baseIndex < fr.scoreForBase.length ; baseIndex++ ) {
+ final byte base = SequenceUtil.VALID_BASES_UPPER[baseIndex];
+ boolean incl = base != fr.readBases[0];
+ if ( incl ) {
+ fr.scoreForBase[baseIndex] = scoreFeature(fr, base);
+ } else {
+ fr.scoreForBase[baseIndex] = Double.NaN;
+ }
+ }
+ }
+
// emit feature if filters in
if ( filterFeature(fr) ) {
featureQueue.add(fr);
@@ -413,14 +483,17 @@ private void enrichFeature(final MappedFeature fr) {
}
private double scoreFeature(final MappedFeature fr) {
+ return scoreFeature(fr, (byte)0);
+ }
+ private double scoreFeature(final MappedFeature fr, byte altBase) {
// build haplotypes
final FlowBasedReadUtils.ReadGroupInfo rgInfo = FlowBasedReadUtils.getReadGroupInfo(getHeaderForReads(), fr.read);
- final FlowBasedHaplotype[] haplotypes = buildHaplotypes(fr, rgInfo.flowOrder);
+ final FlowBasedHaplotype[] haplotypes = buildHaplotypes(fr, rgInfo.flowOrder, altBase);
// create flow read
- final FlowBasedRead flowRead = new FlowBasedRead(fr.read, rgInfo.flowOrder,
- rgInfo.maxClass, fbargs);
+ final FlowBasedRead flowRead = new FlowBasedRead(fr.read, rgInfo.flowOrder, rgInfo.maxClass, fbargs);
+
final int diffLeft = haplotypes[0].getStart() - flowRead.getStart() + fr.offsetDelta;
final int diffRight = flowRead.getEnd() - haplotypes[0].getEnd();
flowRead.applyBaseClipping(Math.max(0, diffLeft), Math.max(diffRight, 0), false);
@@ -524,16 +597,24 @@ public static double computeLikelihoodLocal(final FlowBasedRead read, final Flow
return result;
}
- private FlowBasedHaplotype[] buildHaplotypes(final MappedFeature fr, final String flowOrder) {
+ private FlowBasedHaplotype[] buildHaplotypes(final MappedFeature fr, final String flowOrder, byte altBase) {
// build bases for flow haplotypes
// NOTE!!!: this code assumes length of feature on read and reference is the same
// this is true for SNP but not for INDELs - it will have to be re-written!
// TODO: write for INDEL
- final byte[] bases = fr.read.getBasesNoCopy();
+ byte[] bases = fr.read.getBasesNoCopy();
int offset = fr.readBasesOffset;
int refStart = fr.start;
int refModOfs = 0;
+
+ // install alt base?
+ byte orgBase = 0;
+ if ( altBase != 0 ) {
+ orgBase = fr.refBases[0];
+ fr.refBases[0] = altBase;
+ }
+
if ( offset > 0 ) {
// reach into hmer before
offset--;
@@ -569,6 +650,11 @@ private FlowBasedHaplotype[] buildHaplotypes(final MappedFeature fr, final Strin
new FlowBasedHaplotype(refHaplotype, flowOrder)
};
+ // restore changes
+ if ( altBase != 0 ) {
+ fr.refBases[0] = orgBase;
+ }
+
// return
return result;
}
@@ -590,7 +676,11 @@ private void emitFeature(final MappedFeature fr) {
// create alleles
final Collection alleles = new LinkedList<>();
- alleles.add(Allele.create(fr.readBases, false));
+ if ( fmArgs.reportAllAlts && Arrays.equals(fr.readBases, fr.refBases) ) {
+ alleles.add(Allele.create("*".getBytes(), false));
+ } else {
+ alleles.add(Allele.create(fr.readBases, false));
+ }
alleles.add(Allele.create(fr.refBases, true));
// create variant context builder
@@ -625,11 +715,34 @@ private void emitFeature(final MappedFeature fr) {
vcb.attribute(VCF_SMQ_RIGHT_MEAN, fr.smqRightMean);
}
- for ( String name : fmArgs.copyAttr ) {
- if ( fr.read.hasAttribute(name) ) {
- vcb.attribute(fmArgs.copyAttrPrefix + name, fr.read.getAttributeAsString(name));
+ for ( CopyAttrInfo info : copyAttrInfo ) {
+ if ( fr.read.hasAttribute(info.name) ) {
+ final String attrName = fmArgs.copyAttrPrefix + info.name;
+ if ( info.type == VCFHeaderLineType.Integer ) {
+ vcb.attribute(attrName, fr.read.getAttributeAsInteger(info.name));
+ } else if ( info.type == VCFHeaderLineType.Float ) {
+ vcb.attribute(attrName, fr.read.getAttributeAsFloat(info.name));
+ } else {
+ vcb.attribute(attrName, fr.read.getAttributeAsString(info.name));
+ }
+ }
+ }
+
+ // validation mode?
+ if ( fmArgs.reportAllAlts ) {
+ if ( fr.scoreForBase != null ) {
+ for (int baseIndex = 0; baseIndex < SequenceUtil.VALID_BASES_UPPER.length; baseIndex++) {
+ if (!Double.isNaN(fr.scoreForBase[baseIndex])) {
+ vcb.attribute(scoreNameForBase(baseIndex), String.format("%.5f", fr.scoreForBase[baseIndex]));
+ }
+ }
}
}
+ if ( fr.adjacentRefDiff ) {
+ vcb.attribute(VCF_ADJACENT_REF_DIFF, true);
+ }
+
+ // build it!
final VariantContext vc = vcb.make();
// write to file
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/FlowFeatureMapperArgumentCollection.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/FlowFeatureMapperArgumentCollection.java
index 22c13b00878..03660746d2b 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/FlowFeatureMapperArgumentCollection.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/FlowFeatureMapperArgumentCollection.java
@@ -32,7 +32,7 @@ enum MappingFeatureEnum {
/**
* attributes to copy from bam
**/
- @Argument(fullName = "copy-attr", doc = "attributes to copy from bam", optional = true)
+ @Argument(fullName = "copy-attr", doc = "attributes to copy from bam. format ,,. types: Integer, Float, String, Character, Flag", optional = true)
public List copyAttr = new LinkedList<>();
/**
@@ -116,4 +116,18 @@ enum MappingFeatureEnum {
@Hidden
@Argument(fullName = "surrounding-mean-quality-size", doc = "number of bases around the feature to calculate surrounding mean quality", optional = true)
public Integer surroundingMeanQualitySize = null;
+
+ /**
+ * validation mode - if not specified, this feature is off
+ **/
+ @Hidden
+ @Argument(fullName = "report-all-alts", doc = "In this mode (aka validation mode), every base of every read in the input CRAM and interval is reported, and an X_SCORE value is calculated for all 3 possible alts", optional = true)
+ public boolean reportAllAlts = false;
+
+ /**
+ * adjacent-ref-diff mode - if not specified, this feature is off
+ **/
+ @Hidden
+ @Argument(fullName = "tag-bases-with-adjacent-ref-diff", doc = "In this mode bases that have an adjacent difference from the reference on the same read are not discarded, and tagged with X_ADJACENT_REF_DIFFm", optional = true)
+ public boolean tagBasesWithAdjacentRefDiff = false;
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/SNVMapper.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/SNVMapper.java
index 05c952c317a..3f23cb5a6e6 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/SNVMapper.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/SNVMapper.java
@@ -1,7 +1,7 @@
package org.broadinstitute.hellbender.tools.walkers.featuremapping;
import htsjdk.samtools.CigarElement;
-import org.apache.commons.lang.ArrayUtils;
+import org.apache.commons.lang3.ArrayUtils;
import org.apache.commons.text.similarity.LevenshteinDistance;
import org.broadinstitute.hellbender.engine.ReferenceContext;
import org.broadinstitute.hellbender.exceptions.GATKException;
@@ -21,22 +21,34 @@
public class SNVMapper implements FeatureMapper {
- final int identBefore;
- final int identAfter;
+ final int surroundBefore;
+ final int surroundAfter;
final int minCigarElementLength;
final LevenshteinDistance levDistance = new LevenshteinDistance();
final Integer smqSize;
final Integer smqSizeMean;
+ final boolean ignoreSurround;
+ final int spanBefore;
+ final int spanAfter;
+
+ final FlowFeatureMapperArgumentCollection fmArgs;
+
public SNVMapper(FlowFeatureMapperArgumentCollection fmArgs) {
- identBefore = fmArgs.snvIdenticalBases;
- identAfter = (fmArgs.snvIdenticalBasesAfter != 0) ? fmArgs.snvIdenticalBasesAfter : identBefore;
- minCigarElementLength = identBefore + 1 + identAfter;
+ surroundBefore = fmArgs.snvIdenticalBases;
+ surroundAfter = (fmArgs.snvIdenticalBasesAfter != 0) ? fmArgs.snvIdenticalBasesAfter : surroundBefore;
smqSize = fmArgs.surroundingMediaQualitySize;
smqSizeMean = fmArgs.surroundingMeanQualitySize;
+ this.fmArgs = fmArgs;
+
+ // ignore surround
+ ignoreSurround = fmArgs.reportAllAlts || fmArgs.tagBasesWithAdjacentRefDiff;
+ spanBefore = ignoreSurround ? 0 : surroundBefore;
+ spanAfter = ignoreSurround ? 0 : surroundAfter;
+ minCigarElementLength = spanBefore + 1 + spanAfter;
// adjust minimal read length
- FlowBasedRead.setMinimalReadLength(1 + 1 + identAfter);
+ FlowBasedRead.setMinimalReadLength(1 + 1 + spanAfter);
}
@Override
@@ -93,30 +105,44 @@ public void forEachOnRead(GATKRead read, ReferenceContext referenceContext, Cons
if ( length >= minCigarElementLength &&
cigarElement.getOperator().consumesReadBases() &&
cigarElement.getOperator().consumesReferenceBases() ) {
- readOfs += identBefore;
- refOfs += identBefore;
- for ( int ofs = identBefore ; ofs < length - identAfter ; ofs++, readOfs++, refOfs++ ) {
+ readOfs += spanBefore;
+ refOfs += spanBefore;
+ for ( int ofs = spanBefore ; ofs < length - spanAfter ; ofs++, readOfs++, refOfs++ ) {
- if ( ref[refOfs] != 'N' && bases[readOfs] != ref[refOfs] ) {
+ if ( ref[refOfs] != 'N' && (fmArgs.reportAllAlts || (bases[readOfs] != ref[refOfs])) ) {
// check that this is really a SNV (must be surrounded by identical ref)
boolean surrounded = true;
- for ( int i = 0 ; i < identBefore && surrounded ; i++ ) {
- if ( bases[readOfs-1-i] != ref[refOfs-1-i] ) {
+ for ( int i = 0 ; i < surroundBefore && surrounded ; i++ ) {
+ final int bIndex = readOfs-1-i;
+ final int rIndex = refOfs-1-i;
+ if ( bIndex < 0 || bIndex >= bases.length || rIndex < 0 || rIndex >= ref.length ) {
+ surrounded = false;
+ continue;
+ }
+ if ( bases[bIndex] != ref[rIndex] ) {
surrounded = false;
}
}
- for ( int i = 0 ; i < identAfter && surrounded ; i++ ) {
- if ( bases[readOfs+1+i] != ref[refOfs+1+i] ) {
+ for (int i = 0; i < surroundAfter && surrounded ; i++ ) {
+ final int bIndex = readOfs+1+i;
+ final int rIndex = refOfs+1+i;
+ if ( bIndex < 0 || bIndex >= bases.length || rIndex < 0 || rIndex >= ref.length ) {
+ surrounded = false;
+ continue;
+ }
+ if ( bases[bIndex] != ref[rIndex] ) {
surrounded = false;
}
}
- if ( !surrounded ) {
+ if ( (!fmArgs.reportAllAlts && !fmArgs.tagBasesWithAdjacentRefDiff) && !surrounded ) {
continue;
}
// add this feature
FlowFeatureMapper.MappedFeature feature = FlowFeatureMapper.MappedFeature.makeSNV(read, readOfs, ref[refOfs], referenceContext.getStart() + refOfs, readOfs - refOfs);
+ if ( (fmArgs.reportAllAlts || fmArgs.tagBasesWithAdjacentRefDiff) && !surrounded )
+ feature.adjacentRefDiff = true;
feature.nonIdentMBasesOnRead = nonIdentMBases;
feature.refEditDistance = refEditDistance;
if ( !read.isReverseStrand() )
@@ -166,8 +192,8 @@ public void forEachOnRead(GATKRead read, ReferenceContext referenceContext, Cons
features.add(feature);
}
}
- readOfs += identAfter;
- refOfs += identAfter;
+ readOfs += spanAfter;
+ refOfs += spanAfter;
} else {
@@ -243,8 +269,8 @@ public FilterStatus noFeatureButFilterAt(GATKRead read, ReferenceContext referen
cigarElement.getOperator().consumesReferenceBases() ) {
// break out if not enough clearing
- if ( (start < referenceContext.getStart() + refOfs + identBefore) ||
- (start >= referenceContext.getStart() + refOfs + length - identAfter) )
+ if ( (start < referenceContext.getStart() + refOfs + spanBefore) ||
+ (start >= referenceContext.getStart() + refOfs + length - spanAfter) )
return FilterStatus.Filtered;
int delta = start - (referenceContext.getStart() + refOfs);
@@ -255,13 +281,25 @@ public FilterStatus noFeatureButFilterAt(GATKRead read, ReferenceContext referen
// check that this is really a SNV (must be surrounded by identical ref)
boolean surrounded = true;
- for ( int i = 0 ; i < identBefore && surrounded ; i++ ) {
- if ( bases[readOfs-1-i] != ref[refOfs-1-i] ) {
+ for ( int i = 0 ; i < surroundBefore && surrounded ; i++ ) {
+ final int bIndex = readOfs-1-i;
+ final int rIndex = refOfs-1-i;
+ if ( bIndex < 0 || bIndex >= bases.length || rIndex < 0 || rIndex >= ref.length ) {
+ surrounded = false;
+ continue;
+ }
+ if ( bases[bIndex] != ref[rIndex] ) {
surrounded = false;
}
}
- for ( int i = 0 ; i < identAfter && surrounded ; i++ ) {
- if ( bases[readOfs+1+i] != ref[refOfs+1+i] ) {
+ for (int i = 0; i < surroundAfter && surrounded ; i++ ) {
+ final int bIndex = readOfs+1+i;
+ final int rIndex = refOfs+1+i;
+ if ( bIndex < 0 || bIndex >= bases.length || rIndex < 0 || rIndex >= ref.length ) {
+ surrounded = false;
+ continue;
+ }
+ if ( bases[bIndex] != ref[rIndex] ) {
surrounded = false;
}
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/filters/VariantFiltration.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/filters/VariantFiltration.java
index dd6a2e93838..a2642f20150 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/filters/VariantFiltration.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/filters/VariantFiltration.java
@@ -114,6 +114,7 @@ public final class VariantFiltration extends VariantWalker {
public static final String CLUSTER_WINDOW_SIZE_LONG_NAME = "cluster-window-size";
public static final String MASK_EXTENSION_LONG_NAME = "mask-extension";
public static final String MASK_NAME_LONG_NAME = "mask-name";
+ public static final String MASK_DESCRIPTION_LONG_NAME = "mask-description";
public static final String FILTER_NOT_IN_MASK_LONG_NAME = "filter-not-in-mask";
public static final String MISSING_VAL_LONG_NAME = "missing-values-evaluate-as-failing";
public static final String INVERT_LONG_NAME = "invert-filter-expression";
@@ -238,6 +239,14 @@ public final class VariantFiltration extends VariantWalker {
@Argument(fullName=ALLELE_SPECIFIC_LONG_NAME, optional=true, doc="Set mask at the allele level. This option is not compatible with clustering.")
public boolean applyForAllele = false;
+ /**
+ * If a mask interval list is provided, then set the description of the filter in the VCF header to this String.
+ * Note that if spaces are needed, then the entire description should be enclosed in quotes. Also note that if
+ * --filter-not-in-mask is used, the description should be adapted to reflect the reverse logic.
+ */
+ @Argument(fullName=MASK_DESCRIPTION_LONG_NAME, optional=true, doc="Description to add to the FILTER field in VCF header for the mask filter.")
+ public String maskDescription;
+
// JEXL expressions for the filters
private List filterExps;
private List genotypeFilterExps;
@@ -305,7 +314,9 @@ private void initializeVcfWriter() {
}
if ( mask != null ) {
- if (filterRecordsNotInMask) {
+ if (maskDescription != null) {
+ hInfo.add(new VCFFilterHeaderLine(maskName, maskDescription));
+ } else if (filterRecordsNotInMask) {
hInfo.add(new VCFFilterHeaderLine(maskName, "Doesn't overlap a user-input mask"));
} else {
hInfo.add(new VCFFilterHeaderLine(maskName, "Overlaps a user-input mask"));
@@ -331,6 +342,9 @@ public void onTraversalStart() {
if (filterRecordsNotInMask && mask == null) {
throw new CommandLineException.BadArgumentValue(FILTER_NOT_IN_MASK_LONG_NAME, "argument not allowed if mask argument is not provided");
}
+ if (maskDescription != null && mask == null) {
+ throw new CommandLineException.BadArgumentValue(MASK_DESCRIPTION_LONG_NAME, "argument not allowed if mask argument is not provided");
+ }
filterExps = VariantContextUtils.initializeMatchExps(filterNames, filterExpressions);
genotypeFilterExps = VariantContextUtils.initializeMatchExps(genotypeFilterNames, genotypeFilterExpressions);
howToTreatMissingValues = failMissingValues ? JexlMissingValueTreatment.TREAT_AS_MATCH : JexlMissingValueTreatment.TREAT_AS_MISMATCH;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java
index 60f1ef7b205..aede6c33fc2 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java
@@ -130,7 +130,7 @@ public static GenotypesContext subsetAlleles(final GenotypesContext originalGs,
}
gb.PL(newLikelihoods);
- GATKVariantContextUtils.makeGenotypeCall(g.getPloidy(), gb, assignmentMethod, newLikelihoods, allelesToKeep, g.getAlleles(), gpc);
+ GATKVariantContextUtils.makeGenotypeCall(g.getPloidy(), gb, assignmentMethod, newLikelihoods, allelesToKeep, g, gpc);
// restrict SAC to the new allele subset
if (g.hasExtendedAttribute(GATKVCFConstants.STRAND_COUNT_BY_SAMPLE_KEY)) {
@@ -139,7 +139,7 @@ public static GenotypesContext subsetAlleles(final GenotypesContext originalGs,
}
// restrict AD to the new allele subset
- if(g.hasAD()) {
+ if(g.hasAD() && gb.makeWithShallowCopy().hasAD()) {
final int[] newAD = getNewAlleleBasedReadCountAnnotation(allelesToKeep, allelePermutation, g.getAD());
gb.AD(newAD);
// if we have recalculated AD and the original genotype had AF but was then removed, then recalculate AF based on AD counts
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypeAssignmentMethod.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypeAssignmentMethod.java
index bbc9cb888aa..eaa9bad7033 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypeAssignmentMethod.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/GenotypeAssignmentMethod.java
@@ -56,6 +56,7 @@ public enum GenotypeAssignmentMethod {
/**
* Use PLs unless they are unavailable, in which case use best match to original
+ * GQ0 hom-refs will be converted to no-calls
*/
PREFER_PLS
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/gnarlyGenotyper/GnarlyGenotyperEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/gnarlyGenotyper/GnarlyGenotyperEngine.java
index 7efb5687ac3..9b6c5030497 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/gnarlyGenotyper/GnarlyGenotyperEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/gnarlyGenotyper/GnarlyGenotyperEngine.java
@@ -4,7 +4,7 @@
import com.google.common.primitives.Ints;
import htsjdk.variant.variantcontext.*;
import htsjdk.variant.vcf.VCFConstants;
-import org.apache.commons.lang.StringUtils;
+import org.apache.commons.lang3.StringUtils;
import org.apache.commons.lang3.tuple.Pair;
import org.broadinstitute.hellbender.exceptions.UserException;
import org.broadinstitute.hellbender.tools.walkers.annotator.*;
@@ -349,6 +349,8 @@ protected GenotypesContext iterateOnGenotypes(final VariantContext vc, final Lis
//If GenomicsDB returns no-call genotypes like CombineGVCFs (depending on the GenomicsDBExportConfiguration),
// then we need to actually find the GT from PLs
makeGenotypeCall(g, genotypeBuilder, GenotypeLikelihoods.fromPLs(PLs).getAsVector(), targetAlleles);
+ } else if (g.hasGQ() && g.getGQ() == 0) {
+ makeGenotypeCall(g, genotypeBuilder, null, targetAlleles); //null likelihoods for reblocked hom-ref that we want to no-call
}
final Map attrs = new HashMap<>(g.getExtendedAttributes());
attrs.remove(GATKVCFConstants.MIN_DP_FORMAT_KEY);
@@ -391,14 +393,8 @@ protected GenotypesContext iterateOnGenotypes(final VariantContext vc, final Lis
protected void makeGenotypeCall(final Genotype g, final GenotypeBuilder gb,
final double[] genotypeLikelihoods,
final List allelesToUse) {
- if ( genotypeLikelihoods == null || !GATKVariantContextUtils.isInformative(genotypeLikelihoods) ) {
- //gb.alleles(GATKVariantContextUtils.noCallAlleles(g.getAlleles().size())).noGQ();
- GATKVariantContextUtils.makeGenotypeCall(g.getPloidy(), gb, GenotypeAssignmentMethod.SET_TO_NO_CALL,
- genotypeLikelihoods, allelesToUse, null);
- } else {
- GATKVariantContextUtils.makeGenotypeCall(g.getPloidy(), gb, GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN,
- genotypeLikelihoods, allelesToUse, null);
- }
+ GATKVariantContextUtils.makeGenotypeCall(g.getPloidy(), gb, GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN,
+ genotypeLikelihoods, allelesToUse, g, null);
}
/**
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/AddFlowBaseQuality.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/AddFlowBaseQuality.java
index a0f804ffc11..bd8cf62e851 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/AddFlowBaseQuality.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/AddFlowBaseQuality.java
@@ -239,6 +239,25 @@ private static double[][] extractErrorProbBands(final FlowBasedRead flowRead, fi
return result;
}
+ /**
+ * The following functions estimate the error probability for an hmer. specifically two error
+ * probability values are generated: one for the first base of the hmer and another for the
+ * rest of its bases.
+ *
+ * The computation itself is performed in a subsequent function: generateSidedHmerBaseErrorProbability
+ * It iterates over the possible valid combinations of errors and sums them up.
+ *
+ * @param key - key (hmer length) in flow space
+ * @param errorProbBands - for each flow (position in the key) three error probabilities are provided:
+ * [0] - for the hmer being one base shorter
+ * [1] - for the hmer to be at its length
+ * [2] - for the hmer to be one base longer
+ * @param flow - the flow (index) for which to generate the probabilities (0 <= flow < key.length)
+ * @param flowOrderLength - the cycle length of of the flow order (usually 4)
+ * @return an array of two probabilities:
+ * [0] - probability for the first base of the hmer
+ * [1] - probability for the rest of the bases of the hmer
+ */
@VisibleForTesting
protected static double[] generateHmerBaseErrorProbabilities(final int[] key, final double[][] errorProbBands, final int flow, final int flowOrderLength) {
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/GroundTruthReadsBuilder.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/GroundTruthReadsBuilder.java
index cee262d277c..74434c89cbe 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/GroundTruthReadsBuilder.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/GroundTruthReadsBuilder.java
@@ -5,7 +5,7 @@
import htsjdk.samtools.util.Locatable;
import htsjdk.samtools.util.SequenceUtil;
import htsjdk.samtools.util.Tuple;
-import org.apache.commons.lang.ArrayUtils;
+import org.apache.commons.lang3.ArrayUtils;
import org.apache.commons.lang3.StringUtils;
import org.apache.logging.log4j.LogManager;
import org.apache.logging.log4j.Logger;
@@ -858,9 +858,9 @@ else if ( hasQ || hasZ ) {
cols.put("ReadName", read.getName());
// haplotypes and reference scores
- cols.put("PaternalHaplotypeScore", paternal.score);
- cols.put("MaternalHaplotypeScore", maternal.score);
- cols.put("RefHaplotypeScore", refScore);
+ cols.put("PaternalHaplotypeScore", String.format("%.6f", paternal.score));
+ cols.put("MaternalHaplotypeScore", String.format("%.6f", maternal.score));
+ cols.put("RefHaplotypeScore", String.format("%.6f", refScore));
// build haplotype keys
final FlowBasedReadUtils.ReadGroupInfo rgInfo = FlowBasedReadUtils.getReadGroupInfo(getHeaderForReads(), read);
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/GroundTruthScorer.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/GroundTruthScorer.java
new file mode 100644
index 00000000000..3cee7b801bb
--- /dev/null
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/GroundTruthScorer.java
@@ -0,0 +1,926 @@
+package org.broadinstitute.hellbender.tools.walkers.groundtruth;
+
+import com.opencsv.CSVReader;
+import htsjdk.samtools.CigarOperator;
+import htsjdk.samtools.SAMFileHeader;
+import htsjdk.samtools.SAMReadGroupRecord;
+import htsjdk.variant.variantcontext.VariantContext;
+import org.apache.commons.lang3.StringUtils;
+import org.apache.commons.math3.util.Precision;
+import org.apache.logging.log4j.LogManager;
+import org.apache.logging.log4j.Logger;
+import org.broadinstitute.barclay.argparser.Argument;
+import org.broadinstitute.barclay.argparser.ArgumentCollection;
+import org.broadinstitute.barclay.argparser.CommandLineProgramProperties;
+import org.broadinstitute.barclay.argparser.ExperimentalFeature;
+import org.broadinstitute.barclay.help.DocumentedFeature;
+import org.broadinstitute.hellbender.cmdline.programgroups.FlowBasedProgramGroup;
+import org.broadinstitute.hellbender.engine.*;
+import org.broadinstitute.hellbender.exceptions.GATKException;
+import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.FlowBasedAlignmentArgumentCollection;
+import org.broadinstitute.hellbender.tools.walkers.featuremapping.FlowFeatureMapper;
+import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.AssemblyBasedCallerUtils;
+import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.FlowBasedAlignmentLikelihoodEngine;
+import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.LikelihoodEngineArgumentCollection;
+import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.ReadLikelihoodCalculationEngine;
+import org.broadinstitute.hellbender.utils.BaseUtils;
+import org.broadinstitute.hellbender.utils.SimpleInterval;
+import org.broadinstitute.hellbender.utils.clipping.ReadClipper;
+import org.broadinstitute.hellbender.utils.haplotype.FlowBasedHaplotype;
+import org.broadinstitute.hellbender.utils.haplotype.Haplotype;
+import org.broadinstitute.hellbender.utils.read.*;
+import org.broadinstitute.hellbender.utils.report.GATKReport;
+import org.broadinstitute.hellbender.utils.report.GATKReportTable;
+
+import java.io.*;
+import java.text.DecimalFormat;
+import java.util.*;
+import java.util.zip.GZIPOutputStream;
+
+/**
+ * Converts Ultima reads into flow-based annotation, and provides some general statistics regarding
+ * quality and errors relative to the reference. Ultima data is flow-based, and thus the original computed
+ * quality refers to each flow, rather than each base. In the Ultima cram/bam, there is a per-base representation
+ * of the original flow qualities, where the original quality is distributed along each flow (homopolymer).
+ * In order to reconstitute the original flow information, the tool incorporates the information encoded
+ * in the Ultima cram/bam, and outputs both the read in flow space, as well as a conversion of the aligned
+ * reference portion into flow space, and an alignment score.
+ *
+ * Input
+ *
+ * - Ultima aligned SAM/BAM/CRAM
+ *
+ *
+ * Output
+ *
+ * - Per read ground truth information CSV and a ground truth scoring quality report, in GATK report format
+ *
+ *
+ * CSV Output Description
+ * csv with the read representation in flow space. The csv includes the following columns:
+ * ReadName
+ * ReadKey : The signal of the read at each flow according to the flow order
+ * ReadIsReversed : Whether the read is reversed in the alignment
+ * ReadMQ : The mapping quality of the read
+ * ReadRQ : The read rq value
+ * GroundTruthKey : The aligned reference section, translated into per-flow signals
+ * ReadSequence
+ * Score : A flow-based alignment score. Since the alignment is per-flow, in the case that there’s a cycle skip, the read and reference flow signals will not be aligned, and therefore the score will be inaccurate.
+ * NormalizedScore: A flow-based normalized alignment score
+ * ErrorProbability : The error of each flow (corresponds to the signals in ReadKey)
+ * ReadKeyLength
+ * GroundTruthKeyLength
+ * CycleSkipStatus : One of NS (Non Skip), PCS (Possible Cycle Skip), or CS (Cycle Skip)
+ * Cigar
+ * LowestQBaseTP
+ *
+ * GATK Report Description
+ * In the quality report (optional), the following tables are included:
+ *
+ * qualReport:error rate per qual : The error rate for each quality. Columns:
+ *
+ * - qual: The encoded quality
+ *
- count: The number of times the quality was observed
+ *
- error: The error rate of the flows with this qual
+ *
- phred: the error translated into a phred score
+ *
+ *
+ * qual_hmerReport:error rate per qual by hmer. The error rate for each quality and hmer combination. Columns:
+ *
+ * - qual: The encoded quality
+ *
- hmer: The hmer length
+ *
- count: The number of times the quality was observed
+ *
- error: The error rate of the flows with this qual
+ *
+ *
+ * qual_hmer_deviation_base_Report:error rate per qual by hmer and deviation. The count of errors for each qual, hmer, deviation and base
+ *
+ * - qual: The encoded quality
+ *
- hmer: The hmer length
+ *
- deviation: The deviation (difference in signal, relative to the reference)
+ *
- base: The base
+ *
- count: The number of times the deviation was observed
+ *
+ *
+ * Phred/qual statistics per flow position report. Various statistics for each flow position in relationship to the found quality value. Columns:
+ *
+ * - flow - flow position
+ * - count - count of observations
+ * - min - minimal observed quality
+ * - max - maximal observed quality
+ * - mean - mean observed value
+ * - median - median observed value
+ * - std - standard deviation
+ * - p1...Pn - percentil columns, accotding to the --quality-percentiles parameter
+ *
+ *
+ * Usage examples
+ *
+ * gatk GroundTruthScorer \
+ * -I input.bam \
+ * -R reference.fasta.gz
+ * -L chr20 \
+ * --output-csv output.csv \
+ * --report-file report.txt \
+ * --omit-zeros-from-report \ (optional)
+ * --features-file dbsnp.chr9.vcf.gz \ (optional)
+ * --genome-prior genome_prior.csv (optional)
+ *
+ *
+ * {@GATK.walkertype ReadWalker}
+ */
+
+@CommandLineProgramProperties(
+ summary = "Ground Truth Scorer",
+ oneLineSummary = "Score reads against a reference/ground truth",
+ programGroup = FlowBasedProgramGroup.class
+)
+@DocumentedFeature
+@ExperimentalFeature
+public class GroundTruthScorer extends ReadWalker {
+ private static final Logger logger = LogManager.getLogger(GroundTruthScorer.class);
+
+ public static final String OUTPUT_CSV_LONG_NAME = "output-csv";
+ public static final String REPORT_FILE_LONG_NAME = "report-file";
+ public static final String USE_SOFTCLIPPED_BASES_LONG_NAME = "use-softclipped-bases";
+ public static final String GENOME_PRIOR_LONG_NAME = "genome-prior";
+ public static final String FEATURES_FILE_LONG_NAME = "features-file";
+ public static final String NORMALIZED_SCORE_THRESHOLD_LONG_NAME = "normalized-score-threshold";
+ public static final String ADD_MEAN_CALL_LONG_NAME = "add-mean-call";
+ public static final String GT_NO_OUTPUT_LONG_NAME = "gt-no-output";
+ public static final String OMIT_ZEROS_FROM_REPORT = "omit-zeros-from-report";
+ public static final String QUALITY_PERCENTILES = "quality-percentiles";
+ public static final String EXCLUDE_ZERO_FLOWS = "exclude-zero-flows";
+
+ private static final int QUAL_VALUE_MAX = 60;
+ private static final int HMER_VALUE_MAX = 100; //TODO: This should become a parameter
+ private static final int BASE_VALUE_MAX = FlowBasedRead.DEFAULT_FLOW_ORDER.length() - 1;
+
+ private static final double NORMALIZED_SCORE_THRESHOLD_DEFAULT = -0.1;
+ private static final double DEFAULT_RATIO_THRESHOLD = 0.003;
+
+ /*
+ Private accumulator class for counting false/true observations (hence Boolean).
+
+ Observations are counted at a top level and are also optionally classified into a set of bins (the
+ number of which is fixed upon construction). The bins themselves are also BooleanAccumulator objects,
+ resulting in a tree like multi-level accumulator.
+
+
+ GroundTruthScores builds a four level deep accumulation tree, which can support observations of a
+ boolean event with 3-deep context (bin1,bin2,bin3).
+
+ Once accumulation is done, the instance is able to generate a suitable GATKReportTable for any given
+ bin depth (1, 2 or 3).
+
+ */
+ private static class BooleanAccumulator {
+ long falseCount;
+ long trueCount;
+ BooleanAccumulator[] bins;
+
+ // add an observation to this accumulator
+ void add(final boolean b) {
+ if (b) {
+ trueCount++;
+ } else {
+ falseCount++;
+ }
+ }
+
+ // add an observation to this accumulator and to one of the bins in the level under it (1 deep)
+ void add(final boolean b, final int bin) {
+ add(b);
+ if ( bins != null && bin >= 0 && bin < bins.length ) {
+ bins[bin].add(b);
+ } else {
+ logger.warn("bin out of range; " + bin + ", range: [0," + bins.length + "), clipped");
+ bins[Math.max(0, Math.min(bin, bins.length - 1))].add(b);
+ }
+ }
+
+ // add an observation to this accumulator and to two levels of bins under it (2 deep)
+ void add(final boolean b, final int bin, final int bin2) {
+ add(b);
+ if ( bins != null && bin >= 0 && bin < bins.length ) {
+ bins[bin].add(b, bin2);
+ } else {
+ logger.warn("bin out of range; " + bin + ", range: [0," + bins.length + "), clipped");
+ bins[Math.max(0, Math.min(bin, bins.length - 1))].add(b, bin2);
+ }
+ }
+
+ // add an observation to this accumulator and to three levels of bins under it (3 deep)
+ void add(final boolean b, final int bin, final int bin2, final int bin3) {
+ add(b);
+ if ( bins != null && bin >= 0 && bin < bins.length ) {
+ bins[bin].add(b, bin2, bin3);
+ } else {
+ logger.warn("bin out of range; " + bin + ", range: [0," + bins.length + "), clipped");
+ bins[Math.max(0, Math.min(bin, bins.length - 1))].add(b, bin2, bin3);
+ }
+ }
+
+ // get observation count of this accumulator
+ long getCount() {
+ return falseCount + trueCount;
+ }
+
+ // get the false rate/ration for this accumulator
+ double getFalseRate() {
+ return (getCount() == 0) ? 0.0 : ((double)falseCount / getCount());
+ }
+
+ // create a set of accumulators with 3-deep bin nesting
+ static BooleanAccumulator[] newReport(final int size, final int binCount, final int binCount2, final int binCount3) {
+ BooleanAccumulator[] report = new BooleanAccumulator[size];
+ for ( byte i = 0 ; i < report.length ; i++ ) {
+ report[i] = new BooleanAccumulator();
+ if ( binCount != 0 ) {
+ report[i].bins = new BooleanAccumulator[binCount];
+ for ( int j = 0 ; j < report[i].bins.length ; j++ ) {
+ report[i].bins[j] = new BooleanAccumulator();
+ if ( binCount2 != 0 ) {
+ report[i].bins[j].bins = new BooleanAccumulator[binCount2];
+ for ( int k = 0 ; k < report[i].bins[j].bins.length ; k++ ) {
+ report[i].bins[j].bins[k] = new BooleanAccumulator();
+ if (binCount3 != 0) {
+ report[i].bins[j].bins[k].bins = new BooleanAccumulator[binCount3];
+ for (int m = 0; m < report[i].bins[j].bins[k].bins.length; m++) {
+ report[i].bins[j].bins[k].bins[m] = new BooleanAccumulator();
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ return report;
+ }
+
+ // create a GATK report from a set of accumulators without nesting into their bins
+ static GATKReportTable newReportTable(final BooleanAccumulator[] report, final String name, final double probThreshold, final boolean omitZeros) {
+ final GATKReportTable table = new GATKReportTable(name + "Report", "error rate per " + name, 4);
+ table.addColumn(name, "%d");
+ table.addColumn("count", "%d");
+ table.addColumn("error", "%f");
+ table.addColumn("phred", "%d");
+ int rowIndex = 0;
+ for (int i = 0; i < report.length; i++) {
+ if ( omitZeros && i != 0 && report[i].getCount() == 0 )
+ continue;
+ else {
+ final double rate = report[i].getFalseRate();
+ final double phredRate = (rate == 0.0 && report[i].getCount() != 0 && probThreshold != 0.0) ? probThreshold : rate;
+
+ table.set(rowIndex, 0, i);
+ table.set(rowIndex, 1, report[i].getCount());
+ table.set(rowIndex, 2, rate);
+ table.set(rowIndex, 3, phredRate != 0 ? (int) Math.ceil(-10.0 * Math.log10(phredRate)) : 0);
+ rowIndex++;
+ }
+ }
+ return table;
+ }
+
+ // create a GATK report from a set of accumulators while nesting into one level of their bins (1 deep)
+ static GATKReportTable newReportTable(final BooleanAccumulator[] report, final String name1, final String name2, final boolean omitZeros) {
+ final GATKReportTable table = new GATKReportTable(name1 + "_" + name2 + "Report", "error rate per " + name1 + " by " + name2, 4);
+ table.addColumn(name1, "%d");
+ table.addColumn(name2, "%d");
+ table.addColumn("count", "%d");
+ table.addColumn("error", "%f");
+ int rowIndex = 0;
+ for (int i = 0; i < report.length; i++) {
+ for ( int j = 0; j < report[i].bins.length ; j++ ) {
+ if ( omitZeros && (i != 0 || j != 0) && report[i].bins[j].getCount() == 0 )
+ continue;
+ else {
+ table.set(rowIndex, 0, i);
+ table.set(rowIndex, 1, j);
+ table.set(rowIndex, 2, report[i].bins[j].getCount());
+ table.set(rowIndex, 3, report[i].bins[j].getFalseRate());
+ rowIndex++;
+ }
+ }
+ }
+ return table;
+ }
+
+ // create a GATK report from a set of accumulators while nesting into two levels of their bins (2 deep)
+ static GATKReportTable newReportTable(final BooleanAccumulator[] report, final String name1, final String name2, final String name3, final String name4, final boolean omitZeros) {
+ final GATKReportTable table = new GATKReportTable(name1 + "_" + name2 + "_" + name3 + "_" + name4 + "_Report", "error rate per " + name1 + " by " + name2 + " and " + name3, 5);
+ table.addColumn(name1, "%d");
+ table.addColumn(name2, "%d");
+ table.addColumn(name3, "%s");
+ table.addColumn(name4, "%s");
+ table.addColumn("count", "%d");
+ int rowIndex = 0;
+ for (int i = 0; i < report.length; i++) {
+ for ( int j = 0; j < report[i].bins.length ; j++ ) {
+ for ( int k = 0; k < report[i].bins[j].bins.length ; k++ ) {
+ for ( int m = 0; m < report[i].bins[j].bins[k].bins.length ; m++ ) {
+ if ( omitZeros && (i != 0 || j != 0 || k != 0 || m != 0) && report[i].bins[j].bins[k].bins[m].getCount() == 0 )
+ continue;
+ else {
+ table.set(rowIndex, 0, i);
+ table.set(rowIndex, 1, j);
+ table.set(rowIndex, 2, binToDeviation(k));
+ table.set(rowIndex, 3, String.format("%c", binToBase(m)));
+ table.set(rowIndex, 4, report[i].bins[j].bins[k].bins[m].getCount());
+ rowIndex++;
+ }
+ }
+ }
+ }
+ }
+ return table;
+ }
+ }
+
+ private static class PercentileReport extends SeriesStats {
+
+ static GATKReportTable newReportTable(final Vector report, String qualityPercentiles) {
+ String[] qp = qualityPercentiles.split(",");
+ final GATKReportTable table = new GATKReportTable("PhredBinAccumulator", "PhredBinAccumulator", 8 + qp.length);
+ table.addColumn("flow", "%d");
+ table.addColumn("count", "%d");
+ table.addColumn("min", "%f");
+ table.addColumn("max", "%f");
+ table.addColumn("mean", "%f");
+ table.addColumn("median", "%f");
+ table.addColumn("std", "%f");
+ for ( final String p : qp ) {
+ table.addColumn("p" + p, "%f");
+ }
+ int rowIndex = 0;
+ for ( final PercentileReport r : report ) {
+ int col = 0;
+ table.set(rowIndex, col++, rowIndex);
+ table.set(rowIndex, col++, r.getCount());
+ table.set(rowIndex, col++, r.getMin());
+ table.set(rowIndex, col++, r.getMax());
+ table.set(rowIndex, col++, r.getMean());
+ table.set(rowIndex, col++, r.getMedian());
+ table.set(rowIndex, col++, r.getStd());
+ for ( String p : qp ) {
+ table.set(rowIndex, col++, r.getPercentile(Double.parseDouble(p)));
+ }
+ rowIndex++;
+ }
+ return table;
+ }
+
+ void addProb(double p) {
+ super.add(-10 * Math.log10(p));
+ }
+ }
+
+ @Argument(fullName = OUTPUT_CSV_LONG_NAME, doc="main CSV output file. supported file extensions: .csv, .csv.gz.")
+ public GATKPath outputCsvPath = null;
+
+ @Argument(fullName = REPORT_FILE_LONG_NAME, doc="report output file.", optional = true)
+ public GATKPath reportFilePath = null;
+
+ @ArgumentCollection
+ public LikelihoodEngineArgumentCollection likelihoodArgs = new LikelihoodEngineArgumentCollection();
+
+ @ArgumentCollection
+ public FlowBasedAlignmentArgumentCollection fbargs = new FlowBasedAlignmentArgumentCollection();
+
+ @Argument(fullName = USE_SOFTCLIPPED_BASES_LONG_NAME, doc="", optional = true)
+ public boolean useSoftclippedBases;
+
+ @Argument(fullName = GENOME_PRIOR_LONG_NAME, doc="CSV input file containing genome-prior (one line per base with hmer frequencies).", optional = true)
+ public GATKPath genomePriorPath;
+
+ @Argument(fullName = FEATURES_FILE_LONG_NAME, doc="A VCF file containing features to be used as a use for filtering reads.", optional = true)
+ public FeatureDataSource features;
+
+ @Argument(fullName = NORMALIZED_SCORE_THRESHOLD_LONG_NAME, doc="threshold for normalized score, below which reads are ignored", optional = true)
+ public double normalizedScoreThreshold = NORMALIZED_SCORE_THRESHOLD_DEFAULT;
+
+ @Argument(fullName = ADD_MEAN_CALL_LONG_NAME, doc="Add ReadMeanCall and ReadProbs columns to output", optional = true)
+ public boolean addMeanCalll;
+
+ @Argument(fullName = GT_NO_OUTPUT_LONG_NAME, doc = "do not generate output records", optional = true)
+ public boolean noOutput = false;
+
+ @Argument(fullName = OMIT_ZEROS_FROM_REPORT, doc = "omit zero values from output report", optional = true)
+ public boolean omitZerosFromReport = false;
+
+ @Argument(fullName = QUALITY_PERCENTILES, doc = "list of quality percentiles, defaults to 10,25,50,75,90", optional = true)
+ public String qualityPercentiles = "10,25,50,75,90";
+
+ @Argument(fullName = EXCLUDE_ZERO_FLOWS, doc = "should flows with a call of zero be included in the percentile report?", optional = true)
+ public boolean excludeZeroFlows = false;
+
+ // locals
+ private FlowBasedAlignmentLikelihoodEngine likelihoodCalculationEngine;
+ private PrintWriter outputCsv;
+ private DecimalFormat doubleFormat = new DecimalFormat("0.0#####");
+ private GenomePriorDB genomePriorDB;
+ private BooleanAccumulator[] qualReport;
+ private String[] csvFieldOrder;
+ private Vector percentileReports;
+
+ // static/const
+ static final private String[] CSV_FIELD_ORDER_BASIC = {
+ "ReadName", "ReadKey", "ReadIsReversed", "ReadMQ", "ReadRQ", "GroundTruthKey", "ReadSequence",
+ "Score", "NormalizedScore", "ErrorProbability",
+ "ReadKeyLength", "GroundTruthKeyLength", "CycleSkipStatus", "Cigar", "LowestQBaseTP"
+ };
+ static final private String[] CSV_FIELD_ORDER_MEAN_CALL = {
+ "ReadProbs", "ReadMeanCall"
+ };
+
+ @Override
+ public void onTraversalStart() {
+ super.onTraversalStart();
+
+ // establish csv fields
+ List order = new LinkedList<>(Arrays.asList(CSV_FIELD_ORDER_BASIC));
+ if ( addMeanCalll ) {
+ order.addAll(Arrays.asList(CSV_FIELD_ORDER_MEAN_CALL));
+ }
+ csvFieldOrder = order.toArray(new String[0]);
+
+ // create likelihood engine
+ final ReadLikelihoodCalculationEngine engine = AssemblyBasedCallerUtils.createLikelihoodCalculationEngine(likelihoodArgs, false);
+ if ( engine instanceof FlowBasedAlignmentLikelihoodEngine ) {
+ likelihoodCalculationEngine = (FlowBasedAlignmentLikelihoodEngine)engine;
+ } else {
+ throw new GATKException("must use a flow based likelihood calculation engine");
+ }
+
+ // open genome prior if provided
+ if ( genomePriorPath != null ) {
+ try {
+ genomePriorDB = new GenomePriorDB(genomePriorPath);
+ } catch (IOException e) {
+ throw new GATKException("failed to open genome-prior file: " + genomePriorPath);
+ }
+ }
+
+ // open output, write header
+ try {
+ if (outputCsvPath.toPath().toString().endsWith(".gz")) {
+ outputCsv = new PrintWriter(new GZIPOutputStream(outputCsvPath.getOutputStream()));
+ } else {
+ outputCsv = new PrintWriter(outputCsvPath.getOutputStream());
+ }
+ } catch (IOException e) {
+ throw new GATKException("failed to open csv output: " + outputCsvPath, e);
+ }
+ emitCsvHeaders();
+
+ // initialize reports
+ if ( reportFilePath != null ) {
+ qualReport = BooleanAccumulator.newReport(QUAL_VALUE_MAX + 1, HMER_VALUE_MAX + 1, deviationToBin(HMER_VALUE_MAX + 1), BASE_VALUE_MAX + 1);
+
+ // establish max hmer size of flow input
+ int maxClass = FlowBasedRead.MAX_CLASS;
+ final SAMFileHeader header = getHeaderForReads();
+ if ( header != null ) {
+ for ( SAMReadGroupRecord rg : header.getReadGroups() ) {
+ if ( rg.getAttribute(FlowBasedRead.MAX_CLASS_READ_GROUP_TAG) != null ) {
+ maxClass = Math.max(maxClass, Integer.parseInt(rg.getAttribute(FlowBasedRead.MAX_CLASS_READ_GROUP_TAG)));
+ }
+ }
+ }
+ percentileReports = new Vector<>();
+ }
+ }
+
+ @Override
+ public void closeTool() {
+
+ // close main output csv
+ if ( outputCsv != null ) {
+ outputCsv.close();
+ }
+
+ // write reports
+ if ( reportFilePath != null ) {
+ final GATKReport report = new GATKReport(
+ BooleanAccumulator.newReportTable(qualReport, "qual", DEFAULT_RATIO_THRESHOLD, omitZerosFromReport),
+ BooleanAccumulator.newReportTable(qualReport, "qual", "hmer", omitZerosFromReport),
+ BooleanAccumulator.newReportTable(qualReport, "qual", "hmer", "deviation", "base", omitZerosFromReport),
+ PercentileReport.newReportTable(percentileReports, qualityPercentiles)
+ );
+ try ( final PrintStream ps = new PrintStream(reportFilePath.getOutputStream()) ) {
+ report.print(ps);
+ }
+ }
+
+ super.closeTool();
+ }
+
+ @Override
+ public void apply(final GATKRead read, final ReferenceContext referenceContext, final FeatureContext featureContext) {
+
+ // working with unmapped reads is really not practical, as we need the reference to work out ground truth
+ if ( read.isUnmapped() )
+ return;
+ if ( referenceContext.getWindow() == null )
+ return;
+
+ // handle read clipping
+ final GATKRead clippedRead;
+ if (isSoftClipped(read) ) {
+ if (useSoftclippedBases) {
+ referenceContext.setWindow(read.getStart() - read.getUnclippedStart(), read.getUnclippedEnd() - read.getEnd());;
+ clippedRead = ReadClipper.revertSoftClippedBases(read);
+ } else {
+ clippedRead = ReadClipper.hardClipSoftClippedBases(read);
+ }
+ } else {
+ clippedRead = read;
+ }
+
+ // filter?
+ if ( (features != null) && !filter(clippedRead, referenceContext) ) {
+ return;
+ }
+
+ // create flow read/haplotype
+ final FlowBasedReadUtils.ReadGroupInfo rgInfo = FlowBasedReadUtils.getReadGroupInfo(getHeaderForReads(), clippedRead);
+ final FlowBasedRead flowRead = new FlowBasedRead(clippedRead, rgInfo.flowOrder, rgInfo.maxClass, fbargs);
+ final Haplotype haplotype = new Haplotype(referenceContext.getBases(), true);
+ final FlowBasedHaplotype flowHaplotype = new FlowBasedHaplotype(haplotype, rgInfo.flowOrder);
+
+ // is this really needed?
+ if ( !flowRead.isValid() ) {
+ return;
+ }
+
+ // compute score
+ final int hapKeyLength = flowHaplotype.getKeyLength();
+ final double score = FlowFeatureMapper.computeLikelihoodLocal(flowRead, flowHaplotype, hapKeyLength, false);
+ final double normalizedScore = score / flowRead.getKeyLength();
+ if ( normalizedScore < normalizedScoreThreshold )
+ return;
+
+ // compute error probability
+ final double[] errorProb = computeErrorProb(flowRead, genomePriorDB);
+
+ // cycle skip
+ final FlowBasedReadUtils.CycleSkipStatus cycleSkipStatus = FlowBasedReadUtils.getCycleSkipStatus(flowRead, referenceContext);
+
+ // accumulate reports
+ if ( cycleSkipStatus != FlowBasedReadUtils.CycleSkipStatus.CS && qualReport != null ) {
+ addToQualReport(flowRead, referenceContext, errorProb);
+ }
+
+ // emit
+ try {
+ emit(flowRead, flowHaplotype, score, normalizedScore, errorProb, read, cycleSkipStatus);
+ } catch (IOException e) {
+ throw new GATKException("failed to write output record", e);
+ }
+ }
+
+ private boolean filter(final GATKRead read, final ReferenceContext referenceContext) {
+
+ // loop on features contained in the read, check that they are in agreement with read data
+ Iterator iter = features.query(new SimpleInterval(read));
+ byte[] ref = referenceContext.getBases();
+ while ( iter.hasNext() ) {
+ final VariantContext vc = iter.next();
+ for ( int refCoord = vc.getStart() ; refCoord <= vc.getEnd() ; refCoord++ ) {
+
+ // get byte from read
+ Optional readByte = ReadUtils.getReadBaseAtReferenceCoordinate(read, refCoord);
+ if ( !readByte.isPresent() ) {
+ return false;
+ }
+
+ // get byte from reference
+ byte refByte = ref[refCoord - referenceContext.getWindow().getStart()];
+
+ // compare
+ if ( refByte != readByte.get() ) {
+ return false;
+ }
+ }
+ }
+
+ // if here, no interference detected
+ return true;
+ }
+
+ /*
+ * compute error probability vector for a read
+ *
+ * The vector has one element for each flow key, representing the probability complementing the call-probability to 1
+ * This is further complicated by the optional presence of a genome-prior database, which provides factoring for
+ * each hmer length (on a base basis)
+ */
+ private double[] computeErrorProb(final FlowBasedRead flowRead, final GenomePriorDB genomePriorDB) {
+
+ final int[] key = flowRead.getKey();
+ final byte[] flowOrder = flowRead.getFlowOrderArray();
+
+ final double[] probCol = new double[flowRead.getMaxHmer() + 1];
+ double[] result = new double[key.length];
+
+ for ( int i = 0 ; i < key.length ; i++ ) {
+
+ // step 1 - extract column & sum
+ double sum = 0;
+ for ( int j = 0 ; j < probCol.length ; j++ ) {
+ sum += (probCol[j] = flowRead.getProb(i, j));
+ }
+
+ // step 2 - normalize column
+ if ( sum != 0.0 ) {
+ for (int j = 0; j < probCol.length; j++) {
+ probCol[j] /= sum;
+ }
+ }
+
+ // step 3 - scale according to prior genome?
+ if ( genomePriorDB != null ) {
+
+ long[] prior = genomePriorDB.getPriorForBase(flowOrder[i]);
+ sum = 0;
+ if ( prior != null ) {
+ for (int j = 0; j < probCol.length; j++) {
+ sum += (probCol[j] *= prior[j]);
+ }
+ }
+
+ // assign normalized result
+ if ( sum != 0.0 ) {
+ result[i] = 1 - (probCol[Math.min(probCol.length - 1, Math.min(key[i], flowRead.getMaxHmer()))] / sum);
+ } else {
+ // revert to non-prior normalization
+ result[i] = 1 - probCol[Math.min(key[i], flowRead.getMaxHmer())];
+ }
+ } else {
+
+ // assign normalized result
+ result[i] = 1 - probCol[Math.min(key[i], flowRead.getMaxHmer())];
+ }
+
+ // accumulate error probabilities
+ if ( percentileReports != null ) {
+ if ( key[i] != 0 || !excludeZeroFlows ) {
+ while ( percentileReports.size() < (i + 1) ) {
+ percentileReports.add(new PercentileReport());
+ }
+ percentileReports.get(i).addProb(result[i]);
+ }
+ }
+ }
+
+ return result;
+ }
+
+ /*
+ * compute lowest-quality-base-tp-value vector for a read
+ *
+ * The vector has one element for each flow key, representing the value of the tp for the hmer base
+ * which has the lowest quality value
+ *
+ * example:
+ * Bases: TTTTT
+ * Qs: ABCBA
+ * Tp: -1 1 0 1 -1
+ * Output: -1
+ */
+ private byte[] computeLowestQBaseTP(final FlowBasedRead flowRead) {
+
+ final int[] key = flowRead.getKey();
+ byte[] result = new byte[key.length];
+ final byte[] tp = flowRead.getAttributeAsByteArray("tp");
+ final byte[] qual = flowRead.getBaseQualitiesNoCopy();
+
+ // loop om key
+ int seq_i = 0;
+ for ( int i = 0 ; i < key.length ; i++ ) {
+
+ // extract hmer length, zero is easy
+ int hmer = key[i];
+ if ( hmer == 0 ) {
+ result[i] = 0;
+ continue;
+ }
+
+ // scan qualities for the lowest value, start with first
+ // as qualities and tp are symetric, we can scan up to the middle
+ // when finding the middle (offset) account for even/odd hmers
+ result[i] = tp[seq_i];
+ byte lowestQ = qual[seq_i];
+ int hmer_scan_length = (hmer + 1) / 2;
+ for ( int j = 1 ; j < hmer_scan_length ; j++ ) {
+ if ( qual[seq_i + j] < lowestQ ) {
+ result[i] = tp[seq_i + j];
+ lowestQ = qual[seq_i + j];
+ }
+ }
+
+ // advance
+ seq_i += hmer;
+ }
+
+ return result;
+ }
+
+ private void emitCsvHeaders() {
+
+ outputCsv.println(StringUtils.join(csvFieldOrder, ","));
+ }
+
+ private void emit(final FlowBasedRead flowRead, final FlowBasedHaplotype refHaplotype, double score, final double normalizedScore, final double[] errorProb,
+ GATKRead read,
+ FlowBasedReadUtils.CycleSkipStatus cycleSkipStatus) throws IOException {
+
+ // build line columns
+ final Map cols = new LinkedHashMap<>();
+
+ // read info
+ cols.put("ReadName", flowRead.getName());
+ cols.put("ReadIsReversed", flowRead.isReverseStrand() ? 1 : 0);
+ cols.put("ReadMQ", flowRead.getMappingQuality());
+ cols.put("ReadRQ", flowRead.getAttributeAsFloat("rq"));
+ cols.put("CycleSkipStatus", cycleSkipStatus);
+ cols.put("Cigar", read.getCigar().toString());
+
+
+ // keys, seq, etc
+ cols.put("ReadKey", "\"" + StringUtils.join(flowRead.getKey(), ',') + "\"");
+ cols.put("GroundTruthKey", "\"" + StringUtils.join(refHaplotype.getKey(), ',') + "\"");
+ cols.put("ReadSequence", flowRead.getBasesString());
+ cols.put("ReadKeyLength", flowRead.getKeyLength());
+ cols.put("GroundTruthKeyLength", refHaplotype.getKeyLength());
+
+ // scores
+ cols.put("Score", score);
+ cols.put("NormalizedScore", normalizedScore);
+ cols.put("ErrorProbability", "\"" + StringUtils.join(
+ Arrays.stream(errorProb).mapToObj(v -> doubleFormat.format(v)).toArray(),
+ ',') + "\"");
+
+ // lowest q base tp
+ final byte[] lowestQBaseTP = computeLowestQBaseTP(flowRead);
+ cols.put("LowestQBaseTP", "\"" + StringUtils.join(lowestQBaseTP, ',') + "\"");
+
+ // add read probabilities
+ if ( addMeanCalll ) {
+ double[][] readProbsAndMeanCall = collectReadProbs(flowRead);
+ cols.put("ReadProbs", "\"" + StringUtils.join(
+ Arrays.stream(readProbsAndMeanCall[0]).mapToObj(v -> doubleFormat.format(v)).toArray(),
+ ',') + "\"");
+ cols.put("ReadMeanCall", "\"" + StringUtils.join(
+ Arrays.stream(readProbsAndMeanCall[1]).mapToObj(v -> doubleFormat.format(v)).toArray(),
+ ',') + "\"");
+ }
+
+ // construct line
+ StringBuilder sb = new StringBuilder();
+ int colIndex = 0;
+ for ( String field : csvFieldOrder ) {
+ if ( colIndex++ > 0 ) {
+ sb.append(',');
+ }
+ if ( !cols.containsKey(field) ) {
+ throw new GATKException("column missing from csv line: " + field);
+ }
+ sb.append(cols.get(field));
+ cols.remove(field);
+ }
+ if ( cols.size() > 0 ) {
+ throw new GATKException("invalid columns on csv line: " + cols.keySet());
+ }
+
+ // output line
+ if ( !noOutput ) {
+ outputCsv.println(sb);
+ }
+ }
+
+ private double[][] collectReadProbs(final FlowBasedRead read) {
+
+ final int keyLength = read.getKeyLength();
+ final int maxHmer = read.getMaxHmer();
+ final double[] probs = new double[keyLength * (maxHmer + 1)];
+ final double[] meanCall = new double[keyLength];
+
+ // retrieve probs
+ int pos = 0;
+ for ( int flow = 0 ; flow < keyLength ; flow++ ) {
+ double mc = 0;
+ int mc_sum = 0;
+ for ( int hmer = 0 ; hmer <= maxHmer ; hmer++ ) {
+ final double p = read.getProb(flow, hmer);
+ probs[pos++] = p;
+ mc += (p * hmer);
+ mc_sum += hmer;
+ }
+ meanCall[flow] = mc / mc_sum;
+ }
+
+ return new double[][] {probs, meanCall};
+ }
+
+ static private class GenomePriorDB {
+
+ final private Map db = new LinkedHashMap<>();
+
+ GenomePriorDB(GATKPath path) throws IOException {
+
+ final CSVReader csvReader = new CSVReader(new InputStreamReader(path.getInputStream()));
+ String[] line;
+ while ( (line = csvReader.readNext()) != null ) {
+ long[] prior = new long[HMER_VALUE_MAX + 1];
+ Byte base = line[0].getBytes()[0];
+ for ( int i = 0 ; i < prior.length ; i++ ) {
+ if ( i == 0 ){
+ prior[i] = Long.parseLong(line[i+1]);
+ } else {
+ prior[i] = Long.parseLong(line[i]);
+ }
+ }
+ db.put(base, prior);
+ }
+ }
+
+ long[] getPriorForBase(byte base) {
+ return db.get(base);
+ }
+ }
+
+ private void addToQualReport(FlowBasedRead flowRead, ReferenceContext referenceContext, final double[] errorProb) {
+
+ // convert reference to key space
+ final Haplotype haplotype = new Haplotype(referenceContext.getBases(), true);
+ final FlowBasedHaplotype flowHaplotype = new FlowBasedHaplotype(haplotype, flowRead.getFlowOrder());
+
+ // access keys & flow order
+ final int[] readKey = flowRead.getKey();
+ final int[] hapKey = flowHaplotype.getKey();
+ final byte[] flowOrder = flowRead.getFlowOrder().getBytes();
+
+ // loop on key positions
+ if ( readKey.length != hapKey.length ) {
+ return;
+ }
+ for ( int flow = 0 ; flow < readKey.length ; flow++ ) {
+
+ // determine quality
+ final double prob = Precision.round(errorProb[flow], (QUAL_VALUE_MAX / 10) + 1);
+ final int qual = (int)Math.ceil(-10 * Math.log10(prob));
+
+ // determine if matches reference
+ final int deviation = readKey[flow] - hapKey[flow];
+ final boolean same = (deviation == 0);
+
+ // accumulate
+ if ( qual < qualReport.length ) {
+ int baseBin = baseToBin(flowOrder[flow % flowOrder.length], flowRead.isReverseStrand());
+ qualReport[qual].add(same, readKey[flow], deviationToBin(deviation), baseBin);
+ }
+ }
+ }
+
+ static private int baseToBin(byte base, boolean isReverseStrand) {
+ final byte trueBase = !isReverseStrand ? base : BaseUtils.simpleComplement(base);
+ return FlowBasedRead.DEFAULT_FLOW_ORDER.indexOf(trueBase);
+ }
+
+ static private byte binToBase(int bin) {
+ return (byte)FlowBasedRead.DEFAULT_FLOW_ORDER.charAt(bin);
+ }
+
+ // 0,-1,1,-2,2... -> 0,1,2,3,4...
+ static private int deviationToBin(final int deviation) {
+ if ( deviation >= 0 ) {
+ return deviation * 2;
+ } else {
+ return (-deviation * 2) - 1;
+ }
+ }
+
+ // 0,1,2,3,4... -> 0,-1,1,-2,2...
+ static private String binToDeviation(final int bin) {
+ if ( bin == 0 ) {
+ return "0";
+ } else if ( (bin % 2) == 0 ) {
+ return String.format("+%d", bin / 2);
+ } else {
+ return String.format("%d", -((bin+1) / 2));
+ }
+ }
+
+ static private boolean isSoftClipped( final GATKRead read ) {
+ if ( read.isUnmapped() )
+ return false;
+ if ( read.getCigar().getFirstCigarElement() == null )
+ return false;
+ final CigarOperator firstOperator = read.getCigar().getFirstCigarElement().getOperator();
+ final CigarOperator lastOperator = read.getCigar().getLastCigarElement().getOperator();
+ return (firstOperator == CigarOperator.SOFT_CLIP && lastOperator != CigarOperator.SOFT_CLIP) ||
+ (firstOperator != CigarOperator.SOFT_CLIP && lastOperator == CigarOperator.SOFT_CLIP);
+ }
+}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/SeriesStats.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/SeriesStats.java
new file mode 100644
index 00000000000..2ba594b7058
--- /dev/null
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/SeriesStats.java
@@ -0,0 +1,160 @@
+package org.broadinstitute.hellbender.tools.walkers.groundtruth;
+
+import org.apache.logging.log4j.LogManager;
+import org.apache.logging.log4j.Logger;
+
+import java.io.IOException;
+import java.io.PrintWriter;
+import java.util.LinkedHashMap;
+import java.util.Map;
+import java.util.SortedMap;
+import java.util.TreeMap;
+import java.util.concurrent.atomic.AtomicInteger;
+
+public class SeriesStats {
+
+ private static final Logger logger = LogManager.getLogger(SeriesStats.class);
+
+ // local state
+ private double last = Double.NaN;
+ private int count = 0;
+ private double sum = 0;
+ private double min = Double.NaN;
+ private double max = Double.NaN;
+ private SortedMap bins = new TreeMap<>();
+ private int intCount = 0;
+ private Map auxBins = new LinkedHashMap<>();
+
+ public void csvWrite(final String path) throws IOException {
+ logger.info("Writing SeriesStats " + toDigest() + " into " + path);
+ PrintWriter pw = new PrintWriter(path);
+ pw.println("value,count");
+ boolean intKeys = isIntKeys();
+ for (Map.Entry entry : bins.entrySet() ) {
+ if ( intKeys ) {
+ pw.println(String.format("%d,%d", entry.getKey().intValue(), entry.getValue().get()));
+ } else {
+ pw.println(String.format("%f,%d", entry.getKey(), entry.getValue().get()));
+ }
+ }
+ pw.close();
+ }
+
+ public void add(int v) {
+ add((double)v);
+ intCount++;
+ }
+
+ public void add(double v) {
+
+ // save in simple values
+ last = v;
+ sum += v;
+ if ( count > 0 ) {
+ min = Math.min(min, v);
+ max = Math.max(max, v);
+ } else {
+ min = max = v;
+ }
+ count++;
+
+ // save in bins
+ final Double key = v;
+ if ( bins.containsKey(key) ) {
+ bins.get(key).incrementAndGet();
+ } else {
+ bins.put(key, new AtomicInteger(1));
+ }
+ }
+
+ public double getLast() {
+ return last;
+ }
+
+ public int getCount() {
+ return count;
+ }
+
+ public double getMin() {
+ return (count != 0) ? min : Double.NaN;
+ }
+
+ public double getMax() {
+ return (count != 0) ? max : Double.NaN;
+ }
+
+ public int getUniq() {
+ return bins.size();
+ }
+
+ public double getMean() {
+ return (count != 0) ? (sum / count) : Double.NaN;
+ }
+
+ public double getMedian() {
+ return getPercentile(50);
+ }
+
+ public double getPercentile(double precentile) {
+ if ( count == 0 ) {
+ return Double.NaN;
+ } else if ( count == 1 ) {
+ return last;
+ } else {
+
+ int percentileIndex = (int)(count * precentile / 100);
+ int index = 0;
+ for (Map.Entry entry : bins.entrySet() ) {
+ int binSize = entry.getValue().get();
+ if ( percentileIndex >= index && (percentileIndex < (index + binSize)) ) {
+ return entry.getKey();
+ }
+ index += binSize;
+ }
+
+ // if here, we need the highest entry
+ return bins.lastKey();
+ }
+ }
+
+ public double getStd() {
+
+ if (count == 0) {
+ return Double.NaN;
+ }
+
+ // calculate mean
+ double mean = getMean();
+
+ // calculate sum of sqr(deviations)
+ double vsum = 0;
+ for (Map.Entry entry : bins.entrySet()) {
+ int binSize = entry.getValue().get();
+ vsum += (Math.pow(entry.getKey() - mean, 2) * binSize);
+ }
+
+ // calculate variance and std
+ double variance = vsum / count;
+ return Math.sqrt(variance);
+ }
+
+ public Map getBins() {
+ return this.bins;
+ }
+
+ public Map getAuxBins() {
+ return this.auxBins;
+ }
+
+ public String toDigest() {
+ if ( isIntKeys() ) {
+ return String.format("count=%d, min=%d, max=%d, median=%d, bin.count=%d", getCount(), (int)getMin(), (int)getMax(), (int)getMedian(), getBins().size());
+ } else {
+ return String.format("count=%d, min=%f, max=%f, median=%f, bin.count=%d", getCount(), getMin(), getMax(), getMedian(), getBins().size());
+ }
+ }
+
+ private boolean isIntKeys() {
+ return (count == intCount);
+ }
+}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AlleleFiltering.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AlleleFiltering.java
index 17728fa6035..3a350228128 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AlleleFiltering.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/AlleleFiltering.java
@@ -29,6 +29,8 @@
import java.util.stream.Collectors;
import java.util.stream.IntStream;
+import com.google.common.annotations.VisibleForTesting;
+
/**
* Filtering haplotypes that contribute weak alleles to the genotyping.
*
@@ -278,7 +280,8 @@ private AlleleLikelihoods subsetHaplotypesByAlleles(final A
* @param sorThreshold only variants with SOR above threshold will be considered
* @return list of alleles that can be removed
*/
- private List identifyBadAlleles(final List collectedRPLs, final List collectedSORs,
+ @VisibleForTesting
+ List identifyBadAlleles(final List collectedRPLs, final List collectedSORs,
final List alleles,
final double qualThreshold,
final double sorThreshold) {
@@ -303,9 +306,11 @@ private List identifyBadAlleles(final List collectedRPLs, final
//we then add alleles with high SOR. Note that amongh all allleles with the SOR higher than the SOR_THRESHOLD
//we will first filter the one with the lowest QUAL.
logger.debug(() -> String.format("SHA:: Have %d candidates with low QUAL", rplCount));
- for (int i = sorIndices.length-1 ; (i >= 0) && (collectedSORs.get(sorIndices[i])>SOR_THRESHOLD) ; i--) {
- if (!result.contains(alleles.get(sorIndices[i]))) {
- result.add(alleles.get(sorIndices[i]));
+ for (int i = sorIndices.length-1 ; (i >= 0) ; i--) {
+ if (collectedSORs.get(sorIndices[i])>SOR_THRESHOLD){
+ if (!result.contains(alleles.get(sorIndices[i]))) {
+ result.add(alleles.get(sorIndices[i]));
+ }
}
}
logger.debug(() -> String.format("SHA:: Have %d candidates with high SOR", result.size() - rplCount));
@@ -340,7 +345,7 @@ private AlleleLikelihoods