Cap likelihoods symmetrically #6196
Merged
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
ldgauthier
commented
Oct 3, 2019
| @@ -206,7 +206,7 @@ public static CachingIndexedFastaSequenceFile createReferenceReader(final String | |||
| * @return never {@code null}. | |||
| */ | |||
| public static ReadLikelihoodCalculationEngine createLikelihoodCalculationEngine(final LikelihoodEngineArgumentCollection likelihoodArgs) { | |||
| final double log10GlobalReadMismappingRate = likelihoodArgs.phredScaledGlobalReadMismappingRate < 0 ? - Double.MAX_VALUE | |||
| final double log10GlobalReadMismappingRate = likelihoodArgs.phredScaledGlobalReadMismappingRate < 0 ? Double.NEGATIVE_INFINITY | |||
There was a problem hiding this comment.
The docs say you can turn off the cap by setting this value negative, but it wasn't working before because subsequent methods were looking for -Inf.
davidbenjamin
requested changes
Oct 3, 2019
| @@ -352,7 +352,8 @@ public void normalizeLikelihoods(final double maximumLikelihoodDifferenceCap) { | |||
| private void normalizeLikelihoodsPerEvidence(final double maximumBestAltLikelihoodDifference, | |||
| final double[][] sampleValues, final int sampleIndex, final int evidenceIndex) { | |||
|
|
|||
| final BestAllele bestAlternativeAllele = searchBestAllele(sampleIndex,evidenceIndex,false); | |||
| //allow the best "alternative" allele to be the reference because asymmetry leads to strange artifacts like het calls with >90% alt reads | |||
| final BestAllele bestAlternativeAllele = searchBestAllele(sampleIndex,evidenceIndex,true); | |||
There was a problem hiding this comment.
Can we then just call this bestAllele?
| }; | ||
| runCommandLine(argsNoMaxAlternateAlleles); | ||
|
|
||
| List<VariantContext> vcs = VariantContextTestUtils.readEntireVCFIntoMemory(outputContaminatedHomVarDeletions.getAbsolutePath()).getRight(); |
There was a problem hiding this comment.
VariantContextTestUtils.getVariantContexts is all you need.
|
|
||
| // Run both with and without --max-alternate-alleles over our interval, so that we can | ||
| // prove that the argument is working as intended. | ||
| final String[] argsNoMaxAlternateAlleles = { |
There was a problem hiding this comment.
I think an ArgumentBuilder would be cleaner:
final ArgumentBuilder args = new ArgumentBuilder().addInput(bam)
.addReference(hg38Reference)
.addInterval(new SimpleInterval(intervals))
.addOutput(outputContaminatedHomVarDeletions)
.addNumericArgument(IntervalArgumentCollection.INTERVAL_PADDING_LONG_NAME, 50);
runCommandLine(args);
|
|
||
| final File outputContaminatedHomVarDeletions = createTempFile("testContaminatedHomVarDeletions", ".vcf"); | ||
|
|
||
| // Run both with and without --max-alternate-alleles over our interval, so that we can |
ldgauthier
force-pushed
the
ldg_capLikelihoodsSymmetrically
branch
from
4b4d7ef to
36caaa9
Compare
|
@davidbenjamin take another look? |
ldgauthier
force-pushed
the
ldg_capLikelihoodsSymmetrically
branch
from
36caaa9 to
c5b0ba6
Compare
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Big events or haplotypes with a lot of events can lead to very disparate likelihoods between the reference allele and the alt allele. We put a cap on how unlikely an alt read can be, but the behavior was not symmetric with respect to the reference. As a result, some sites with large indels or many SNPs could get genotyped as heterozygous with >90% alt alleles, where the one or two reference reads are likely derived from cross-sample contamination. This genotyping artifact is now fixed, though the cap is still set by the "global mismapping penalty", or the probability that a read is derived from another locus in the genome and not either of the haplotypes represented in this ActiveRegion.