This change doesn't affect the performance of the Indel Realigner at all (as per tests).
This is just a request from the Picard side (where further testing is happening).
Note that this patch involves ignoring supplementary alignments. Ideally we would want
to fix their mates properly but that would require a major refactoring of this soon-to-be
deprecated tool.
Find out about a dev-bug and added TODOs (reported in #1096).
Addresses issue #1095.
Conflicts:
protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/haplotypecaller/HaplotypeCaller.java
Changed a division by -10.0 to a multiplication by -.1 in QualUtils (typically multiplication is faster than division).
Addresses performance issue #1081.
Now that Ron updated the GATK so that we use star to represent spanning
deletions, we need to catch those cases in the code that remaps alleles.
Otherwise, we try to pad the stars and that's just bad.
Added test from actual failing data.
When a sample has multiple spanning deletions and we are asked to assign
likelihoods to the spanning deletion allele, we currently choose the first
deletion. Valentin pointed out that this isn't desired behavior. I
promised Valentin that I would address this issue, so here it is.
I do not believe that the correct thing to do is to sum the likelihoods
over all spanning deletions (I came up with problematic cases where this
breaks down).
So instead I'm using a simple heuristic approach: using the hom alt PLs, find
the most likely spanning deletion for this position and use its likelihoods.
In the 10K-sample VCF from Monkol there were only 2 cases that this problem
popped up. In both cases the heuristic approach works well.
Add oxoG read count annotation and add as default annotation
Add ##SAMPLE VCF header line in accordance with TCGA VCF spec, specifying "File" line in sample header with BAM file name and "SampleName" with BAM sample name (Don't print sample file path if --no_cmdline_in_header is specified to help with test consistency)
Turn on active region assembly-based physical phasing for M2
Clean up M2-related annotations so UG doesn't crash if M2 annotations are called
-We now pull htsjdk and picard from maven central.
-Updated the GATK codebase as necessary to adapt to changes in the Feature
interface.
-Since VCFHeader now requires that all header lines have unique keys, uniquified
the keys of GVCFBlock header lines by including the min/max GQ in the key.
Updated MD5s accordingly.
-Other MD5s changed as a result of an htsjdk fix to eliminate "-0" in VCF output.
In the case where there's a low quality SNP under a spanning deletion in the gvcfs:
if the SNP is not genotyped by GenotypeGVCFs (because it's just noise) we were still
emitting a record with just the symbolic DEL allele (because that allele is high quality).
We no longer do that.
Previously, if a SNP occurred in sample A at a position that was in the middle of a deletion for sample B,
sample B would be genotyped as homozygous reference there (but it's NOT reference - there's a deletion).
Now, sample B is genotyped as having a symbolic DEL allele.
Minor cleanup added. Note that I also removed Laura's previous fix for this problem.
Existing integration tests change because I've added a new header line to the VCF being output.
I also added several tests for the new functionality showing:
1. genotyping from separate and already combined gvcfs give the same output
2. genotyping over multiple spanning deletions works
3. combining works too
Existing unit tests also cover this case.
New unit test for deprecated mergeVariantsViaLD
Update HaplotypeCallerIntegrationTest.java
Delete duplicate testHaplotypeCallerMergeVariantsViaLDException test.
Exclude MQ0BySample
Move SD and TRA to new StandardUGAnnotation interface
There is now annotation interface (StandardUGAnnotation) holding annots that are standard in UG but should't be used as they are now with HC. This allows us to not have to exclude these annotations explicitly in HC, but still be able to use them for development purposes.
Fairly minor if plentiful fixes to various gatkdocs. Merging this without formal review since all tests pass, the gatkdocs build, and no one really wants to review corrections to grammar, typos and layout for 120+ documents. Review will be done by users in production ;-)
fix blasted license blurbs
updates based on PR comments (abstractify HaplotypeCallerArgumentCollection into AssemblyBasedCallerArgumentCollection)
comments on comments from PR review
Build a ReferenceContext in ActiveRegionWalkers to pass in to annotation engine so we can call the TandemRepeatAnnotator from M2
Make TandemRepeatAnnotator default annotation for M2.
Setup (but don't use yet) HC-style contamination downsampling.
New HC integration test with TandemRepeatAnnotator
- ASEReadCounter (public tool) replce Tuuli's script to produce the input to Manny's tool.
It count the number of reads that support the ref allele and the alt allele, filtereing low qual reads and bases and keep only properPaired reads
- ASECaller (private tool) take both RNA and DNA, and produce ontingencyTables ** still under development **
minor changes in other tools:
- update RNA HC variant calling scala script
- expose FS method pValueForContingencyTable to be able to call it from ASEcaller
In ASEReadCounter:
- allow different option to deal with overlaping read from the same fragment
- add option to ignore or include indels in the pileups
- add option to disabled DuplicateRead
add ASEReadCounterIntegrationTest.java and files for the test
fixed NPE when normal contains no reads
first integration test (micro) and unit tests, also rename of MuTectHC -> M2
adding in standard GATK license terms
incorporated HOSTILE mode to PCR Error Correction
removed tumor and normal name parameters and cleaned up internal name handling
changes to allow for calling without a matched normal (technically, not true 'tumor-only' calling). Used for panel-of-normals creation
additional regression tests, based on DREAM data. Removed accidental addition of TandemRepeatAnnotator to default annotations
updated MD5 based on run from GSA4 to fix bamboo issue
reverted unneeded visibility changes
Now, instead of stripping out the GQs for mono sites, we transfer them to the RGQ.
This is extremely useful for people who want to know how confident the hom ref genotype calls are.
Perhaps this is just what CRSP needs for pertinent negatives.
Note that I also changed the tool to no longer use the GenotypeSummaries annotation by default since
it was adding some seemingly unnecessary annotations (like mean GQ now that we keep the GQ around and
number of no-calls). Let me know if this was a mistake (although Laura gave me a thumbs up).
Using --breakBandsAtMultiplesOf N will ensure that no reference blocks span across
genomic positions that are multiples of N. This is especially important in the
case of scatter-gather where you don't want your scatter intervals to start in the
middle of blocks (because of a limitation in the way -L works in the GATK for VCF
records with the END tag).
For example, running with --breakBandsAtMultiplesOf 5 on this record:
1 69491 . G <NON_REF> . . END=69523 GT:DP:GQ:MIN_DP:MIN_GQ:PL ./.:94:99:82:99:0,120,1800
Will produce the following records:
1 69491 . G <NON_REF> . . END=69494 GT:DP:GQ:MIN_DP:MIN_GQ:PL ./.:94:99:82:99:0,120,1800
1 69495 . C <NON_REF> . . END=69499 GT:DP:GQ:MIN_DP:MIN_GQ:PL ./.:94:99:82:99:0,120,1800
1 69500 . T <NON_REF> . . END=69504 GT:DP:GQ:MIN_DP:MIN_GQ:PL ./.:94:99:82:99:0,120,1800
etc.
Added docs and a new test.
GenotypeGVCFs now has the ability to unique-ify samples so I can genotype together two different datasets containing the same sample
Modify InbreedingCoeff so that it works when genotyping uniquified samples
* The value of this element (default true) determines whether Queue will explicitly run this walker over unmapped reads
* This patch fixes a runtime error when FindCoveredIntervals was used with Queue
* PT 81777160
* TextCigarCodec.decode() is now static, and the getSingleton() method is gone
* MergingSamRecordIterator now wants a Collection<SamReader> rather than Collection<SAMFileReader> in the constructor
* SeekableBufferedStream now correctly reads the requested number of bytes, removed workaround in GATKBAMIndex
* Removed unused annotations (CCC and HWP)
* Renamed one of the two GC annotations to "IGC" (for Interval GC)
* Revved picard & htsjdk (GATK constants are now removed from htsjdk)
* PT 82046038
-- Active Region Traversal was using per sample limits on the number of reads that were too low, especially now that we are running one sample at a time. This caused issues with high confidence variants being dropped in high coverage data.
-- HaplotypeCallerGVCFIntegrationTest PL/annotation changes due to using more reads in those tests
-- Removed a CountReadsInActiveRegionsIntegrationTest test for excessive coverage because the read coverage no longer goes over the limits in ART
Story:
=====
- https://www.pivotaltracker.com/story/show/83803796
Changes:
=======
- From a fix maximum ploidy indel RCM likelihood cache to a
dynamically resizable one.
- Used the occassion to removed an unused and deprecated method from ReferenceConfidenceModel
Testing:
=======
- Added integration test to check on ploidies larger than the previous limit of 20.
Add multi-allele test for info field annotations
Fix to process all types of INFO annotations
roll back to previous version, removes INFO and FORMAT
Correct @return for VariantAnnotatorEngine.getNonReferenceAlleles()
Enhance comments and clean up multi-allelic logic, handle header info number = R
only parse counts of A & R
Add INFO for AC
update MD5
Performance enhancement, only parse multiallelic with a count A or R
Make argument final in getNonReferenceAlleles()
Code cleanup, add exceptions for bad expression/allele size mismatch and missing header info for an expression
Change exception to warning for expression value/number of alleles check
remove adevertised exceptions
-- Ignore SNP matches that lie outside the clipped read window
-- This fixes an issue where GATK would skip the entire read if a SNP is entirely
contained within a sequencing adapter.
Story:
=====
- https://www.pivotaltracker.com/story/show/83259038
Changes:
=======
- Done minimal changes to make the fix after an arduous attempt to understand
CombineGVCFs code.
Test:
====
- Added a integration test to explicitly test for the bug.
- Updated a md5 changes as the bug was actually affecting one of the existing
integration tests.
* PT 84242218
* Note that FORMAT fields behave the same as INFO fields - if the annotation has a count of A (one entry per Alt Allele), it is split across the multiple output lines. Otherwise, the entire list is output with each field
Story:
-----
- https://www.pivotaltracker.com/story/show/83800586
Changes:
-------
- In GVCFWriter GQ is now recalculated out of the fianl PL array for the block.
Testing:
-------
- Updated affected integration test md5s
Add more logging to annotators, change loggers from info to warn
Add comments to testStrandBiasBySample()
Clarify comments in testStrandBiasBySample
remove logic for not prcossing an indel if strand bias (SB) was not computed
remove per variant warnings in annotate()
Log warnings if using the wrong annotator or missing a pedgree file
Log test failures once in annotate(), because HaplotypeCaller does not call initialize(). Avoid using exceptions
Fix so only log once in annotate(), Hardey-Weinberg does not require pedigree files, fix test MD5s so pass
Check if founderIds == null
Update MD5s from HaplotypeCaller integrations tests and clean up code
Change logic so SnpEff does not throw excpetions, change engine to utils in imports
Update test MD5s, return immediately if cannot annotate in SnpEff.initialization()
Post peer review, add more logging warnings
Update MD5 for testHaplotypeCallerMultiSampleComplex1, return null if PossibleDeNovo.annotate() is not called by VariantAnnotator
Story:
-----
https://www.pivotaltracker.com/story/show/80684230
Changes:
-------
- Corrected the bug: AlignmentUtils#createReadAlignedToRef was
not realigning against the reference but the best haplotype for
the read.
Test:
----
- Added integration test in HaplotypeCallerIntegrationTest to check
that the bug has been fixed.
- Fixed md5s modified by this change; these are cause due to small
changes in the state of the random-number generator and read vs
variant site overlapping.
CombineGVCFs now outputs ref conf for the duration of deletions so that SNPs occuring in other samples aligned with those deletions will be genotyped correctly
Reading the multiple GATKText files as a single stream, especially with new top level target executable jar files pointing to a lib folder.
Don't dirty the build with a new GATKText.properties if input files are unmodified.
Stop warning on undocumented abstract classes.
Fixed ClassNotFoundException/NoClassDefFoundError by fixing ResourceBundleExtractorDoclet artifact.
Excluding Exceptions from documentation.
Removed custom log4j dependency from ResourceBundleExtractorDoclet.
Stop generating the dependency reduced pom during shade.
Stop regenerating gsalib when the files are already up to date.
Disabled mvn site generation from external-example.
Moved top level target symlinks to package jar files to under target/package.
Executable jar files are placed under target/executable with the new target[/lib] directories.
Under top level target, symlinks to *either* the package *or* the executable jars replace what was a symlink to the package jar path.
Allow disabling of the shade package.
ant-bridge.sh by default only builds executable jars, and doesn't package by default, as did the old ant build.xml.
Added a new package_path.sh utility script for other scripts to use instead of anything in the target folder.
remove final keyword before refMap and altMap, constructHaplotype() changes their values
return ArtificialHaplotype from constructHaplotype instaed of passing as an argument
Add logic so arraycopy does not throw an IndexOutOfBoundsException, add test for a long insert
* This argument is intended to be used in conjunction with -bamout, and disable early-exit optimizations to allow reference regions to be contained in the output bam
* Also forcibly includes the reference haplotype in the set of haplotypes given to the BAMWriter
* Made -dontTrimActiveRegions visible, as it is likely also desirable in this use case
* Addresses PT 77731660
remove TODO comment after activeProbThreshold
recover static ACTIVE_PROB_THRESHOLD for unit tests
Add min/max values for active_probability_threshold parameter
Move activeProbThreshold parameter to GATKArguemtnCollection
define ACTIVE_PROB_THRESHOLD in unit tests
add construction of argCollection in in ctor
Move arguments from GATKArgumentCollection to ActiveRegionWalker
Throw exception if threshold < 0 or > 1 in ActivityProfile ctor
max propogation distance parameter to ActiveRegionWalker for AcrtivityProfile
Use polymorphic getMaxProbPropagationDistance() so BandPassActivityProfile computes the crrect region size cutoff
Get the maxProbPropagationDistance from the super class's method, instead of directly, this is safer
Removed extraneous command line imports and make maxProbPropagationDistance a hidden argument
remove limit check for activeProbThreshold, not necessary because the check is made when imput as a command line arg
Remove extra 'region' in the doxygen param description for maxProbPropagationDistance
Rename parameters using camel case and add to integration test
Correct documentation for maxReadsInRegionPerSample and minReadsPerAlignmentStart
Change the argument--minReadsPerAlignmentStart in the integration test from 50 to 5
'each genomic location' only pertains to minReadsPerAlignmentStart, not maxReadsInRegionPerSample
The QUAL value calculated by this Exact AF Calculator is very underestimated when
there are more than one alternative allele (non-biallelic sites). The reason is
that the QUAL was roughly calculated by adding the QUALs resulting of each alternative
alleles vs all other alleles, reference and alts, collapsed. This is ok for MLEAC
calculations but not for QUAL.
Now, for calculating the QUAL we collapse all the alternatives as only one. This change
improves sensitivy with a cost of additional false positives, but this is naturally expected.
The resulting QUAL column is much closer to the one returned by the reference implementation.
Story:
https://www.pivotaltracker.com/story/show/75926368.
Changes:
Changed the QUAL calculation as described above.
Updated MD5s.
Fixed MD5s
The problem whas that the MLE table calculation aborted "unlikely"
genotype combinations to aggresively.
This also uncovered another bug where GeneralPloidyExactAFCalculation
makes a slightly different use of StateTracker
as compared to DiploidExactAFCalculation. We have changed StateTracker
generalizing it to be able to work with both using code behaviors.
Story:
-----
* https://www.pivotaltracker.com/story/show/78920568
Changes:
-------
* Fixes in GeneralPloidyExactAFCalculator.
* Needed changes in StateTracker API and its consequences in DiploidExactAFCalculation.
* Updated affected integrated tests' MD5s after fixing the GeneralPloidyExactAF.
Changes:
-------
* Updated current unit and integration test to use the new API components.
* Added unit tests for new classes AFPriorProvider and AFCalculatorProviders.
* Added integration test for mixed ploidy GenotypeGVCFs and CombineGVCFs
Changes:
-------
* GenotypingEngine uses now a AFCalc provider instead of
its own thread-local with one-time initialized and fixed
AF calculator.
* All walkers that use a GenotypingEngine now are passing
the appropiate AF calculator provider. For now most
just use a fix calculator (FixedAFCalculatorProvider)
except GenotypeGVCFs as this one now can cope with
mixture of ploidies failing-over to a general-ploidy
calculator when the preferred implementation is not
capable to handle a site's analysis.
to the total-ploidy (added ploidy accross samples).
Changes:
--------
* Instead of calculate a fixed log10 prior array with a fix
total likelihood we use a new component, the AFPriorProvider
to generate the priors for different total plodies on
demand; these are cached however so there is no unecessary
recompute involved.
with mixed ploidies and max-alt-allele number changes dynamically.
Changes:
--------
* Moved the AFCalcFactory.Calculation enum in a top level class
AFCalculatorImplementation.
* Given more reponsabilities to the enum like resolving the constructor
method once per implementation and the best-model selection algorithm.
* Removed test-code only fields and methods from AFCalc; just used to perform
unit-testing and not any actual functionality of this component.
* Removed the fixed ploidy constraint of GeneralPloidyExactAFCalc
implementation... now can deal with mixed ploidies that may change
per site and sample.
* Removed the fixed maxAltAllele restriction by allowing resizing of
the stateTracker structures.
* Due to previous two points now call the the AFCalc object are passed
the default-ploidy to assume in case some genotype in the input
VC does not have it and the max-alt-allele.
* Also due to those changes, removed the now totally useless 3 int
parameters from all AFCalc constructors.
* Cleaned the code a bit from no further used components and methods.
Dangling head merging (like with tails) in now enabled by default.
The --recoverDanglingHeads argument is now deprecated so that users know not to use it anymore.
We now also allow the user to set the minimum branch length for merging. This will be different
for exomes and RNA (see below).
The other changes in the code itself:
1. We no longer allow an arbitrarily large number of mismatches in the dangling head for merging
2. The max number of mismatches allowed in a dangling head is proportional to the kmer size
There will be a difference in the RNA calling pipeline. Instead of invoking '--recoverDanglingHeads'
the user will instead want to use '--minDanglingBranchLength 0'.
Below are the knowledgebase results of the master branch vs. this one.
For NA12878 DNA Exome:
master SNPS TRUE_POSITIVE 36722
master SNPS CALLED_NOT_IN_DB_AT_ALL 2699
master SNPS REASONABLE_FILTERS_WOULD_FILTER_FP_SITE 292
master SNPS FALSE_POSITIVE_SITE_IS_FP 70
branch SNPS TRUE_POSITIVE 36867
branch SNPS CALLED_NOT_IN_DB_AT_ALL 2952
branch SNPS REASONABLE_FILTERS_WOULD_FILTER_FP_SITE 387
branch SNPS FALSE_POSITIVE_SITE_IS_FP 94
As I discussed with Ryan in person, there are a good number of FPs that are called in the new
code, but they nearly all have bad strand bias and should be easily filtered by VQSR.
Note that there is no change for indels.
For NA12878 RNA from Ami:
master SNPS TRUE_POSITIVE 11055
master SNPS CALLED_NOT_IN_DB_AT_ALL 831
master SNPS REASONABLE_FILTERS_WOULD_FILTER_FP_SITE 44
master SNPS FALSE_POSITIVE_SITE_IS_FP 96
branch SNPS TRUE_POSITIVE 11113
branch SNPS CALLED_NOT_IN_DB_AT_ALL 874
branch SNPS REASONABLE_FILTERS_WOULD_FILTER_FP_SITE 47
branch SNPS FALSE_POSITIVE_SITE_IS_FP 92
Again, there's basically no change for indels.
* Arguments involved are --no_cmdline_in_header, --sites_only, and --bcf for VCF files and --bam_compression, --simplifyBAM, --disable_bam_indexing, and --generate_md5 for BAM files
* PT 52740563
* Removed ReadUtils.createSAMFileWriterWithCompression(), replaced with ReadUtils.createSAMFileWriter(), which applies all appropriate engine-level arguments
* Replaced hard-coded field names in ArgumentDefinitionField (Queue extension generator) with a Reflections-based lookup that will fail noisily during extension generation if there's an error
Explicitly including gatk/queue test-jar artifacts in package test classpaths.
SelectVariantsIntegrationTest#testInvalidJexl now resets the JexlEngine silent flag that VariantFiltration.initialize() toggles.
External example no longer tries to unpack nonexistent gatk artifact jars during package tests.
Eric Banks
ldgauthier
Ron Levine
Samuel Lee
Mark Fleharty
Bertrand Haas
Geraldine Van der Auwera
David Benjamin
vruano
Valentin Ruano Rubio
meganshand
Louis Bergelson
Joseph White
Sheila Chandran
melonistic
David Roazen
Kristian Cibulskis
Yossi Farjoun
Phillip Dexheimer
Ami Levy-Moonshine
rpoplin
Khalid Shakir
Ron Levine
Menachem Fromer
Ron Levine