\
+-qsub -jobQueue gsa -jobResReq virtual_free=5G -sc 50
+```
+
+The HapMap bams get run as tumors without normals because we're not interested in specificity here, so we don't need the normals to filter out noise
+
+###Inputs
+Bam lists:
+5- and 10-plex have 3 replicates, 20-plex has 9
+
+Interval files:
+If we're only interested in sensitivity, then we only need to run the caller around known true positive sites, which we take from the truth VCFs
+This workaround repeats the truth filename for the number of bams -- in theory each could have a separate truth VCF, but they are the same titration mixture
+
+SNP/INDEL counts:
+This is the number of events in the truth VCFs so we can find the sensitivity across all samples
+TODO: this could be generalized
+
+###Outputs
+Each run outputs its own SNP and INDEL sensitivity combined across all samples:
+```
+Sensitivity across all samples:
+SNPs: 0.95156
+INDELs: 0.7328859
+```
+
+Note that these are not filtered for depth as described in the CRSP documentation
+
+###Resources
+Truth file preparation for 5-plex:
+Start with /cga/tcga-gsc/benchmark/data/crsp-truth/1kg_5plex_wgs_hc_calls.codingIndelSnp.db135.recode.vcf
+Select out allele fraction greater than 20% using "vc.isBiallelic() ? AF >= 0.2 : vc.hasGenotypes() && vc.getCalledChrCount(vc.getAltAlleleWithHighestAlleleCount())/(1.0*vc.getCalledChrCount()) >= 0.2"
+
+Similar for 10-plex source:
+/cga/tcga-gsc/benchmark/data/crsp-truth/1kg_10plex_wgs_hc_calls.codingIndelSnp.db135.recode.vcf
+
+And 20-plex source:
+/cga/tcga-gsc/benchmark/data/crsp-truth/1kg_20plex_wgs_hc_calls.codingIndelSnp.db135.recode.vcf
+
+both also using AF filter of 0.2
\ No newline at end of file
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/MuTect2.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/MuTect2.java
index 3c9bc8402..5e04f39d1 100644
--- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/MuTect2.java
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/MuTect2.java
@@ -71,6 +71,7 @@ import org.broadinstitute.gatk.tools.walkers.haplotypecaller.*;
import org.broadinstitute.gatk.tools.walkers.haplotypecaller.readthreading.ReadThreadingAssembler;
import org.broadinstitute.gatk.utils.GenomeLoc;
import org.broadinstitute.gatk.utils.GenomeLocParser;
+import org.broadinstitute.gatk.utils.MathUtils;
import org.broadinstitute.gatk.utils.QualityUtils;
import org.broadinstitute.gatk.utils.activeregion.ActiveRegion;
import org.broadinstitute.gatk.utils.activeregion.ActiveRegionReadState;
@@ -108,26 +109,19 @@ import static java.lang.Math.pow;
/**
* Call somatic SNPs and indels via local re-assembly of haplotypes
*
- * MuTect2 is a somatic SNP and indel caller that combines the DREAM challenge-winning somatic genotyping engine of the original MuTect (Cibulskis et al., 2013) with the assembly-based machinery of HaplotypeCaller.
- *
- * The basic operation of MuTect2 proceeds similarly to that of the HaplotypeCaller
+ * MuTect2 is a somatic SNP and indel caller that combines the DREAM challenge-winning somatic genotyping engine of the original MuTect (Cibulskis et al., 2013) with the assembly-based machinery of HaplotypeCaller. The basic operation of MuTect2 proceeds similarly to that of the HaplotypeCaller.
*
* Differences from HaplotypeCaller
- * While the HaplotypeCaller relies on a ploidy assumption (diploid by default) to inform its genotype likelihood and
- * variant quality calculations, MuTect2 allows for a varying allelic fraction for each variant, as is often seen in tumors with purity less
- * than 100%, multiple subclones, and/or copy number variation (either local or aneuploidy). MuTect2 also differs from the HaplotypeCaller in that it does apply some hard filters
- * to variants before producing output.
+ * While the HaplotypeCaller relies on a ploidy assumption (diploid by default) to inform its genotype likelihood and variant quality calculations, MuTect2 allows for a varying allelic fraction for each variant, as is often seen in tumors with purity less than 100%, multiple subclones, and/or copy number variation (either local or aneuploidy). MuTect2 also differs from the HaplotypeCaller in that it does apply some hard filters to variants before producing output.
+ * Note that the GVCF generation capabilities of HaplotypeCaller are NOT available in MuTect2, even though some of the relevant arguments are listed below. There are currently no plans to make GVCF calling available in MuTect2.
*
* Usage examples
- * These are example commands that show how to run MuTect2 for typical use cases. Square brackets ("[ ]")
- * indicate optional arguments. Note that parameter values shown here may not be the latest recommended; see the
- * Best Practices documentation for detailed recommendations.
+ * These are example commands that show how to run MuTect2 for typical use cases. Square brackets ("[ ]") indicate optional arguments. Note that parameter values shown here may not be the latest recommended; see the Best Practices documentation for detailed recommendations.
*
*
* Tumor/Normal variant calling
*
- * java
- * -jar GenomeAnalysisTK.jar \
+ * java -jar GenomeAnalysisTK.jar \
* -T MuTect2 \
* -R reference.fasta \
* -I:tumor tumor.bam \
@@ -140,10 +134,9 @@ import static java.lang.Math.pow;
*
* Normal-only calling for panel of normals creation
*
- * java
- * -jar GenomeAnalysisTK.jar
- * -T MuTect2
- * -R reference.fasta
+ * java -jar GenomeAnalysisTK.jar \
+ * -T MuTect2 \
+ * -R reference.fasta \
* -I:tumor normal1.bam \
* [--dbsnp dbSNP.vcf] \
* [--cosmic COSMIC.vcf] \
@@ -152,12 +145,11 @@ import static java.lang.Math.pow;
* -o output.normal1.vcf
*
*
- * For full PON creation, call each of your normals separately in artifact detection mode. Then use CombineVariants to
- * output only sites where a variant was seen in at least two samples:
+ * For full PON creation, call each of your normals separately in artifact detection mode as shown above. Then use CombineVariants to output only sites where a variant was seen in at least two samples:
*
- * java -jar GenomeAnalysisTK.jar
- * -T CombineVariants
- * -R reference.fasta
+ * java -jar GenomeAnalysisTK.jar \
+ * -T CombineVariants \
+ * -R reference.fasta \
* -V output.normal1.vcf -V output.normal2.vcf [-V output.normal2.vcf ...] \
* -minN 2 \
* --setKey "null" \
@@ -169,7 +161,10 @@ import static java.lang.Math.pow;
*
* Caveats
*
- * - MuTect2 currently only supports the calling of a single tumor-normal pair at a time
+ * - As noted in several places in the documentation, MuTect2 is currently released under BETA status; it is NOT recommended for production work and is NOT available for commercial/for-profit licensing.
+ * - MuTect2 currently only supports the calling of a single tumor-normal pair at a time.
+ * - Tumor-only variant calling is possible but it is NOT supported and we will not answer any questions about it until it becomes a supported feature.
+ * - Some of the arguments listed below are not functional; they are exclusive to HaplotypeCaller and are listed here due to technical entanglements in the code. This will be resolved in the upcoming GATK 4 release.
*
*
*/
@@ -189,23 +184,16 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
protected SampleList samplesList;
protected boolean printTCGAsampleHeader = false;
- // fasta reference reader to supplement the edges of the reference sequence
protected CachingIndexedFastaSequenceFile referenceReader;
-
- // the assembly engine
protected LocalAssemblyEngine assemblyEngine = null;
-
- // the likelihoods engine
protected ReadLikelihoodCalculationEngine likelihoodCalculationEngine = null;
-
- // the genotyping engine
- protected HaplotypeCallerGenotypingEngine genotypingEngine = null;
-
+ protected SomaticGenotypingEngine genotypingEngine = null;
+ private HaplotypeBAMWriter haplotypeBAMWriter;
private byte MIN_TAIL_QUALITY;
private double log10GlobalReadMismappingRate;
-
-
+ private static final int REFERENCE_PADDING = 500;
+ private static final byte MIN_TAIL_QUALITY_WITH_ERROR_CORRECTION = 6;
@ArgumentCollection
protected M2ArgumentCollection MTAC = new M2ArgumentCollection();
@@ -216,34 +204,98 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
@ArgumentCollection
protected LikelihoodEngineArgumentCollection LEAC = new LikelihoodEngineArgumentCollection();
-
@Argument(fullName = "debug_read_name", required = false, doc="trace this read name through the calling process")
- public String DEBUG_READ_NAME = null;
+ protected String DEBUG_READ_NAME = null;
@Hidden
@Advanced
@Argument(fullName = "MQ_filtering_level", shortName = "MQthreshold", required = false, doc="Set an alternate MQ threshold for debugging")
final public int MQthreshold = 20;
-
- /***************************************/
- // Reference Metadata inputs
- /***************************************/
/**
- * MuTect2 has the ability to use COSMIC data in conjunction with dbSNP to adjust the threshold for evidence of a variant
- * in the normal. If a variant is present in dbSNP, but not in COSMIC, then more evidence is required from the normal
- * sample to prove the variant is not present in germline.
+ * If a call overlaps with a record from the provided comp track, the INFO field will be annotated
+ * as such in the output with the track name (e.g. -comp:FOO will have 'FOO' in the INFO field).
+ * Records that are filtered in the comp track will be ignored.
+ * Note that 'dbSNP' has been special-cased (see the --dbsnp argument).
*/
- @Input(fullName="cosmic", shortName = "cosmic", doc="VCF file of COSMIC sites", required=false)
- public List> cosmicRod = Collections.emptyList();
+ @Advanced
+ @Input(fullName="comp", shortName = "comp", doc="comparison VCF file", required=false)
+ public List> comps = Collections.emptyList();
+ public List> getCompRodBindings() { return comps; }
/**
- * A panel of normals can be a useful (optional) input to help filter out commonly seen sequencing noise that may appear as low allele-fraction somatic variants.
+ * Which annotations to add to the output VCF file. See the VariantAnnotator -list argument to view available annotations.
*/
- @Input(fullName="normal_panel", shortName = "PON", doc="VCF file of sites observed in normal", required=false)
- public List> normalPanelRod = Collections.emptyList();
+ @Advanced
+ @Argument(fullName="annotation", shortName="A", doc="One or more specific annotations to apply to variant calls", required=false)
+ protected List annotationsToUse = new ArrayList<>(Arrays.asList(new String[]{"DepthPerAlleleBySample", "BaseQualitySumPerAlleleBySample", "TandemRepeatAnnotator", "OxoGReadCounts"}));
- private HaplotypeBAMWriter haplotypeBAMWriter;
+ /**
+ * Which annotations to exclude from output in the VCF file. Note that this argument has higher priority than the -A or -G arguments,
+ * so annotations will be excluded even if they are explicitly included with the other options.
+ */
+ @Advanced
+ @Argument(fullName="excludeAnnotation", shortName="XA", doc="One or more specific annotations to exclude", required=false)
+ protected List annotationsToExclude = new ArrayList<>(Arrays.asList(new String[]{"SpanningDeletions"}));
+
+ /**
+ * Which groups of annotations to add to the output VCF file. See the VariantAnnotator -list argument to view available groups.
+ */
+ @Argument(fullName="group", shortName="G", doc="One or more classes/groups of annotations to apply to variant calls", required=false)
+ protected String[] annotationClassesToUse = { };
+
+ @Output(doc="File to which variants should be written")
+ protected VariantContextWriter vcfWriter = null;
+
+ /**
+ * Active region trimmer reference.
+ */
+ @ArgumentCollection
+ protected ActiveRegionTrimmer trimmer = new ActiveRegionTrimmer();
+
+ @Hidden
+ @Argument(fullName="keepRG", shortName="keepRG", doc="Only use read from this read group when making calls (but use all reads to build the assembly)", required = false)
+ protected String keepRG = null;
+
+ @Argument(fullName = "min_base_quality_score", shortName = "mbq", doc = "Minimum base quality required to consider a base for calling", required = false)
+ public byte MIN_BASE_QUALTY_SCORE = 10;
+
+ public PairHMMLikelihoodCalculationEngine.PCR_ERROR_MODEL pcrErrorModel = PairHMMLikelihoodCalculationEngine.PCR_ERROR_MODEL.HOSTILE;
+
+ @Hidden
+ @Argument(fullName="errorCorrectReads", shortName="errorCorrectReads", doc = "Use an exploratory algorithm to error correct the kmers used during assembly. May cause fundamental problems with the assembly graph itself", required=false)
+ protected boolean errorCorrectReads = false;
+
+ @Hidden
+ @Argument(fullName="captureAssemblyFailureBAM", shortName="captureAssemblyFailureBAM", doc="If specified, we will write a BAM called assemblyFailure.bam capturing all of the reads that were in the active region when the assembler failed for any reason", required = false)
+ protected boolean captureAssemblyFailureBAM = false;
+
+ @Advanced
+ @Argument(fullName="dontUseSoftClippedBases", shortName="dontUseSoftClippedBases", doc="If specified, we will not analyze soft clipped bases in the reads", required = false)
+ protected boolean dontUseSoftClippedBases = false;
+
+ @Hidden
+ @Argument(fullName="justDetermineActiveRegions", shortName="justDetermineActiveRegions", doc = "If specified, the HC won't actually do any assembly or calling, it'll just run the upfront active region determination code. Useful for benchmarking and scalability testing", required=false)
+ protected boolean justDetermineActiveRegions = false;
+
+ /**
+ * As of GATK 3.3, HaplotypeCaller outputs physical (read-based) information (see version 3.3 release notes and documentation for details). This argument disables that behavior.
+ */
+ @Advanced
+ @Argument(fullName="doNotRunPhysicalPhasing", shortName="doNotRunPhysicalPhasing", doc="Disable physical phasing", required = false)
+ protected boolean doNotRunPhysicalPhasing = false;
+
+ /**
+ * When downsampling, level the coverage of the reads in each sample to no more than maxReadsInRegionPerSample reads,
+ * not reducing coverage at any read start to less than minReadsPerAlignmentStart
+ */
+ @Argument(fullName = "maxReadsInRegionPerSample", shortName = "maxReadsInRegionPerSample", doc="Maximum reads in an active region", required = false)
+ protected int maxReadsInRegionPerSample = 1000;
+
+ @Argument(fullName = "minReadsPerAlignmentStart", shortName = "minReadsPerAlignStart", doc="Minimum number of reads sharing the same alignment start for each genomic location in an active region", required = false)
+ protected int minReadsPerAlignmentStart = 5;
+
+ public RodBinding getDbsnpRodBinding() { return MTAC.dbsnp.dbsnp; }
@Override
public void initialize() {
@@ -252,15 +304,15 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
samplesList = new IndexedSampleList(new ArrayList<>(ReadUtils.getSAMFileSamples(getToolkit().getSAMFileHeader())));
// MUTECT: check that we have at least one tumor bam
- for(SAMReaderID id : getToolkit().getReadsDataSource().getReaderIDs()) {
+ for(final SAMReaderID id : getToolkit().getReadsDataSource().getReaderIDs()) {
if (id.getTags().getPositionalTags().size() == 0) {
throw new RuntimeException("BAMs must be tagged as either 'tumor' or 'normal'");
}
// only supports single-sample BAMs (ie first read group is representative)
- String bamSampleName = getToolkit().getReadsDataSource().getHeader(id).getReadGroups().get(0).getSample();
+ final String bamSampleName = getToolkit().getReadsDataSource().getHeader(id).getReadGroups().get(0).getSample();
- for(String tag : id.getTags().getPositionalTags()) {
+ for(final String tag : id.getTags().getPositionalTags()) {
if (BAM_TAG_TUMOR.equalsIgnoreCase(tag)) {
tumorSAMReaderIDs.add(id);
if (tumorSampleName == null) {
@@ -332,10 +384,9 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
final MergeVariantsAcrossHaplotypes variantMerger = new MergeVariantsAcrossHaplotypes();
final GenomeAnalysisEngine toolkit = getToolkit();
- final GenomeLocParser genomeLocParser = toolkit.getGenomeLocParser();
-
- genotypingEngine = new SomaticGenotypingEngine( MTAC, samplesList, genomeLocParser, FixedAFCalculatorProvider.createThreadSafeProvider(getToolkit(), MTAC, logger), !doNotRunPhysicalPhasing, MTAC);
+ genotypingEngine = new SomaticGenotypingEngine( MTAC, samplesList, toolkit.getGenomeLocParser(), !doNotRunPhysicalPhasing, MTAC,
+ tumorSampleName, normalSampleName, DEBUG_READ_NAME);
genotypingEngine.setCrossHaplotypeEventMerger(variantMerger);
genotypingEngine.setAnnotationEngine(annotationEngine);
@@ -355,7 +406,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
trimmer.snpPadding = 50;
samplesList = toolkit.getReadSampleList();
- Set sampleSet = SampleListUtils.asSet(samplesList);
+ final Set sampleSet = SampleListUtils.asSet(samplesList);
if( MTAC.CONTAMINATION_FRACTION_FILE != null )
MTAC.setSampleContamination(AlleleBiasedDownsamplingUtils.loadContaminationFile(MTAC.CONTAMINATION_FRACTION_FILE, MTAC.CONTAMINATION_FRACTION, sampleSet, logger));
@@ -364,9 +415,12 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
private VariantAnnotatorEngine initializeVCFOutput() {
// initialize the output VCF header
+ if (MTAC.ENABLE_CLUSTERED_READ_POSITION_FILTER) {
+ annotationsToUse.add("ClusteredReadPosition");
+ }
final VariantAnnotatorEngine annotationEngine = new VariantAnnotatorEngine(Arrays.asList(annotationClassesToUse), annotationsToUse, annotationsToExclude, this, getToolkit());
- Set headerInfo = new HashSet<>();
+ final Set headerInfo = new HashSet<>();
// all annotation fields from VariantAnnotatorEngine
headerInfo.addAll(annotationEngine.getVCFAnnotationDescriptions());
@@ -382,7 +436,8 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
headerInfo.addAll(getM2HeaderLines());
headerInfo.addAll(getSampleHeaderLines());
- List outputSampleNames = getOutputSampleNames();
+ // if printTCGAsampleHeader, we already checked for exactly 1 tumor and 1 normal in printTCGAsampleHeader assignment in initialize()
+ final List outputSampleNames = printTCGAsampleHeader ? Arrays.asList("TUMOR", "NORMAL") : SampleListUtils.asList(samplesList);
vcfWriter.writeHeader(new VCFHeader(headerInfo, outputSampleNames));
@@ -390,7 +445,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
private Set getM2HeaderLines(){
- Set headerInfo = new HashSet<>();
+ final Set headerInfo = new HashSet<>();
headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.NORMAL_LOD_KEY));
headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.TUMOR_LOD_KEY));
headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.PANEL_OF_NORMALS_COUNT_KEY));
@@ -399,6 +454,13 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.EVENT_DISTANCE_MIN_KEY));
headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.EVENT_DISTANCE_MAX_KEY));
+ if (MTAC.ENABLE_STRAND_ARTIFACT_FILTER){
+ headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.TLOD_FWD_KEY));
+ headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.TLOD_REV_KEY));
+ headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.TUMOR_SB_POWER_FWD_KEY));
+ headerInfo.add(GATKVCFHeaderLines.getInfoLine(GATKVCFConstants.TUMOR_SB_POWER_REV_KEY));
+ }
+
headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.ALLELE_FRACTION_KEY));
headerInfo.add(GATKVCFHeaderLines.getFilterLine(GATKVCFConstants.STR_CONTRACTION_FILTER_NAME));
@@ -410,6 +472,9 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
headerInfo.add(GATKVCFHeaderLines.getFilterLine(GATKVCFConstants.TUMOR_LOD_FILTER_NAME));
headerInfo.add(GATKVCFHeaderLines.getFilterLine(GATKVCFConstants.GERMLINE_RISK_FILTER_NAME));
headerInfo.add(GATKVCFHeaderLines.getFilterLine(GATKVCFConstants.TRIALLELIC_SITE_FILTER_NAME));
+ headerInfo.add(GATKVCFHeaderLines.getFilterLine(GATKVCFConstants.STRAND_ARTIFACT_FILTER_NAME));
+ headerInfo.add(GATKVCFHeaderLines.getFilterLine(GATKVCFConstants.CLUSTERED_READ_POSITION_FILTER_NAME));
+
if ( ! doNotRunPhysicalPhasing ) {
headerInfo.add(GATKVCFHeaderLines.getFormatLine(GATKVCFConstants.HAPLOTYPE_CALLER_PHASING_ID_KEY));
@@ -419,21 +484,21 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
private Set getSampleHeaderLines(){
- Set sampleLines = new HashSet<>();
+ final Set sampleLines = new HashSet<>();
if (printTCGAsampleHeader) {
//NOTE: This will only list the first bam file for each tumor/normal sample if there is more than one
- Map normalSampleHeaderAttributes = new HashMap<>();
+ final Map normalSampleHeaderAttributes = new HashMap<>();
normalSampleHeaderAttributes.put("ID", "NORMAL");
normalSampleHeaderAttributes.put("SampleName", normalSampleName);
if (normalSAMReaderIDs.iterator().hasNext() && !getToolkit().getArguments().disableCommandLineInVCF)
normalSampleHeaderAttributes.put("File", normalSAMReaderIDs.iterator().next().getSamFilePath());
- VCFSimpleHeaderLine normalSampleHeader = new VCFSimpleHeaderLine("SAMPLE", normalSampleHeaderAttributes);
- Map tumorSampleHeaderAttributes = new HashMap<>();
+ final VCFSimpleHeaderLine normalSampleHeader = new VCFSimpleHeaderLine("SAMPLE", normalSampleHeaderAttributes);
+ final Map tumorSampleHeaderAttributes = new HashMap<>();
tumorSampleHeaderAttributes.put("ID", "TUMOR");
tumorSampleHeaderAttributes.put("SampleName", tumorSampleName);
if (tumorSAMReaderIDs.iterator().hasNext() && !getToolkit().getArguments().disableCommandLineInVCF)
tumorSampleHeaderAttributes.put("File", tumorSAMReaderIDs.iterator().next().getSamFilePath());
- VCFSimpleHeaderLine tumorSampleHeader = new VCFSimpleHeaderLine("SAMPLE", tumorSampleHeaderAttributes);
+ final VCFSimpleHeaderLine tumorSampleHeader = new VCFSimpleHeaderLine("SAMPLE", tumorSampleHeaderAttributes);
sampleLines.add(normalSampleHeader);
sampleLines.add(tumorSampleHeader);
@@ -441,28 +506,15 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
return sampleLines;
}
- private List getOutputSampleNames(){
- if (printTCGAsampleHeader) {
- //Already checked for exactly 1 tumor and 1 normal in printTCGAsampleHeader assignment in initialize()
- List sampleNamePlaceholders = new ArrayList<>(2);
- sampleNamePlaceholders.add("TUMOR");
- sampleNamePlaceholders.add("NORMAL");
- return sampleNamePlaceholders;
- }
- else {
- return SampleListUtils.asList(samplesList);
- }
- }
-
@Override
- public ActivityProfileState isActive(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
+ public ActivityProfileState isActive(final RefMetaDataTracker tracker, final ReferenceContext ref, final AlignmentContext context) {
if( context == null || context.getBasePileup().isEmpty() )
// if we don't have any data, just abort early
return new ActivityProfileState(ref.getLocus(), 0.0);
final Map splitContexts = AlignmentContextUtils.splitContextBySampleName(context);
- AlignmentContext tumorContext = splitContexts.get(tumorSampleName);
- AlignmentContext normalContext = splitContexts.get(normalSampleName);
+ final AlignmentContext tumorContext = splitContexts.get(tumorSampleName);
+ final AlignmentContext normalContext = splitContexts.get(normalSampleName);
// if there are no tumor reads... there is no activity!
if (tumorContext == null) {
@@ -470,7 +522,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
// KCIBUL -- this method was inlined and modified from ReferenceConfidenceModel
- ReadBackedPileup tumorPileup = tumorContext.getBasePileup().getMappingFilteredPileup(MQthreshold);
+ final ReadBackedPileup tumorPileup = tumorContext.getBasePileup().getMappingFilteredPileup(MQthreshold);
final double[] tumorGLs = calcGenotypeLikelihoodsOfRefVsAny(tumorPileup, ref.getBase(), MIN_BASE_QUALTY_SCORE);
final double tumorLod = tumorGLs[1] - tumorGLs[0];
@@ -484,7 +536,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
// in any case, we have to handle the case where there is no normal (and thus no normal context) which is
// different than having a normal but having no reads (where we should not enter the active region)
if (normalSampleName != null && normalContext != null) {
- int nonRefInNormal = getCountOfNonRefEvents(normalContext.getBasePileup(), ref.getBase(), MIN_BASE_QUALTY_SCORE);
+ final int nonRefInNormal = getCountOfNonRefEvents(normalContext.getBasePileup(), ref.getBase(), MIN_BASE_QUALTY_SCORE);
final double[] normalGLs = calcGenotypeLikelihoodsOfRefVsAny(normalContext.getBasePileup(), ref.getBase(), MIN_BASE_QUALTY_SCORE, 0.5f);
final double normalLod = normalGLs[0] - normalGLs[1];
@@ -498,7 +550,6 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
prob = 1;
logger.debug("At " + ref.getLocus().toString() + " tlod: " + tumorLod + " and no-normal calling");
}
-
}
return new ActivityProfileState( ref.getLocus(), prob, ActivityProfileState.Type.NONE, null);
@@ -507,58 +558,42 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
private final static List NO_CALLS = Collections.emptyList();
@Override
public List map( final ActiveRegion originalActiveRegion, final RefMetaDataTracker metaDataTracker ) {
- if ( justDetermineActiveRegions )
- // we're benchmarking ART and/or the active region determination code in the HC, just leave without doing any work
+ if ( justDetermineActiveRegions ) {
return NO_CALLS;
-
- if( !originalActiveRegion.isActive() )
- // Not active so nothing to do!
+ } else if( !originalActiveRegion.isActive() ) {
return referenceModelForNoVariation(originalActiveRegion, true);
-
- // No reads here so nothing to do!
- if( originalActiveRegion.size() == 0 ) { return referenceModelForNoVariation(originalActiveRegion, true); }
-
+ } else if( originalActiveRegion.size() == 0 ) {
+ return referenceModelForNoVariation(originalActiveRegion, true);
+ }
logReadInfo(DEBUG_READ_NAME, originalActiveRegion.getReads(), "Present in original active region");
// create the assembly using just high quality reads (Q20 or higher). We want to use lower
// quality reads in the PairHMM (and especially in the normal) later, so we can't use a ReadFilter
- ActiveRegion assemblyActiveRegion = new ActiveRegion(originalActiveRegion.getLocation(), originalActiveRegion.getSupportingStates(),originalActiveRegion.isActive(), getToolkit().getGenomeLocParser(), originalActiveRegion.getExtension());
- for (GATKSAMRecord rec : originalActiveRegion.getReads()) {
+ final ActiveRegion assemblyActiveRegion = new ActiveRegion(originalActiveRegion.getLocation(), originalActiveRegion.getSupportingStates(),originalActiveRegion.isActive(), getToolkit().getGenomeLocParser(), originalActiveRegion.getExtension());
+ for (final GATKSAMRecord rec : originalActiveRegion.getReads()) {
if (rec.getMappingQuality() >= MQthreshold ) {
assemblyActiveRegion.add(rec);
}
}
-
logReadInfo(DEBUG_READ_NAME, assemblyActiveRegion.getReads(), "Present in assembly active region");
// run the local assembler, getting back a collection of information on how we should proceed
- final List givenAlleles = new ArrayList<>();
- final AssemblyResultSet untrimmedAssemblyResult = assembleReads(assemblyActiveRegion, givenAlleles);
-
-
+ final AssemblyResultSet untrimmedAssemblyResult = assembleReads(assemblyActiveRegion, Collections.EMPTY_LIST);
final TreeSet allVariationEvents = untrimmedAssemblyResult.getVariationEvents();
- // TODO - line bellow might be unecessary : it might be that assemblyResult will always have those alleles anyway
- // TODO - so check and remove if that is the case:
- allVariationEvents.addAll(givenAlleles);
-
final ActiveRegionTrimmer.Result trimmingResult = trimmer.trim(originalActiveRegion,allVariationEvents);
-
-
- // Stop the trimming madness!!!
- if (!trimmingResult.isVariationPresent())
- return referenceModelForNoVariation(originalActiveRegion,false);
-
+ if (!trimmingResult.isVariationPresent()) {
+ return referenceModelForNoVariation(originalActiveRegion, false);
+ }
logReadInfo(DEBUG_READ_NAME, trimmingResult.getCallableRegion().getReads(), "Present in trimming result");
final AssemblyResultSet assemblyResult =
trimmingResult.needsTrimming() ? untrimmedAssemblyResult.trimTo(trimmingResult.getCallableRegion()) : untrimmedAssemblyResult;
-// final AssemblyResultSet assemblyResult = untrimmedAssemblyResult;
+ // it is conceivable that the region is active yet has no events upon assembling only the well-mapped reads
+ if( ! assemblyResult.isVariationPresent() ) {
+ return referenceModelForNoVariation(originalActiveRegion, false);
+ }
- // after talking to Ryan -- they grab the reads out of the assembly (and trim then) to pass into the PairHMM
- // because at one point they were trying error correcting of the reads based on the haplotypes.. but that is not
- // working out, so it's safe for us just to take the reads
-//
final ActiveRegion regionForGenotyping = assemblyResult.getRegionForGenotyping();
logReadInfo(DEBUG_READ_NAME, regionForGenotyping.getReads(), "Present in region for genotyping");
@@ -566,11 +601,6 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
haplotypeBAMWriter.addDroppedReadsFromDelta(DroppedReadsTracker.Reason.TRIMMMED, originalActiveRegion.getReads(), regionForGenotyping.getReads());
}
-//
-// final ActiveRegion regionForGenotyping = trimmingResult.getCallableRegion();
-
-// final ActiveRegion regionForGenotyping = originalActiveRegion;
-
// filter out reads from genotyping which fail mapping quality based criteria
//TODO - why don't do this before any assembly is done? Why not just once at the beginning of this method
//TODO - on the originalActiveRegion?
@@ -583,49 +613,36 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
final Map> perSampleFilteredReadList = splitReadsBySample(filteredReads);
-
logReadInfo(DEBUG_READ_NAME, regionForGenotyping.getReads(), "Present in region for genotyping after filtering reads");
- // abort early if something is out of the acceptable range
- // TODO is this ever true at this point??? perhaps GGA. Need to check.
- if( ! assemblyResult.isVariationPresent() )
- return referenceModelForNoVariation(originalActiveRegion, false);
-
- // TODO is this ever true at this point??? perhaps GGA. Need to check.
- if( regionForGenotyping.size() == 0 ) {
- // no reads remain after filtering so nothing else to do!
- return referenceModelForNoVariation(originalActiveRegion, false);
- }
-
// evaluate each sample's reads against all haplotypes
final List haplotypes = assemblyResult.getHaplotypeList();
final Map> reads = splitReadsBySample( regionForGenotyping.getReads() );
- for (List rec : reads.values()) {
+ for (final List rec : reads.values()) {
logReadInfo(DEBUG_READ_NAME, rec, "Present after splitting assemblyResult by sample");
}
+ //TODO: this should be written as
+ //TODO final Map ARreads_origNormalMQ = regionForGenotyping.getReads().stream()
+ //TODO .collect(Collectors.toMap(GATKSAMRecord::getReadName, GATKSAMRecord::getMappingQuality));
+ //TODO: but for some reason sometimes streamifying Mutect2 code causes a MalformedWalkerArgumentsException
+ //TODO: this will probably not occur after the port to GATK4
final HashMap ARreads_origNormalMQ = new HashMap<>();
- for (GATKSAMRecord read : regionForGenotyping.getReads()) {
+ for (final GATKSAMRecord read : regionForGenotyping.getReads()) {
ARreads_origNormalMQ.put(read.getReadName(), read.getMappingQuality());
}
// modify MAPQ scores in normal to be high so that we don't do any base quality score capping
- for(GATKSAMRecord rec : regionForGenotyping.getReads()) {
+ for(final GATKSAMRecord rec : regionForGenotyping.getReads()) {
if (isReadFromNormal(rec)) {
rec.setMappingQuality(60);
}
}
logger.debug("Computing read likelihoods with " + regionForGenotyping.getReads().size() + " reads against " + haplotypes.size() + " haplotypes across region " + assemblyResult.getRegionForGenotyping().toString());
+ final ReadLikelihoods readLikelihoods = likelihoodCalculationEngine.computeReadLikelihoods(assemblyResult,samplesList,reads);
-
- // Calculate the likelihoods: CPU intensive part.
- final ReadLikelihoods readLikelihoods =
- likelihoodCalculationEngine.computeReadLikelihoods(assemblyResult,samplesList,reads);
-
- // Realign reads to their best haplotype.
- // KCIBUL: this is new stuff -- review it!
final Map readRealignments = realignReadsToTheirBestHaplotype(readLikelihoods, assemblyResult.getReferenceHaplotype(), assemblyResult.getPaddedReferenceLoc());
if ( MTAC.bamWriter != null && MTAC.emitDroppedReads ) {
@@ -636,34 +653,18 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
readLikelihoods.changeReads(readRealignments);
-
- for (GATKSAMRecord rec : readRealignments.keySet()) {
+ for (final GATKSAMRecord rec : readRealignments.keySet()) {
logReadInfo(DEBUG_READ_NAME, rec, "Present after computing read likelihoods");
}
- // Note: we used to subset down at this point to only the "best" haplotypes in all samples for genotyping, but there
- // was a bad interaction between that selection and the marginalization that happens over each event when computing
- // GLs. In particular, for samples that are heterozygous non-reference (B/C) the marginalization for B treats the
- // haplotype containing C as reference (and vice versa). Now this is fine if all possible haplotypes are included
- // in the genotyping, but we lose information if we select down to a few haplotypes. [EB]
-
- final HaplotypeCallerGenotypingEngine.CalledHaplotypes calledHaplotypes = ((SomaticGenotypingEngine)genotypingEngine).callMutations(
- haplotypes,
+ final HaplotypeCallerGenotypingEngine.CalledHaplotypes calledHaplotypes = genotypingEngine.callMutations(
readLikelihoods,
ARreads_origNormalMQ,
perSampleFilteredReadList,
assemblyResult.getFullReferenceWithPadding(),
assemblyResult.getPaddedReferenceLoc(),
regionForGenotyping.getLocation(),
- getToolkit().getGenomeLocParser(),
- metaDataTracker,
- givenAlleles, false ,
- tumorSampleName,
- normalSampleName,
- dbsnp.dbsnp,
- cosmicRod,
- DEBUG_READ_NAME
- );
+ metaDataTracker);
if ( MTAC.bamWriter != null ) {
final Set calledHaplotypeSet = new HashSet<>(calledHaplotypes.getCalledHaplotypes());
@@ -682,90 +683,18 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
if( MTAC.DEBUG ) { logger.info("----------------------------------------------------------------------------------"); }
-
-
- List annotatedCalls = new ArrayList<>();
- int eventCount = calledHaplotypes.getCalls().size();
- Integer minEventDistance = null;
- Integer maxEventDistance = null;
- Integer lastPosition = null;
- for (VariantContext vc : calledHaplotypes.getCalls()) {
- if (lastPosition == null) {
- lastPosition = vc.getStart();
- } else {
- int dist = Math.abs(vc.getStart() - lastPosition);
- if (maxEventDistance == null || dist > maxEventDistance) {
- maxEventDistance = dist;
- }
- if (minEventDistance == null || dist < minEventDistance) {
- minEventDistance = dist;
- }
- }
- }
- Map eventDistanceAttributes = new HashMap<>();
- eventDistanceAttributes.put(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY, eventCount);
- eventDistanceAttributes.put(GATKVCFConstants.EVENT_DISTANCE_MIN_KEY, minEventDistance);
- eventDistanceAttributes.put(GATKVCFConstants.EVENT_DISTANCE_MAX_KEY, maxEventDistance);
-
-
- // can we do this with the Annotation classes instead?
- for (VariantContext originalVC : calledHaplotypes.getCalls()) {
- VariantContextBuilder vcb = new VariantContextBuilder(originalVC);
-
- Map attributes = new HashMap<>(originalVC.getAttributes());
- attributes.putAll(eventDistanceAttributes);
- vcb.attributes(attributes);
-
- Set filters = new HashSet<>(originalVC.getFilters());
-
- double tumorLod = originalVC.getAttributeAsDouble(GATKVCFConstants.TUMOR_LOD_KEY, -1);
- if (tumorLod < MTAC.TUMOR_LOD_THRESHOLD) {
- filters.add(GATKVCFConstants.TUMOR_LOD_FILTER_NAME);
- }
-
- // if we are in artifact detection mode, apply the thresholds for the LOD scores
- if (!MTAC.ARTIFACT_DETECTION_MODE) {
- filters.addAll(calculateFilters(metaDataTracker, originalVC, eventDistanceAttributes));
- }
-
- if (filters.size() > 0) {
- vcb.filters(filters);
- } else {
- vcb.passFilters();
- }
-
- if (printTCGAsampleHeader) {
- GenotypesContext genotypesWithBamSampleNames = originalVC.getGenotypes();
- List renamedGenotypes = new ArrayList<>();
- GenotypeBuilder GTbuilder = new GenotypeBuilder(genotypesWithBamSampleNames.get(tumorSampleName));
- GTbuilder.name("TUMOR");
- renamedGenotypes.add(GTbuilder.make());
- GTbuilder = new GenotypeBuilder(genotypesWithBamSampleNames.get(normalSampleName));
- GTbuilder.name("NORMAL");
- renamedGenotypes.add(GTbuilder.make());
- vcb.genotypes(renamedGenotypes);
- }
-
- annotatedCalls.add(vcb.make());
- }
-
-
-
-
-
-
- return annotatedCalls;
+ return annotateVCs(calledHaplotypes, metaDataTracker);
}
- private Set calculateFilters(RefMetaDataTracker metaDataTracker, VariantContext vc, Map eventDistanceAttributes) {
- Set filters = new HashSet<>();
+ private Set calculateFilters(final RefMetaDataTracker metaDataTracker, final VariantContext vc, final Map eventDistanceAttributes) {
+ final Set filters = new HashSet<>();
- Integer eventCount = (Integer) eventDistanceAttributes.get(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY);
- Integer maxEventDistance = (Integer) eventDistanceAttributes.get(GATKVCFConstants.EVENT_DISTANCE_MAX_KEY);
+ final Integer eventCount = (Integer) eventDistanceAttributes.get(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY);
+ final Integer maxEventDistance = (Integer) eventDistanceAttributes.get(GATKVCFConstants.EVENT_DISTANCE_MAX_KEY);
- Collection panelOfNormalsVC = metaDataTracker.getValues(normalPanelRod,
+ final Collection panelOfNormalsVC = metaDataTracker.getValues(MTAC.normalPanelRod,
getToolkit().getGenomeLocParser().createGenomeLoc(vc.getChr(), vc.getStart()));
- VariantContext ponVc = panelOfNormalsVC.isEmpty()?null:panelOfNormalsVC.iterator().next();
+ final VariantContext ponVc = panelOfNormalsVC.isEmpty()?null:panelOfNormalsVC.iterator().next();
if (ponVc != null) {
filters.add(GATKVCFConstants.PON_FILTER_NAME);
@@ -780,13 +709,13 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
double normalF = 0;
int normalAltQualityScoreSum = 0;
if (hasNormal()) {
- Genotype normalGenotype = vc.getGenotype(normalSampleName);
+ final Genotype normalGenotype = vc.getGenotype(normalSampleName);
// NOTE: how do we get the non-ref depth here?
normalAltCounts = normalGenotype.getAD()[1];
normalF = (Double) normalGenotype.getExtendedAttribute(GATKVCFConstants.ALLELE_FRACTION_KEY);
- Object qss = normalGenotype.getExtendedAttribute(GATKVCFConstants.QUALITY_SCORE_SUM_KEY);
+ final Object qss = normalGenotype.getExtendedAttribute(GATKVCFConstants.QUALITY_SCORE_SUM_KEY);
if (qss != null) {
normalAltQualityScoreSum = (Integer) ((Object[]) qss)[1];
} else {
@@ -796,26 +725,21 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
if ( (normalAltCounts > MTAC.MAX_ALT_ALLELES_IN_NORMAL_COUNT || normalF > MTAC.MAX_ALT_ALLELE_IN_NORMAL_FRACTION ) && normalAltQualityScoreSum > MTAC.MAX_ALT_ALLELES_IN_NORMAL_QSCORE_SUM) {
filters.add(GATKVCFConstants.ALT_ALLELE_IN_NORMAL_FILTER_NAME);
- } else {
-
- // NOTE: does normal alt counts presume the normal had all these events in CIS?
- if ( eventCount > 1 && normalAltCounts >= 1) {
- filters.add(GATKVCFConstants.MULTI_EVENT_ALT_ALLELE_IN_NORMAL_FILTER_NAME);
- } else if (eventCount >= 3) {
- filters.add(GATKVCFConstants.HOMOLOGOUS_MAPPING_EVENT_FILTER_NAME);
- }
-
+ } else if ( eventCount > 1 && normalAltCounts >= 1) {
+ filters.add(GATKVCFConstants.MULTI_EVENT_ALT_ALLELE_IN_NORMAL_FILTER_NAME);
+ } else if (eventCount >= 3) {
+ filters.add(GATKVCFConstants.HOMOLOGOUS_MAPPING_EVENT_FILTER_NAME);
}
// STR contractions, that is the deletion of one repeat unit of a short repeat (>1bp repeat unit)
// such as ACTACTACT -> ACTACT, are overwhelmingly false positives so we
// hard filter them out by default
if (vc.isIndel()) {
- ArrayList rpa = (ArrayList) vc.getAttribute(GATKVCFConstants.REPEATS_PER_ALLELE_KEY);
- String ru = vc.getAttributeAsString(GATKVCFConstants.REPEAT_UNIT_KEY, "");
+ final ArrayList rpa = (ArrayList) vc.getAttribute(GATKVCFConstants.REPEATS_PER_ALLELE_KEY);
+ final String ru = vc.getAttributeAsString(GATKVCFConstants.REPEAT_UNIT_KEY, "");
if (rpa != null && rpa.size() > 1 && ru.length() > 1) {
- int refCount = (Integer) rpa.get(0);
- int altCount = (Integer) rpa.get(1);
+ final int refCount = (Integer) rpa.get(0);
+ final int altCount = (Integer) rpa.get(1);
if (refCount - altCount == 1) {
filters.add(GATKVCFConstants.STR_CONTRACTION_FILTER_NAME);
@@ -829,10 +753,22 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
filters.add(GATKVCFConstants.CLUSTERED_EVENTS_FILTER_NAME);
}
+ // clustered read position filter
+ if (MTAC.ENABLE_CLUSTERED_READ_POSITION_FILTER){
+ final Double tumorFwdPosMedian = (Double) vc.getAttribute(GATKVCFConstants.MEDIAN_LEFT_OFFSET_KEY);
+ final Double tumorRevPosMedian = (Double) vc.getAttribute(GATKVCFConstants.MEDIAN_RIGHT_OFFSET_KEY);
+ final Double tumorFwdPosMAD = (Double) vc.getAttribute(GATKVCFConstants.MAD_MEDIAN_LEFT_OFFSET_KEY);
+ final Double tumorRevPosMAD = (Double) vc.getAttribute(GATKVCFConstants.MAD_MEDIAN_RIGHT_OFFSET_KEY);
+ //If the variant is near the read end (median threshold) and the positions are very similar (MAD threshold) then filter
+ if ( (tumorFwdPosMedian != null && tumorFwdPosMedian <= MTAC.PIR_MEDIAN_THRESHOLD && tumorFwdPosMAD != null && tumorFwdPosMAD <= MTAC.PIR_MAD_THRESHOLD) ||
+ (tumorRevPosMedian != null && tumorRevPosMedian <= MTAC.PIR_MEDIAN_THRESHOLD && tumorRevPosMAD != null && tumorRevPosMAD <= MTAC.PIR_MAD_THRESHOLD))
+ filters.add(GATKVCFConstants.CLUSTERED_READ_POSITION_FILTER_NAME);
+ }
+ // TODO: Move strand bias filter here
+
return filters;
}
-
private final static byte REF_MODEL_DELETION_QUAL = (byte) 30;
/**
* Calculate the genotype likelihoods for the sample in pileup for being hom-ref contrasted with being ref vs. alt
@@ -844,14 +780,15 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
*/
protected double[] calcGenotypeLikelihoodsOfRefVsAny(final ReadBackedPileup pileup, final byte refBase, final byte minBaseQual, final double f) {
final double[] genotypeLikelihoods = new double[2];
- int AA = 0, AB=1;
+ final int AA = 0;
+ final int AB=1;
for( final PileupElement p : pileup ) {
final byte qual = (p.isDeletion() ? REF_MODEL_DELETION_QUAL : p.getQual());
if( p.isDeletion() || qual > minBaseQual ) {
// TODO: why not use base qualities here?
//double pobs = QualityUtils.qualToErrorProbLog10(qual);
- double pobs = 1.0d - pow(10, (30 / -10.0));
+ final double pobs = 1.0d - pow(10, (30 / -10.0));
if( isNonRef(refBase, p)) {
genotypeLikelihoods[AB] += Math.log10(f*pobs + (1-f)*pobs/3.0d);
genotypeLikelihoods[AA] += Math.log10((1-pobs)/3);
@@ -869,6 +806,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
return (normalSampleName != null);
}
+ //TODO: streamify in GATK4
protected int getCountOfNonRefEvents(final ReadBackedPileup pileup, final byte refBase, final byte minBaseQual) {
int i=0;
for( final PileupElement p : pileup ) {
@@ -883,7 +821,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
protected double[] calcGenotypeLikelihoodsOfRefVsAny(final ReadBackedPileup pileup, final byte refBase, final byte minBaseQual) {
- double f = calculateF(pileup, refBase, minBaseQual);
+ final double f = calculateF(pileup, refBase, minBaseQual);
return calcGenotypeLikelihoodsOfRefVsAny(pileup, refBase, minBaseQual, f);
}
@@ -901,11 +839,10 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
}
}
- double f = (double) altCount / ((double) refCount + (double) altCount);
- return f;
+ return (double) altCount / (refCount + altCount);
}
- private boolean isNonRef(byte refBase, PileupElement p) {
+ private boolean isNonRef(final byte refBase, final PileupElement p) {
return p.getBase() != refBase || p.isDeletion() || p.isBeforeDeletionStart() || p.isAfterDeletionEnd() || p.isBeforeInsertion() || p.isAfterInsertion() || p.isNextToSoftClip();
}
@@ -914,37 +851,21 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
protected Set filterNonPassingReads( final ActiveRegion activeRegion) {
final Set readsToRemove = new LinkedHashSet<>();
for( final GATKSAMRecord rec : activeRegion.getReads() ) {
-
+ //TODO: Takuto points out that this is questionable. Let's think hard abut it.
// KCIBUL: only perform read quality filtering on tumor reads...
if (isReadFromNormal(rec)) {
-
if( rec.getReadLength() < MIN_READ_LENGTH ) {
readsToRemove.add(rec);
}
-
- } else {
-
-
- if( rec.getReadLength() < MIN_READ_LENGTH ||
- rec.getMappingQuality() < MQthreshold ||
- BadMateFilter.hasBadMate(rec) ||
-
+ } else if( rec.getReadLength() < MIN_READ_LENGTH || rec.getMappingQuality() < MQthreshold || BadMateFilter.hasBadMate(rec) ||
(keepRG != null && !rec.getReadGroup().getId().equals(keepRG)) ) {
- readsToRemove.add(rec);
- }
+ readsToRemove.add(rec);
}
}
activeRegion.removeAll(readsToRemove);
return readsToRemove;
}
- private static GATKSAMRecord findReadByName(Collection reads, String name) {
- for(GATKSAMRecord read : reads) {
- if (name.equals(read.getReadName())) return read;
- }
- return null;
- }
-
/**
* Instantiates the appropriate likelihood calculation engine.
*
@@ -979,13 +900,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
// enable non primary and extended reads in the active region
@Override
public EnumSet desiredReadStates() {
-// if ( includeUnmappedReads )
-// throw new UserException.BadArgumentValue("includeUnmappedReads", "is not yet functional");
-// else
- return EnumSet.of(
- ActiveRegionReadState.PRIMARY,
- ActiveRegionReadState.NONPRIMARY,
- ActiveRegionReadState.EXTENDED);
+ return EnumSet.of(ActiveRegionReadState.PRIMARY, ActiveRegionReadState.NONPRIMARY, ActiveRegionReadState.EXTENDED);
}
//---------------------------------------------------------------------------------------------------------------
@@ -1000,7 +915,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
@Override
- public Integer reduce(List callsInRegion, Integer numCalledRegions) {
+ public Integer reduce(final List callsInRegion, final Integer numCalledRegions) {
for( final VariantContext call : callsInRegion ) {
vcfWriter.add( call );
}
@@ -1008,7 +923,7 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
}
@Override
- public void onTraversalDone(Integer result) {
+ public void onTraversalDone(final Integer result) {
// if ( SCAC.emitReferenceConfidence == ReferenceConfidenceMode.GVCF ) ((GVCFWriter)vcfWriter).close(false); // GROSS -- engine forces us to close our own VCF writer since we wrapped it
// referenceConfidenceModel.close();
//TODO remove the need to call close here for debugging, the likelihood output stream should be managed
@@ -1017,120 +932,11 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
logger.info("Ran local assembly on " + result + " active regions");
}
-
// The following are not used but are required by the AnnotatorCompatible interface
public RodBinding getSnpEffRodBinding() { return null; }
public List> getResourceRodBindings() { return Collections.emptyList(); }
public boolean alwaysAppendDbsnpId() { return false; }
- /**
- * rsIDs from this file are used to populate the ID column of the output. Also, the DB INFO flag will be set when appropriate.
- * dbSNP overlap is only used to require more evidence of absence in the normal if the variant in question has been seen before in germline.
- */
- @ArgumentCollection
- protected DbsnpArgumentCollection dbsnp = new DbsnpArgumentCollection();
- public RodBinding getDbsnpRodBinding() { return dbsnp.dbsnp; }
-
- /**
- * If a call overlaps with a record from the provided comp track, the INFO field will be annotated
- * as such in the output with the track name (e.g. -comp:FOO will have 'FOO' in the INFO field).
- * Records that are filtered in the comp track will be ignored.
- * Note that 'dbSNP' has been special-cased (see the --dbsnp argument).
- */
- @Advanced
- @Input(fullName="comp", shortName = "comp", doc="comparison VCF file", required=false)
- public List> comps = Collections.emptyList();
- public List> getCompRodBindings() { return comps; }
-
-
-
- /**
- * Which annotations to add to the output VCF file. See the VariantAnnotator -list argument to view available annotations.
- */
- @Advanced
- @Argument(fullName="annotation", shortName="A", doc="One or more specific annotations to apply to variant calls", required=false)
-// protected List annotationsToUse = new ArrayList<>(Arrays.asList(new String[]{"ClippingRankSumTest", "DepthPerSampleHC"}));
- //protected List annotationsToUse = new ArrayList<>(Arrays.asList(new String[]{"DepthPerAlleleBySample", "BaseQualitySumPerAlleleBySample", "TandemRepeatAnnotator",
- // "RMSMappingQuality","MappingQualityRankSumTest","FisherStrand","StrandOddsRatio","ReadPosRankSumTest","QualByDepth", "Coverage"}));
- protected List annotationsToUse = new ArrayList<>(Arrays.asList(new String[]{"DepthPerAlleleBySample", "BaseQualitySumPerAlleleBySample", "TandemRepeatAnnotator", "OxoGReadCounts"}));
-
- /**
- * Which annotations to exclude from output in the VCF file. Note that this argument has higher priority than the -A or -G arguments,
- * so annotations will be excluded even if they are explicitly included with the other options.
- */
- @Advanced
- @Argument(fullName="excludeAnnotation", shortName="XA", doc="One or more specific annotations to exclude", required=false)
- protected List annotationsToExclude = new ArrayList<>(Arrays.asList(new String[]{"SpanningDeletions"}));
-
- /**
- * Which groups of annotations to add to the output VCF file. See the VariantAnnotator -list argument to view available groups.
- */
- @Argument(fullName="group", shortName="G", doc="One or more classes/groups of annotations to apply to variant calls", required=false)
- //protected String[] annotationGroupsToUse = { StandardAnnotation.class.getSimpleName() };
- protected String[] annotationClassesToUse = { };
-
- /**
- * A raw, unfiltered, highly sensitive callset in VCF format.
- */
- @Output(doc="File to which variants should be written")
- protected VariantContextWriter vcfWriter = null;
-
- /**
- * Active region trimmer reference.
- */
- @ArgumentCollection
- protected ActiveRegionTrimmer trimmer = new ActiveRegionTrimmer();
-
- @Hidden
- @Argument(fullName="keepRG", shortName="keepRG", doc="Only use read from this read group when making calls (but use all reads to build the assembly)", required = false)
- protected String keepRG = null;
-
-
-
-
- /**
- * The minimum confidence needed for a given base for it to be used in variant calling.
- */
- @Argument(fullName = "min_base_quality_score", shortName = "mbq", doc = "Minimum base quality required to consider a base for calling", required = false)
- public byte MIN_BASE_QUALTY_SCORE = 10;
-
-
-
-// PAIR-HMM-Related Goodness
-
-// public PairHMMLikelihoodCalculationEngine.PCR_ERROR_MODEL pcrErrorModel = PairHMMLikelihoodCalculationEngine.PCR_ERROR_MODEL.CONSERVATIVE;
-// public PairHMMLikelihoodCalculationEngine.PCR_ERROR_MODEL pcrErrorModel = PairHMMLikelihoodCalculationEngine.PCR_ERROR_MODEL.AGGRESSIVE;
- public PairHMMLikelihoodCalculationEngine.PCR_ERROR_MODEL pcrErrorModel = PairHMMLikelihoodCalculationEngine.PCR_ERROR_MODEL.HOSTILE;
-
- // Parameters to control read error correction
- @Hidden
- @Argument(fullName="errorCorrectReads", shortName="errorCorrectReads", doc = "Use an exploratory algorithm to error correct the kmers used during assembly. May cause fundamental problems with the assembly graph itself", required=false)
- protected boolean errorCorrectReads = false;
-
- @Hidden
- @Argument(fullName="captureAssemblyFailureBAM", shortName="captureAssemblyFailureBAM", doc="If specified, we will write a BAM called assemblyFailure.bam capturing all of the reads that were in the active region when the assembler failed for any reason", required = false)
- protected boolean captureAssemblyFailureBAM = false;
-
- @Advanced
- @Argument(fullName="dontUseSoftClippedBases", shortName="dontUseSoftClippedBases", doc="If specified, we will not analyze soft clipped bases in the reads", required = false)
- protected boolean dontUseSoftClippedBases = false;
-
- @Hidden
- @Argument(fullName="justDetermineActiveRegions", shortName="justDetermineActiveRegions", doc = "If specified, the HC won't actually do any assembly or calling, it'll just run the upfront active region determination code. Useful for benchmarking and scalability testing", required=false)
- protected boolean justDetermineActiveRegions = false;
-
-
-
-
- // reference base padding size
- private static final int REFERENCE_PADDING = 500;
-
- private static final byte MIN_TAIL_QUALITY_WITH_ERROR_CORRECTION = 6;
- private final static int maxReadsInRegionPerSample = 1000; // TODO -- should be an argument
- private final static int minReadsPerAlignmentStart = 5; // TODO -- should be an argument
-
-
-
/**
* High-level function that runs the assembler on the active region reads,
* returning a data structure with the resulting information needed
@@ -1149,15 +955,12 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
final Haplotype referenceHaplotype = createReferenceHaplotype(activeRegion, paddedReferenceLoc);
// Create ReadErrorCorrector object if requested - will be used within assembly engine.
- ReadErrorCorrector readErrorCorrector = null;
- if (errorCorrectReads)
- readErrorCorrector = new ReadErrorCorrector(RTAC.kmerLengthForReadErrorCorrection, MIN_TAIL_QUALITY_WITH_ERROR_CORRECTION, RTAC.minObservationsForKmerToBeSolid, MTAC.DEBUG, fullReferenceWithPadding);
-
+ final ReadErrorCorrector readErrorCorrector = errorCorrectReads ? new ReadErrorCorrector(RTAC.kmerLengthForReadErrorCorrection, MIN_TAIL_QUALITY_WITH_ERROR_CORRECTION,
+ RTAC.minObservationsForKmerToBeSolid, MTAC.DEBUG, fullReferenceWithPadding) : null;
try {
final AssemblyResultSet assemblyResultSet = assemblyEngine.runLocalAssembly( activeRegion, referenceHaplotype, fullReferenceWithPadding, paddedReferenceLoc, giveAlleles,readErrorCorrector );
assemblyResultSet.debugDump(logger);
return assemblyResultSet;
-
} catch ( final Exception e ) {
// Capture any exception that might be thrown, and write out the assembly failure BAM if requested
if ( captureAssemblyFailureBAM ) {
@@ -1202,24 +1005,17 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
if( !clippedRead.isEmpty() && clippedRead.getCigar().getReadLength() > 0 ) {
clippedRead = ReadClipper.hardClipToRegion(clippedRead, activeRegion.getExtendedLoc().getStart(), activeRegion.getExtendedLoc().getStop());
if( activeRegion.readOverlapsRegion(clippedRead) && clippedRead.getReadLength() > 0 ) {
- //logger.info("Keeping read " + clippedRead + " start " + clippedRead.getAlignmentStart() + " end " + clippedRead.getAlignmentEnd());
readsToUse.add(clippedRead);
}
}
}
- // TODO -- Performance optimization: we partition the reads by sample 4 times right now; let's unify that code.
-
final List downsampledReads = DownsamplingUtils.levelCoverageByPosition(ReadUtils.sortReadsByCoordinate(readsToUse), maxReadsInRegionPerSample, minReadsPerAlignmentStart);
if ( MTAC.bamWriter != null && MTAC.emitDroppedReads ) {
haplotypeBAMWriter.addDroppedReadsFromDelta(DroppedReadsTracker.Reason.DOWNSAMPLED, activeRegion.getReads(), downsampledReads);
}
- // handle overlapping read pairs from the same fragment
- // KC: commented out as we handle overlapping read pairs in a different way...
- //cleanOverlappingReadPairs(downsampledReads, normalSampleNames);
-
activeRegion.clearReads();
activeRegion.addAll(downsampledReads);
activeRegion.setFinalized(true);
@@ -1246,32 +1042,27 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
*
* @param reads the list of reads to consider
*/
- private void cleanOverlappingReadPairs(final List reads, Set normalSampleNames) {
- Map> data = splitReadsBySample(reads);
- for ( String sampleName : data.keySet() ) {
+ private void cleanOverlappingReadPairs(final List reads, final Set normalSampleNames) {
+ final Map> data = splitReadsBySample(reads);
+ for ( final String sampleName : data.keySet() ) {
final boolean isTumor = !normalSampleNames.contains(sampleName);
final List perSampleReadList = data.get(sampleName);
final FragmentCollection fragmentCollection = FragmentUtils.create(perSampleReadList);
for ( final List overlappingPair : fragmentCollection.getOverlappingPairs() )
-
- // in MuTect -- right now we compare the
FragmentUtils.adjustQualsOfOverlappingPairedFragments(overlappingPair);
-
-
}
}
- public static void logReadInfo(String readName, Collection records, String message) {
+ public static void logReadInfo(final String readName, final Collection records, final String message) {
if (readName != null) {
- for (GATKSAMRecord rec : records) {
+ for (final GATKSAMRecord rec : records) {
logReadInfo(readName, rec, message);
}
-
}
}
- public static void logReadInfo(String readName, GATKSAMRecord rec, String message) {
+ public static void logReadInfo(final String readName, final GATKSAMRecord rec, final String message) {
if (readName != null && rec != null && readName.equals(rec.getReadName())) {
logger.info("Found " + rec.toString() + " - " + message);
}
@@ -1300,18 +1091,63 @@ public class MuTect2 extends ActiveRegionWalker, Integer> i
return result;
}
- private boolean isReadFromNormal(GATKSAMRecord rec) {
+ private boolean isReadFromNormal(final GATKSAMRecord rec) {
return normalSampleName != null && normalSampleName.equals(rec.getReadGroup().getSample());
-
}
- // KCIBUL: new stuff -- read up on this!!
- /**
- * As of GATK 3.3, HaplotypeCaller outputs physical (read-based) information (see version 3.3 release notes and documentation for details). This argument disables that behavior.
- */
- @Advanced
- @Argument(fullName="doNotRunPhysicalPhasing", shortName="doNotRunPhysicalPhasing", doc="Disable physical phasing", required = false)
- protected boolean doNotRunPhysicalPhasing = false;
+ public String getTumorSampleName(){
+ return tumorSampleName;
+ }
+
+ final List annotateVCs(final HaplotypeCallerGenotypingEngine.CalledHaplotypes calledHaplotypes, final RefMetaDataTracker metaDataTracker) {
+ final int eventCount = calledHaplotypes.getCalls().size();
+ final Map eventDistanceAttributes = new HashMap<>(); //TODO: should be Map -- see TODO below
+ eventDistanceAttributes.put(GATKVCFConstants.EVENT_COUNT_IN_HAPLOTYPE_KEY, eventCount);
+ if (eventCount > 1) {
+ final int lastPosition = calledHaplotypes.getCalls().get(0).getStart();
+ final int[] eventDistances = new int[calledHaplotypes.getCalls().size() - 1];
+ for (int n = 0; n < eventDistances.length; n++) {
+ eventDistances[n] = Math.abs(calledHaplotypes.getCalls().get(n + 1).getStart() - lastPosition);
+ }
+ eventDistanceAttributes.put(GATKVCFConstants.EVENT_DISTANCE_MIN_KEY, MathUtils.arrayMin(eventDistances));
+ eventDistanceAttributes.put(GATKVCFConstants.EVENT_DISTANCE_MAX_KEY, MathUtils.arrayMax(eventDistances));
+ } else { //TODO: putting null is a hack -- we should remove this and update the integration test md5s
+ eventDistanceAttributes.put(GATKVCFConstants.EVENT_DISTANCE_MIN_KEY, null);
+ eventDistanceAttributes.put(GATKVCFConstants.EVENT_DISTANCE_MAX_KEY, null);
+ }
+
+ final List annotatedCalls = new ArrayList<>();
+ // can we do this with the Annotation classes instead?
+ for (final VariantContext originalVC : calledHaplotypes.getCalls()) {
+ final VariantContextBuilder vcb = new VariantContextBuilder(originalVC);
+
+ final Map attributes = new HashMap<>(originalVC.getAttributes());
+ attributes.putAll(eventDistanceAttributes);
+ vcb.attributes(attributes);
+
+ final Set filters = new HashSet<>(originalVC.getFilters());
+
+ final double tumorLod = originalVC.getAttributeAsDouble(GATKVCFConstants.TUMOR_LOD_KEY, -1);
+ if (tumorLod < MTAC.TUMOR_LOD_THRESHOLD) {
+ filters.add(GATKVCFConstants.TUMOR_LOD_FILTER_NAME);
+ }
+
+ // if we are in artifact detection mode, apply the thresholds for the LOD scores
+ if (!MTAC.ARTIFACT_DETECTION_MODE) {
+ filters.addAll(calculateFilters(metaDataTracker, originalVC, eventDistanceAttributes));
+ }
+
+ vcb.filters(filters.isEmpty() ? VariantContext.PASSES_FILTERS : filters);
+
+ if (printTCGAsampleHeader) {
+ final Genotype tumorGenotype = new GenotypeBuilder(originalVC.getGenotype(tumorSampleName)).name("TUMOR").make();
+ final Genotype normalGenotype = new GenotypeBuilder(originalVC.getGenotype(normalSampleName)).name("NORMAL").make();
+ vcb.genotypes(Arrays.asList(tumorGenotype, normalGenotype));
+ }
+ annotatedCalls.add(vcb.make());
+ }
+ return annotatedCalls;
+ }
}
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/NA12878_Evaluations.md b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/NA12878_Evaluations.md
index b7d0c58f2..e5ed5572c 100644
--- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/NA12878_Evaluations.md
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/NA12878_Evaluations.md
@@ -22,16 +22,18 @@ TODO: write a simple tool to do this more easily
To calculate per pair-counts, run:
```
# for SNPs
-for vcf in *.bam.vcf
-do
- cat $vcf | grep PASS | awk '{ if ( length($4) + length($5) == 2) print $0 }' | wc -l
-done
+for vcf in *.vcf
+do
+ cat $vcf | grep PASS | awk '{ if ( length($4) + length($5) == 2) print $0 }' | wc -l
+done > snp-fps.txt
+cat snp-fps.txt | awk '{ sum += $1 } END { print sum }'
# for INDELs
-for vcf in *.bam.vcf
-do
- cat $vcf | grep PASS | awk '{ if ( length($4) + length($5) != 2) print $0 }' | wc -l
-done
+for vcf in *.vcf
+do
+ cat $vcf | grep PASS | awk '{ if ( length($4) + length($5) != 2) print $0 }' | wc -l
+done > indel-fps.txt
+cat indel-fps.txt | awk '{ sum += $1 } END { print sum }'
```
### Current M1 and Indelocator Performance
@@ -72,7 +74,7 @@ and then run the following Queue command
java \
-Djava.io.tmpdir=$TEMPDIR \
-jar $QUEUE_JAR \
- -S $GSA_UNSTABLE_HOME/private/gatk-tools-private/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/run_M2_ICE_NN.scala \
+ -S $GSA_UNSTABLE_HOME/private/gatk-queue-extensions-internal/src/main/qscripts/org/broadinstitute/gatk/queue/qscripts/m2/run_M2_ICE_NN.scala \
-sc 50 \
--job_queue gsa -qsub -jobResReq virtual_free=5G -startFromScratch \
--allbams /humgen/gsa-hpprojects/NA12878Collection/bams/crsp_ice_validation//NA12878.intra.flowcell.replicate.bam_list \
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/PerAlleleCollection.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/PerAlleleCollection.java
new file mode 100644
index 000000000..d5ac3094d
--- /dev/null
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/PerAlleleCollection.java
@@ -0,0 +1,195 @@
+/*
+* By downloading the PROGRAM you agree to the following terms of use:
+*
+* BROAD INSTITUTE
+* SOFTWARE LICENSE AGREEMENT
+* FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
+*
+* This Agreement is made between the Broad Institute, Inc. with a principal address at 415 Main Street, Cambridge, MA 02142 ("BROAD") and the LICENSEE and is effective at the date the downloading is completed ("EFFECTIVE DATE").
+*
+* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
+* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
+* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
+*
+* 1. DEFINITIONS
+* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK3 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute.org/gatk on the EFFECTIVE DATE.
+*
+* 2. LICENSE
+* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM. LICENSEE hereby automatically grants to BROAD a non-exclusive, royalty-free, irrevocable license to any LICENSEE bug fixes or modifications to the PROGRAM with unlimited rights to sublicense and/or distribute. LICENSEE agrees to provide any such modifications and bug fixes to BROAD promptly upon their creation.
+* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
+* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
+* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
+*
+* 3. PHONE-HOME FEATURE
+* LICENSEE expressly acknowledges that the PROGRAM contains an embedded automatic reporting system ("PHONE-HOME") which is enabled by default upon download. Unless LICENSEE requests disablement of PHONE-HOME, LICENSEE agrees that BROAD may collect limited information transmitted by PHONE-HOME regarding LICENSEE and its use of the PROGRAM. Such information shall include LICENSEE'S user identification, version number of the PROGRAM and tools being run, mode of analysis employed, and any error reports generated during run-time. Collection of such information is used by BROAD solely to monitor usage rates, fulfill reporting requirements to BROAD funding agencies, drive improvements to the PROGRAM, and facilitate adjustments to PROGRAM-related documentation.
+*
+* 4. OWNERSHIP OF INTELLECTUAL PROPERTY
+* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
+* Copyright 2012-2016 Broad Institute, Inc.
+* Notice of attribution: The GATK3 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
+* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
+*
+* 5. INDEMNIFICATION
+* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
+*
+* 6. NO REPRESENTATIONS OR WARRANTIES
+* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
+* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
+*
+* 7. ASSIGNMENT
+* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
+*
+* 8. MISCELLANEOUS
+* 8.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
+* 8.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
+* 8.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
+* 8.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
+* 8.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
+* 8.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
+* 8.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
+*/
+
+package org.broadinstitute.gatk.tools.walkers.cancer.m2;
+
+import htsjdk.variant.variantcontext.Allele;
+
+import java.util.*;
+
+/**
+ * A container for allele to value mapping.
+ *
+ * Each PerAlleleCollection may hold a value for each ALT allele and, optionally, a value for the REF allele.
+ * For example,
+ *
+ * PerAlleleCollection alleleFractions = PerAlleleCollection.createPerAltAlleleCollection()
+ *
+ * may be a container for allele fractions for ALT alleles in a variant context. While
+ *
+ * PerAlleleCollection alleleCount = PerAlleleCollection.createPerRefAndAltAlleleCollection()
+ *
+ * may hold the allele counts for the REF allele and all ALT alleles in a variant context.
+ *
+ *
+ **/
+public class PerAlleleCollection {
+ // TODO: consider using Optional for ref allele
+ private Optional refAllele;
+ private Optional refValue;
+ private Map altAlleleValueMap;
+ private boolean altOnly;
+
+ private PerAlleleCollection(boolean altOnly){
+ this.altOnly = altOnly;
+ this.altAlleleValueMap = new HashMap();
+ this.refAllele = Optional.empty();
+
+ }
+
+ public static PerAlleleCollection createPerAltAlleleCollection(){
+ return new PerAlleleCollection(true);
+ }
+
+ public static PerAlleleCollection createPerRefAndAltAlleleCollection(){
+ return new PerAlleleCollection(false);
+ }
+
+ /**
+ * Take an allele, REF or ALT, and update its value appropriately
+ *
+ * @param allele : REF or ALT allele
+ * @param newValue :
+ */
+ public void set(Allele allele, X newValue){
+ if (allele == null || newValue == null){
+ throw new IllegalArgumentException("allele or newValue is null");
+ }
+ if (allele.isReference() && altOnly){
+ throw new IllegalArgumentException("Collection stores values for alternate alleles only");
+ }
+ if (allele.isReference()){
+ this.setRef(allele, newValue);
+ } else {
+ this.setAlt(allele, newValue);
+ }
+ }
+
+ public void setRef(Allele refAllele, X newValue){
+ if (refAllele == null || newValue == null){
+ throw new IllegalArgumentException("refAllele or newValue is null");
+ }
+ if (refAllele.isNonReference()){
+ throw new IllegalArgumentException("Setting Non-reference allele as reference");
+ }
+
+ if (this.refAllele.isPresent()){
+ throw new IllegalArgumentException("Resetting the reference allele not permitted");
+ }
+
+ this.refAllele = Optional.of(refAllele);
+ this.refValue = Optional.of(newValue);
+ }
+
+ public void setAlt(Allele altAllele, X newValue){
+ if (altAllele == null || newValue == null){
+ throw new IllegalArgumentException("altAllele or newValue is null");
+ }
+ if (altAllele.isReference()){
+ throw new IllegalArgumentException("Setting reference allele as alt");
+ }
+
+ altAlleleValueMap.put(altAllele, newValue);
+ }
+
+ /**
+ * Get the value for an allele, REF or ALT
+ * @param allele
+ */
+ public X get(Allele allele){
+ if (allele == null){
+ throw new IllegalArgumentException("allele is null");
+ }
+
+ if (allele.isReference()){
+ if (allele.equals(this.refAllele.get())){
+ return(getRef());
+ } else {
+ throw new IllegalArgumentException("Requested ref allele does not match the stored ref allele");
+ }
+ } else {
+ return(getAlt(allele));
+ }
+ }
+
+ public X getRef(){
+ if (altOnly) {
+ throw new IllegalStateException("Collection does not hold the REF allele");
+ }
+
+ if (this.refAllele.isPresent()){
+ return(refValue.get());
+ } else {
+ throw new IllegalStateException("Collection's ref allele has not been set yet");
+ }
+ }
+
+ public X getAlt(Allele allele){
+ if (allele == null){
+ throw new IllegalArgumentException("allele is null");
+ }
+ if (allele.isReference()){
+ throw new IllegalArgumentException("allele is not an alt allele");
+ }
+
+ if (altAlleleValueMap.containsKey(allele)) {
+ return(altAlleleValueMap.get(allele));
+ } else {
+ throw new IllegalArgumentException("Requested alt allele is not in the collection");
+ }
+
+ }
+
+
+ public Set getAltAlleles(){
+ return(altAlleleValueMap.keySet());
+ }
+}
\ No newline at end of file
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/SomaticGenotypingEngine.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/SomaticGenotypingEngine.java
index 0a0d42b06..a0b259e1e 100644
--- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/SomaticGenotypingEngine.java
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/SomaticGenotypingEngine.java
@@ -51,16 +51,15 @@
package org.broadinstitute.gatk.tools.walkers.cancer.m2;
-import com.google.java.contract.Ensures;
-import htsjdk.samtools.util.StringUtil;
import htsjdk.variant.variantcontext.*;
+import org.apache.commons.lang.mutable.MutableDouble;
+import org.apache.commons.lang.mutable.MutableInt;
import org.apache.log4j.Logger;
-import org.broadinstitute.gatk.tools.walkers.genotyper.GenotypeLikelihoodsCalculationModel;
+import org.broadinstitute.gatk.tools.walkers.genotyper.afcalc.AFCalculator;
import org.broadinstitute.gatk.tools.walkers.genotyper.afcalc.AFCalculatorProvider;
import org.broadinstitute.gatk.tools.walkers.haplotypecaller.HaplotypeCallerGenotypingEngine;
import org.broadinstitute.gatk.utils.GenomeLoc;
import org.broadinstitute.gatk.utils.GenomeLocParser;
-import org.broadinstitute.gatk.utils.commandline.RodBinding;
import org.broadinstitute.gatk.utils.contexts.ReferenceContext;
import org.broadinstitute.gatk.utils.genotyper.MostLikelyAllele;
import org.broadinstitute.gatk.utils.genotyper.PerReadAlleleLikelihoodMap;
@@ -77,13 +76,41 @@ import java.util.*;
public class SomaticGenotypingEngine extends HaplotypeCallerGenotypingEngine {
- protected M2ArgumentCollection MTAC;
+ private final M2ArgumentCollection MTAC;
+ private final TumorPowerCalculator strandArtifactPowerCalculator;
+
+ private final String tumorSampleName;
+ private final String matchedNormalSampleName;
+ private final String DEBUG_READ_NAME;
+
+ //Mutect2 does not run in GGA mode
+ private static final List NO_GIVEN_ALLELES = Collections.EMPTY_LIST;
+
+ // {@link GenotypingEngine} requires a non-null {@link AFCalculatorProvider} but this class doesn't need it. Thus we make a dummy
+ private static AFCalculatorProvider DUMMY_AF_CALCULATOR_PROVIDER = new AFCalculatorProvider() {
+ public AFCalculator getInstance(final int ploidy, final int maximumAltAlleles) { return null; }
+ };
private final static Logger logger = Logger.getLogger(SomaticGenotypingEngine.class);
- public SomaticGenotypingEngine(final M2ArgumentCollection configuration, final SampleList samples, final GenomeLocParser genomeLocParser, final AFCalculatorProvider afCalculatorProvider, final boolean doPhysicalPhasing, final M2ArgumentCollection MTAC) {
- super(configuration, samples, genomeLocParser, afCalculatorProvider, doPhysicalPhasing);
+ public SomaticGenotypingEngine(final M2ArgumentCollection configuration,
+ final SampleList samples,
+ final GenomeLocParser genomeLocParser,
+ final boolean doPhysicalPhasing,
+ final M2ArgumentCollection MTAC,
+ final String tumorSampleName,
+ final String matchedNormalSampleName,
+ final String DEBUG_READ_NAME) {
+ super(configuration, samples, genomeLocParser, DUMMY_AF_CALCULATOR_PROVIDER, doPhysicalPhasing);
this.MTAC = MTAC;
+ this.tumorSampleName = tumorSampleName;
+ this.matchedNormalSampleName = matchedNormalSampleName;
+ this.DEBUG_READ_NAME = DEBUG_READ_NAME;
+
+ // coverage related initialization
+ //TODO: in GATK4, use a QualityUtils method
+ final double errorProbability = Math.pow(10, -MTAC.POWER_CONSTANT_QSCORE/10);
+ strandArtifactPowerCalculator = new TumorPowerCalculator(errorProbability, MTAC.STRAND_ARTIFACT_LOD_THRESHOLD, 0.0f);
}
/**
@@ -91,373 +118,392 @@ public class SomaticGenotypingEngine extends HaplotypeCallerGenotypingEngine {
* genotype likelihoods and assemble into a list of variant contexts and genomic events ready for calling
*
* The list of samples we're working with is obtained from the readLikelihoods
- *
- * @param haplotypes Haplotypes to assign likelihoods to
* @param readLikelihoods Map from reads->(haplotypes,likelihoods)
* @param perSampleFilteredReadList Map from sample to reads that were filtered after assembly and before calculating per-read likelihoods.
* @param ref Reference bytes at active region
* @param refLoc Corresponding active region genome location
* @param activeRegionWindow Active window
- * @param genomeLocParser GenomeLocParser
- * @param activeAllelesToGenotype Alleles to genotype
- * @param emitReferenceConfidence whether we should add a <NON_REF> alternative allele to the result variation contexts.
*
* @return A CalledHaplotypes object containing a list of VC's with genotyped events and called haplotypes
*
*/
-// @Requires({"refLoc.containsP(activeRegionWindow)", "haplotypes.size() > 0"})
- @Ensures("result != null")
- // TODO - can this be refactored? this is hard to follow!
- public HaplotypeCallerGenotypingEngine.CalledHaplotypes callMutations (
- final List haplotypes,
- //final Map haplotypeReadMap,
- final ReadLikelihoods readLikelihoods,
- final Map originalNormalReadQualities,
- final Map> perSampleFilteredReadList,
- final byte[] ref,
- final GenomeLoc refLoc,
- final GenomeLoc activeRegionWindow,
- final GenomeLocParser genomeLocParser,
- final RefMetaDataTracker tracker,
- final List activeAllelesToGenotype,
- final boolean emitReferenceConfidence,
- final String tumorSampleName,
- final String matchedNormalSampleName,
- final RodBinding dbsnpRod,
- final List> cosmicRod,
- final String DEBUG_READ_NAME
-
- ) {
-
- // sanity check input arguments
- if (haplotypes == null || haplotypes.isEmpty()) throw new IllegalArgumentException("haplotypes input should be non-empty and non-null, got "+haplotypes);
+ public CalledHaplotypes callMutations (
+ final ReadLikelihoods readLikelihoods,
+ final Map originalNormalReadQualities,
+ final Map> perSampleFilteredReadList,
+ final byte[] ref,
+ final GenomeLoc refLoc,
+ final GenomeLoc activeRegionWindow,
+ final RefMetaDataTracker tracker) {
+ //TODO: in GATK4 use Utils.nonNull
if (readLikelihoods == null || readLikelihoods.sampleCount() == 0) throw new IllegalArgumentException("readLikelihoods input should be non-empty and non-null, got "+readLikelihoods);
if (ref == null || ref.length == 0 ) throw new IllegalArgumentException("ref bytes input should be non-empty and non-null, got "+ref);
if (refLoc == null || refLoc.size() != ref.length) throw new IllegalArgumentException(" refLoc must be non-null and length must match ref bytes, got "+refLoc);
if (activeRegionWindow == null ) throw new IllegalArgumentException("activeRegionWindow must be non-null, got "+activeRegionWindow);
- if (activeAllelesToGenotype == null ) throw new IllegalArgumentException("activeAllelesToGenotype must be non-null, got "+activeAllelesToGenotype);
- if (genomeLocParser == null ) throw new IllegalArgumentException("genomeLocParser must be non-null, got "+genomeLocParser);
+ final List haplotypes = readLikelihoods.alleles();
// Somatic Tumor/Normal Sample Handling
- verifySamplePresence(tumorSampleName, readLikelihoods.samples());
- final boolean hasNormal = (matchedNormalSampleName != null);
+ if (!readLikelihoods.samples().contains(tumorSampleName)) {
+ throw new IllegalArgumentException("readLikelihoods does not contain the tumor sample " + tumorSampleName);
+ }
+ final boolean hasNormal = matchedNormalSampleName != null;
// update the haplotypes so we're ready to call, getting the ordered list of positions on the reference
// that carry events among the haplotypes
- final TreeSet startPosKeySet = decomposeHaplotypesIntoVariantContexts(haplotypes, readLikelihoods, ref, refLoc, activeAllelesToGenotype);
+ final TreeSet startPosKeySet = decomposeHaplotypesIntoVariantContexts(haplotypes, readLikelihoods, ref, refLoc, NO_GIVEN_ALLELES);
// Walk along each position in the key set and create each event to be outputted
final Set calledHaplotypes = new HashSet<>();
final List returnCalls = new ArrayList<>();
for( final int loc : startPosKeySet ) {
- if( loc >= activeRegionWindow.getStart() && loc <= activeRegionWindow.getStop() ) { // genotyping an event inside this active region
- final List eventsAtThisLoc = getVCsAtThisLocation(haplotypes, loc, activeAllelesToGenotype);
-
- if( eventsAtThisLoc.isEmpty() ) { continue; }
-
- // Create the event mapping object which maps the original haplotype events to the events present at just this locus
- final Map> eventMapper = createEventMapper(loc, eventsAtThisLoc, haplotypes);
-
- // Sanity check the priority list for mistakes
- final List priorityList = makePriorityList(eventsAtThisLoc);
-
- // Merge the event to find a common reference representation
-
- VariantContext mergedVC = GATKVariantContextUtils.simpleMerge(eventsAtThisLoc, priorityList,
- GATKVariantContextUtils.FilteredRecordMergeType.KEEP_IF_ANY_UNFILTERED,
- GATKVariantContextUtils.GenotypeMergeType.PRIORITIZE, false, false, null, false, false);
-
- if( mergedVC == null ) { continue; }
-
- final int numAlts = mergedVC.getNAlleles()-1;
-
-// final VariantContextBuilder vcb = new VariantContextBuilder(mergedVC);
-
- final GenotypeLikelihoodsCalculationModel.Model calculationModel = mergedVC.isSNP()
- ? GenotypeLikelihoodsCalculationModel.Model.SNP : GenotypeLikelihoodsCalculationModel.Model.INDEL;
-
- if (emitReferenceConfidence)
- mergedVC = addNonRefSymbolicAllele(mergedVC);
-
- final Map mergeMap = new LinkedHashMap<>();
- mergeMap.put(null, mergedVC.getReference()); // the reference event (null) --> the reference allele
- for(int iii = 0; iii < eventsAtThisLoc.size(); iii++) {
- mergeMap.put(eventsAtThisLoc.get(iii), mergedVC.getAlternateAllele(iii)); // BUGBUG: This is assuming that the order of alleles is the same as the priority list given to simpleMerge function
- }
-
- final Map> alleleMapper = createAlleleMapper(mergeMap, eventMapper);
-
- if( configuration.DEBUG && logger != null ) {
- if (logger != null) logger.info("Genotyping event at " + loc + " with alleles = " + mergedVC.getAlleles());
- }
-
- ReadLikelihoods readAlleleLikelihoods = readLikelihoods.marginalize(alleleMapper, genomeLocParser.createPaddedGenomeLoc(genomeLocParser.createGenomeLoc(mergedVC), ALLELE_EXTENSION));
-
- //LDG: do we want to do this before or after pulling out overlapping reads?
- if (MTAC.isSampleContaminationPresent())
- readAlleleLikelihoods.contaminationDownsampling(MTAC.getSampleContamination());
-
- //if (!mergedVC.isBiallelic()) {
- // logger.info("[UNSUPPORTED] Detected non-Biallelic VC" + mergedVC.toString());
- // continue;
- //}
-
- // TODO: once tests are passing, refactor to use the new data structure (not the deprecated one)
- // handle overlapping fragments
- // TODO: CONFIRM WITH GSA IF IT IS OK TO REMOVE READS FROM THE PRALM (should be... they do it in filterPoorlyModeledReads!)
- PerReadAlleleLikelihoodMap tumorPRALM = readAlleleLikelihoods.toPerReadAlleleLikelihoodMap(readAlleleLikelihoods.sampleIndex(tumorSampleName));
- filterPRALMForOverlappingReads(tumorPRALM, mergedVC.getReference(), loc, false);
- MuTect2.logReadInfo(DEBUG_READ_NAME, tumorPRALM.getLikelihoodReadMap().keySet(), "Present after filtering for overlapping reads");
- // extend to multiple samples
-
- //handle existence of secondary alts
- double[] afs = estimateAlleleFraction(mergedVC, tumorPRALM);
-
- if( configuration.DEBUG && logger != null ) {
- String output = "Calculated allelic fraction at " + loc + " = ";
- for (int i = 0; i altAlleleFractions = estimateAlleleFraction(mergedVC, tumorPRALM, false);
+ final PerAlleleCollection tumorHetGenotypeLLs = getHetGenotypeLogLikelihoods(mergedVC, tumorPRALM, originalNormalReadQualities, altAlleleFractions);
+
+ final PerAlleleCollection tumorLods = PerAlleleCollection.createPerAltAlleleCollection();
+ for (final Allele altAllele : mergedVC.getAlternateAlleles()){
+ tumorLods.set(altAllele, tumorHetGenotypeLLs.get(altAllele) - tumorHetGenotypeLLs.getRef());
+ }
+
+ // TODO: another good breakpoint e.g. compute normal LOD/set thresholds
+ // TODO: anything related to normal should be encapsulated in Optional
+
+ // A variant candidate whose normal LOD is below this threshold will be filtered as 'germline_risk'
+ // This is a more stringent threshold than normalLodThresholdForVCF
+ double normalLodFilterThreshold = -Double.MAX_VALUE;
+ PerReadAlleleLikelihoodMap normalPRALM = null;
+ final PerAlleleCollection normalLods = PerAlleleCollection.createPerAltAlleleCollection();
+
+ // if normal bam is available, compute normal LOD
+ // TODO: this if statement should be a standalone method for computing normal LOD
+ // TODO: then we can do something like normalLodThreshold = hasNormal ? thisMethod() : Optional.empty()
+ if (hasNormal) {
+ normalPRALM = readAlleleLikelihoods.toPerReadAlleleLikelihoodMap(readAlleleLikelihoods.sampleIndex(matchedNormalSampleName));
+ filterPRALMForOverlappingReads(normalPRALM, mergedVC.getReference(), loc, true);
+ MuTect2.logReadInfo(DEBUG_READ_NAME, normalPRALM.getLikelihoodReadMap().keySet(), "Present after in Nomral PRALM filtering for overlapping reads");
+
+ final GenomeLoc eventGenomeLoc = genomeLocParser.createGenomeLoc(activeRegionWindow.getContig(), loc);
+ final Collection cosmicVC = tracker.getValues(MTAC.cosmicRod, eventGenomeLoc);
+ final Collection dbsnpVC = tracker.getValues(MTAC.dbsnp.dbsnp, eventGenomeLoc);
+ final boolean germlineAtRisk = !dbsnpVC.isEmpty() && cosmicVC.isEmpty();
+
+ normalLodFilterThreshold = germlineAtRisk ? MTAC.NORMAL_DBSNP_LOD_THRESHOLD : MTAC.NORMAL_LOD_THRESHOLD;
+
+ // compute normal LOD = LL(X|REF)/LL(X|ALT) where REF is the diploid HET with AF = 0.5
+ // note normal LOD is REF over ALT, the reciprocal of the tumor LOD
+ final PerAlleleCollection diploidHetAlleleFractions = PerAlleleCollection.createPerRefAndAltAlleleCollection();
+ for (final Allele allele : mergedVC.getAlternateAlleles()){
+ diploidHetAlleleFractions.setAlt(allele, 0.5);
+ }
+
+ final PerAlleleCollection normalGenotypeLLs = getHetGenotypeLogLikelihoods(mergedVC, normalPRALM, originalNormalReadQualities, diploidHetAlleleFractions);
+
+ for (final Allele altAllele : mergedVC.getAlternateAlleles()){
+ normalLods.setAlt(altAllele, normalGenotypeLLs.getRef() - normalGenotypeLLs.getAlt(altAllele));
+ }
+ }
+
+ int numPassingAlts = 0;
+ final Set allelesThatPassThreshold = new HashSet<>();
+ Allele alleleWithHighestTumorLOD = null;
+
+ for (final Allele altAllele : mergedVC.getAlternateAlleles()) {
+ final boolean passesTumorLodThreshold = tumorLods.getAlt(altAllele) >= MTAC.INITIAL_TUMOR_LOD_THRESHOLD;
+ final boolean passesNormalLodThreshold = hasNormal ? normalLods.getAlt(altAllele) >= MTAC.INITIAL_NORMAL_LOD_THRESHOLD : true;
+ if (passesTumorLodThreshold && passesNormalLodThreshold) {
+ numPassingAlts++;
+ allelesThatPassThreshold.add(altAllele);
+ if (alleleWithHighestTumorLOD == null || tumorLods.getAlt(altAllele) > tumorLods.getAlt(alleleWithHighestTumorLOD)){
+ alleleWithHighestTumorLOD = altAllele;
+ }
+ }
+ }
+
+ if (numPassingAlts == 0) {
+ continue;
+ }
+
+ final VariantContextBuilder callVcb = new VariantContextBuilder(mergedVC);
+ final int haplotypeCount = alleleMapper.get(alleleWithHighestTumorLOD).size();
+ callVcb.attribute(GATKVCFConstants.HAPLOTYPE_COUNT_KEY, haplotypeCount);
+ callVcb.attribute(GATKVCFConstants.TUMOR_LOD_KEY, tumorLods.getAlt(alleleWithHighestTumorLOD));
+
+ if (hasNormal) {
+ callVcb.attribute(GATKVCFConstants.NORMAL_LOD_KEY, normalLods.getAlt(alleleWithHighestTumorLOD));
+ if (normalLods.getAlt(alleleWithHighestTumorLOD) < normalLodFilterThreshold) {
+ callVcb.filter(GATKVCFConstants.GERMLINE_RISK_FILTER_NAME);
+ }
+ }
+
+ // TODO: this should be a separate method
+ // TODO: move code to MuTect2::calculateFilters()
+ if (MTAC.ENABLE_STRAND_ARTIFACT_FILTER && numPassingAlts == 1) {
+ final PerReadAlleleLikelihoodMap forwardPRALM = new PerReadAlleleLikelihoodMap();
+ final PerReadAlleleLikelihoodMap reversePRALM = new PerReadAlleleLikelihoodMap();
+ splitPRALMintoForwardAndReverseReads(tumorPRALM, forwardPRALM, reversePRALM);
+
+ MuTect2.logReadInfo(DEBUG_READ_NAME, tumorPRALM.getLikelihoodReadMap().keySet(), "Present in tumor PRALM after PRALM is split");
+ MuTect2.logReadInfo(DEBUG_READ_NAME, forwardPRALM.getLikelihoodReadMap().keySet(), "Present in forward PRALM after PRALM is split");
+ MuTect2.logReadInfo(DEBUG_READ_NAME, reversePRALM.getLikelihoodReadMap().keySet(), "Present in reverse PRALM after PRALM is split");
+
+ // TODO: build a new type for probability, likelihood, and log_likelihood. e.g. f_fwd :: probability[], tumorGLs_fwd :: likelihood[]
+ // TODO: don't want to call getHetGenotypeLogLikelihoods on more than one alternate alelle. May need to overload it to take a scalar f_fwd.
+ final PerAlleleCollection alleleFractionsForward = estimateAlleleFraction(mergedVC, forwardPRALM, true);
+ final PerAlleleCollection tumorGenotypeLLForward = getHetGenotypeLogLikelihoods(mergedVC, forwardPRALM, originalNormalReadQualities, alleleFractionsForward);
+
+ final PerAlleleCollection alleleFractionsReverse = estimateAlleleFraction(mergedVC, reversePRALM, true);
+ final PerAlleleCollection tumorGenotypeLLReverse = getHetGenotypeLogLikelihoods(mergedVC, reversePRALM, originalNormalReadQualities, alleleFractionsReverse);
+
+ final double tumorLod_fwd = tumorGenotypeLLForward.getAlt(alleleWithHighestTumorLOD) - tumorGenotypeLLForward.getRef();
+ final double tumorLod_rev = tumorGenotypeLLReverse.getAlt(alleleWithHighestTumorLOD) - tumorGenotypeLLReverse.getRef();
+
+ // Note that we use the observed combined (+ and -) allele fraction for power calculation in either direction
+ final double tumorSBpower_fwd = strandArtifactPowerCalculator.cachedPowerCalculation(forwardPRALM.getNumberOfStoredElements(), altAlleleFractions.getAlt(alleleWithHighestTumorLOD));
+ final double tumorSBpower_rev = strandArtifactPowerCalculator.cachedPowerCalculation(reversePRALM.getNumberOfStoredElements(), altAlleleFractions.getAlt(alleleWithHighestTumorLOD));
+
+ callVcb.attribute(GATKVCFConstants.TLOD_FWD_KEY, tumorLod_fwd);
+ callVcb.attribute(GATKVCFConstants.TLOD_REV_KEY, tumorLod_rev);
+ callVcb.attribute(GATKVCFConstants.TUMOR_SB_POWER_FWD_KEY, tumorSBpower_fwd);
+ callVcb.attribute(GATKVCFConstants.TUMOR_SB_POWER_REV_KEY, tumorSBpower_rev);
+
+ if ((tumorSBpower_fwd > MTAC.STRAND_ARTIFACT_POWER_THRESHOLD && tumorLod_fwd < MTAC.STRAND_ARTIFACT_LOD_THRESHOLD) ||
+ (tumorSBpower_rev > MTAC.STRAND_ARTIFACT_POWER_THRESHOLD && tumorLod_rev < MTAC.STRAND_ARTIFACT_LOD_THRESHOLD))
+ callVcb.filter(GATKVCFConstants.STRAND_ARTIFACT_FILTER_NAME);
+ }
+
+ // TODO: this probably belongs in M2::calculateFilters()
+ if (numPassingAlts > 1) {
+ callVcb.filter(GATKVCFConstants.TRIALLELIC_SITE_FILTER_NAME);
+ }
+
+ // build genotypes TODO: this part needs review and refactor
+ final List tumorAlleles = Arrays.asList(mergedVC.getReference(), alleleWithHighestTumorLOD);
+ // TODO: estimateAlleleFraction should not repeat counting allele depths
+ final PerAlleleCollection tumorAlleleDepths = getRefAltCount(mergedVC, tumorPRALM, false);
+ final int tumorRefAlleleDepth = tumorAlleleDepths.getRef();
+ final int tumorAltAlleleDepth = tumorAlleleDepths.getAlt(alleleWithHighestTumorLOD);
+ final Genotype tumorGenotype = new GenotypeBuilder(tumorSampleName, tumorAlleles)
+ .AD(new int[] { tumorRefAlleleDepth, tumorAltAlleleDepth })
+ .attribute(GATKVCFConstants.ALLELE_FRACTION_KEY, altAlleleFractions.getAlt(alleleWithHighestTumorLOD))
+ .make();
+
+ final List genotypes = new ArrayList<>();
+ genotypes.add(tumorGenotype);
+
+ // We assume that the genotype in the normal is 0/0
+ // TODO: is normal always homozygous reference?
+ final List homRefAllelesforNormalGenotype = Collections.nCopies(2, mergedVC.getReference());
+
+ // if we are calling with a normal, build the genotype for the sample to appear in vcf
+ if (hasNormal) {
+ final PerAlleleCollection normalAlleleDepths = getRefAltCount(mergedVC, normalPRALM, false);
+ final int normalRefAlleleDepth = normalAlleleDepths.getRef();
+ final int normalAltAlleleDepth = normalAlleleDepths.getAlt(alleleWithHighestTumorLOD);
+ final double normalAlleleFraction = (double) normalAltAlleleDepth / ( normalRefAlleleDepth + normalAltAlleleDepth);
+
+ final Genotype normalGenotype = new GenotypeBuilder(matchedNormalSampleName, homRefAllelesforNormalGenotype)
+ .AD(new int[] { normalRefAlleleDepth, normalAltAlleleDepth })
+ .attribute(GATKVCFConstants.ALLELE_FRACTION_KEY, normalAlleleFraction)
+ .make();
+ genotypes.add(normalGenotype);
+ }
+
+ final VariantContext call = new VariantContextBuilder(callVcb).alleles(tumorAlleles).genotypes(genotypes).make();
+
+ // how should we be making use of _perSampleFilteredReadList_?
+ readAlleleLikelihoods = prepareReadAlleleLikelihoodsForAnnotation(readLikelihoods, perSampleFilteredReadList,
+ genomeLocParser, false, alleleMapper, readAlleleLikelihoods, call);
+
+ final ReferenceContext referenceContext = new ReferenceContext(genomeLocParser, genomeLocParser.createGenomeLoc(mergedVC.getChr(), mergedVC.getStart(), mergedVC.getEnd()), refLoc, ref);
+ VariantContext annotatedCall = annotationEngine.annotateContextForActiveRegion(referenceContext, tracker, readAlleleLikelihoods, call, false);
+
+ if( call.getAlleles().size() != mergedVC.getAlleles().size() )
+ annotatedCall = GATKVariantContextUtils.reverseTrimAlleles(annotatedCall);
+
+ // maintain the set of all called haplotypes
+ call.getAlleles().stream().map(alleleMapper::get).filter(Objects::nonNull).forEach(calledHaplotypes::addAll);
+ returnCalls.add( annotatedCall );
}
// TODO: understand effect of enabling this for somatic calling...
- final List phasedCalls = doPhysicalPhasing ? phaseCalls(returnCalls, calledHaplotypes) : returnCalls;
- return new CalledHaplotypes(phasedCalls, calledHaplotypes);
- //return new CalledHaplotypes(returnCalls, calledHaplotypes);
+ final List outputCalls = doPhysicalPhasing ? phaseCalls(returnCalls, calledHaplotypes) : returnCalls;
+ return new CalledHaplotypes(outputCalls, calledHaplotypes);
}
- private void verifySamplePresence(String sampleName, List samples) {
- if (!samples.contains(sampleName)) {
- throw new IllegalArgumentException("Unable to find sample name "+sampleName+"in sample list of " + StringUtil.join(",", samples));
- }
- }
-
- /** Calculate the genotype likelihoods for variable allele fraction
+ /** Calculate the likelihoods of hom ref and each het genotype of the form ref/alt
*
* @param mergedVC input VC
* @param tumorPRALM read likelihoods
* @param originalNormalMQs original MQs, before boosting normals to avoid qual capping
- * @param afs allele fraction(s) for alternate allele(s)
+ * @param alleleFractions allele fraction(s) for alternate allele(s)
*
- * @return genotype likelihoods for homRef (index 0) and het for each alternate allele
- */
- private double[] getVariableGenotypeLikelihoods(final VariantContext mergedVC, final PerReadAlleleLikelihoodMap tumorPRALM,
- final Map originalNormalMQs, double[] afs) {
- double[] genotypeLikelihoods = new double[mergedVC.getNAlleles()];
- for(Map.Entry> e : tumorPRALM.getLikelihoodReadMap().entrySet()) {
- Map m = e.getValue();
- Double refLL = m.get(mergedVC.getReference());
- if (originalNormalMQs.get(e.getKey().getReadName()) != 0) {
- genotypeLikelihoods[0] += Math.log10(Math.pow(10, refLL));
-
- for (int altInd = 0; altInd < mergedVC.getNAlleles()-1; altInd++) {
- Double altLL = m.get(mergedVC.getAlternateAllele(altInd));
- genotypeLikelihoods[altInd+1] += Math.log10(Math.pow(10, refLL) * (1 - afs[altInd]) + Math.pow(10, altLL) * afs[altInd]);
- }
- }
+ * @return genotype likelihoods for homRef and het for each alternate allele
+ */
+ private PerAlleleCollection getHetGenotypeLogLikelihoods(final VariantContext mergedVC,
+ final PerReadAlleleLikelihoodMap tumorPRALM,
+ final Map originalNormalMQs,
+ final PerAlleleCollection alleleFractions) {
+ // make sure that alleles in alleleFraction are a subset of alleles in the variant context
+ if (! mergedVC.getAlternateAlleles().containsAll(alleleFractions.getAltAlleles()) ){
+ throw new IllegalArgumentException("alleleFractions has alleles that are not in the variant context");
}
- return genotypeLikelihoods;
+
+ final PerAlleleCollection genotypeLogLikelihoods = PerAlleleCollection.createPerRefAndAltAlleleCollection();
+ mergedVC.getAlleles().forEach(a -> genotypeLogLikelihoods.set(a, new MutableDouble(0)));
+
+ final Allele refAllele = mergedVC.getReference();
+ for(Map.Entry> readAlleleLikelihoodMap : tumorPRALM.getLikelihoodReadMap().entrySet()) {
+ final Map alleleLikelihoodMap = readAlleleLikelihoodMap.getValue();
+ if (originalNormalMQs.get(readAlleleLikelihoodMap.getKey().getReadName()) == 0) {
+ continue;
+ }
+
+ final double readRefLogLikelihood = alleleLikelihoodMap.get(refAllele);
+ genotypeLogLikelihoods.getRef().add(readRefLogLikelihood);
+
+ for (final Allele altAllele : alleleFractions.getAltAlleles()) {
+ final double readAltLogLikelihood = alleleLikelihoodMap.get(altAllele);
+ final double adjustedReadAltLL = Math.log10(
+ Math.pow(10, readRefLogLikelihood) * (1 - alleleFractions.getAlt(altAllele)) +
+ Math.pow(10, readAltLogLikelihood) * alleleFractions.getAlt(altAllele)
+ );
+ genotypeLogLikelihoods.get(altAllele).add(adjustedReadAltLL);
+ }
+
+ }
+
+ final PerAlleleCollection result = PerAlleleCollection.createPerRefAndAltAlleleCollection();
+ mergedVC.getAlleles().stream().forEach(a -> result.set(a,genotypeLogLikelihoods.get(a).toDouble()));
+
+ return result;
}
/**
* Find the allele fractions for each alternate allele
*
* @param vc input VC, for alleles
- * @param map read likelihoods
+ * @param pralm read likelihoods
* @return estimated AF for each alt
*/
// FIXME: calculate using the uncertainty rather than this cheap approach
- private double[] estimateAlleleFraction(VariantContext vc, PerReadAlleleLikelihoodMap map) {
- int[] counts = getRefAltCount(vc, map);
- int numAlts = vc.getNAlleles()-1;
- double[] afs = new double[numAlts];
- int refCount = counts[0];
- int altCount;
+ private PerAlleleCollection estimateAlleleFraction(final VariantContext vc,
+ final PerReadAlleleLikelihoodMap pralm,
+ final boolean oneStrandOnly) {
+ final PerAlleleCollection alleleCounts = getRefAltCount(vc, pralm, oneStrandOnly);
+ final PerAlleleCollection alleleFractions = PerAlleleCollection.createPerAltAlleleCollection();
- for(int altInd = 0; altInd < numAlts; altInd++) {
- altCount = counts[altInd+1];
- afs[altInd] = (double) altCount / ((double) refCount + (double) altCount);
- //logger.info("Counted " + refCount + " ref and " + altCount + " alt " );
+ final int refCount = alleleCounts.getRef();
+ for ( final Allele altAllele : vc.getAlternateAlleles() ) {
+ final int altCount = alleleCounts.getAlt(altAllele);
+ double alleleFraction = (double) altCount / (refCount + altCount);
+ // weird case, but I've seen it happen in one strand cases
+ if (refCount == 0 && altCount == refCount ) {
+ alleleFraction = 0;
+ }
+ alleleFractions.setAlt(altAllele, alleleFraction);
+ // logger.info("Counted " + refCount + " ref and " + altCount + " alt " );
}
- return afs;
+ return alleleFractions;
}
/**
- * Evalutate the most likely allele for each read, if it is in fact informative
+ * Go through the PRALM and tally the most likely allele in each read. Only count informative reads.
*
- * @param mergedVC input VC, for alleles
- * @param afMap read likelihoods
+ * @param vc input VC, for alleles
+ * @param pralm read likelihoods
* @return an array giving the read counts for the ref and each alt allele
*/
- // TODO: ensure there are only two alleles in the VC
- private int[] getRefAltCount(VariantContext mergedVC, PerReadAlleleLikelihoodMap afMap) {
- int counts[] = new int[mergedVC.getNAlleles()];
- int REF = 0;
+ private PerAlleleCollection getRefAltCount(final VariantContext vc,
+ final PerReadAlleleLikelihoodMap pralm,
+ final boolean oneStrandOnly) {
+ // Check that the alleles in Variant Context are in PRALM
+ // Skip the check for strand-conscious PRALM; + reads may not have alleles in - reads, for example.
+ final Set vcAlleles = new HashSet<>(vc.getAlleles());
+ if ( ! oneStrandOnly && ! pralm.getAllelesSet().containsAll( vcAlleles ) ) {
+ StringBuilder message = new StringBuilder();
+ message.append("At Locus chr" + vc.getContig() + ":" + vc.getStart() + ", we detected that variant context had alleles that not in PRALM. ");
+ message.append("VC alleles = " + vcAlleles + ", PRALM alleles = " + pralm.getAllelesSet());
+ logger.warn(message);
+ }
- for(Map.Entry> e : afMap.getLikelihoodReadMap().entrySet()) {
- Map m = e.getValue();
- Double rl = m.get(mergedVC.getReference());
- for(int altInd=0; altInd alleleCounts = PerAlleleCollection.createPerRefAndAltAlleleCollection();
+ vcAlleles.stream().forEach(a -> alleleCounts.set(a, new MutableInt(0)));
+
+ for (final Map.Entry> readAlleleLikelihoodMap : pralm.getLikelihoodReadMap().entrySet()) {
+ final GATKSAMRecord read = readAlleleLikelihoodMap.getKey();
+ final Map alleleLikelihoodMap = readAlleleLikelihoodMap.getValue();
+ final MostLikelyAllele mostLikelyAllele = PerReadAlleleLikelihoodMap.getMostLikelyAllele(alleleLikelihoodMap, vcAlleles);
+
+ if (read.getMappingQuality() > 0 && mostLikelyAllele.isInformative()) {
+ alleleCounts.get(mostLikelyAllele.getMostLikelyAllele()).increment();
}
-// if (al >= rl) logger.info("Alt found in " + e.getKey().getReadName());
}
- return counts;
- }
+ final PerAlleleCollection result = PerAlleleCollection.createPerRefAndAltAlleleCollection();
+ vc.getAlleles().stream().forEach(a -> result.set(a, alleleCounts.get(a).toInteger()));
+
+ return(result);
+ }
private void logM2Debug(String s) {
if (MTAC.M2_DEBUG) {
@@ -465,26 +511,16 @@ public class SomaticGenotypingEngine extends HaplotypeCallerGenotypingEngine {
}
}
- // would have used org.broadinstitute.sting.utils.genotyper.PerReadAlleleLikelihoodMap.getMostLikelyAllele but we have this case where
- // there is a read that doesn't overlap the variant site, and thus supports both alleles equally.
- private boolean arePairHMMLikelihoodsInformative(double l1, double l2) {
- // TODO: should this be parameterized, or simply encoded
- double EPSILON = 0.1;
- return (Math.abs(l1 - l2) >= EPSILON);
- }
-
- private void filterPRALMForOverlappingReads(PerReadAlleleLikelihoodMap pralm, Allele ref, int location, boolean retainMismatches) {
-
- Map> m = pralm.getLikelihoodReadMap();
-
+ private void filterPRALMForOverlappingReads(final PerReadAlleleLikelihoodMap pralm, final Allele ref, final int location, final boolean retainMismatches) {
+ final Map> m = pralm.getLikelihoodReadMap();
// iterate through the reads, if the name has been seen before we have overlapping (potentially) fragments, so handle them
- Map nameToRead = new HashMap<>();
- Set readsToKeep = new HashSet<>();
+ final Map nameToRead = new HashMap<>();
+ final Set readsToKeep = new HashSet<>();
- for(GATKSAMRecord rec : m.keySet()) {
+ for(final GATKSAMRecord rec : m.keySet()) {
// if we haven't seen it... just record the name and add it to the list of reads to keep
- GATKSAMRecord existing = nameToRead.get(rec.getReadName());
+ final GATKSAMRecord existing = nameToRead.get(rec.getReadName());
if (existing == null) {
nameToRead.put(rec.getReadName(), rec);
readsToKeep.add(rec);
@@ -496,11 +532,11 @@ public class SomaticGenotypingEngine extends HaplotypeCallerGenotypingEngine {
// TODO: CHECK IF THE READS BOTH OVERLAP THE POSITION!!!!
if ( ReadUtils.isInsideRead(existing, location) && ReadUtils.isInsideRead(rec, location) ) {
- MostLikelyAllele existingMLA = pralm.getMostLikelyAllele(pralm.getLikelihoodReadMap().get(existing));
- Allele existingAllele = existingMLA.getMostLikelyAllele();
+ final MostLikelyAllele existingMLA = PerReadAlleleLikelihoodMap.getMostLikelyAllele(pralm.getLikelihoodReadMap().get(existing));
+ final Allele existingAllele = existingMLA.getMostLikelyAllele();
- MostLikelyAllele recMLA = pralm.getMostLikelyAllele(pralm.getLikelihoodReadMap().get(rec));
- Allele recAllele = recMLA.getMostLikelyAllele();
+ final MostLikelyAllele recMLA = PerReadAlleleLikelihoodMap.getMostLikelyAllele(pralm.getLikelihoodReadMap().get(rec));
+ final Allele recAllele = recMLA.getMostLikelyAllele();
// if the reads disagree at this position...
if (!existingAllele.equals(recAllele)) {
@@ -545,15 +581,20 @@ public class SomaticGenotypingEngine extends HaplotypeCallerGenotypingEngine {
}
}
- // Move to utility class so we can use one shared with HaplotypeCallerGenotypingEngine
- private VariantContext addNonRefSymbolicAllele(final VariantContext mergedVC) {
- final VariantContextBuilder vcb = new VariantContextBuilder(mergedVC);
- final List originalList = mergedVC.getAlleles();
- final List alleleList = new ArrayList<>(originalList.size() + 1);
- alleleList.addAll(mergedVC.getAlleles());
- alleleList.add(GATKVCFConstants.NON_REF_SYMBOLIC_ALLELE);
- vcb.alleles(alleleList);
- return vcb.make();
- }
+ private void splitPRALMintoForwardAndReverseReads(final PerReadAlleleLikelihoodMap originalPRALM, final PerReadAlleleLikelihoodMap forwardPRALM, final PerReadAlleleLikelihoodMap reversePRALM) {
+ final Map> origReadAlleleLikelihoodMap = originalPRALM.getLikelihoodReadMap();
+ for (final GATKSAMRecord read : origReadAlleleLikelihoodMap.keySet()) {
+ if (read.isStrandless())
+ continue;
+ for (final Map.Entry alleleLikelihoodMap : origReadAlleleLikelihoodMap.get(read).entrySet()) {
+ final Allele allele = alleleLikelihoodMap.getKey();
+ final Double likelihood = alleleLikelihoodMap.getValue();
+ if (read.getReadNegativeStrandFlag())
+ reversePRALM.add(read, allele, likelihood);
+ else
+ forwardPRALM.add(read, allele, likelihood);
+ }
+ }
+ }
}
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/TumorPowerCalculator.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/TumorPowerCalculator.java
new file mode 100644
index 000000000..8c09160ff
--- /dev/null
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/cancer/m2/TumorPowerCalculator.java
@@ -0,0 +1,201 @@
+/*
+* By downloading the PROGRAM you agree to the following terms of use:
+*
+* BROAD INSTITUTE
+* SOFTWARE LICENSE AGREEMENT
+* FOR ACADEMIC NON-COMMERCIAL RESEARCH PURPOSES ONLY
+*
+* This Agreement is made between the Broad Institute, Inc. with a principal address at 415 Main Street, Cambridge, MA 02142 ("BROAD") and the LICENSEE and is effective at the date the downloading is completed ("EFFECTIVE DATE").
+*
+* WHEREAS, LICENSEE desires to license the PROGRAM, as defined hereinafter, and BROAD wishes to have this PROGRAM utilized in the public interest, subject only to the royalty-free, nonexclusive, nontransferable license rights of the United States Government pursuant to 48 CFR 52.227-14; and
+* WHEREAS, LICENSEE desires to license the PROGRAM and BROAD desires to grant a license on the following terms and conditions.
+* NOW, THEREFORE, in consideration of the promises and covenants made herein, the parties hereto agree as follows:
+*
+* 1. DEFINITIONS
+* 1.1 PROGRAM shall mean copyright in the object code and source code known as GATK3 and related documentation, if any, as they exist on the EFFECTIVE DATE and can be downloaded from http://www.broadinstitute.org/gatk on the EFFECTIVE DATE.
+*
+* 2. LICENSE
+* 2.1 Grant. Subject to the terms of this Agreement, BROAD hereby grants to LICENSEE, solely for academic non-commercial research purposes, a non-exclusive, non-transferable license to: (a) download, execute and display the PROGRAM and (b) create bug fixes and modify the PROGRAM. LICENSEE hereby automatically grants to BROAD a non-exclusive, royalty-free, irrevocable license to any LICENSEE bug fixes or modifications to the PROGRAM with unlimited rights to sublicense and/or distribute. LICENSEE agrees to provide any such modifications and bug fixes to BROAD promptly upon their creation.
+* The LICENSEE may apply the PROGRAM in a pipeline to data owned by users other than the LICENSEE and provide these users the results of the PROGRAM provided LICENSEE does so for academic non-commercial purposes only. For clarification purposes, academic sponsored research is not a commercial use under the terms of this Agreement.
+* 2.2 No Sublicensing or Additional Rights. LICENSEE shall not sublicense or distribute the PROGRAM, in whole or in part, without prior written permission from BROAD. LICENSEE shall ensure that all of its users agree to the terms of this Agreement. LICENSEE further agrees that it shall not put the PROGRAM on a network, server, or other similar technology that may be accessed by anyone other than the LICENSEE and its employees and users who have agreed to the terms of this agreement.
+* 2.3 License Limitations. Nothing in this Agreement shall be construed to confer any rights upon LICENSEE by implication, estoppel, or otherwise to any computer software, trademark, intellectual property, or patent rights of BROAD, or of any other entity, except as expressly granted herein. LICENSEE agrees that the PROGRAM, in whole or part, shall not be used for any commercial purpose, including without limitation, as the basis of a commercial software or hardware product or to provide services. LICENSEE further agrees that the PROGRAM shall not be copied or otherwise adapted in order to circumvent the need for obtaining a license for use of the PROGRAM.
+*
+* 3. PHONE-HOME FEATURE
+* LICENSEE expressly acknowledges that the PROGRAM contains an embedded automatic reporting system ("PHONE-HOME") which is enabled by default upon download. Unless LICENSEE requests disablement of PHONE-HOME, LICENSEE agrees that BROAD may collect limited information transmitted by PHONE-HOME regarding LICENSEE and its use of the PROGRAM. Such information shall include LICENSEE'S user identification, version number of the PROGRAM and tools being run, mode of analysis employed, and any error reports generated during run-time. Collection of such information is used by BROAD solely to monitor usage rates, fulfill reporting requirements to BROAD funding agencies, drive improvements to the PROGRAM, and facilitate adjustments to PROGRAM-related documentation.
+*
+* 4. OWNERSHIP OF INTELLECTUAL PROPERTY
+* LICENSEE acknowledges that title to the PROGRAM shall remain with BROAD. The PROGRAM is marked with the following BROAD copyright notice and notice of attribution to contributors. LICENSEE shall retain such notice on all copies. LICENSEE agrees to include appropriate attribution if any results obtained from use of the PROGRAM are included in any publication.
+* Copyright 2012-2016 Broad Institute, Inc.
+* Notice of attribution: The GATK3 program was made available through the generosity of Medical and Population Genetics program at the Broad Institute, Inc.
+* LICENSEE shall not use any trademark or trade name of BROAD, or any variation, adaptation, or abbreviation, of such marks or trade names, or any names of officers, faculty, students, employees, or agents of BROAD except as states above for attribution purposes.
+*
+* 5. INDEMNIFICATION
+* LICENSEE shall indemnify, defend, and hold harmless BROAD, and their respective officers, faculty, students, employees, associated investigators and agents, and their respective successors, heirs and assigns, (Indemnitees), against any liability, damage, loss, or expense (including reasonable attorneys fees and expenses) incurred by or imposed upon any of the Indemnitees in connection with any claims, suits, actions, demands or judgments arising out of any theory of liability (including, without limitation, actions in the form of tort, warranty, or strict liability and regardless of whether such action has any factual basis) pursuant to any right or license granted under this Agreement.
+*
+* 6. NO REPRESENTATIONS OR WARRANTIES
+* THE PROGRAM IS DELIVERED AS IS. BROAD MAKES NO REPRESENTATIONS OR WARRANTIES OF ANY KIND CONCERNING THE PROGRAM OR THE COPYRIGHT, EXPRESS OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, NONINFRINGEMENT, OR THE ABSENCE OF LATENT OR OTHER DEFECTS, WHETHER OR NOT DISCOVERABLE. BROAD EXTENDS NO WARRANTIES OF ANY KIND AS TO PROGRAM CONFORMITY WITH WHATEVER USER MANUALS OR OTHER LITERATURE MAY BE ISSUED FROM TIME TO TIME.
+* IN NO EVENT SHALL BROAD OR ITS RESPECTIVE DIRECTORS, OFFICERS, EMPLOYEES, AFFILIATED INVESTIGATORS AND AFFILIATES BE LIABLE FOR INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, INCLUDING, WITHOUT LIMITATION, ECONOMIC DAMAGES OR INJURY TO PROPERTY AND LOST PROFITS, REGARDLESS OF WHETHER BROAD SHALL BE ADVISED, SHALL HAVE OTHER REASON TO KNOW, OR IN FACT SHALL KNOW OF THE POSSIBILITY OF THE FOREGOING.
+*
+* 7. ASSIGNMENT
+* This Agreement is personal to LICENSEE and any rights or obligations assigned by LICENSEE without the prior written consent of BROAD shall be null and void.
+*
+* 8. MISCELLANEOUS
+* 8.1 Export Control. LICENSEE gives assurance that it will comply with all United States export control laws and regulations controlling the export of the PROGRAM, including, without limitation, all Export Administration Regulations of the United States Department of Commerce. Among other things, these laws and regulations prohibit, or require a license for, the export of certain types of software to specified countries.
+* 8.2 Termination. LICENSEE shall have the right to terminate this Agreement for any reason upon prior written notice to BROAD. If LICENSEE breaches any provision hereunder, and fails to cure such breach within thirty (30) days, BROAD may terminate this Agreement immediately. Upon termination, LICENSEE shall provide BROAD with written assurance that the original and all copies of the PROGRAM have been destroyed, except that, upon prior written authorization from BROAD, LICENSEE may retain a copy for archive purposes.
+* 8.3 Survival. The following provisions shall survive the expiration or termination of this Agreement: Articles 1, 3, 4, 5 and Sections 2.2, 2.3, 7.3, and 7.4.
+* 8.4 Notice. Any notices under this Agreement shall be in writing, shall specifically refer to this Agreement, and shall be sent by hand, recognized national overnight courier, confirmed facsimile transmission, confirmed electronic mail, or registered or certified mail, postage prepaid, return receipt requested. All notices under this Agreement shall be deemed effective upon receipt.
+* 8.5 Amendment and Waiver; Entire Agreement. This Agreement may be amended, supplemented, or otherwise modified only by means of a written instrument signed by all parties. Any waiver of any rights or failure to act in a specific instance shall relate only to such instance and shall not be construed as an agreement to waive any rights or fail to act in any other instance, whether or not similar. This Agreement constitutes the entire agreement among the parties with respect to its subject matter and supersedes prior agreements or understandings between the parties relating to its subject matter.
+* 8.6 Binding Effect; Headings. This Agreement shall be binding upon and inure to the benefit of the parties and their respective permitted successors and assigns. All headings are for convenience only and shall not affect the meaning of any provision of this Agreement.
+* 8.7 Governing Law. This Agreement shall be construed, governed, interpreted and applied in accordance with the internal laws of the Commonwealth of Massachusetts, U.S.A., without regard to conflict of laws principles.
+*/
+
+package org.broadinstitute.gatk.tools.walkers.cancer.m2;
+
+import org.apache.commons.math.MathException;
+import org.apache.commons.math.distribution.BinomialDistribution;
+import org.apache.commons.math.distribution.BinomialDistributionImpl;
+import org.apache.commons.math3.util.Pair;
+import org.broadinstitute.gatk.utils.exceptions.GATKException;
+
+import java.util.Arrays;
+import java.util.HashMap;
+import java.util.OptionalInt;
+import java.util.stream.IntStream;
+
+/**
+ * We store a memo to avoid repeated computation of statistical power to detect a variant.
+ * The key of the memo is a pair of numbers: number of reads and estimated allele fraction
+ */
+public class TumorPowerCalculator {
+ private final double errorProbability;
+ private final double tumorLODThreshold;
+ private final double contamination;
+ private final boolean enableSmoothing;
+
+ private final HashMap cache = new HashMap();
+
+ public TumorPowerCalculator(double errorProbability, double constantLodThreshold, double contamination) {
+ this(errorProbability, constantLodThreshold, contamination, true);
+ }
+
+ public TumorPowerCalculator(double errorProbability, double tumorLODThreshold, double contamination, boolean enableSmoothing) {
+ this.errorProbability = errorProbability;
+ this.tumorLODThreshold = tumorLODThreshold;
+ this.contamination = contamination;
+ this.enableSmoothing = enableSmoothing;
+ }
+
+ /**
+ * A helper class that acts as the key to the memo of pre-computed power
+ *
+ * TODO: Not ideal to use double as a key. Refactor such that we use as keys numAlts and numReads, which are integers. Then calculate numAlts/numReads when we need allele fraction.
+ *
+ */
+ private static class PowerCacheKey extends Pair {
+ private final Double alleleFraction;
+ private final Integer numReads;
+
+ public PowerCacheKey(final int numReads, final double alleleFraction) {
+ super(numReads, alleleFraction);
+ this.alleleFraction = alleleFraction;
+ this.numReads = numReads;
+ }
+
+ private boolean closeEnough(final double x, final double y, final double epsilon){
+ return(Math.abs(x - y) < epsilon);
+ }
+
+ @Override
+ public boolean equals(Object o) {
+ if (this == o) return true;
+ if (o == null || getClass() != o.getClass()) return false;
+
+ PowerCacheKey that = (PowerCacheKey) o;
+ return (closeEnough(alleleFraction, that.alleleFraction, 0.001) && numReads != that.numReads);
+ }
+
+ @Override
+ public int hashCode() {
+ int result;
+ long temp;
+ result = numReads;
+ temp = alleleFraction != +0.0d ? Double.doubleToLongBits(alleleFraction) : 0L;
+ result = 31 * result + (int) (temp ^ (temp >>> 32));
+ return result;
+ }
+ }
+
+ /**
+ *
+ * @param numReads total number of reads, REF and ALT combined, in + or - strand
+ * @param alleleFraction the true allele fraction estimated as the combined allele fraction from + and - reads
+ * @return probability of correctly calling the variant (i.e. power) given the above estimated allele fraction and number of reads.
+ * we compute power separately for each strand (+ and -)
+ * @throws MathException
+ *
+ */
+ public double cachedPowerCalculation(final int numReads, final double alleleFraction) {
+ PowerCacheKey key = new PowerCacheKey(numReads, alleleFraction);
+ // we first look up if power for given number of read and allele fraction has already been computed and stored in the cache.
+ // if not we compute it and store it in teh cache.
+ Double power = cache.get(key);
+ if (power == null) {
+ try {
+ power = calculatePower(numReads, alleleFraction);
+ } catch (final Exception ex) {
+ throw new GATKException("Power calculation failed", ex);
+ }
+ cache.put(key, power);
+ }
+ return power;
+ }
+
+ /* helper function for calculateTumorLod */
+ private double calculateLogLikelihood(final int numReads, final int numAlts, final double alleleFraction) {
+ return((numReads-numAlts) * Math.log10( alleleFraction * errorProbability + (1 - alleleFraction)*(1 - errorProbability) ) +
+ numAlts * Math.log10(alleleFraction * (1 - errorProbability) + (1 - alleleFraction) * errorProbability));
+ }
+
+ private double calculateTumorLod(final int numReads, final int numAlts) {
+ final double alleleFraction = (double) numAlts / (double) numReads;
+ final double altLikelihod = calculateLogLikelihood(numReads, numAlts, alleleFraction);
+ final double refLikelihood = calculateLogLikelihood(numReads, numAlts, contamination);
+ return(altLikelihod - refLikelihood);
+}
+
+ private double calculatePower(final int numReads, final double alleleFraction) throws MathException {
+ if (numReads==0) return 0;
+
+ // TODO: add the factor of 1/3
+ final double probAltRead = alleleFraction*(1 - errorProbability) + (1/3)*(1 - alleleFraction) * errorProbability;
+ final BinomialDistribution binom = new BinomialDistributionImpl(numReads, probAltRead);
+ final double[] binomialProbabilities = IntStream.range(0, numReads + 1).mapToDouble(binom::probability).toArray();
+
+ // find the smallest number of ALT reads k such that tumorLOD(k) > tumorLODThreshold
+ final OptionalInt smallestKAboveLogThreshold = IntStream.range(0, numReads + 1)
+ .filter(k -> calculateTumorLod(numReads, k) > tumorLODThreshold)
+ .findFirst();
+
+ if (! smallestKAboveLogThreshold.isPresent()){
+ return 0;
+ }
+
+ if (smallestKAboveLogThreshold.getAsInt() <= 0){
+ throw new IllegalStateException("smallest k that meets the tumor LOD threshold is less than or equal to 0");
+ }
+
+ double power = Arrays.stream(binomialProbabilities, smallestKAboveLogThreshold.getAsInt(), binomialProbabilities.length).sum();
+
+ // here we correct for the fact that the exact lod threshold is likely somewhere between
+ // the k and k-1 bin, so we prorate the power from that bin
+ if ( enableSmoothing ){
+ final double tumorLODAtK = calculateTumorLod(numReads, smallestKAboveLogThreshold.getAsInt());
+ final double tumorLODAtKMinusOne = calculateTumorLod(numReads, smallestKAboveLogThreshold.getAsInt()-1);
+ final double weight = 1 - (tumorLODThreshold - tumorLODAtKMinusOne ) / (tumorLODAtK - tumorLODAtKMinusOne);
+ power += weight * binomialProbabilities[smallestKAboveLogThreshold.getAsInt() - 1];
+ }
+
+ return(power);
+ }
+
+}
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeAlleleCounts.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeAlleleCounts.java
index f1edc7eae..daf7e653c 100644
--- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeAlleleCounts.java
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeAlleleCounts.java
@@ -52,6 +52,7 @@
package org.broadinstitute.gatk.tools.walkers.genotyper;
import htsjdk.variant.variantcontext.Allele;
+import org.broadinstitute.gatk.utils.IndexRange;
import org.broadinstitute.gatk.utils.MathUtils;
import java.util.Arrays;
@@ -108,7 +109,18 @@ import java.util.List;
*/
public class GenotypeAlleleCounts implements Comparable, Cloneable {
- private double log10CombinationCount;
+ private static final double UNCOMPUTED_LOG_10_COMBINATION_COUNT = -1;
+
+ /**
+ * The log10 number of phased genotypes corresponding to this unphased genotype. For example,
+ * [0, 1, 1, 1] = AB: log10(2)
+ * [0, 2] = AA: log10(1)
+ * [0, 1, 1, 1, 2, 1] = ABC: log10(6)
+ * [0, 2, 1, 2] = AABB: log10(4!/(2!2!))
+ * This is evaluated lazily i.e. it is initialized to {@link GenotypeAlleleCounts::UNCOMPUTED_LOG_10_COMBINATION_COUNT}
+ * and only calculated if its getter is invoked.
+ */
+ private double log10CombinationCount = UNCOMPUTED_LOG_10_COMBINATION_COUNT;
/**
* The ploidy of the genotype.
@@ -156,38 +168,30 @@ public class GenotypeAlleleCounts implements Comparable, C
* @param index the genotype index.
*/
private GenotypeAlleleCounts(final int ploidy, final int index, final int... sortedAlleleCounts) {
+ this(ploidy, index, sortedAlleleCounts, sortedAlleleCounts.length >> 1);
+ }
+
+ private GenotypeAlleleCounts(final int ploidy, final int index, final int[] sortedAlleleCounts, final int distinctAlleleCount){
this.ploidy = ploidy;
- this.sortedAlleleCounts = sortedAlleleCounts;
- distinctAlleleCount = sortedAlleleCounts.length >> 1;
- log10CombinationCount = -1;
this.index = index;
+ this.sortedAlleleCounts = sortedAlleleCounts;
+ this.distinctAlleleCount = distinctAlleleCount;
}
/**
- * Returns the log10 of the number of possible allele combinations that would give raise to this allele count.
- * @return 0 or less.
+ * Gets the log10 combination count, computing it if uninitialized. Note that the invoked MathUtils method uses fast cached
+ * log10 values of integers for any reasonable ploidy.
+ *
+ * This method should be invoked on instances of {@link GenotypeAlleleCounts} cached in {@link GenotypeLikelihoodCalculators::genotypeTableByPloidy}.
+ * Such usage allows the result of this computation to be cached once for an entire run of HaplotypeCaller.
+ * @return
*/
public double log10CombinationCount() {
- if (log10CombinationCount == -1)
- return log10CombinationCount = calculateLog10CombinationCount();
- else
- return log10CombinationCount;
- }
-
- /**
- * Calculates log10 combination count.
- *
- * @return 0 or less.
- */
- private double calculateLog10CombinationCount() {
- if (ploidy <= 1)
- return 0;
- else {
- final int[] counts = new int[distinctAlleleCount];
- for (int i = 0, j = 1; i < distinctAlleleCount; i++, j+=2)
- counts[i] = sortedAlleleCounts[j];
- return MathUtils.log10MultinomialCoefficient(ploidy, counts);
+ if (log10CombinationCount == UNCOMPUTED_LOG_10_COMBINATION_COUNT) {
+ log10CombinationCount = MathUtils.log10Factorial(ploidy)
+ - new IndexRange(0, distinctAlleleCount).sum(n -> MathUtils.log10Factorial(sortedAlleleCounts[2*n+1]));
}
+ return log10CombinationCount;
}
/**
@@ -785,4 +789,22 @@ public class GenotypeAlleleCounts implements Comparable, C
}
return result;
}
+
+ @FunctionalInterface
+ public interface IntBiConsumer {
+ void accept(final int alleleIndex, final int alleleCount);
+ }
+
+ @FunctionalInterface
+ public interface IntToDoubleBiFunction {
+ double apply(final int alleleIndex, final int alleleCount);
+ }
+
+ public void forEachAlleleIndexAndCount(final IntBiConsumer action) {
+ new IndexRange(0, distinctAlleleCount).forEach(n -> action.accept(sortedAlleleCounts[2*n], sortedAlleleCounts[2*n+1]));
+ }
+
+ public double sumOverAlleleIndicesAndCounts(final IntToDoubleBiFunction func) {
+ return new IndexRange(0, distinctAlleleCount).sum(n -> func.apply(sortedAlleleCounts[2*n], sortedAlleleCounts[2*n+1]));
+ }
}
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculator.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculator.java
index 1a3cdf2b8..aead2ab15 100644
--- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculator.java
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculator.java
@@ -163,7 +163,7 @@ public class GenotypeLikelihoodCalculator {
* This is in fact a shallow copy if {@link GenotypeLikelihoodCalculators#ploidyLog10}
and is not meant to be modified by
* this class.
*/
- private final double[] log10;
+ private final double[] ploidyLog10;
/**
* Buffer field use as a temporal container for sorted allele counts when calculating the likelihood of a
@@ -202,24 +202,22 @@ public class GenotypeLikelihoodCalculator {
* Creates a new calculator providing its ploidy and number of genotyping alleles.
*/
protected GenotypeLikelihoodCalculator(final int ploidy, final int alleleCount,
- final int[][] alleleFirstGenotypeOffsetByPloidy,
- final GenotypeAlleleCounts[][] genotypeTableByPloidy,
- final double[] ploidyLog10) {
+ final int[][] alleleFirstGenotypeOffsetByPloidy,
+ final GenotypeAlleleCounts[][] genotypeTableByPloidy,
+ final double[] ploidyLog10) {
+
this.alleleFirstGenotypeOffsetByPloidy = alleleFirstGenotypeOffsetByPloidy;
genotypeAlleleCounts = genotypeTableByPloidy[ploidy];
this.alleleCount = alleleCount;
this.ploidy = ploidy;
genotypeCount = this.alleleFirstGenotypeOffsetByPloidy[ploidy][alleleCount];
- if (genotypeCount == GenotypeLikelihoodCalculators.GENOTYPE_COUNT_OVERFLOW)
- throw new IllegalArgumentException(
- String.format("the combination of ploidy (%s) and number of alleles (%s) results in a very large number of genotypes (> %s). You need to limit ploidy or the number of alternative alleles to analyze this locus",
- ploidy,alleleCount,Integer.MAX_VALUE));
+
alleleHeap = new IntMaxHeap(ploidy);
readLikelihoodsByGenotypeIndex = new double[genotypeCount][];
- log10 = ploidyLog10;
+ this.ploidyLog10 = ploidyLog10;
// The number of possible components is limited by distinct allele count and ploidy.
maximumDistinctAllelesInGenotype = Math.min(ploidy, alleleCount);
- genotypeAllelesAndCounts = new int[maximumDistinctAllelesInGenotype << 1];
+ genotypeAllelesAndCounts = new int[maximumDistinctAllelesInGenotype*2];
}
/**
@@ -349,7 +347,7 @@ public class GenotypeLikelihoodCalculator {
*/
private double[] genotypeLikelihoods(final double[][] readLikelihoodsByGenotypeIndex, final int readCount) {
final double[] result = new double[genotypeCount];
- final double denominator = readCount * log10[ploidy]; // instead of dividing each read likelihood by ploidy
+ final double denominator = readCount * ploidyLog10[ploidy]; // instead of dividing each read likelihood by ploidy
// ( so subtract log10(ploidy) ) we multiply them all and the divide by ploidy^readCount (so substract readCount * log10(ploidy) )
for (int g = 0; g < genotypeCount; g++) {
final double[] likelihoodsByRead = readLikelihoodsByGenotypeIndex[g];
@@ -464,7 +462,9 @@ public class GenotypeLikelihoodCalculator {
* exactly one allele present in the genotype.
*/
private void singleComponentGenotypeLikelihoodByRead(final GenotypeAlleleCounts genotypeAlleleCounts,
- final double[] likelihoodByRead, final double[] readLikelihoodComponentsByAlleleCount, final int readCount) {
+ final double[] likelihoodByRead,
+ final double[] readLikelihoodComponentsByAlleleCount,
+ final int readCount) {
final int allele = genotypeAlleleCounts.alleleIndexAt(0);
// the count of the only component must be = ploidy.
int offset = (allele * (ploidy + 1) + ploidy) * readCount;
@@ -493,7 +493,7 @@ public class GenotypeLikelihoodCalculator {
// p = 2 because the frequency == 1 we already have it.
for (int frequency = 2, destinationOffset = frequency1Offset + readCount; frequency <= ploidy; frequency++) {
- final double log10frequency = log10[frequency];
+ final double log10frequency = ploidyLog10[frequency];
for (int r = 0, sourceOffset = frequency1Offset; r < readCount; r++)
readAlleleLikelihoodByAlleleCount[destinationOffset++] =
readAlleleLikelihoodByAlleleCount[sourceOffset++] + log10frequency;
@@ -620,7 +620,11 @@ public class GenotypeLikelihoodCalculator {
* @param destination where to store the new genotype index mapping to old.
* @param sortedAlleleCountsBuffer a buffer to re-use to get the genotype-allele-count's sorted allele counts.
*/
- private void genotypeIndexMapPerGenotypeIndex(final int newGenotypeIndex, final GenotypeAlleleCounts alleleCounts, final int[] oldToNewAlleleIndexMap, final int[] destination, final int[] sortedAlleleCountsBuffer) {
+ private void genotypeIndexMapPerGenotypeIndex(final int newGenotypeIndex,
+ final GenotypeAlleleCounts alleleCounts,
+ final int[] oldToNewAlleleIndexMap,
+ final int[] destination,
+ final int[] sortedAlleleCountsBuffer) {
final int distinctAlleleCount = alleleCounts.distinctAlleleCount();
alleleCounts.copyAlleleCounts(sortedAlleleCountsBuffer,0);
for (int j = 0, jj = 0; j < distinctAlleleCount; j++) {
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculators.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculators.java
index 83b280bff..da0ff0975 100644
--- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculators.java
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculators.java
@@ -51,7 +51,11 @@
package org.broadinstitute.gatk.tools.walkers.genotyper;
+import org.broadinstitute.gatk.utils.MathUtils;
+import org.broadinstitute.gatk.utils.exceptions.GATKException;
+
import java.util.Arrays;
+import java.util.stream.IntStream;
/**
* Genotype likelihood calculator utility.
@@ -116,7 +120,10 @@ public class GenotypeLikelihoodCalculators {
private volatile static GenotypeAlleleCounts[][] genotypeTableByPloidy =
buildGenotypeAlleleCountsTable(maximumPloidy,maximumAllele,alleleFirstGenotypeOffsetByPloidy);
-
+ /**
+ * Cached log10 values for the first integer up to the maximum ploidy requested thus far.
+ */
+ private volatile static double[] ploidyLog10 = IntStream.range(0, maximumPloidy + 1).mapToDouble(Math::log10).toArray();
/**
* Build the table with the genotype offsets based on ploidy and the maximum allele index with representation
@@ -291,40 +298,29 @@ public class GenotypeLikelihoodCalculators {
return result;
}
- /**
- * Cached log10 values for the first integer up to the maximum ploidy requested thus far.
- */
- private volatile static double[] ploidyLog10;
-
- // Initialize {@link #ploidyLog10}.
- static {
- ploidyLog10 = new double[maximumPloidy + 1];
- for (int i = 0; i <= maximumPloidy; i++)
- ploidyLog10[i] = Math.log10(i);
- }
-
/**
* Returns an instance given its ploidy and the number of alleles.
*
* @param alleleCount the required allele-count.
* @param ploidy the required ploidy-count.
*
- * @throws IllegalArgumentException if either {@code ploidy} or {@code alleleCount} is {@code null}, or
- * the resulting number of genotypes is too large.
+ * @throws IllegalArgumentException if either {@code ploidy} or {@code alleleCount} is negative, or the resulting number of genotypes is too large.
*
* @return never {@code null}.
*/
- public static GenotypeLikelihoodCalculator getInstance(final int ploidy,
- final int alleleCount) {
+ public static GenotypeLikelihoodCalculator getInstance(final int ploidy, final int alleleCount) {
checkPloidyAndMaximumAllele(ploidy, alleleCount);
// Non-thread safe (fast) check on tables capacities,
- // if not enough capacity we expand the tables in a thread-safe manner:
- if (alleleCount > maximumAllele || ploidy > maximumPloidy)
- ensureCapacity(alleleCount, ploidy);
+ // if not enough capacity we expand the tables in a thread-safe manner
+ // also checks if the requested ploidy and allele count result in a genotype count too large to deal with
+ if(calculateGenotypeCountUsingTables(ploidy, alleleCount) == GENOTYPE_COUNT_OVERFLOW){
+ final double largeGenotypeCount = MathUtils.binomialCoefficient(ploidy + alleleCount - 1, alleleCount - 1);
+ throw new IllegalArgumentException(String.format("the number of genotypes is too large for ploidy %d and allele %d: approx. %.0f", ploidy, alleleCount, largeGenotypeCount));
+ }
// At this point the tables must have at least the requested capacity, likely to be much more.
- return new GenotypeLikelihoodCalculator(ploidy,alleleCount,alleleFirstGenotypeOffsetByPloidy,genotypeTableByPloidy,ploidyLog10);
+ return new GenotypeLikelihoodCalculator(ploidy, alleleCount, alleleFirstGenotypeOffsetByPloidy, genotypeTableByPloidy, ploidyLog10);
}
/**
@@ -413,14 +409,59 @@ public class GenotypeLikelihoodCalculators {
* @param ploidy the requested ploidy.
* @param alleleCount the requested number of alleles.
*
- * @throws IllegalArgumentException if {@code ploidy} or {@code alleleCount} is negative.
+ * @throws IllegalArgumentException if {@code ploidy} or {@code alleleCount} is negative or
+ * the number of genotypes is too large (more than {@link Integer#MAX_VALUE}).
*
- * @return 0 or greater.
+ * @return the number of genotypes given ploidy and allele count (0 or greater).
*/
public final static int genotypeCount(final int ploidy, final int alleleCount) {
+
+ final int result = calculateGenotypeCountUsingTables(ploidy, alleleCount);
+ if (result == GENOTYPE_COUNT_OVERFLOW) {
+ final double largeGenotypeCount = MathUtils.binomialCoefficient(ploidy + alleleCount - 1, alleleCount - 1);
+ throw new IllegalArgumentException(String.format("the number of genotypes is too large for ploidy %d and allele %d: approx. %.0f", ploidy, alleleCount, largeGenotypeCount));
+ }
+ return result;
+ }
+
+ /**
+ * Compute the maximally acceptable allele count (ref allele included) given the maximally acceptable genotype count.
+ * @param ploidy sample ploidy
+ * @param maxGenotypeCount maximum number of genotype count used to calculate upper bound on number of alleles given ploidy
+ * @throws IllegalArgumentException if {@code ploidy} or {@code alleleCount} is negative.
+ * @return the maximally acceptable allele count given ploidy and maximum number of genotypes acceptable
+ */
+ public static int computeMaxAcceptableAlleleCount(final int ploidy, final int maxGenotypeCount){
+
+ checkPloidyAndMaximumAllele(ploidy, ploidy); // a hack to check ploidy makes sense (could duplicate code but choice must be made)
+
+ final double log10MaxGenotypeCount = Math.log10(maxGenotypeCount);
+
+ // Math explanation: genotype count is determined by ${P+A-1 \choose A-1}$, this leads to constraint
+ // $\log(\frac{(P+A-1)!}{(A-1)!}) \le \log(P!G)$,
+ // where $P$ is ploidy, $A$ is allele count, and $G$ is maxGenotypeCount
+ // The upper and lower bounds of the left hand side of the constraint are $P \log(A-1+P)$ and $P \log(A)$
+ // which require $A$ to be searched in interval $[10^{\log(P!G)/P} - (P-1), 10^{\log(P!G)/P}]$
+ // Denote $[10^{\log(P!G)/P}$ as $x$ in the code.
+
+ final double x = Math.pow(10, (MathUtils.log10Factorial(ploidy) + log10MaxGenotypeCount)/ploidy );
+ final int lower = (int)Math.floor(x) - ploidy - 1;
+ final int upper = (int)Math.ceil(x);
+ for(int a=upper; a>=lower; --a){// check one by one
+
+ final double log10GTCnt = MathUtils.log10BinomialCoefficient(ploidy+a-1, a-1);
+ if(log10MaxGenotypeCount >= log10GTCnt) {
+ return a;
+ }
+ }
+ throw new GATKException("Code should never reach here.");
+ }
+
+ private static int calculateGenotypeCountUsingTables(int ploidy, int alleleCount) {
checkPloidyAndMaximumAllele(ploidy, alleleCount);
- if (ploidy > maximumPloidy || alleleCount > maximumAllele)
- ensureCapacity(alleleCount,ploidy);
+ if (ploidy > maximumPloidy || alleleCount > maximumAllele) {
+ ensureCapacity(alleleCount, ploidy);
+ }
return alleleFirstGenotypeOffsetByPloidy[ploidy][alleleCount];
}
}
\ No newline at end of file
diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypingEngine.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypingEngine.java
index 83c7ed533..2635fea66 100644
--- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypingEngine.java
+++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypingEngine.java
@@ -60,6 +60,7 @@ import org.broadinstitute.gatk.tools.walkers.annotator.VariantAnnotatorEngine;
import org.broadinstitute.gatk.tools.walkers.genotyper.afcalc.AFCalculationResult;
import org.broadinstitute.gatk.tools.walkers.genotyper.afcalc.AFCalculator;
import org.broadinstitute.gatk.tools.walkers.genotyper.afcalc.AFCalculatorProvider;
+import org.broadinstitute.gatk.tools.walkers.genotyper.afcalc.AlleleFrequencyCalculator;
import org.broadinstitute.gatk.utils.GenomeLoc;
import org.broadinstitute.gatk.utils.GenomeLocParser;
import org.broadinstitute.gatk.utils.MathUtils;
@@ -68,6 +69,7 @@ import org.broadinstitute.gatk.utils.contexts.AlignmentContext;
import org.broadinstitute.gatk.utils.contexts.AlignmentContextUtils;
import org.broadinstitute.gatk.utils.contexts.ReferenceContext;
import org.broadinstitute.gatk.utils.exceptions.UserException;
+import org.broadinstitute.gatk.utils.genotyper.PerReadAlleleLikelihoodMap;
import org.broadinstitute.gatk.utils.genotyper.SampleList;
import org.broadinstitute.gatk.utils.gga.GenotypingGivenAllelesUtils;
import org.broadinstitute.gatk.utils.pileup.ReadBackedPileup;
@@ -104,8 +106,9 @@ public abstract class GenotypingEngine upstreamDeletionsLoc = new LinkedList<>();
+
+ protected final AFCalculator newAFCalculator;
/**
* Construct a new genotyper engine, on a specific subset of samples.
@@ -139,6 +142,11 @@ public abstract class GenotypingEngine stratifiedContexts,
final VariantContext vc, final GenotypeLikelihoodsCalculationModel.Model model,
final boolean inheritAttributesFromInputVC,
- final Map