From 9a1d0d7076853670dcace2dcd4e37a3d55d0566c Mon Sep 17 00:00:00 2001 From: depristo Date: Wed, 11 May 2011 14:17:14 +0000 Subject: [PATCH] Simple bug fix to allow multiple records at same site when genotyping given alleles. Takes only the first record (respecting filters, SNP type, etc), and issues a warning if there is more than one valid record at a site git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@5789 348d0f76-0448-11de-a6fe-93d51630548a --- ...NPGenotypeLikelihoodsCalculationModel.java | 26 +++++++++++++------ 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java index 54ce35bbd..ae5c81642 100755 --- a/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java +++ b/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/SNPGenotypeLikelihoodsCalculationModel.java @@ -73,18 +73,28 @@ public class SNPGenotypeLikelihoodsCalculationModel extends GenotypeLikelihoodsC if ( alternateAlleleToUse != null ) { bestAlternateAllele = alternateAlleleToUse.getBases()[0]; } else if ( useAlleleFromVCF ) { - final VariantContext vcInput = tracker.getVariantContext(ref, "alleles", null, ref.getLocus(), true); - if ( vcInput == null || vcInput.isFiltered() || !vcInput.isSNP() ) - return null; - if ( !vcInput.isSNP() ) { - logger.info("Record at position " + ref.getLocus() + " is not a SNP; skipping..."); - return null; + VariantContext vc = null; + + // search for usable record + for( final VariantContext vc_input : tracker.getVariantContexts(ref, "alleles", null, ref.getLocus(), true, false) ) { + if ( vc_input != null && ! vc_input.isFiltered() && vc_input.isSNP() ) { + if ( vc == null ) { + vc = vc_input; + } else { + logger.warn("Multiple valid VCF records detected at site " + ref.getLocus() + ", only considering alleles from first record only"); + } + } } - if ( !vcInput.isBiallelic() ) { + + // ignore places where we don't have a variant + if ( vc == null ) + return null; + + if ( !vc.isBiallelic() ) { // for multi-allelic sites go back to the reads and find the most likely alternate allele initializeBestAlternateAllele(refBase, contexts); } else { - bestAlternateAllele = vcInput.getAlternateAllele(0).getBases()[0]; + bestAlternateAllele = vc.getAlternateAllele(0).getBases()[0]; } } else { initializeBestAlternateAllele(refBase, contexts);