+ * This exists to mimic GATK3 interval traversal patterns
+ */
+ @Advanced
+ @Argument(fullName = IGNORE_VARIANTS_THAT_START_OUTSIDE_INTERVAL,
+ doc = "Restrict variant output to sites that start within provided intervals (only applies when an interval is specified)",
+ optional = true)
+ protected boolean ignoreIntervalsOutsideStart = false;
+
+ @Advanced
+ @Argument(fullName = COMBINE_VARIANTS_DISTANCE, doc = "Maximum distance for variants to be grouped together", optional = true)
+ protected int distanceToCombineVariants = defaultDistanceToGroupVariants();
+
+ @Advanced
+ @Argument(fullName = MAX_COMBINED_DISTANCE, doc = "Maximum distance for variants to be grouped together", optional = true)
+ protected int maxCombinedDistance = defaultMaxGroupedSpan();
+
+ @Advanced
+ @Argument(fullName = REFERENCE_WINDOW_PADDING, doc = "Number of bases on either side to expand spanning reference window", optional = true)
+ protected int referenceWindowPadding = defaultReferenceWindowPadding();
+
+ @Advanced
+ @Argument(fullName = "ignore-reference-blocks", shortName = "ignore-ref-blocks", optional = true)
+ protected boolean ignoreReferenceBlocks = false;
+
+ // override to group variants that start nearby but not at the same locus
+ protected int defaultDistanceToGroupVariants() {
+ return 10000;
+ }
+
+ // override to change reference padding
+ protected int defaultReferenceWindowPadding() {
+ return 1;
+ }
+
+ // override to cap the size span of combined variants
+ protected int defaultMaxGroupedSpan() {
+ return Integer.MAX_VALUE;
+ }
+
+ @Override
+ public boolean requiresReference() {
+ return true;
+ }
+
+ /**
+ * This method keeps track of all the variants it is passed and will feed all the variants that start at the same
+ * site to the reduce method.
+ *
+ * {@inheritDoc}
+ */
+ @Override
+ public final void apply(VariantContext variant, ReadsContext readsContext, ReferenceContext referenceContext, FeatureContext featureContext) {
+
+ // Filtering out variants that start outside of the specified intervals
+ if (ignoreIntervalsOutsideStart && !isWithinInterval(new SimpleInterval(variant.getContig(), variant.getStart(), variant.getStart()))) {
+ return;
+ }
+
+ if (ignoreReferenceBlocks && variant.getAlternateAlleles().size() ==1 && variant.getAlternateAllele(0).equals(Allele.NON_REF_ALLELE)) {
+ return;
+ }
+
+ // Collecting all the reads that start at a particular base into one.
+ if (currentVariants.isEmpty()) {
+ lastCurrentVariantContig = variant.getContig();
+ } else if (!currentVariants.get(0).contigsMatch(variant)
+ || variant.getStart() > lastCurrentVariantEnd) {
+ // Emptying any sites which should emit a new VC since the last one
+ apply(new ArrayList<>(currentVariants), currentReadsContexts);
+ currentVariants.clear();
+ currentReadsContexts.clear();
+ }
+
+ currentVariants.add(variant);
+ currentReadsContexts.add(readsContext);
+ if (variant.getEnd() > lastCurrentVariantEnd || !lastCurrentVariantContig.equals(variant.getContig())) {
+ lastCurrentVariantEnd = variant.getEnd();
+ lastCurrentVariantContig = variant.getContig();
+ }
+ }
+
+ /**
+ * This method must be implemented by tool authors.
+ *
+ * This is the primary traversal for any MultiVariantWalkerGroupedOnStart walkers. Will traverse over input variant contexts
+ * and call #apply() exactly once for each unique reference start position. All variants starting at each locus
+ * across source files will be grouped and passed as a list of VariantContext objects.
+ * @param variantContexts VariantContexts from driving variants with matching start position
+ * NOTE: This will never be empty
+ * @param referenceContext ReferenceContext object covering the reference of the longest spanning VariantContext
+ * @param readsContexts
+ */
+ public abstract void apply(final List This issue affects GATK versions 4.3.0.0 through 4.5.0.0, and is fixed in GATK 4.6.0.0. This issue also affects Picard versions 2.27.3 through 3.1.1, and is fixed in Picard 3.2.0. The bug is triggered when writing a CRAM file using one of the affected GATK/Picard versions,
+ * and both of the following conditions are met: When both of these conditions are met, the resulting CRAM file may have corrupt containers containing reads
+ * with an incorrect sequence. This tool writes a report to an output text file indicating whether the CRAM file appears to have read base corruption caused by issue 8768,
+ * and listing the affected containers. By default, the output report will have a summary of the average mismatch rate for all suspected bad containers
+ * and a few presumed good containers in order to determine if there is a large difference in the base mismatch rate. Optionally, a TSV file with the same information as the textual report, but in tabular form, can be written
+ * using the "--output-tsv" argument. To analyze the base mismatch rate for ALL containers, use the "verbose" option. Works on files ending in .cram. OpenMP and MKL parallelism can be controlled by setting the Advanced users may wish to set the Advanced users may wish to set the OpenMP and MKL parallelism can be controlled by setting the Advanced users may wish to set the Advanced users may wish to set the
+ *
+ *
+ *
+ *
+ * Sample Usage:
+ *
+ * gatk CRAMIssue8768Detector \
+ * -I input.cram \
+ * -O output_report.txt \
+ * -R reference.fasta
+ *
+ *
+ * gatk CRAMIssue8768Detector \
+ * -I input.cram \
+ * -O output_report.txt \
+ * -R reference.fasta \
+ * --output-tsv output_report_as_table.tsv
+ *
+ */
+@ExperimentalFeature
+@WorkflowProperties
+@CommandLineProgramProperties(
+ summary = "Analyze a CRAM file to check for base corruption caused by GATK issue 8768",
+ oneLineSummary = "Analyze a CRAM file to check for base corruption caused by GATK issue 8768",
+ programGroup = OtherProgramGroup.class
+)
+public class CRAMIssue8768Detector extends CommandLineProgram {
+ // default average mismatch rate threshold above which we consider the file to be corrupt
+ private static final double DEFAULT_MISMATCH_RATE_THRESHOLD = 0.05;
+
+ @Argument(fullName = StandardArgumentDefinitions.INPUT_LONG_NAME,
+ shortName = StandardArgumentDefinitions.INPUT_SHORT_NAME,
+ doc = "Input path of CRAM file to analyze",
+ common = true)
+ @WorkflowInput
+ public GATKPath inputPath;
+
+ @Argument(fullName = StandardArgumentDefinitions.OUTPUT_LONG_NAME,
+ shortName = StandardArgumentDefinitions.OUTPUT_SHORT_NAME,
+ doc = "Output diagnostics text file",
+ common = true)
+ @WorkflowOutput
+ public GATKPath textOutputPath;
+
+ public static final String OUTPUT_TSV__ARG_NAME = "output-tsv";
+ @Argument(fullName = OUTPUT_TSV__ARG_NAME,
+ shortName = OUTPUT_TSV__ARG_NAME,
+ doc = "Output diagnostics tsv file",
+ optional = true)
+ @WorkflowOutput
+ public GATKPath tsvOutputPath;
+
+ @Argument(fullName = StandardArgumentDefinitions.REFERENCE_LONG_NAME,
+ shortName = StandardArgumentDefinitions.REFERENCE_SHORT_NAME,
+ doc = "Reference for the CRAM file",
+ common = true)
+ @WorkflowOutput
+ public GATKPath referencePath;
+
+ public static final String MISMATCH_RATE_THRESHOLD_ARG_NAME = "mismatch-rate-threshold";
+ @Argument(fullName = MISMATCH_RATE_THRESHOLD_ARG_NAME,
+ shortName = MISMATCH_RATE_THRESHOLD_ARG_NAME,
+ doc = "Mismatch rate threshold above which we consider the file to be corrupt",
+ optional = true)
+ public double mismatchRateThreshold = DEFAULT_MISMATCH_RATE_THRESHOLD;
+
+ public static final String VERBOSE_ARG_NAME = "verbose";
+ @Argument(fullName = VERBOSE_ARG_NAME,
+ shortName= VERBOSE_ARG_NAME,
+ doc="Calculate and print the mismatch rate for all containers",
+ optional=true)
+ public boolean verbose = false;
+
+ public static final String ECHO_ARG_NAME = "echo-to-stdout";
+ @Argument(fullName = ECHO_ARG_NAME,
+ shortName= ECHO_ARG_NAME,
+ doc="Echo text output to stdout",
+ optional=true)
+ public boolean echoToStdout = false;
+
+ private CRAMIssue8768Analyzer cramAnalyzer;
+
+ @Override
+ protected Object doWork() {
+ cramAnalyzer = new CRAMIssue8768Analyzer(
+ inputPath,
+ textOutputPath,
+ tsvOutputPath,
+ referencePath,
+ mismatchRateThreshold,
+ verbose,
+ echoToStdout);
+ cramAnalyzer.doAnalysis();
+ return cramAnalyzer.getRetCode();
+ }
+
+ @Override
+ protected void onShutdown() {
+ if ( cramAnalyzer != null ) {
+ try {
+ cramAnalyzer.close();
+ } catch (Exception e) {
+ throw new RuntimeException(e);
+ }
+ }
+ }
+}
+
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/FlowBasedArgumentCollection.java b/src/main/java/org/broadinstitute/hellbender/tools/FlowBasedArgumentCollection.java
index 3fa41b46a04..695d1385e4c 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/FlowBasedArgumentCollection.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/FlowBasedArgumentCollection.java
@@ -10,7 +10,6 @@ public class FlowBasedArgumentCollection implements Serializable {
private static final long serialVersionUID = 0;
public static final String FLOW_USE_T0_TAG = "flow-use-t0-tag";
- public static final String PROBABILITY_RATIO_THRESHOLD_LONG_NAME = "flow-probability-threshold";
public static final String REMOVE_LONGER_THAN_ONE_INDELS_LONG_NAME = "flow-remove-non-single-base-pair-indels";
public static final String REMOVE_ONE_TO_ZERO_PROBS_LONG_NAME = "flow-remove-one-zero-probs";
public static final String NUMBER_OF_POSSIBLE_PROBS_LONG_NAME = "flow-quantization-bins";
@@ -27,8 +26,7 @@ public class FlowBasedArgumentCollection implements Serializable {
- private static final double DEFAULT_RATIO_THRESHOLD = 0.003;
- private static final double DEFAULT_FILLING_VALUE = 0.001;
+ public static final double DEFAULT_FILLING_VALUE = 0.001;
private static final boolean DEFAULT_REMOVE_LONGER_INDELS = false;
private static final boolean DEFAULT_REMOVE_ONE_TO_ZERO = false;
private static final boolean DEFAULT_SYMMETRIC_INDELS = false;
@@ -45,10 +43,6 @@ public class FlowBasedArgumentCollection implements Serializable {
@Argument(fullName = FLOW_USE_T0_TAG, doc = "Use t0 tag if exists in the read to create flow matrix", optional = true)
public boolean useT0Tag = DEFAULT_FLOW_USE_T0_TAG;
- @Advanced
- @Argument(fullName = PROBABILITY_RATIO_THRESHOLD_LONG_NAME, doc = "Lowest probability ratio to be used as an option", optional = true)
- public double probabilityRatioThreshold = DEFAULT_RATIO_THRESHOLD;
-
@Advanced
@Argument(fullName = REMOVE_LONGER_THAN_ONE_INDELS_LONG_NAME, doc = "Should the probabilities of more then 1 indel be used", optional = true)
public boolean removeLongerThanOneIndels = DEFAULT_REMOVE_LONGER_INDELS;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/GetSampleName.java b/src/main/java/org/broadinstitute/hellbender/tools/GetSampleName.java
index f66556f539d..6aa20a9dc5c 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/GetSampleName.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/GetSampleName.java
@@ -86,18 +86,16 @@ public void onTraversalStart() {
}
final ListOMP_NUM_THREADS and MKL_NUM_THREADS
* environment variables, respectively.THEANO_FLAGS environment variable to override the GATK theano
+ * PYTENSOR_FLAGS environment variable to override the GATK PyTensor
* configuration. For example, by running
- * THEANO_FLAGS="base_compiledir=PATH/TO/BASE_COMPILEDIR" gatk DetermineGermlineContigPloidy ..., users can specify
- * the theano compilation directory (which is set to $HOME/.theano by default). See theano documentation
- * at
- * https://theano-pymc.readthedocs.io/en/latest/library/config.html.
+ * PYTENSOR_FLAGS="base_compiledir=PATH/TO/BASE_COMPILEDIR" gatk DetermineGermlineContigPloidy ..., users can specify
+ * the PyTensor compilation directory (which is set to $HOME/.pytensor by default). See PyTensor documentation
+ * at
+ * https://pytensor.readthedocs.io/en/latest/library/config.html.
* Tool run modes
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/copynumber/GermlineCNVCaller.java b/src/main/java/org/broadinstitute/hellbender/tools/copynumber/GermlineCNVCaller.java
index 095ffb76850..aac440912bd 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/copynumber/GermlineCNVCaller.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/copynumber/GermlineCNVCaller.java
@@ -92,12 +92,12 @@
* OMP_NUM_THREADS and MKL_NUM_THREADS
* environment variables, respectively.THEANO_FLAGS environment variable to override the GATK theano
+ * PYTENSOR_FLAGS environment variable to override the GATK PyTensor
* configuration. For example, by running
- * THEANO_FLAGS="base_compiledir=PATH/TO/BASE_COMPILEDIR" gatk GermlineCNVCaller ..., users can specify
- * the theano compilation directory (which is set to $HOME/.theano by default). See theano documentation
- * at
- * https://theano-pymc.readthedocs.io/en/latest/library/config.html.
+ * PYTENSOR_FLAGS="base_compiledir=PATH/TO/BASE_COMPILEDIR" gatk DetermineGermlineContigPloidy ..., users can specify
+ * the PyTensor compilation directory (which is set to $HOME/.pytensor by default). See PyTensor documentation
+ * at
+ * https://pytensor.readthedocs.io/en/latest/library/config.html.
* Resource usage
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/copynumber/PostprocessGermlineCNVCalls.java b/src/main/java/org/broadinstitute/hellbender/tools/copynumber/PostprocessGermlineCNVCalls.java
index fa68bb4f34d..b9ac8561e59 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/copynumber/PostprocessGermlineCNVCalls.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/copynumber/PostprocessGermlineCNVCalls.java
@@ -89,12 +89,12 @@
* the python environment is already set up. Otherwise, the environment must be created and activated as described in the
* main GATK README.md file.
Advanced users may wish to set the THEANO_FLAGS environment variable to override the GATK theano
+ *
Advanced users may wish to set the PYTENSOR_FLAGS environment variable to override the GATK PyTensor
* configuration. For example, by running
- * THEANO_FLAGS="base_compiledir=PATH/TO/BASE_COMPILEDIR" gatk PostprocessGermlineCNVCalls ..., users can specify
- * the theano compilation directory (which is set to $HOME/.theano by default). See theano documentation
- * at
- * https://theano-pymc.readthedocs.io/en/latest/library/config.html.
+ * PYTENSOR_FLAGS="base_compiledir=PATH/TO/BASE_COMPILEDIR" gatk DetermineGermlineContigPloidy ..., users can specify
+ * the PyTensor compilation directory (which is set to $HOME/.pytensor by default). See PyTensor documentation
+ * at
+ * https://pytensor.readthedocs.io/en/latest/library/config.html.
*
Main class for SV clustering. Two items are clustered together if they:
*This tool outputs no annotations by default, all annotations/groups must be specified explicitly.
+ *This tool outputs no annotations by default, all annotations/groups must be specified explicitly. This tool accepts VCF format files only. Using GVCF files as input may result in unexpected behavior.
* *RankSumAnnotations produced by this tool are not the same as those produced by the HaplotypeCaller. Without the
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
+ * 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 ListThe {@code GtfInfo.Type} enum specifies whether the features is a + * {@code GENE} or a {@code TRANSCRIPT}.
+ * + *The interval specifies the feature's contig (chromosome), start position, + * and end position, providing the precise location of the gene or transcript + * on the genome.
+ * + *Example usage:
+ *
+ * Interval interval = new Interval("chr1", 1000, 2000);
+ * GtfInfo gtfInfo = new GtfInfo(interval, GtfInfo.Type.GENE, "MAPK1");
+ *
+ **/
+
+public class GtfInfo {
+
+ public enum Type {
+ GENE,
+ TRANSCRIPT
+ }
+
+ private final Type type;
+ private final String geneName;
+ private final Interval interval;
+
+ public GtfInfo(Interval interval, Type type, String geneName) {
+ this.interval = interval;
+ this.type = type;
+ this.geneName = geneName;
+ }
+
+ public Type getType() {
+ return type;
+ }
+
+ public String getGeneName() {
+ return geneName;
+ }
+
+ public Interval getInterval() {
+ return interval;
+ }
+
+ public Integer getStart() {
+ return interval.getStart();
+ }
+
+ public Integer getEnd() {
+ return interval.getEnd();
+ }
+
+ @Override
+ public String toString() {
+ return "GtfInfo{ " + "type = " + type + " geneName = " + geneName + "interval = " + interval;
+ }
+
+}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/conversion/GtfToBed.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/conversion/GtfToBed.java
new file mode 100644
index 00000000000..f1610405d5b
--- /dev/null
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/conversion/GtfToBed.java
@@ -0,0 +1,322 @@
+package org.broadinstitute.hellbender.tools.walkers.conversion;
+
+import htsjdk.samtools.SAMSequenceDictionary;
+import htsjdk.samtools.util.*;
+import htsjdk.tribble.Feature;
+import org.apache.log4j.LogManager;
+import org.apache.log4j.Logger;
+import org.broadinstitute.barclay.argparser.Argument;
+import org.broadinstitute.barclay.argparser.CommandLineProgramProperties;
+import org.broadinstitute.barclay.argparser.WorkflowProperties;
+import org.broadinstitute.barclay.help.DocumentedFeature;
+import org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions;
+import org.broadinstitute.hellbender.cmdline.programgroups.ShortVariantDiscoveryProgramGroup;
+import org.broadinstitute.hellbender.engine.*;
+import org.broadinstitute.hellbender.exceptions.GATKException;
+import org.broadinstitute.hellbender.exceptions.UserException;
+import org.broadinstitute.hellbender.utils.Utils;
+import org.broadinstitute.hellbender.utils.codecs.gtf.GencodeGtfFeature;
+
+import java.io.IOException;
+import java.io.OutputStream;
+import java.nio.file.Files;
+import java.util.*;
+
+/**
+ * Convert Gencode GTF files to BED format with options for gene and transcript level processing. + * This tool allows for the extraction of gene and transcript information from Gencode GTF files and + * outputs the data in BED format.
+ * + * + *The conversion process includes sorting entries + * by karyotype, providing flexibility in the selection of either gene or transcript level + * data, and an option to only use basic tags. It ensures that the BED output is sorted and formatted correctly for subsequent use. + * Note that it has been tested for both human and mouse Gencode GTFs.
+ * + *Example commands to run GtfToBed for typical scenarios:
+ * + *This mode extracts and converts gene data from the input GTF file to BED format:
+ * + *+ * java -jar GtfToBed.jar \ + * -gtf-path input.gtf \ + * -sequence-dictionary dictionary.dict \ += * -output output.bed \ + *+ * + *
This mode extracts and converts transcript data from the input GTF file to BED format:
+ * + *+ * java -jar GtfToBed.jar \ + * -gtf-path input.gtf \ + * -sequence-dictionary dictionary.dict \ + * -sort-by-transcript \ + * -output output.bed \ + *+ * + *
This mode extracts and converts basic transcript data from the input GTF file to BED format:
+ * + *+ * java -jar GtfToBed.jar \ + * -gtf-path input.gtf \ + * -sequence-dictionary dictionary.dict \ + * -sort-by-transcript \ + * -use-basic-transcript \ + * -output output.bed \ + *+ * + *
This mode extracts and converts transcript data from the input GTF file to BED format using a FASTA file:
+ * + *+ * java -jar GtfToBed.jar \ + * -gtf-path input.gtf \ + * -reference ref.fa \ + * -sort-by-transcript \ + * -output output.bed \ + *+ */ + +@CommandLineProgramProperties( + summary = "Converts Gencode GTF files to Bed file format with each row of bed file being either a gene or a transcript.", + oneLineSummary = "Gencode GTF to BED", + programGroup = ShortVariantDiscoveryProgramGroup.class +) + +@DocumentedFeature +@WorkflowProperties +public class GtfToBed extends FeatureWalker
Mutect2 automatically sets parameters appropriately for calling on mitochondria with the
Mutect2 automatically sets parameters appropriately for calling on mitochondria with the
* gatk Mutect2 \
@@ -262,7 +265,8 @@ public boolean shouldTrackPileupsForAssemblyRegions() {
@Override
public void onTraversalStart() {
VariantAnnotatorEngine annotatorEngine = new VariantAnnotatorEngine(makeVariantAnnotations(), null, Collections.emptyList(), false, false);
- m2Engine = new Mutect2Engine(MTAC, assemblyRegionArgs, createOutputBamIndex, createOutputBamMD5, getHeaderForReads(), referenceArguments.getReferenceSpecifier(), annotatorEngine);
+ m2Engine = new Mutect2Engine(MTAC, assemblyRegionArgs, createOutputBamIndex, createOutputBamMD5, getHeaderForReads(),
+ getBestAvailableSequenceDictionary(), referenceArguments.getReferenceSpecifier(), annotatorEngine);
vcfWriter = createVCFWriter(outputVCF);
if (m2Engine.emitReferenceConfidence()) {
logger.warn("Note that the Mutect2 reference confidence mode is in BETA -- the likelihoods model and output format are subject to change in subsequent versions.");
@@ -281,17 +285,6 @@ public void onTraversalStart() {
m2Engine.writeHeader(vcfWriter, getDefaultToolVCFHeaderLines());
}
- @Override
- public Collection makeVariantAnnotations(){
- final Collection annotations = super.makeVariantAnnotations();
-
- if (MTAC.mitochondria) {
- annotations.add(new OriginalAlignment());
- }
-
- return annotations;
- }
-
@Override
public Object onTraversalSuccess() {
m2Engine.writeExtraOutputs(new File(outputVCF + DEFAULT_STATS_EXTENSION));
@@ -320,6 +313,12 @@ public void closeTool() {
*/
@Override
protected String[] customCommandLineValidation() {
+ if (MTAC.mitochondria) {
+ ModeArgumentUtils.setArgValues(
+ getCommandLineParser(),
+ MTAC.getMitochondriaModeNameValuePairs(),
+ M2ArgumentCollection.MITOCHONDRIA_MODE_LONG_NAME);
+ }
if (MTAC.flowMode != M2ArgumentCollection.FlowMode.NONE) {
ModeArgumentUtils.setArgValues(
getCommandLineParser(),
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java
index 1d2096f3c2b..c831a4169e0 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java
@@ -1,6 +1,7 @@
package org.broadinstitute.hellbender.tools.walkers.mutect;
import htsjdk.samtools.SAMFileHeader;
+import htsjdk.samtools.SAMSequenceDictionary;
import htsjdk.samtools.util.Locatable;
import htsjdk.variant.variantcontext.Allele;
import htsjdk.variant.variantcontext.Genotype;
@@ -72,7 +73,7 @@
public final class Mutect2Engine implements AssemblyRegionEvaluator, AutoCloseable {
private static final List STANDARD_MUTECT_INFO_FIELDS = Arrays.asList(GATKVCFConstants.NORMAL_LOG_10_ODDS_KEY, GATKVCFConstants.TUMOR_LOG_10_ODDS_KEY, GATKVCFConstants.NORMAL_ARTIFACT_LOG_10_ODDS_KEY,
- GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY, GATKVCFConstants.IN_PON_KEY, GATKVCFConstants.POPULATION_AF_KEY,
+ GATKVCFConstants.EVENT_COUNT_IN_REGION_KEY, GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY, GATKVCFConstants.IN_PON_KEY, GATKVCFConstants.POPULATION_AF_KEY,
GATKVCFConstants.GERMLINE_QUAL_KEY, GATKVCFConstants.CONTAMINATION_QUAL_KEY, GATKVCFConstants.SEQUENCING_QUAL_KEY,
GATKVCFConstants.POLYMERASE_SLIPPAGE_QUAL_KEY, GATKVCFConstants.READ_ORIENTATION_QUAL_KEY,
GATKVCFConstants.STRAND_QUAL_KEY, GATKVCFConstants.ORIGINAL_CONTIG_MISMATCH_KEY, GATKVCFConstants.N_COUNT_KEY, GATKVCFConstants.AS_UNIQUE_ALT_READ_SET_COUNT_KEY);
@@ -94,6 +95,7 @@ public final class Mutect2Engine implements AssemblyRegionEvaluator, AutoCloseab
final private M2ArgumentCollection MTAC;
private SAMFileHeader header;
+ private SAMSequenceDictionary sequenceDictionary;
private final int minCallableDepth;
public static final String CALLABLE_SITES_NAME = "callable";
@@ -136,9 +138,12 @@ public final class Mutect2Engine implements AssemblyRegionEvaluator, AutoCloseab
* @param referenceSpec reference specifier for the reference
* @param annotatorEngine annotator engine built with desired annotations
*/
- public Mutect2Engine(final M2ArgumentCollection MTAC, AssemblyRegionArgumentCollection assemblyRegionArgs, final boolean createBamOutIndex, final boolean createBamOutMD5, final SAMFileHeader header, final GATKPath referenceSpec, final VariantAnnotatorEngine annotatorEngine) {
+ public Mutect2Engine(final M2ArgumentCollection MTAC, AssemblyRegionArgumentCollection assemblyRegionArgs,
+ final boolean createBamOutIndex, final boolean createBamOutMD5, final SAMFileHeader header,
+ final SAMSequenceDictionary sequenceDictionary, final GATKPath referenceSpec, final VariantAnnotatorEngine annotatorEngine) {
this.MTAC = Utils.nonNull(MTAC);
this.header = Utils.nonNull(header);
+ this.sequenceDictionary = sequenceDictionary;
minCallableDepth = MTAC.callableDepth;
referenceReader = ReferenceUtils.createReferenceReader(Utils.nonNull(referenceSpec));
aligner = SmithWatermanAligner.getAligner(MTAC.smithWatermanImplementation);
@@ -162,7 +167,7 @@ public Mutect2Engine(final M2ArgumentCollection MTAC, AssemblyRegionArgumentColl
annotationEngine = Utils.nonNull(annotatorEngine);
assemblyEngine = MTAC.createReadThreadingAssembler();
likelihoodCalculationEngine = AssemblyBasedCallerUtils.createLikelihoodCalculationEngine(MTAC.likelihoodArgs, MTAC.fbargs, true, MTAC.likelihoodArgs.likelihoodEngineImplementation);
- genotypingEngine = new SomaticGenotypingEngine(MTAC, normalSamples, annotationEngine);
+ genotypingEngine = new SomaticGenotypingEngine(MTAC, normalSamples, annotationEngine, header, sequenceDictionary);
haplotypeBAMWriter = AssemblyBasedCallerUtils.createBamWriter(MTAC, createBamOutIndex, createBamOutMD5, header);
trimmer = new AssemblyRegionTrimmer(assemblyRegionArgs, header.getSequenceDictionary());
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect3DatasetEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect3DatasetEngine.java
index 942c2de15b2..ad708c8dac8 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect3DatasetEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect3DatasetEngine.java
@@ -1,5 +1,9 @@
package org.broadinstitute.hellbender.tools.walkers.mutect;
+import htsjdk.samtools.SAMFileHeader;
+import htsjdk.samtools.SAMReadGroupRecord;
+import htsjdk.samtools.SAMSequenceDictionary;
+import htsjdk.samtools.SAMSequenceRecord;
import htsjdk.variant.variantcontext.Allele;
import htsjdk.variant.variantcontext.Genotype;
import htsjdk.variant.variantcontext.VariantContext;
@@ -43,6 +47,10 @@ private enum Label {
ARTIFACT, VARIANT, UNLABELED, IGNORE
}
+ private final SAMSequenceDictionary sequenceDictionary;
+
+ private final Map readGroupIndices = new HashMap<>();
+
// number of additional variant features for assembly complexity, reference context
private static final int NUM_EXTRA_FEATURES = 9;
@@ -65,6 +73,8 @@ private enum Label {
private static final int MIN_REF = 5;
private final PrintWriter printWriter;
+ private final PrintWriter contigPrintWriter;
+ private final PrintWriter readGroupPrintWriter;
// number of nonartifact data to keep for each artifact datum
private final int nonArtifactPerArtifact;
@@ -79,9 +89,15 @@ private enum Label {
private final EnumMap> unmatchedArtifactAltCounts;
- public Mutect3DatasetEngine(final File datasetFile, final boolean trainingMode, final int maxRefCount, final int maxAltCount, final int nonArtifactPerArtifact, final Set normalSamples) {
+ public Mutect3DatasetEngine(final File datasetFile, final boolean trainingMode, final int maxRefCount,
+ final int maxAltCount, final int nonArtifactPerArtifact, final Set normalSamples,
+ final SAMFileHeader header, final SAMSequenceDictionary sequenceDictionary) {
try {
printWriter = new PrintWriter(new FileWriter(Utils.nonNull(datasetFile)));
+ final File contigTableFile = datasetFile.toPath().resolveSibling("contigs.table").toFile();
+ final File readGroupTableFile = datasetFile.toPath().resolveSibling("read-groups.table").toFile();
+ contigPrintWriter = new PrintWriter(new FileWriter(contigTableFile));
+ readGroupPrintWriter = new PrintWriter(new FileWriter(readGroupTableFile));
} catch (IOException ex) {
throw new UserException.BadInput("Could not create dataset file writer");
}
@@ -92,6 +108,12 @@ public Mutect3DatasetEngine(final File datasetFile, final boolean trainingMode,
this.maxRefCount = maxRefCount;
this.maxAltCount = maxAltCount;
+ this.sequenceDictionary = sequenceDictionary;
+ final List readGroups = header.getReadGroups();
+ for (int n = 0; n < readGroups.size(); n++) {
+ readGroupIndices.put(readGroups.get(n).getReadGroupId(), n);
+ }
+
unmatchedArtifactAltCounts = new EnumMap<>(VariantType.class);
for (final VariantType type : VariantType.values()) {
unmatchedArtifactAltCounts.put(type, new ArrayBlockingQueue<>(CAPACITY));
@@ -106,7 +128,7 @@ public void addData(final ReferenceContext ref, final VariantContext vc, Optiona
final M2ArgumentCollection.Mutect3DatasetMode mutect3DatasetMode) {
final String refBases = ReferenceBases.annotate(ref, vc);
final String refAllele = vc.getReference().getBaseString();
- final String contig = vc.getContig();
+ final int contigIndex = sequenceDictionary.getSequenceIndex(vc.getContig());
final int position = vc.getStart();
final Set tumorSamples = likelihoods.samples().stream().filter(sample -> !normalSamples.contains(sample)).collect(Collectors.toSet());
final int numAlt = vc.getNAlleles() - 1;
@@ -153,7 +175,7 @@ public void addData(final ReferenceContext ref, final VariantContext vc, Optiona
final int diff = altAlleleString.length() - refAllele.length();
final VariantType type = diff == 0 ? VariantType.SNV : ( diff > 0 ? VariantType.INSERTION : VariantType.DELETION);
- if (trainingMode) {
+ if (trainingMode) { // training mode -- collecting tensors to train the Permutect artifact model
final ArrayBlockingQueue unmatchedQueue = unmatchedArtifactAltCounts.get(type);
final boolean likelySeqError = tumorLods[n] < TLOD_THRESHOLD;
@@ -182,8 +204,13 @@ public void addData(final ReferenceContext ref, final VariantContext vc, Optiona
} else {
labels.add(Label.IGNORE);
}
- } else {
- labels.add(Label.UNLABELED);
+ } else { // not training mode -- we are generating tensors in order to apply the Permutect artifact model to a callset
+ if (truthVCs.isPresent()) {
+ // here, for the purposes of test data, both sequencing errors and technical artifacts get the "ARTIFACT" label
+ labels.add(truthAlleles.contains(remappedAltAlelle) ? Label.VARIANT : Label.ARTIFACT);
+ } else {
+ labels.add(Label.UNLABELED);
+ }
}
}
@@ -199,9 +226,9 @@ public void addData(final ReferenceContext ref, final VariantContext vc, Optiona
// TODO: for now we don't really need normal reads
// note that the following use the VC's allele order, not necessarily the likelihoods' allele order
final List>> normalReadVectorsByAllele = FeaturizedReadSets.getReadVectors(vc, normalSamples,
- likelihoods, logFragmentLikelihoods, maxRefCount, maxAltCount, mutect3DatasetMode);
+ likelihoods, logFragmentLikelihoods, maxRefCount, maxAltCount, mutect3DatasetMode, readGroupIndices);
final List>> tumorReadVectorsByAllele = FeaturizedReadSets.getReadVectors(vc, tumorSamples,
- likelihoods, logFragmentLikelihoods, maxRefCount, maxAltCount, altDownsampleMap, mutect3DatasetMode);
+ likelihoods, logFragmentLikelihoods, maxRefCount, maxAltCount, altDownsampleMap, mutect3DatasetMode, readGroupIndices);
// ref and alt reads have already been downsampled by the read featurizer
final List> tumorRefReads = tumorReadVectorsByAllele.get(0);
@@ -222,7 +249,7 @@ public void addData(final ReferenceContext ref, final VariantContext vc, Optiona
final List> normalAltReads = normalReadVectorsByAllele.get(n+1);
printWriter.println(labels.get(n).toString());
- printWriter.printf("%s:%d,%s->%s%n", contig, position, refAllele, altAllele);
+ printWriter.printf("%d:%d,%s->%s%n", contigIndex, position, refAllele, altAllele);
printWriter.println(refBases);
printWriter.println(numberString(variantFeatureVector, "%.2f", " "));
//printWriter.printf("%d %d %d %d%n", tumorRefReads.size(), tumorAltReads.size(), normalRefReads.size(), normalAltReads.size());
@@ -322,5 +349,16 @@ private int[] sumADsOverSamples(final VariantContext vc, final Set sampl
@Override
public void close() {
printWriter.close();
+
+ for (final SAMSequenceRecord contigRecord : sequenceDictionary.getSequences()) {
+ contigPrintWriter.println(String.format("%s\t%d", contigRecord.getContig(), contigRecord.getSequenceIndex()));
+ }
+
+ for (final Map.Entry entry : readGroupIndices.entrySet()) {
+ readGroupPrintWriter.println(String.format("%s\t%d", entry.getKey(), entry.getValue()));
+ }
+
+ contigPrintWriter.close();
+ readGroupPrintWriter.close();
}
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/SomaticGenotypingEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/SomaticGenotypingEngine.java
index e4504398998..55bef4dc806 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/SomaticGenotypingEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/SomaticGenotypingEngine.java
@@ -4,10 +4,12 @@
import com.google.common.collect.ImmutableMap;
import com.google.common.primitives.Doubles;
import htsjdk.samtools.SAMFileHeader;
+import htsjdk.samtools.SAMSequenceDictionary;
import htsjdk.samtools.util.Locatable;
import htsjdk.variant.variantcontext.*;
import htsjdk.variant.vcf.VCFConstants;
import org.apache.commons.collections4.ListUtils;
+import org.apache.commons.lang3.mutable.MutableInt;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor;
import org.apache.commons.math3.linear.RealMatrix;
@@ -51,7 +53,9 @@ public class SomaticGenotypingEngine implements AutoCloseable {
private final double refPseudocount = 1;
private final double altPseudocount;
- public SomaticGenotypingEngine(final M2ArgumentCollection MTAC, final Set normalSamples, final VariantAnnotatorEngine annotationEngine) {
+ public SomaticGenotypingEngine(final M2ArgumentCollection MTAC, final Set normalSamples,
+ final VariantAnnotatorEngine annotationEngine,
+ final SAMFileHeader header, final SAMSequenceDictionary sequenceDictionary) {
this.MTAC = MTAC;
altPseudocount = MTAC.minAF == 0.0 ? 1 : 1 - Math.log(2)/Math.log(MTAC.minAF);
@@ -61,7 +65,7 @@ public SomaticGenotypingEngine(final M2ArgumentCollection MTAC, final Set logFragmentLikelihoods = logReadLikelihoods.groupEvidence(MTAC.independentMates ? read -> read : GATKRead::getName, Fragment::createAndAvoidFailure);
+ final Set potentialSomaticEventsInRegion = new HashSet<>();
for( final int loc : eventStarts ) {
final List eventsAtThisLoc = AssemblyBasedCallerUtils.getVariantsFromActiveHaplotypes(loc, haplotypes, false);
- VariantContext mergedVC = AssemblyBasedCallerUtils.makeMergedVariantContext(eventsAtThisLoc);
- if( mergedVC == null ) {
+ VariantContext merged = AssemblyBasedCallerUtils.makeMergedVariantContext(eventsAtThisLoc);
+ if( merged == null ) {
continue;
}
+ final VariantContext mergedVC = emitRefConf ? ReferenceConfidenceUtils.addNonRefSymbolicAllele(merged) : merged;
// converting haplotype likelihoods to allele likelihoods
final Map> alleleMapper = AssemblyBasedCallerUtils.createAlleleMapper(mergedVC, loc, haplotypes, true);
@@ -127,7 +133,6 @@ public CalledHaplotypes callMutations(
logLikelihoods.retainEvidence(variantCallingRelevantFragmentOverlap::overlaps);
if (emitRefConf) {
- mergedVC = ReferenceConfidenceUtils.addNonRefSymbolicAllele(mergedVC);
logLikelihoods.addNonReferenceAllele(Allele.NON_REF_ALLELE);
}
final List> tumorMatrices = IntStream.range(0, logLikelihoods.numberOfSamples())
@@ -152,13 +157,21 @@ public CalledHaplotypes callMutations(
.filter(allele -> forcedAlleles.contains(allele) || tumorLogOdds.getAlt(allele) > MTAC.getEmissionLogOdds())
.collect(Collectors.toList());
- final long somaticAltCount = tumorAltAlleles.stream()
+ final List allelesToGenotype = tumorAltAlleles.stream()
.filter(allele -> forcedAlleles.contains(allele) || !hasNormal || MTAC.genotypeGermlineSites || normalLogOdds.getAlt(allele) > MathUtils.log10ToLog(MTAC.normalLog10Odds))
- .count();
+ .toList();
+
+ // record somatic alleles for later use in the Event Count annotation
+ // note that in tumor-only calling it does not attempt to detect germline events
+ mergedVC.getAlternateAlleles().stream()
+ .filter(allele -> tumorLogOdds.getAlt(allele) > MTAC.getEmissionLogOdds())
+ .filter(allele -> !hasNormal || normalLogOdds.getAlt(allele) > MathUtils.log10ToLog(MTAC.normalLog10Odds))
+ .map(allele -> new Event(mergedVC.getContig(), mergedVC.getStart(), mergedVC.getReference(), allele))
+ .forEach(potentialSomaticEventsInRegion::add);
// if every alt allele is germline, skip this variant. However, if some alt alleles are germline and others
// are not we emit them all so that the filtering engine can see them
- if (somaticAltCount == 0) {
+ if (allelesToGenotype.isEmpty()) {
continue;
}
@@ -222,8 +235,41 @@ public CalledHaplotypes callMutations(
final List outputCalls = AssemblyBasedCallerUtils.phaseCalls(returnCalls, calledHaplotypes);
final int eventCount = outputCalls.size();
+
+ // calculate the number of somatic events in the best haplotype of each called variant
+ final Map haplotypeSupportCounts = logReadLikelihoods.alleles().stream()
+ .collect(Collectors.toMap(hap -> hap, label -> new MutableInt(0)));
+ logReadLikelihoods.bestAllelesBreakingTies()
+ .forEach(bestHaplotype -> haplotypeSupportCounts.get(bestHaplotype.allele).increment());
+
+ final Map> haplotypesByEvent= new HashMap<>();
+ for (final Haplotype haplotype : logReadLikelihoods.alleles()) {
+ for (final Event event : haplotype.getEventMap().getEvents()) {
+ haplotypesByEvent.computeIfAbsent(event, e -> new ArrayList<>()).add(haplotype);
+ }
+ }
+ final Map> eventCountAnnotations = new HashMap<>();
+ for (final VariantContext outputCall : outputCalls) {
+ for (final Allele allele : outputCall.getAlternateAlleles()) {
+ // note: this creates the minimal representation behind the scenes
+ final Event event = new Event(outputCall.getContig(), outputCall.getStart(), outputCall.getReference(), allele);
+ // haplotypesByEvent contains every *assembled* event, including events injected into the original assembly,
+ // but there are some modes where we genotype events that were never in an assembly graph, in which case
+ // this annotation is irrelevant
+ if (haplotypesByEvent.containsKey(event)) {
+ final Haplotype bestHaplotype = haplotypesByEvent.get(event).stream()
+ .sorted(Comparator.comparingInt(h -> haplotypeSupportCounts.getOrDefault(h, new MutableInt(0)).intValue()).reversed())
+ .findFirst().get();
+
+ eventCountAnnotations.computeIfAbsent(outputCall, vc -> new ArrayList<>())
+ .add((int) bestHaplotype.getEventMap().getEvents().stream().filter(potentialSomaticEventsInRegion::contains).count());
+ }
+ }
+ }
final List outputCallsWithEventCountAnnotation = outputCalls.stream()
- .map(vc -> new VariantContextBuilder(vc).attribute(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY, eventCount).make())
+ .map(vc -> new VariantContextBuilder(vc)
+ .attribute(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY, eventCountAnnotations.get(vc))
+ .attribute(GATKVCFConstants.EVENT_COUNT_IN_REGION_KEY, potentialSomaticEventsInRegion.size()).make())
.collect(Collectors.toList());
return new CalledHaplotypes(outputCallsWithEventCountAnnotation, calledHaplotypes);
}
@@ -361,9 +407,7 @@ private Optional getForNormal(final Supplier supplier) {
private static Map getNegativeLogPopulationAFAnnotation(List germlineResourceVariants,
final List allAlleles,
final double afOfAllelesNotInGermlineResource) {
- final Optional germlineVC = germlineResourceVariants.isEmpty() ? Optional.empty()
- : Optional.of(germlineResourceVariants.get(0)); // assume only one VC per site
- final double[] populationAlleleFrequencies = getGermlineAltAlleleFrequencies(allAlleles, germlineVC, afOfAllelesNotInGermlineResource);
+ final double[] populationAlleleFrequencies = getGermlineAltAlleleFrequencies(allAlleles, germlineResourceVariants, afOfAllelesNotInGermlineResource);
return ImmutableMap.of(GATKVCFConstants.POPULATION_AF_KEY, MathUtils.applyToArray(populationAlleleFrequencies, x -> - Math.log10(x)));
}
@@ -373,27 +417,35 @@ private static Map getNegativeLogPopulationAFAnnotation(List allAlleles, final Optional germlineVC, final double afOfAllelesNotInGermlineResource) {
+ static double[] getGermlineAltAlleleFrequencies(final List allAlleles, final List germlineVCs, final double afOfAllelesNotInGermlineResource) {
Utils.validateArg(!allAlleles.isEmpty(), "allAlleles are empty -- there is not even a reference allele.");
- if (germlineVC.isPresent()) {
- if (! germlineVC.get().hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY)) {
- logger.warn("Germline resource variant at " + germlineVC.get().getContig() + ":" + germlineVC.get().getStart() +" missing AF attribute");
- return Doubles.toArray(Collections.nCopies(allAlleles.size() - 1, afOfAllelesNotInGermlineResource));
+ final Map alleleFrequencies = new HashMap<>();
+ allAlleles.forEach(a -> alleleFrequencies.put(a, afOfAllelesNotInGermlineResource)); // initialize everything to the default
+
+ // look through every germline resource variant context at this locus and fill in the AFs
+ for (final VariantContext germlineVC : germlineVCs) {
+ if (! germlineVC.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY)) {
+ logger.warn("Germline resource variant at " + germlineVC.getContig() + ":" + germlineVC.getStart() +" missing AF attribute");
+ }
+
+ List germlineIndices = GATKVariantContextUtils.alleleIndices(allAlleles, germlineVC.getAlleles());
+ final List germlineAltAFs = Mutect2Engine.getAttributeAsDoubleList(germlineVC, VCFConstants.ALLELE_FREQUENCY_KEY, afOfAllelesNotInGermlineResource);
+
+ if (germlineAltAFs.size() == (germlineVC.getNAlleles() - 1)) { // skip VCs with a bad AF field that got parsed as a wrong-length list
+ for (int alleleIndex = 1; alleleIndex < allAlleles.size(); alleleIndex++) { // start at 1 to skip the reference, which doesn't have an AF annotation
+ final Allele allele = allAlleles.get(alleleIndex);
+ // note the -1 since germlineAltAFs do not include ref
+ germlineIndices.get(alleleIndex).ifPresent(germlineIndex -> alleleFrequencies.put(allele, germlineAltAFs.get(germlineIndex - 1)));
+ }
}
- List germlineIndices = GATKVariantContextUtils.alleleIndices(allAlleles, germlineVC.get().getAlleles());
- final List germlineAltAFs = Mutect2Engine.getAttributeAsDoubleList(germlineVC.get(), VCFConstants.ALLELE_FREQUENCY_KEY, afOfAllelesNotInGermlineResource);
-
- return germlineIndices.stream().skip(1) // skip the reference allele
- .mapToDouble(idx -> idx.isPresent() ? germlineAltAFs.get(idx.getAsInt() - 1) : afOfAllelesNotInGermlineResource) // note the -1 since germlineAltAFs do not include ref
- .toArray();
- } else {
- return Doubles.toArray(Collections.nCopies(allAlleles.size() - 1, afOfAllelesNotInGermlineResource));
}
+
+ return allAlleles.stream().skip(1).mapToDouble(alleleFrequencies::get).toArray(); // skip the reference allele
}
/**
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/ClusteredEventsFilter.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/ClusteredEventsFilter.java
index b7212562bb2..88cb589ef29 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/ClusteredEventsFilter.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/ClusteredEventsFilter.java
@@ -8,9 +8,11 @@
public class ClusteredEventsFilter extends HardFilter {
private final int maxEventsInRegion;
+ private final int maxEventsInHaplotype;
- public ClusteredEventsFilter(final int maxEventsInRegion) {
+ public ClusteredEventsFilter(final int maxEventsInRegion, final int maxEventsInHaplotype) {
this.maxEventsInRegion = maxEventsInRegion;
+ this.maxEventsInHaplotype = maxEventsInHaplotype;
}
@Override
@@ -18,8 +20,9 @@ public ClusteredEventsFilter(final int maxEventsInRegion) {
@Override
public boolean isArtifact(final VariantContext vc, final Mutect2FilteringEngine filteringEngine) {
- final Integer eventCount = vc.getAttributeAsInt(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY, -1);
- return eventCount > maxEventsInRegion;
+ final List haplotypeEventCounts = vc.getAttributeAsIntList(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY, 0);
+ final int regionEventCounts = vc.getAttributeAsInt(GATKVCFConstants.EVENT_COUNT_IN_REGION_KEY, 0);
+ return haplotypeEventCounts.stream().mapToInt(n -> n).max().getAsInt() > maxEventsInHaplotype || regionEventCounts > maxEventsInRegion;
}
@Override
@@ -28,5 +31,5 @@ public String filterName() {
}
@Override
- protected List requiredInfoAnnotations() { return Collections.singletonList(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY); }
+ protected List requiredInfoAnnotations() { return List.of(GATKVCFConstants.EVENT_COUNT_IN_REGION_KEY, GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY); }
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/FilteredHaplotypeFilter.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/FilteredHaplotypeFilter.java
index 9a6375931ae..b725ee2587e 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/FilteredHaplotypeFilter.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/FilteredHaplotypeFilter.java
@@ -13,6 +13,7 @@
import java.util.*;
public class FilteredHaplotypeFilter extends Mutect2VariantFilter {
+ private static final double GERMLINE_PROBABILITY_TO_IGNORE_NORMAL_ARTIFACT = 0.25;
private final double maxIntraHaplotypeDistance;
// for each pgt + pid phasing string, a list of loci-error probability pairs
@@ -54,10 +55,21 @@ public double calculateErrorProbability(final VariantContext vc, final Mutect2Fi
@Override
protected void accumulateDataForLearning(final VariantContext vc, final ErrorProbabilities errorProbabilities, final Mutect2FilteringEngine filteringEngine) {
- // we record the maximum non-sequencing artifact that is not this filter itself
- final double artifactProbability = errorProbabilities.getProbabilitiesByFilter().entrySet().stream()
- .filter(e -> e.getKey().errorType() != ErrorType.SEQUENCING)
- .filter(e -> !e.getKey().filterName().equals(filterName()))
+ // we record the maximum non-sequencing, non-germline, artifact probability that is not from this filter itself
+ final Map> probabilitiesByFilter = errorProbabilities.getProbabilitiesByFilter();
+
+ final double germlineProbability = probabilitiesByFilter.entrySet().stream()
+ .filter(e -> e.getKey().filterName() == GATKVCFConstants.GERMLINE_RISK_FILTER_NAME)
+ .flatMap(e -> e.getValue().stream()) // the value is a list of double, we need the max of all the lists
+ .max(Double::compareTo).orElse(0.0);
+
+ // the normal artifact filter often lights up when there's a non-artifactual germline event, which we don't want here
+ final boolean ignoreNormalArtifact = germlineProbability > GERMLINE_PROBABILITY_TO_IGNORE_NORMAL_ARTIFACT;
+
+ final double artifactProbability = probabilitiesByFilter.entrySet().stream()
+ .filter(e -> e.getKey().errorType() != ErrorType.NON_SOMATIC)
+ .filter(e -> !(ignoreNormalArtifact && e.getKey().filterName() == GATKVCFConstants.ARTIFACT_IN_NORMAL_FILTER_NAME))
+ .filter(e -> !e.getKey().filterName().equals(filterName())) // exclude the haplotype filter itself, which would be circular
.flatMap(e -> e.getValue().stream()) // the value is a list of double, we need the max of all the lists
.max(Double::compareTo).orElse(0.0);
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/M2FiltersArgumentCollection.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/M2FiltersArgumentCollection.java
index cec046e4026..26218b5008d 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/M2FiltersArgumentCollection.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/M2FiltersArgumentCollection.java
@@ -55,6 +55,7 @@ public class M2FiltersArgumentCollection {
* Hard filter thresholds
*/
public static final String MAX_EVENTS_IN_REGION_LONG_NAME = "max-events-in-region";
+ public static final String MAX_EVENTS_IN_HAPLOTYPE_LONG_NAME = "max-events-in-haplotype";
public static final String MAX_ALT_ALLELE_COUNT_LONG_NAME = "max-alt-allele-count";
public static final String UNIQUE_ALT_READ_COUNT_LONG_NAME = "unique-alt-read-count";
public static final String MIN_MEDIAN_MAPPING_QUALITY_LONG_NAME = "min-median-mapping-quality";
@@ -65,7 +66,8 @@ public class M2FiltersArgumentCollection {
public static final String MIN_READS_ON_EACH_STRAND_LONG_NAME = "min-reads-per-strand";
public static final String MIN_AF_LONG_NAME = "min-allele-fraction";
- private static final int DEFAULT_MAX_EVENTS_IN_REGION = 2;
+ private static final int DEFAULT_MAX_EVENTS_IN_REGION = 3;
+ private static final int DEFAULT_MAX_EVENTS_IN_HAPLOTYPE = 2;
private static final int DEFAULT_MAX_ALT_ALLELES = 1;
private static final int DEFAULT_MIN_UNIQUE_ALT_READS = 0;
private static final int DEFAULT_MIN_MEDIAN_MAPPING_QUALITY = 30;
@@ -77,9 +79,12 @@ public class M2FiltersArgumentCollection {
private static final int DEFAULT_MIN_READS_ON_EACH_STRAND = 0;
private static final double DEFAULT_MIN_AF = 0;
- @Argument(fullName = MAX_EVENTS_IN_REGION_LONG_NAME, optional = true, doc = "Maximum events in a single assembly region. Filter all variants if exceeded.")
+ @Argument(fullName = MAX_EVENTS_IN_REGION_LONG_NAME, optional = true, doc = "Maximum number of non-germline events in a single assembly region. Filter all variants if exceeded.")
public int maxEventsInRegion = DEFAULT_MAX_EVENTS_IN_REGION;
+ @Argument(fullName = MAX_EVENTS_IN_HAPLOTYPE_LONG_NAME, optional = true, doc = "Maximum number of non-germline events in a variant allele's best haplotype.")
+ public int maxEventsInHaplotype = DEFAULT_MAX_EVENTS_IN_HAPLOTYPE;
+
@Argument(fullName = MAX_ALT_ALLELE_COUNT_LONG_NAME, optional = true, doc = "Maximum alt alleles per site.")
public int numAltAllelesThreshold = DEFAULT_MAX_ALT_ALLELES;
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java
index 159f67df68d..77d2055587c 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java
@@ -312,7 +312,7 @@ private void buildFiltersList(final M2FiltersArgumentCollection MTFAC) {
}
if (!MTFAC.mitochondria && !MTFAC.microbial) {
- filters.add(new ClusteredEventsFilter(MTFAC.maxEventsInRegion));
+ filters.add(new ClusteredEventsFilter(MTFAC.maxEventsInRegion, MTFAC.maxEventsInHaplotype));
filters.add(new MultiallelicFilter(MTFAC.numAltAllelesThreshold));
filters.add(new FragmentLengthFilter(MTFAC.maxMedianFragmentLengthDifference));
filters.add(new PolymeraseSlippageFilter(MTFAC.minSlippageLength, MTFAC.slippageRate));
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/JointGermlineCNVSegmentation.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/JointGermlineCNVSegmentation.java
index 81138035bb1..e447d2c1086 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/JointGermlineCNVSegmentation.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/JointGermlineCNVSegmentation.java
@@ -646,11 +646,11 @@ public void closeTool(){
* @param minQuality drop events with quality lower than this
* @return a new record or null
*/
- public static SVCallRecord createDepthOnlyFromGCNVWithOriginalGenotypes(final VariantContext variant,
- final double minQuality,
- final Set allosomalContigs,
- final int refAutosomalCopyNumber,
- final SampleDB sampleDB) {
+ public SVCallRecord createDepthOnlyFromGCNVWithOriginalGenotypes(final VariantContext variant,
+ final double minQuality,
+ final Set allosomalContigs,
+ final int refAutosomalCopyNumber,
+ final SampleDB sampleDB) {
Utils.nonNull(variant);
if (variant.getGenotypes().size() == 1) {
//only cluster good variants
@@ -672,7 +672,7 @@ public static SVCallRecord createDepthOnlyFromGCNVWithOriginalGenotypes(final Va
.collect(Collectors.toList());
svBuilder.genotypes(genotypesWithECN);
- final SVCallRecord baseRecord = SVCallRecordUtils.create(svBuilder.make(), true);
+ final SVCallRecord baseRecord = SVCallRecordUtils.create(svBuilder.make(), true, dictionary);
final List nonRefGenotypes = baseRecord.getGenotypes().stream()
.filter(g -> !(g.isHomRef() || (g.isNoCall() && !g.hasExtendedAttribute(GATKSVVCFConstants.COPY_NUMBER_FORMAT))))
.collect(Collectors.toList());
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVAnnotateEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVAnnotateEngine.java
index 0e840e4652d..d49a351c39f 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVAnnotateEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVAnnotateEngine.java
@@ -9,10 +9,10 @@
import org.broadinstitute.hellbender.exceptions.UserException;
import org.broadinstitute.hellbender.tools.spark.sv.utils.GATKSVVCFConstants;
import org.broadinstitute.hellbender.tools.spark.sv.utils.SVUtils;
+import org.broadinstitute.hellbender.tools.sv.SVCallRecordUtils;
import org.broadinstitute.hellbender.utils.SVInterval;
import org.broadinstitute.hellbender.utils.SVIntervalTree;
import org.broadinstitute.hellbender.utils.SimpleInterval;
-import org.broadinstitute.hellbender.utils.Utils;
import org.broadinstitute.hellbender.utils.codecs.gtf.GencodeGtfFeature;
import org.broadinstitute.hellbender.utils.codecs.gtf.GencodeGtfTranscriptFeature;
import org.broadinstitute.hellbender.utils.variant.GATKSVVariantContextUtils;
@@ -57,36 +57,6 @@ public class SVAnnotateEngine {
GATKSVVCFConstants.ComplexVariantSubtype.delINVdup,
GATKSVVCFConstants.ComplexVariantSubtype.dDUP_iDEL);
- // Mini class to package SV type and interval into one object
- @VisibleForTesting
- protected static final class SVSegment {
- private final GATKSVVCFConstants.StructuralVariantAnnotationType intervalSVType;
- private final SimpleInterval interval;
- protected SVSegment(final GATKSVVCFConstants.StructuralVariantAnnotationType svType, final SimpleInterval interval) {
- this.intervalSVType = svType;
- this.interval = interval;
- }
- public GATKSVVCFConstants.StructuralVariantAnnotationType getIntervalSVType() {
- return intervalSVType;
- }
- public SimpleInterval getInterval() {
- return interval;
- }
-
- @Override
- public boolean equals(Object o) {
- if (this == o) return true;
- if (o == null || getClass() != o.getClass()) return false;
- final SVSegment svSegment = (SVSegment) o;
- return intervalSVType == svSegment.intervalSVType && interval.equals(svSegment.interval);
- }
-
- @Override
- public int hashCode() {
- return Objects.hash(intervalSVType, interval);
- }
- }
-
// Container class for all SVIntervalTree trees created from the GTF
@VisibleForTesting
public static final class GTFIntervalTreesContainer {
@@ -862,12 +832,7 @@ protected static boolean includesDispersedDuplication(final GATKSVVCFConstants.C
protected Map annotateStructuralVariant(final VariantContext variant) {
final Map> variantConsequenceDict = new HashMap<>();
final GATKSVVCFConstants.StructuralVariantAnnotationType overallSVType = getSVType(variant);
- final String complexTypeString = variant.getAttributeAsString(GATKSVVCFConstants.CPX_TYPE, null);
- GATKSVVCFConstants.ComplexVariantSubtype complexType = null;
- if (complexTypeString != null) {
- // replace / in CTX_PP/QQ and CTX_PQ/QP with _ to match ComplexVariantSubtype constants which cannot contain slashes
- complexType = GATKSVVCFConstants.ComplexVariantSubtype.valueOf(complexTypeString.replace("/", "_"));
- }
+ final GATKSVVCFConstants.ComplexVariantSubtype complexType = SVCallRecordUtils.getComplexSubtype(variant);
final boolean includesDispersedDuplication = includesDispersedDuplication(complexType, COMPLEX_SUBTYPES_WITH_DISPERSED_DUP);
final List svSegmentsForGeneOverlaps = getSVSegments(variant, overallSVType, maxBreakendLen, complexType);
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVCluster.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVCluster.java
index 4e483652d3c..791befb79fb 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVCluster.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVCluster.java
@@ -378,7 +378,7 @@ public void closeTool() {
@Override
public void apply(final VariantContext variant, final ReadsContext readsContext,
final ReferenceContext referenceContext, final FeatureContext featureContext) {
- final SVCallRecord call = SVCallRecordUtils.create(variant);
+ final SVCallRecord call = SVCallRecordUtils.create(variant, dictionary);
final SVCallRecord filteredCall;
if (fastMode && call.getType() != GATKSVVCFConstants.StructuralVariantAnnotationType.CNV) {
// Strip out non-carrier genotypes to save memory and compute
@@ -447,9 +447,9 @@ public VariantContext buildVariantContext(final SVCallRecord call) {
// Build new variant
final SVCallRecord finalCall = new SVCallRecord(newId, call.getContigA(), call.getPositionA(), call.getStrandA(),
- call.getContigB(), call.getPositionB(), call.getStrandB(), call.getType(), call.getComplexSubtype(), call.getLength(),
- call.getAlgorithms(), call.getAlleles(), filledGenotypes, call.getAttributes(), call.getFilters(),
- call.getLog10PError(), dictionary);
+ call.getContigB(), call.getPositionB(), call.getStrandB(), call.getType(), call.getComplexSubtype(),
+ call.getComplexEventIntervals(), call.getLength(), call.getAlgorithms(), call.getAlleles(), filledGenotypes,
+ call.getAttributes(), call.getFilters(), call.getLog10PError(), dictionary);
final VariantContextBuilder builder = SVCallRecordUtils.getVariantBuilder(finalCall);
if (omitMembers) {
builder.rmAttribute(GATKSVVCFConstants.CLUSTER_MEMBER_IDS_KEY);
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVConcordance.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVConcordance.java
index d53c3a22b0b..56c389539c6 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVConcordance.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVConcordance.java
@@ -180,7 +180,7 @@ public void apply(final TruthVersusEval truthVersusEval, final ReadsContext read
}
private void add(final VariantContext variant, final boolean isTruth) {
- SVCallRecord record = SVCallRecordUtils.create(variant);
+ SVCallRecord record = SVCallRecordUtils.create(variant, dictionary);
if (!record.getContigA().equals(currentContig)) {
flushClusters(true);
currentContig = record.getContigA();
@@ -199,8 +199,8 @@ protected SVCallRecord minimizeTruthFootprint(final SVCallRecord item) {
final List genotypes = item.getGenotypes().stream().map(SVConcordance::stripTruthGenotype).collect(Collectors.toList());
return new SVCallRecord(item.getId(), item.getContigA(), item.getPositionA(),
item.getStrandA(), item.getContigB(), item.getPositionB(), item.getStrandB(), item.getType(),
- item.getCpxSubtype(), item.getLength(), item.getAlgorithms(), item.getAlleles(), genotypes,
- item.getAttributes(), item.getFilters(), item.getLog10PError(), dictionary);
+ item.getComplexSubtype(), item.getComplexEventIntervals(), item.getLength(), item.getAlgorithms(),
+ item.getAlleles(), genotypes, item.getAttributes(), item.getFilters(), item.getLog10PError(), dictionary);
}
private static Genotype stripTruthGenotype(final Genotype genotype) {
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVSegment.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVSegment.java
new file mode 100644
index 00000000000..29f4d5464d0
--- /dev/null
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/sv/SVSegment.java
@@ -0,0 +1,56 @@
+package org.broadinstitute.hellbender.tools.walkers.sv;
+
+import htsjdk.samtools.util.Locatable;
+import org.broadinstitute.hellbender.tools.spark.sv.utils.GATKSVVCFConstants;
+import org.broadinstitute.hellbender.utils.SimpleInterval;
+import org.broadinstitute.hellbender.utils.Utils;
+
+import java.util.Objects;
+
+// Mini class to package SV type and interval into one object
+public class SVSegment implements Locatable {
+ protected final GATKSVVCFConstants.StructuralVariantAnnotationType intervalSVType;
+ protected final SimpleInterval interval;
+
+ public SVSegment(final GATKSVVCFConstants.StructuralVariantAnnotationType svType, final SimpleInterval interval) {
+ Utils.nonNull(interval);
+ this.intervalSVType = svType;
+ this.interval = interval;
+ }
+
+ public GATKSVVCFConstants.StructuralVariantAnnotationType getIntervalSVType() {
+ return intervalSVType;
+ }
+
+ @Override
+ public String getContig() {
+ return interval.getContig();
+ }
+
+ @Override
+ public int getStart() {
+ return interval.getStart();
+ }
+
+ @Override
+ public int getEnd() {
+ return interval.getEnd();
+ }
+
+ public SimpleInterval getInterval() {
+ return interval;
+ }
+
+ @Override
+ public boolean equals(Object o) {
+ if (this == o) return true;
+ if (o == null || getClass() != o.getClass()) return false;
+ final SVSegment svSegment = (SVSegment) o;
+ return intervalSVType == svSegment.intervalSVType && interval.equals(svSegment.interval);
+ }
+
+ @Override
+ public int hashCode() {
+ return Objects.hash(intervalSVType, interval);
+ }
+}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantrecalling/HaplotypeBasedVariantRecaller.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantrecalling/HaplotypeBasedVariantRecaller.java
index 6e02a782a15..33108ec5968 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantrecalling/HaplotypeBasedVariantRecaller.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantrecalling/HaplotypeBasedVariantRecaller.java
@@ -135,11 +135,11 @@ public void traverse() {
}
// get reads overlapping haplotypes
- final Map> readsByReader = readsReader.getReads(haplotypeSpan, vcLoc);
+ final Map> readsByReader = readsReader.getReads(haplotypeSpan, vcLoc);
final List variants = new LinkedList<>(Arrays.asList(vc));
if ( logger.isDebugEnabled() ) {
int readCount = 0;
- for ( Collection reads : readsByReader.values() )
+ for ( Collection reads : readsByReader.values() )
readCount += reads.size();
logger.debug(String.format("vcLoc %s, haplotypeSpan: %s, %d haplotypes, %d reads",
vcLoc.toString(), haplotypeSpan.toString(), processedHaplotypes.size(), readCount, variants.size()));
@@ -150,16 +150,16 @@ public void traverse() {
final List>> genotypeLikelihoodsList = new LinkedList<>();
final List assemblyResultList = new LinkedList<>();
final List readsHeaderList = new LinkedList<>();
- for ( Map.Entry> entry : readsByReader.entrySet() ) {
+ for ( Map.Entry> entry : readsByReader.entrySet() ) {
final AssemblyResultSet assemblyResult = new AssemblyResultSet();
processedHaplotypes.forEach(haplotype -> assemblyResult.add(haplotype));
final Map> perSampleReadList = new LinkedHashMap<>();
final SamReader samReader = entry.getKey();
- final Collection reads = entry.getValue();
+ final Collection reads = entry.getValue();
List gtakReads = new LinkedList<>();
- reads.forEach(flowBasedRead -> gtakReads.add(flowBasedRead));
+ reads.forEach(read -> gtakReads.add(read));
perSampleReadList.put(sampleNames[0], gtakReads);
AssemblyRegion regionForGenotyping = new AssemblyRegion(haplotypeSpan, 0, samReader.getFileHeader());
assemblyResult.setPaddedReferenceLoc(haplotypeSpan);
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantrecalling/TrimmedReadsReader.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantrecalling/TrimmedReadsReader.java
index 23f47eec203..94600f97ea1 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantrecalling/TrimmedReadsReader.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantrecalling/TrimmedReadsReader.java
@@ -26,8 +26,6 @@ public class TrimmedReadsReader {
private final List samReaders = new LinkedList<>();
private CountingReadFilter readFilter;
- private final Map readGroupMaxClass = new LinkedHashMap<>();
- private final Map readGroupFlowOrder = new LinkedHashMap<>();
private final FlowBasedArgumentCollection fbArgs = new FlowBasedArgumentCollection();
public TrimmedReadsReader(final List readsFiles, final Path referencePath, final int cloudPrefetchBuffer) {
@@ -44,11 +42,12 @@ public SAMSequenceDictionary getSamSequenceDictionary(final SamReader samReader)
return ((samReader != null) ? samReader : samReaders.get(0)).getFileHeader().getSequenceDictionary();
}
- public Map> getReads(final Locatable span, final Locatable vcLoc) {
- final Map> readsByReader = new LinkedHashMap<>();
+ public Map> getReads(final Locatable span, final Locatable vcLoc) {
+
+ final Map> readsByReader = new LinkedHashMap<>();
for ( SamReader samReader : samReaders ) {
- final List reads = new LinkedList<>();
+ final List reads = new LinkedList<>();
final SAMRecordIterator iter = samReader.query(span.getContig(), span.getStart(), span.getEnd(), false);
while (iter.hasNext()) {
@@ -72,7 +71,10 @@ public Map> getReads(final Locatable span,
gatkRead = ReadClipper.hardClipToRegion(gatkRead, span.getStart(), span.getEnd());
if (gatkRead.isUnmapped() || gatkRead.getCigar().isEmpty())
continue;
-
+ if (!FlowBasedReadUtils.isFlowPlatform(samReader.getFileHeader(), gatkRead)){
+ reads.add(gatkRead);
+ continue;
+ }
// convert to a flow based read
FlowBasedReadUtils.ReadGroupInfo rgInfo = FlowBasedReadUtils.getReadGroupInfo(samReader.getFileHeader(), gatkRead);
final FlowBasedRead fbr = new FlowBasedRead(gatkRead, rgInfo.flowOrder, rgInfo.maxClass, fbArgs);
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/FamilyLikelihoods.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/FamilyLikelihoods.java
index 012b8396799..8f261162d6b 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/FamilyLikelihoods.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/FamilyLikelihoods.java
@@ -123,7 +123,7 @@ private Genotype getUpdatedGenotype(final VariantContext vc, final Genotype geno
final GenotypeBuilder builder = new GenotypeBuilder(genotype);
//update genotype types based on posteriors
- GATKVariantContextUtils.makeGenotypeCall(vc.getMaxPloidy(2), builder, GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN, log10Posteriors, vc.getAlleles(), null);
+ GATKVariantContextUtils.makeGenotypeCall(vc.getMaxPloidy(2), builder, GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN, log10Posteriors, vc.getAlleles(), genotype, null);
builder.attribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY,
Utils.listFromPrimitives(GenotypeLikelihoods.fromLog10Likelihoods(log10Posteriors).getAsPLs()));
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java
index 3b1ef696d53..b919530f23a 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtils.java
@@ -7,6 +7,7 @@
import org.broadinstitute.hellbender.exceptions.UserException;
import org.broadinstitute.hellbender.tools.walkers.ReferenceConfidenceVariantContextMerger;
import org.broadinstitute.hellbender.tools.walkers.genotyper.GenotypeAssignmentMethod;
+import org.broadinstitute.hellbender.tools.walkers.genotyper.GenotypesCache;
import org.broadinstitute.hellbender.utils.MathUtils;
import org.broadinstitute.hellbender.utils.SimpleInterval;
import org.broadinstitute.hellbender.utils.Utils;
@@ -134,9 +135,16 @@ else if (a.isSymbolic()) { //use the greater of the priors for non-ref
builder.phased(vc1.getGenotype(genoIdx).isPhased());
if ( posteriors.get(genoIdx) != null ) {
GATKVariantContextUtils.makeGenotypeCall(vc1.getMaxPloidy(HomoSapiensConstants.DEFAULT_PLOIDY), builder,
- GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN, posteriors.get(genoIdx), vc1.getAlleles(), null);
+ GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN, posteriors.get(genoIdx), vc1.getAlleles(), vc1.getGenotype(genoIdx), null);
+
builder.attribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY,
- Utils.listFromPrimitives(GenotypeLikelihoods.fromLog10Likelihoods(posteriors.get(genoIdx)).getAsPLs()));
+ Utils.listFromPrimitives(GenotypeLikelihoods.fromLog10Likelihoods(posteriors.get(genoIdx)).getAsPLs()));
+ //TODO: refactor me -- this is to make up for the GGVCFs changes fixing GQ0 hom-refs
+ if (!builder.makeWithShallowCopy().hasPL()) {
+ final int maxPosteriorIndex = MathUtils.maxElementIndex(posteriors.get(genoIdx));
+ builder.log10PError(GenotypeLikelihoods.getGQLog10FromLikelihoods(maxPosteriorIndex, posteriors.get(genoIdx)));
+ builder.alleles(GenotypesCache.get(vc1.getGenotype(genoIdx).getPloidy(), maxPosteriorIndex).asAlleleList(vc1.getAlleles()));
+ }
}
newContext.add(builder.make());
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java
index 8dd4c6a7eea..073b43fa0aa 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java
@@ -159,7 +159,7 @@ public final class ReblockGVCF extends MultiVariantWalker {
@Advanced
@Argument(fullName=ANNOTATIONS_TO_REMOVE_LONG_NAME, doc="FORMAT level annotations to remove from all genotypes in final GVCF.", optional = true)
- private List annotationsToRemove = new ArrayList<>();
+ private List formatAnnotationsToRemove = new ArrayList<>();
@Advanced
@Argument(fullName=QUAL_APPROX_LONG_NAME, shortName=QUAL_APPROX_SHORT_NAME, doc="Add necessary INFO field annotation to perform QUAL approximation downstream; required for GnarlyGenotyper", optional = true)
@@ -238,9 +238,6 @@ public void onTraversalStart() {
+ ", but the " + GATKVCFConstants.TREE_SCORE + " annotation is not present in the input GVCF.");
}
- List missingAnnotationsToRemove = annotationsToRemove.stream().filter(a -> inputHeader.getFormatHeaderLine(a)==null).toList();
- missingAnnotationsToRemove.forEach(a -> logger.warn("FORMAT level annotation " + a + ", which was requested to be removed by --" + ANNOTATIONS_TO_REMOVE_LONG_NAME + ", not found in input GVCF header."));
-
final Set inputHeaders = inputHeader.getMetaDataInSortedOrder();
final Set headerLines = new HashSet<>(inputHeaders);
@@ -249,6 +246,11 @@ public void onTraversalStart() {
(vcfHeaderLine.getKey().equals("INFO")) && ((VCFInfoHeaderLine)vcfHeaderLine).getID().equals(GATKVCFConstants.RAW_RMS_MAPPING_QUALITY_DEPRECATED) || //remove old (maybe wrong type) and add new with deprecated note
(vcfHeaderLine.getKey().equals("INFO")) && infoFieldAnnotationKeyNamesToRemove.contains(((VCFInfoHeaderLine)vcfHeaderLine).getID()));
+ List missingAnnotationsToRemove = formatAnnotationsToRemove.stream().filter(a -> inputHeader.getFormatHeaderLine(a)==null).toList();
+ missingAnnotationsToRemove.forEach(a -> logger.warn("FORMAT level annotation " + a + ", which was requested to be removed by --" + ANNOTATIONS_TO_REMOVE_LONG_NAME + ", not found in input GVCF header."));
+ headerLines.removeIf(vcfHeaderLine -> vcfHeaderLine.getKey().equals("FORMAT") &&
+ formatAnnotationsToRemove.contains(((VCFFormatHeaderLine) vcfHeaderLine).getID()));
+
headerLines.addAll(getDefaultToolVCFHeaderLines());
genotypingEngine = createGenotypingEngine(new IndexedSampleList(getSamplesForVariants()));
@@ -338,7 +340,7 @@ public void apply(VariantContext variant, ReadsContext reads, ReferenceContext r
if (!variant.hasAllele(Allele.NON_REF_ALLELE)) {
throw new UserException("Variant Context at " + variant.getContig() + ":" + variant.getStart() + " does not contain a allele. This tool is only intended for use with GVCFs.");
}
- VariantContext newVC = annotationsToRemove.size() > 0 ? removeVCFFormatAnnotations(variant) : variant;
+ VariantContext newVC = formatAnnotationsToRemove.size() > 0 ? removeVCFFormatAnnotations(variant) : variant;
regenotypeVC(newVC);
}
@@ -351,7 +353,12 @@ public void apply(VariantContext variant, ReadsContext reads, ReferenceContext r
private VariantContext removeVCFFormatAnnotations(final VariantContext vc) {
final Genotype genotype = vc.getGenotype(0);
Map extendedAttributes = genotype.getExtendedAttributes();
- for (String annotation : annotationsToRemove) {
+ // If extendedAttributes is empty, we get an unmodifiable empty map, so skip removing any annotations
+ if (extendedAttributes.isEmpty()) {
+ return vc;
+ }
+
+ for (String annotation : formatAnnotationsToRemove) {
extendedAttributes.remove(annotation);
}
final Genotype newGenotype = new GenotypeBuilder(genotype).noAttributes().attributes(extendedAttributes).make();
@@ -480,7 +487,7 @@ private boolean isMonomorphicCallWithAlts(final VariantContext result) {
final Genotype genotype = result.getGenotype(0);
return (hasGenotypeValuesArray(posteriorsKey, genotype)
&& (genotype.isHomRef() || isNoCalledHomRef(posteriorsKey, genotype) || (genotype.hasPL() && MathUtils.minElementIndex(genotype.getPL()) == 0))
- && result.getAlternateAlleles().stream().anyMatch(this::isConcreteAlt));
+ && result.getAlternateAlleles().stream().anyMatch(GATKVariantContextUtils::isConcreteAlt));
}
/**
@@ -519,26 +526,13 @@ boolean shouldBeReblocked(final VariantContext vc) {
final List finalAlleles = alleleCounts.asAlleleList(vc.getAlleles());
return (pls != null && pls[0] < rgqThreshold)
- || !genotypeHasConcreteAlt(finalAlleles)
+ || !GATKVariantContextUtils.genotypeHasConcreteAlt(finalAlleles)
|| finalAlleles.stream().anyMatch(a -> a.equals(Allele.NON_REF_ALLELE))
|| (!genotype.hasPL() && !genotype.hasGQ()
// If TREE_SCORE threshold is set to be > 0 then sites without TREE_SCORE should be reblocked.
|| vc.getAttributeAsDouble(GATKVCFConstants.TREE_SCORE, 0) < treeScoreThreshold);
}
- /**
- * Is there a real ALT allele that's not or * ?
- * @param alleles alleles in the called genotype
- * @return true if the genotype has a called allele that is a "real" alternate
- */
- private boolean genotypeHasConcreteAlt(final List alleles) {
- return alleles.stream().anyMatch(this::isConcreteAlt);
- }
-
- private boolean isConcreteAlt(final Allele a) {
- return !a.isReference() && !a.isSymbolic() && !a.equals(Allele.SPAN_DEL);
- }
-
/**
* "reblock" a variant by converting its genotype to homRef, changing PLs, adding reblock END tags and other attributes
* @param lowQualityVariant a variant already determined to be low quality
@@ -819,7 +813,7 @@ private Genotype getCalledGenotype(final VariantContext variant) {
final GenotypeBuilder builderToCallAlleles = new GenotypeBuilder(noCallGT);
//TODO: update to support DRAGEN posteriors
GATKVariantContextUtils.makeGenotypeCall(noCallGT.getPloidy(), builderToCallAlleles, GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN,
- noCallGT.getLikelihoods().getAsVector(), variant.getAlleles(), null);
+ noCallGT.getLikelihoods().getAsVector(), variant.getAlleles(), noCallGT, null);
return builderToCallAlleles.make();
} else {
return variant.getGenotype(0);
@@ -936,7 +930,10 @@ private static void copyInfoAnnotations(final Map destinationAtt
final boolean allelesNeedSubsetting, final int[] relevantIndices) {
//copy over info annotations
final Map origMap = sourceVC.getAttributes();
- for(final InfoFieldAnnotation annotation : annotationEngine.getInfoAnnotations()) {
+ final List engineAnnotations = new ArrayList<>();
+ engineAnnotations.addAll(annotationEngine.getInfoAnnotations());
+ engineAnnotations.addAll(annotationEngine.getJumboInfoAnnotations());
+ for(final VariantAnnotation annotation : engineAnnotations) {
for (final String key : annotation.getKeyNames()) {
if (infoFieldAnnotationKeyNamesToRemove.contains(key)) {
continue;
@@ -945,17 +942,34 @@ private static void copyInfoAnnotations(final Map destinationAtt
destinationAttrMap.put(key, origMap.get(key));
}
}
- if (annotation instanceof ReducibleAnnotation) {
- for (final String rawKey : ((ReducibleAnnotation)annotation).getRawKeyNames()) {
+ if (annotation instanceof AlleleSpecificAnnotation) {
+ final boolean isReducible = annotation instanceof ReducibleAnnotation;
+ final List keyNames = isReducible ? ((ReducibleAnnotation)annotation).getRawKeyNames() :
+ annotation.getKeyNames();
+ for (final String rawKey : keyNames) {
if (infoFieldAnnotationKeyNamesToRemove.contains(rawKey)) {
continue;
}
if (origMap.containsKey(rawKey)) {
- if (allelesNeedSubsetting && AnnotationUtils.isAlleleSpecific(annotation)) {
- final List alleleSpecificValues = AnnotationUtils.getAlleleLengthListOfString(sourceVC.getAttributeAsString(rawKey, null));
+ if (allelesNeedSubsetting && AnnotationUtils.isAlleleSpecific(annotation) && AnnotationUtils.isAlleleSpecificGatkKey(rawKey)) {
+ final List alleleSpecificValues;
+ if (isReducible) {
+ alleleSpecificValues = AnnotationUtils.getAlleleLengthListOfStringFromRawData(sourceVC.getAttributeAsString(rawKey, null));
+ } else {
+ String value = sourceVC.getAttributeAsString(rawKey, null);
+ if (value == null) {
+ alleleSpecificValues = Collections.emptyList();
+ } else {
+ alleleSpecificValues = AnnotationUtils.getAlleleLengthListOfString(value);
+ }
+ }
final List subsetList;
if (alleleSpecificValues.size() > 0) {
- subsetList = AlleleSubsettingUtils.remapRLengthList(alleleSpecificValues, relevantIndices, "");
+ if (isReducible) {
+ subsetList = AlleleSubsettingUtils.remapRLengthList(alleleSpecificValues, relevantIndices, "");
+ } else {
+ subsetList = AlleleSubsettingUtils.remapALengthList(alleleSpecificValues, relevantIndices, "");
+ }
if (sourceVC.getAlleles().get(relevantIndices[relevantIndices.length - 1]).equals(Allele.NON_REF_ALLELE)) {
//zero out non-ref value, just in case
subsetList.set(subsetList.size() - 1, ((AlleleSpecificAnnotation) annotation).getEmptyRawValue());
@@ -964,7 +978,11 @@ private static void copyInfoAnnotations(final Map destinationAtt
subsetList = Collections.nCopies(relevantIndices.length, "");
}
- destinationAttrMap.put(rawKey, AnnotationUtils.encodeAnyASListWithRawDelim(subsetList));
+ if (isReducible) {
+ destinationAttrMap.put(rawKey, AnnotationUtils.encodeAnyASListWithRawDelim(subsetList));
+ } else {
+ destinationAttrMap.put(rawKey, AnnotationUtils.encodeStringList(subsetList));
+ }
} else {
destinationAttrMap.put(rawKey, origMap.get(rawKey));
}
diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/VCFComparator.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/VCFComparator.java
new file mode 100644
index 00000000000..01395e54f21
--- /dev/null
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/VCFComparator.java
@@ -0,0 +1,981 @@
+package org.broadinstitute.hellbender.tools.walkers.variantutils;
+
+import htsjdk.variant.variantcontext.*;
+import htsjdk.variant.vcf.VCFConstants;
+import org.apache.commons.collections4.CollectionUtils;
+import org.broadinstitute.barclay.argparser.*;
+import org.broadinstitute.hellbender.engine.*;
+import org.broadinstitute.hellbender.exceptions.GATKException;
+import org.broadinstitute.hellbender.exceptions.UserException;
+import org.broadinstitute.hellbender.tools.walkers.annotator.*;
+import org.broadinstitute.hellbender.tools.walkers.annotator.allelespecific.AS_QualByDepth;
+import org.broadinstitute.hellbender.tools.walkers.annotator.allelespecific.AS_StandardAnnotation;
+import org.broadinstitute.hellbender.tools.walkers.genotyper.AlleleSubsettingUtils;
+import org.broadinstitute.hellbender.tools.walkers.genotyper.GenotypeAssignmentMethod;
+import org.broadinstitute.hellbender.utils.MathUtils;
+import org.broadinstitute.hellbender.utils.variant.GATKVCFConstants;
+import org.broadinstitute.hellbender.utils.variant.GATKVariantContextUtils;
+import picard.cmdline.programgroups.VariantEvaluationProgramGroup;
+
+import java.io.PrintStream;
+import java.util.*;
+import java.util.stream.Collectors;
+
+/**
+ * A tool to allow for comparison of two single-sample GVCFs or multi-sample VCFs when some differences are expected
+ * input variants should be tagged "actual" or "expected"
+ *
+ * Example usage:
+ * java -jar gatk.jar VCFComparator -R ref.fasta -V:actual actual.vcf -V:expected expected.vcf
+ *
+ * Caveats:
+ * Reference blocks are not checked
+ * VCFComparator will ignore version-specific NaN versus empty AS_RAW annotation discrepancies
+ * VCFs are both expected to NOT have split multi-alleleics
+ */
+@ExperimentalFeature
+@CommandLineProgramProperties(
+ summary = "Compare two VCFs, as for pipeline test updates",
+ oneLineSummary = "Compare two VCFs",
+ programGroup = VariantEvaluationProgramGroup.class
+)
+public class VCFComparator extends MultiVariantWalkerGroupedByOverlap {
+ public static final String ALLOW_NEW_STARS_LONG_NAME = "allow-new-stars";
+ public static final String IGNORE_QUALS_LONG_NAME = "ignore-quals";
+ public static final String DP_CHANGE_ALLOWED_LONG_NAME = "dp-change-allowed";
+ public static final String RANK_SUM_CHANGE_ALLOWED_LONG_NAME = "ranksum-change-allowed";
+ public static final String POSITIONS_ONLY_LONG_NAME = "positions-only";
+
+ private VariantAnnotatorEngine annotatorEngine;
+ private boolean isSingleSample = true;
+ private boolean alleleNumberIsDifferent = false;
+ private boolean inbreedingCoeffIsDifferent = false;
+ private final List annotationsThatVaryWithNoCalls = Arrays.asList(VCFConstants.ALLELE_NUMBER_KEY, GATKVCFConstants.INBREEDING_COEFFICIENT_KEY,
+ GATKVCFConstants.AS_INBREEDING_COEFFICIENT_KEY, GATKVCFConstants.EXCESS_HET_KEY,
+ GATKVCFConstants.MLE_ALLELE_COUNT_KEY, GATKVCFConstants.MLE_ALLELE_FREQUENCY_KEY, VCFConstants.ALLELE_FREQUENCY_KEY,
+ VCFConstants.DEPTH_KEY, GATKVCFConstants.QUAL_BY_DEPTH_KEY, GATKVCFConstants.AS_QUAL_BY_DEPTH_KEY,
+ GATKVCFConstants.VQS_LOD_KEY, GATKVCFConstants.AS_VQS_LOD_KEY, GATKVCFConstants.AS_FILTER_STATUS_KEY,
+ GATKVCFConstants.CULPRIT_KEY, GATKVCFConstants.AS_CULPRIT_KEY); //AC should be the same
+ private boolean failOnCompletion = false;
+
+ @Argument(fullName = "warn-on-errors",
+ shortName = "warn-on-errors",
+ doc = "just emit warnings on errors instead of terminating the run at the first instance. Won't fail unless --finish-before-failing is set",
+ optional = true)
+ private boolean warnOnError = false;
+
+ @Argument(fullName = "finish-before-failing",
+ doc = "emit warnings as necesarry, but perform the whole comparison before failing",
+ optional= true)
+ private boolean finishBeforeFailing = false;
+
+ @Argument(fullName = "default-ploidy", optional = true, doc = "Overall expected ploidy. Default value is 2.")
+ private int defaultPloidy = 2;
+
+ @Argument(fullName = IGNORE_QUALS_LONG_NAME, optional = true, mutex={POSITIONS_ONLY_LONG_NAME})
+ private boolean ignoreQuals = false;
+
+ @Argument(fullName = "qual-change-allowed", optional = true)
+ private double qualTolerance = 0.001;
+
+ @Argument(fullName = "inbreeding-coeff-change-allowed", optional = true)
+ private double inbreedingCoeffTolerance = 0.001;
+
+ @Argument(fullName = "good-qual-threshold", optional = true, doc = "Variants with QUAL at or above this value are considered high quality and should appear in both files")
+ private double goodQualThreshold = 100.0;
+
+ @Argument(fullName = DP_CHANGE_ALLOWED_LONG_NAME, optional = true, doc = "Note that this is a signed change (actual - expected)", mutex={POSITIONS_ONLY_LONG_NAME})
+ private int dpChange = 0;
+
+ @Argument(fullName = RANK_SUM_CHANGE_ALLOWED_LONG_NAME, optional = true, doc = "Amount of numerical 'noise' allowed in rank sum annotations")
+ private double ranksumTolerance = 0.0;
+
+ @Argument(fullName = "likelihood-change-allowed", optional = true)
+ private int likelihoodTolerance = 0;
+
+ @Argument(fullName = "ignore-non-ref-data", optional = true)
+ private boolean ignoreNonRefData = false;
+
+ @Argument(fullName = "ignore-annotations", optional = true, doc = "Only match on position and alleles, ignoring all INFO annotations", mutex={POSITIONS_ONLY_LONG_NAME})
+ private boolean ignoreAnnotations = false;
+
+ @Argument(fullName = "ignore-genotype-annotations", optional = true, doc = "Only match on genotype call, ignoring all FORMAT annotations", mutex={POSITIONS_ONLY_LONG_NAME})
+ private boolean ignoreGenotypeAnnotations = false;
+
+ @Argument(fullName = "ignore-genotype-phasing", optional = true, doc = "Only check called genotype alleles, not order or phasing status")
+ private boolean ignoreGenotypePhasing = false;
+
+ @Argument(fullName = "ignore-filters", optional = true, doc = "Ignore filter status when comparing variants", mutex={POSITIONS_ONLY_LONG_NAME})
+ private boolean ignoreFilters = false;
+
+ @Argument(fullName = "ignore-attribute", optional = true, doc = "Ignore INFO attributes with this key", mutex={POSITIONS_ONLY_LONG_NAME})
+ private List ignoreAttributes = new ArrayList<>(5);
+
+ @Argument(fullName = POSITIONS_ONLY_LONG_NAME, optional = true, doc = "Only match on position, ignoring alleles and annotations", mutex={})
+ private boolean positionsOnly = false;
+
+ @Argument(fullName = ALLOW_NEW_STARS_LONG_NAME, optional = true, doc = "Allow additional * alleles in actual if there is a corresponding deletion", mutex={POSITIONS_ONLY_LONG_NAME})
+ private boolean allowNewStars = false;
+
+ @Argument(fullName = "allow-extra-alleles", optional = true, doc = "Allow extra alleles in actual provided actual is a superset of expected", mutex={ALLOW_NEW_STARS_LONG_NAME, POSITIONS_ONLY_LONG_NAME})
+ private boolean allowExtraAlleles = false;
+
+ @Argument(fullName = "allow-missing-stars", optional = true, doc = "Allow missing * in actual if actual has no corresponding deletions")
+ private boolean allowMissingStars = false;
+
+ @Argument(fullName = "ignore-star-attributes", optional = true, doc = "Ignore annotation values for spanning deletion (*) alleles")
+ private boolean ignoreStarAttributes = false;
+
+ @Argument(fullName = "allow-nan-mismatch", optional = true, doc = "Allow VCF values of missing (.) and NaN to match")
+ private boolean allowNanMismatch = false;
+
+ @Argument(fullName = "mute-acceptable-diffs", optional = true, doc = "Suppress warnings or exceptions for differences that are consequences of low quality genotypes or AN discrepancies")
+ private boolean muteDiffs = false;
+
+ @Argument(fullName = "ignore-hom-ref-attributes", optional = true, doc = "Skip attribute comparison for homozygous reference genotypes")
+ private boolean ignoreHomRefAttributes = false;
+
+ @Argument(fullName = "ignore-dbsnp-ids", optional = true)
+ private boolean ignoreDbsnp = false;
+
+ @Argument(fullName = "ignore-gq0", optional = true)
+ private boolean ignoreGq0 = false;
+
+ @Argument(fullName = "ignore-some-multi-allelics", optional = true, doc="Ignore sites where the AC length in the actual matches the actual number of alleles, but doesn't match the expected VC.")
+ private boolean ignoreSomeMultiAllelics = false;
+
+ @Advanced
+ @Argument(fullName=ReblockGVCF.ANNOTATIONS_TO_KEEP_LONG_NAME, doc="Annotations that are not recognized by GATK that should be checked, otherwise they are removed if ignore-non-ref-data is set.", optional = true)
+ private List annotationsToKeep = new ArrayList<>();
+
+ @Hidden
+ @Argument(fullName = "output-warnings", optional = true, doc = "Output warnings to a file, useful for testing")
+ private GATKPath warningsOutput = null;
+
+ private PrintStream outputStream = null;
+
+ @Override
+ public boolean doDictionaryCrossValidation() {
+ return false;
+ } //speed things up a bit
+
+ @Override
+ public boolean useVariantAnnotations() { return true;}
+
+ @Override
+ public List> getDefaultVariantAnnotationGroups() {
+ return Arrays.asList(StandardAnnotation.class, AS_StandardAnnotation.class);
+ }
+
+ @Override
+ public void onTraversalStart() {
+ if (getDrivingVariantsFeatureInputs().size() != 2) {
+ throw new UserException.BadInput("VCFComparator expects exactly two inputs -- one actual and one expected.");
+ }
+ final List> expected = getDrivingVariantsFeatureInputs().stream().filter(file -> file.getName().equals("expected")).collect(Collectors.toList());
+ if (expected.size() != 1) {
+ throw new UserException.BadInput("Tool requires exactly one expected input file");
+ }
+
+ annotatorEngine = new VariantAnnotatorEngine(makeVariantAnnotations(), null, Collections.emptyList(), false, false);
+
+ isSingleSample = getSamplesForVariants().size() == 1 ? true : false; //if either expected or actual has more than one sample then we're not single-sample
+
+ if (warningsOutput != null) {
+ outputStream = new PrintStream(warningsOutput.getOutputStream());
+ }
+ }
+
+ @Override
+ public void apply(final List variantContexts, final ReferenceContext referenceContext, final List readsContexts) {
+ if (isSingleSample && !variantContexts.get(0).getAlleles().contains(Allele.NON_REF_ALLELE)) {
+ throw new UserException.BadInput("Single-sample mode expects two GVCFs with data for comparison");
+ }
+
+ if (variantContexts.size() == 1) {
+ final VariantContext vc = variantContexts.get(0);
+ //TODO: if there's a variant that overlaps the requested interval the start may not be in the requested intervals
+ if (!muteDiffs) {
+ if (!ignoreGq0 && vc.getGenotypes().stream().anyMatch(g -> g.hasGQ() && g.getGQ() == 0)) {
+ throwOrWarn(new UserException("Unmatched variant in " + vc.getSource() + " at position " + vc.getContig() + ":" + vc.getStart()));
+ } else {
+ return;
+ }
+ //low coverage sites will have QUAL for singleton hom-vars boosted because of GQ0 hom-refs
+ } else if (hasGoodEvidence(vc)) {
+ if (!ignoreGq0 && vc.getGenotypes().stream().anyMatch(g -> g.hasGQ() && g.getGQ() == 0)) {
+ throwOrWarn(new UserException("Unmatched variant in " + vc.getSource() + " at position " + vc.getContig() + ":" + vc.getStart()));
+ } else {
+ return;
+ }
+ }
+ }
+
+ //here we have all variants overlapping a given position
+ //if there's no actual for the expected (at the same position), see if the expected has an overlapping deletion
+ for (final VariantContext vc : variantContexts) {
+ if (vc.getSource().equals("actual")) { //there may be more expected than actual variants, but only compare once
+ continue; //TODO: this means we don't report cases where actual has multiple reference blocks spanning expected as unmatched
+ }
+ //vc is guaranteed to be from expected
+
+ //variantContexts can have a bunch of overlapping variants, but use the start position to compare. Could still have two alleles with same start.
+ final List matches = variantContexts.stream().filter(v -> v.getStart() == vc.getStart()).collect(Collectors.toList());
+ //matches includes vc (both actual and expected)
+ //TODO: matches.size() can be up to 2 times the number of haplotypes (based on ploidy). So checking only one match won't work for non-diploid organisms
+ if (matches.size() == 1 ) {
+ // if there's only one match, then it's the expected matching itself, so actual is missing the same vc start
+ // if it's a high quality site, throw an error/warning otherwise it's low quality so skip it
+ if (isHighQuality(vc) && hasGoodEvidence(vc)
+ && vc.getGenotypes().stream().noneMatch(g -> g.getAlleles().contains(Allele.SPAN_DEL))){
+ throwOrWarn(new UserException("Apparent unmatched high quality variant in " + vc.getSource() + " at " + vc.getContig() + ":" + vc.getStart()));
+ } else {
+ return;
+ }
+ } else {
+ if ((!isHighQuality(vc) || (isSingleSample && !GATKVariantContextUtils.genotypeHasConcreteAlt(vc.getGenotype(0).getAlleles())))
+ && ignoreHomRefAttributes) {
+ return;
+ }
+ final VariantContext match;
+ // VCFs are expected to NOT have split multi-alleleics
+ if (!matches.get(0).getSource().equals(vc.getSource())) {
+ match = matches.get(0);
+ } else {
+ match = matches.get(1);
+ }
+
+ final List overlappingDels = variantContexts.stream().filter(v -> v.getStart() < vc.getStart() && v.overlaps(vc))
+ .collect(Collectors.toList());
+
+ final VariantContext expectedTrimmed = trimAlleles(vc, overlappingDels);
+ final VariantContext actualTrimmed = trimAlleles(match, overlappingDels);
+
+ final boolean vcShouldBeSkippedDueToNoCalls = variantContexts.stream().anyMatch(v -> v.getGenotypes().stream().anyMatch(g -> g.isNoCall() || (g.hasGQ() && g.getGQ() == 0)));
+ //do single-sample GVCF checks, including deletion trimming and dropping
+ if (isSingleSample) {
+ try {
+ if (!(ignoreGq0 && vcShouldBeSkippedDueToNoCalls)) {
+ validateSingleSampleDeletions(vc, match, expectedTrimmed, actualTrimmed, overlappingDels, vcShouldBeSkippedDueToNoCalls);
+ }
+ } catch (UserException e) {
+ throwOrWarn(e);
+ }
+ }
+
+ if (positionsOnly) {
+ return;
+ }
+
+ //more rigorous checks on annotations and genotypes
+ try {
+ checkVariantContextsAreMatching(actualTrimmed, expectedTrimmed, overlappingDels, vcShouldBeSkippedDueToNoCalls);
+ } catch (UserException e) {
+ final boolean hasLowQualityGenotype = expectedTrimmed.getGenotypes().stream().anyMatch(g -> g.getGQ() < 20);
+ if (!muteDiffs || !(alleleNumberIsDifferent || inbreedingCoeffIsDifferent || hasLowQualityGenotype)) {
+ throwOrWarn(e);
+ }
+ }
+
+ if (alleleNumberIsDifferent && !muteDiffs) {
+ logger.warn("Observed allele number differed at position " + vc.getContig() + ":" + vc.getStart());
+ }
+ if (inbreedingCoeffIsDifferent && expectedTrimmed.getGenotypes().stream().anyMatch(g -> g.getGQ() < 20) && !muteDiffs) {
+ logger.warn("Low quality genotype may have caused inbreeding coeff differences at position " + vc.getContig() + ":" + vc.getStart());
+ }
+ }
+ }
+ }
+
+ private void validateSingleSampleDeletions(final VariantContext vc, final VariantContext match,
+ final VariantContext expectedTrimmed, final VariantContext actualTrimmed,
+ final List overlappingDels, final boolean nearbyGq0) {
+ final Genotype expectedGenotype = expectedTrimmed.getGenotype(0);
+ final List expectedGenotypeAlleles = expectedGenotype.getAlleles();
+ final Genotype actualGenotype = actualTrimmed.getGenotype(0);
+ final List actualGenotypeAlleles = actualGenotype.getAlleles();
+
+ if (ignoreGenotypePhasing && !actualHasNewAlleles(expectedTrimmed, actualTrimmed)) {
+ return;
+ }
+
+ // make sure both genotypes are not haploid and do use order here so we can check phasing
+ if (!expectedGenotypeAlleles.get(0).equals(actualGenotypeAlleles.get(0)) || ( expectedGenotypeAlleles.size() > 1 && actualGenotypeAlleles.size() > 1 && !expectedGenotypeAlleles.get(1).equals(actualGenotypeAlleles.get(1)))) {
+ if (expectedGenotype.isPhased() && expectedGenotypeAlleles.get(1).equals(actualGenotypeAlleles.get(0)) && expectedGenotypeAlleles.get(0).equals(actualGenotypeAlleles.get(1))) {
+ throw wrapWithPosition(vc.getContig(), vc.getStart(), new UserException("phasing is swapped. Actual in phaseset "
+ + actualGenotype.getExtendedAttribute(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY, "") + " has "
+ + actualGenotype.getGenotypeString(false)
+ + " expected in phaseset " + expectedGenotype.getExtendedAttribute(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY, "")
+ +" has " + expectedGenotype.getGenotypeString(false)));
+ }
+ if (overlappingDels.size() == 0) {
+ throw wrapWithPosition(vc.getContig(), vc.getStart(), makeGenotypeExceptionFromDifference("called genotype alleles",
+ actualGenotype.getGenotypeString(false),
+ expectedGenotype.getGenotypeString(false)));
+ }
+
+ //this is okay if a star got corrected by dropping an upstream hom ref del
+ //could be dropped hom ref
+ if (overlappingDels.stream().anyMatch(v -> v.getGenotype(0).isHomRef())) {
+ try {
+ checkAttributes(vc.getAttributes(), match.getAttributes(), vc.getAlleles(), match.getAlleles(), vc.getPhredScaledQual(), getDpForQual(vc), getAsAdSum(vc));
+ } catch (UserException e) {
+ throw wrapWithPosition(vc.getContig(), vc.getStart(), new UserException("INFO attributes do not match at " + vc.getContig() + ":" + vc.getStart()));
+ }
+ try {
+ checkGenotypes(expectedGenotype, actualGenotype, nearbyGq0);
+ } catch (final UserException e) {
+ throw wrapWithPosition(vc.getContig(), vc.getStart(), e);
+ }
+ } else if (overlappingDels.size() == 0 && expectedGenotypeAlleles.contains(Allele.SPAN_DEL) && !actualGenotypeAlleles.contains(Allele.SPAN_DEL)) {
+ return;
+ } else {
+ if (overlappingDels.stream().anyMatch(v -> v.getSource().equals("actual"))) { //there should be one overlapping that got trimmed and the other shouldn't overlap anymore
+ if (!allowMissingStars) {
+ throw wrapWithPosition(vc.getContig(), vc.getStart(), new UserException("genotype alleles do not match at spanning deletion site. "
+ + vc.getSource() + " has " + expectedGenotype.toString() + " and " + match.getSource() + " has "
+ + actualGenotype.toString()));
+ }
+ return;
+ } else {
+ final VariantContext overlapper = overlappingDels.get(0);
+ final VariantContext trimmedOverlapper = trimAlleles(overlapper, overlappingDels);
+ if (overlapper.overlaps(vc) && !trimmedOverlapper.overlaps(vc)) {
+ return;
+ }
+ }
+
+ }
+ }
+ }
+
+ private static int getAsAdSum(VariantContext vc) {
+ List AS_ADs = AS_QualByDepth.getAlleleDepths(vc.getGenotypes());
+ return AS_ADs == null ? 0 : AS_ADs.stream().mapToInt(Integer::intValue).sum();
+ }
+
+ private static int getDpForQual(VariantContext vc) {
+ return QualByDepth.getDepth(vc.getGenotypes(), null);
+ }
+
+ private boolean hasGoodEvidence(final VariantContext vc) {
+ return (vc.getPhredScaledQual() > goodQualThreshold)
+ && (vc.getAttributeAsInt(VCFConstants.DEPTH_KEY,0)/(double)vc.getAttributeAsInt(VCFConstants.ALLELE_NUMBER_KEY,0)) > 5
+ && vc.getGenotypes().stream().filter(g -> !g.isHomRef()).mapToInt(Genotype::getGQ).sum() > goodQualThreshold;
+ }
+
+ private void checkVariantContextsAreMatching(final VariantContext actual, final VariantContext expected,
+ final List overlappingDels, boolean nearbyGq0) {
+ if (!actual.getContig().equals(expected.getContig())) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), new UserException("contigs differ for VCs"));
+ }
+
+ if (actual.getStart() != expected.getStart()) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), new UserException("start positions differ for VCs"));
+ }
+
+ //don't check end in case we're being lenient about alleles
+
+ //check alleles
+ if (!(ignoreGq0 && nearbyGq0)) {
+ if (actualHasNewAlleles(expected, actual)) {
+ try {
+ checkAlleles(actual, expected, overlappingDels);
+ } catch (final UserException e) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), e);
+ }
+ }
+ }
+
+ if (!ignoreDbsnp && !actual.getID().equals(expected.getID())) { //more alleles might mean more dbSNP matches
+ //alleles are checked elsewhere, so if that check was ignored the dbsnp ID check should be ignored too
+ if (actual.getAlleles().size() == expected.getAlleles().size()) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), new UserException("dbsnp IDs differ for VCs"));
+ }
+ }
+
+ if (ignoreAnnotations) {
+ return;
+ }
+
+ if (!ignoreQuals) {
+ if (!(ignoreGq0 && nearbyGq0)) {
+ final double diff = Math.abs(actual.getPhredScaledQual() - expected.getPhredScaledQual());
+ if (diff > qualTolerance) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), new UserException("qual scores differ by " + diff + ", which is more than " + qualTolerance));
+ }
+ }
+ }
+
+ if (!(ignoreGq0 && nearbyGq0)) {
+ try {
+ checkAttributes(actual.getAttributes(), expected.getAttributes(), actual.getAlternateAlleles(), expected.getAlternateAlleles(), expected.getPhredScaledQual(),
+ getDpForQual(expected), getAsAdSum(expected));
+ } catch (final UserException e) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), e);
+ }
+ }
+
+ if (!alleleNumberIsDifferent) {
+ if (!ignoreFilters) {
+ if (actual.filtersWereApplied() != expected.filtersWereApplied()) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), new UserException(" filters were not applied to both variants"));
+ }
+ if (!actual.getFilters().equals(expected.getFilters())) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), new UserException("variants have different filters: expected has "
+ + (expected.getFilters().isEmpty() ? "PASS" : expected.getFilters()) + " and actual has " + actual.getFilters()));
+ }
+ }
+
+ for (int i = 0; i < actual.getGenotypes().size(); i++) {
+ try {
+ checkGenotypes(actual.getGenotype(0), expected.getGenotype(0), nearbyGq0);
+ } catch (final UserException e) {
+ throw wrapWithPosition(expected.getContig(), expected.getStart(), e);
+ }
+ }
+ }
+ }
+
+ private void checkAlleles(final VariantContext actual, final VariantContext expected, final List overlappingDels) {
+ //TODO: this could be more rigorous
+ if (ignoreGenotypePhasing && !actualHasNewAlleles(actual, expected)) {
+ return;
+ }
+ if (!allowExtraAlleles && actualHasNewAlleles(actual, expected)) {
+ throw new UserException("Alleles are mismatched at " + actual.getContig() + ":" + actual.getStart() + ": actual has "
+ + actual.getAlternateAlleles() + " and expected has " + expected.getAlternateAlleles());
+
+ } else if (!allowNewStars && hasNewStar(actual, expected)) {
+ if (overlappingDels.size() == 0 || !overlappingDels.stream().anyMatch(vc -> vc.getSource().equals("actual"))) {
+ throw new UserException("Actual has new unmatched * allele. Alleles are mismatched at " + actual.getContig() + ":" + actual.getStart() + ": actual has "
+ + actual.getAlternateAlleles() + " and expected has " + expected.getAlternateAlleles());
+ }
+ //this is a GenomicsDB/CombineGVCFs bug -- there is an overlapping deletion and * should be output, but those tools don't account for multiple haplotypes and upstream variant "ends" the deletion
+ } else if (!allowMissingStars && hasMissingStar(actual, expected)) {
+ final Set remainder = new LinkedHashSet<>(expected.getAlleles());
+ remainder.removeAll(actual.getAlleles());
+ if (remainder.size() > 1 || !remainder.contains(Allele.SPAN_DEL)) {
+ throw new UserException("Actual missing * allele. Alleles are mismatched at " + actual.getContig() + ":" + actual.getStart() + ": actual has "
+ + actual.getAlternateAlleles() + " and expected has " + expected.getAlternateAlleles());
+ }
+ } else {
+ throw new UserException("Alleles are mismatched at " + actual.getContig() + ":" + actual.getStart() + ": actual has "
+ + actual.getAlternateAlleles() + " and expected has " + expected.getAlternateAlleles());
+ }
+ }
+
+ /**
+ * Compare VC alternate alleles
+ * @param actual
+ * @param expected
+ * @return
+ */
+ private boolean actualHasNewAlleles(final VariantContext actual, final VariantContext expected) {
+ return !actual.getAlternateAlleles().stream().allMatch(a -> GATKVariantContextUtils.isAlleleInList(
+ actual.getReference(), a, expected.getReference(), expected.getAlternateAlleles()));
+ }
+
+ private boolean hasMissingStar(final VariantContext actual, final VariantContext expected) {
+ return expected.getAlleles().contains(Allele.SPAN_DEL) && !actual.getAlleles().contains(Allele.SPAN_DEL);
+ }
+
+ private boolean hasNewStar(final VariantContext actual, final VariantContext expected) {
+ return actual.getAlleles().contains(Allele.SPAN_DEL) && !expected.getAlleles().contains(Allele.SPAN_DEL);
+ }
+
+ private UserException wrapWithPosition(final String contig, final int start, final UserException e) {
+ return new UserException("At position " + contig + ":" +start + " " + e.getMessage());
+ }
+
+ @SuppressWarnings("unchecked")
+ private void checkAttributes(final Map actual, final Map expected,
+ final List actualAlts, final List expectedAlts,
+ final double expectedQual, final int expectedDP, final int expectedAS_AD) {
+ final Set expectedKeys = new LinkedHashSet<>(expected.keySet());
+
+ //do a precheck on AN because then we can't expect the rest of the annotations to match
+ alleleNumberIsDifferent = actual.containsKey(VCFConstants.ALLELE_NUMBER_KEY) && !isAttributeValueEqual(VCFConstants.ALLELE_NUMBER_KEY, actual.get(VCFConstants.ALLELE_NUMBER_KEY),
+ expected.get(VCFConstants.ALLELE_NUMBER_KEY));
+
+ if (ignoreSomeMultiAllelics && actual.containsKey(VCFConstants.ALLELE_COUNT_KEY)) {
+ final Object actualAcValue = actual.get(VCFConstants.ALLELE_COUNT_KEY);
+ final Object expectedAcValue = expected.get(VCFConstants.ALLELE_COUNT_KEY);
+ final List