A couple of fixes for writing VCFs in GenotypeConcordance #810

meganshand · 2017-05-05T17:54:48Z

Description

There are a couple of fixes when running GenotypeConcordance with OUTPUT_VCF=true:

Fix in GenotypeConcordance for writing VCFs with no-call genotypes #785 ensured that the VCF builder used the original site alleles instead of the genotype alleles in order to prevent the writing of no calls. Unfortunately, this meant that the site alleles were not being normalized by normalizeAlleles, which is needed when you are merging two variant contexts with different REF alleles (due to indels). This is now fixed by using the normalized genotype alleles, but removing no calls from the set before making the builder.
When a no call site is normalized the new genotype allele would add the surplus bases after the .. For example if the ref allele was changing from A to ACCC, the no call genotype would change from . to .CCC. This is fixed/tested in the last two commits.

Checklist (never delete this)

Never delete this, it is our record that procedure was followed. If you find that for whatever reason one of the checklist points doesn't apply to your PR, you can leave it unchecked but please add an explanation below.

Content

Added or modified tests to cover changes and any new functionality
Edited the README / documentation (if applicable)
All tests passing on Travis

Review

Final thumbs-up from reviewer
Rebase, squash and reword as applicable

For more detailed guidelines, see https://github.com/broadinstitute/picard/wiki/Guidelines-for-pull-requests

coveralls · 2017-05-05T18:05:36Z

Coverage increased (+0.007%) to 73.986% when pulling 10282e8 on ms_gcWriteMultiRefAlleles into da60c2b on master.

yfarjoun

nice catch @meganshand

a few comments

yfarjoun · 2017-05-05T19:48:01Z

testdata/picard/vcf/normalize_no_calls_truth.vcf

+##contig=<ID=GL000203.1,length=37498>
+##contig=<ID=GL000246.1,length=38154>
+##contig=<ID=GL000249.1,length=38502>
+##contig=<ID=GL000196.1,length=38914>


for berevity, could you remove these GL*** contigs?

yfarjoun · 2017-05-05T19:48:30Z

testdata/picard/vcf/normalize_no_calls_truth.vcf

+##INFO=<ID=HapNoVar,Number=1,Type=Integer,Description="Number of datasets for which HaplotypeCaller called a variant within 35bp and did not call a variant at this location">
+##INFO=<ID=LEN,Number=A,Type=Integer,Description="allele length">
+##INFO=<ID=NoCG,Number=0,Type=Flag,Description="Present if no consensus reached, so looked at all datasets except Complete Genomics since it may have a different representation of complex variants">
+##INFO=<ID=NoPLTot,Number=1,Type=Integer,Description="Number of datasets with likelihood ratio > 20 for a genotype different from the called genotype">


also, could you remove all these unused header lines?

ideally, one would be able to see both VCFs at the same time in the github diff...

yfarjoun · 2017-05-05T19:48:57Z

src/main/java/picard/vcf/GenotypeConcordance.java

@@ -432,14 +432,11 @@ private void writeVcfTuple(final VcfTuple tuple, final VariantContextWriter writ
            final List<Allele> truthAlleles = alleles.truthAlleles();
            final List<Allele> callAlleles  = alleles.callAlleles();

-            // Get the alleles present at this site for both samples to use for the output variant context.
+


no need for extra newline

yfarjoun · 2017-05-05T19:49:56Z

src/main/java/picard/vcf/GenotypeConcordance.java

-                siteAlleles.addAll(callContext.getAlleles());
-            }
+            siteAlleles.addAll(allAlleles);
+            siteAlleles.remove(Allele.NO_CALL);


please add a symbolic allele to your test.

It seems that GenotypeConcordance doesn't support symbolic alleles:

Note that only SNP and INDEL variants are considered, MNP, Symbolic, and Mixed classes of variants are not included.

That said, a symbolic allele does cause it to blow up if you are writing an output VCF. In this case what should the tool do? Skip that site or blow up with an appropriate error message?

yfarjoun · 2017-05-05T19:53:59Z

src/test/java/picard/vcf/GenotypeConcordanceTest.java

+        genotypeConcordance.OUTPUT = new File(OUTPUT_DATA_PATH, "MultipleRefAlleles");
+        genotypeConcordance.OUTPUT_VCF = true;
+
+        Assert.assertEquals(genotypeConcordance.instanceMain(new String[0]), 0);


can you also read in the metrics file and validate the results in it, rather than just checking that the program didn't blowup?

meganshand · 2017-05-08T22:03:52Z

src/main/java/picard/vcf/GenotypeConcordance.java

@@ -327,11 +327,11 @@ private boolean indexExists(final File vcf) {
                final String condition = truthVariantContextType + " " + callVariantContextType;
                final Integer count = unClassifiedStatesMap.getOrDefault(condition, 0) + 1;
                unClassifiedStatesMap.put(condition, count);
+            } else {
+                // write to the output VCF
+                writer.ifPresent(w -> writeVcfTuple(tuple, w, scheme));


Note that if we prevent writing unless stateClassified is true, we will not write sites with no variation in either original VCF. This means that a site that is HOM REF in both the truth and call VCFs will not be output. (This is currently breaking other GenotypeConcordance tests.) Is the desired behavior to output no variation sites?

yfarjoun · 2017-05-28T18:19:08Z

👍 Thanks for this @meganshand

…ymbolic allele

coveralls · 2017-06-01T19:40:35Z

Coverage remained the same at 74.077% when pulling dd4dc88 on ms_gcWriteMultiRefAlleles into 1c54e16 on master.

yfarjoun reviewed May 5, 2017

View reviewed changes

meganshand commented May 8, 2017

View reviewed changes

Fixes GenotypeConcordance VCF context with different Ref alleles or s…

dd4dc88

…ymbolic allele

meganshand force-pushed the ms_gcWriteMultiRefAlleles branch from 2834099 to dd4dc88 Compare June 1, 2017 19:30

meganshand merged commit d4a632a into master Jun 1, 2017

meganshand deleted the ms_gcWriteMultiRefAlleles branch June 1, 2017 19:46

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

A couple of fixes for writing VCFs in GenotypeConcordance #810

A couple of fixes for writing VCFs in GenotypeConcordance #810

meganshand commented May 5, 2017 •

edited

coveralls commented May 5, 2017

yfarjoun left a comment

yfarjoun May 5, 2017

meganshand May 8, 2017

yfarjoun May 5, 2017

yfarjoun May 5, 2017

meganshand May 8, 2017

yfarjoun May 5, 2017

meganshand May 8, 2017

yfarjoun May 5, 2017

meganshand May 8, 2017

yfarjoun May 5, 2017

meganshand May 8, 2017

meganshand May 8, 2017

yfarjoun commented May 28, 2017

coveralls commented Jun 1, 2017 •

edited

A couple of fixes for writing VCFs in GenotypeConcordance #810

A couple of fixes for writing VCFs in GenotypeConcordance #810

Conversation

meganshand commented May 5, 2017 • edited

Description

Checklist (never delete this)

Content

Review

coveralls commented May 5, 2017

yfarjoun left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

yfarjoun commented May 28, 2017

coveralls commented Jun 1, 2017 • edited

meganshand commented May 5, 2017 •

edited

coveralls commented Jun 1, 2017 •

edited