From 32516a2f603cab4d77a9ac36e14dcd5be621ec15 Mon Sep 17 00:00:00 2001 From: Eric Banks Date: Thu, 26 Jul 2012 01:50:39 -0400 Subject: [PATCH] Initial checkpoint commit of VariantContext/Allele refactoring. There were just too many problems associated with the different representation of alleles in VCF (padded) vs. VariantContext (unpadded). We are moving VC to use the VCF representation. No more reference base for indels in VC and no more trimming and padding of alleles. Even reverse trimming has been stopped (the theory being that writers of VCF now know what they are doing and often want the reverse padding if they put it there; this has been requested on GetSatisfaction). Code compiles but presumably pretty much all tests with indels with fail at this point. --- ...olGenotypeLikelihoodsCalculationModel.java | 1 - .../haplotypecaller/GenotypingEngine.java | 28 +- .../GenotypingEngineUnitTest.java | 46 +- .../gatk/refdata/VariantContextAdaptors.java | 72 +-- .../annotator/DepthPerAlleleBySample.java | 57 +-- .../walkers/beagle/BeagleOutputToVCF.java | 2 - .../walkers/beagle/ProduceBeagleInput.java | 2 +- .../beagle/VariantsToBeagleUnphased.java | 2 +- .../genotyper/ConsensusAlleleCounter.java | 12 +- ...elGenotypeLikelihoodsCalculationModel.java | 2 +- .../genotyper/UnifiedGenotyperEngine.java | 14 +- .../walkers/indels/SomaticIndelDetector.java | 28 +- .../validationsiteselector/GenomeEvent.java | 8 +- .../KeepAFSpectrumFrequencySelector.java | 2 +- .../UniformSamplingFrequencySelector.java | 2 +- .../evaluators/ThetaVariantEvaluator.java | 2 +- .../variantutils/LeftAlignVariants.java | 2 +- .../variantutils/LiftoverVariants.java | 8 +- .../variantutils/ValidateVariants.java | 39 +- .../walkers/variantutils/VariantsToTable.java | 11 +- .../walkers/variantutils/VariantsToVCF.java | 14 +- .../broadinstitute/sting/utils/BaseUtils.java | 31 ++ .../sting/utils/codecs/bcf2/BCF2Codec.java | 23 - .../utils/codecs/vcf/AbstractVCFCodec.java | 31 +- .../utils/codecs/vcf/VCFAlleleClipper.java | 434 ------------------ .../sting/utils/variantcontext/Allele.java | 67 +-- .../utils/variantcontext/VariantContext.java | 96 +--- .../variantcontext/VariantContextBuilder.java | 26 +- .../variantcontext/VariantContextUtils.java | 94 +--- .../variantcontext/writer/BCF2Writer.java | 3 - .../variantcontext/writer/VCFWriter.java | 1 - .../utils/variantcontext/AlleleUnitTest.java | 10 +- .../VariantContextUnitTest.java | 1 - 33 files changed, 218 insertions(+), 953 deletions(-) delete mode 100644 public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFAlleleClipper.java diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PoolGenotypeLikelihoodsCalculationModel.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PoolGenotypeLikelihoodsCalculationModel.java index 37b676601..3e0bdd2ea 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PoolGenotypeLikelihoodsCalculationModel.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/PoolGenotypeLikelihoodsCalculationModel.java @@ -90,7 +90,6 @@ public abstract class PoolGenotypeLikelihoodsCalculationModel extends GenotypeLi return new VariantContextBuilder("pc",referenceSampleVC.getChr(), referenceSampleVC.getStart(), referenceSampleVC.getEnd(), referenceSampleVC.getAlleles()) - .referenceBaseForIndel(referenceSampleVC.getReferenceBaseForIndel()) .genotypes(new GenotypeBuilder(UAC.referenceSampleName, referenceAlleles).GQ(referenceGenotype.getGQ()).make()) .make(); } diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java index e2445e926..ad468f657 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngine.java @@ -33,7 +33,6 @@ import org.apache.commons.lang.ArrayUtils; import org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyperEngine; import org.broadinstitute.sting.gatk.walkers.genotyper.VariantCallContext; import org.broadinstitute.sting.utils.*; -import org.broadinstitute.sting.utils.codecs.vcf.VCFAlleleClipper; import org.broadinstitute.sting.utils.collections.Pair; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; import org.broadinstitute.sting.utils.variantcontext.*; @@ -419,8 +418,8 @@ public class GenotypingEngine { protected static VariantContext createMergedVariantContext( final VariantContext thisVC, final VariantContext nextVC, final byte[] ref, final GenomeLoc refLoc ) { final int thisStart = thisVC.getStart(); final int nextStart = nextVC.getStart(); - byte[] refBases = ( thisVC.hasReferenceBaseForIndel() ? new byte[]{ thisVC.getReferenceBaseForIndel() } : new byte[]{} ); - byte[] altBases = ( thisVC.hasReferenceBaseForIndel() ? new byte[]{ thisVC.getReferenceBaseForIndel() } : new byte[]{} ); + byte[] refBases = ( new byte[]{} ); + byte[] altBases = ( new byte[]{} ); refBases = ArrayUtils.addAll(refBases, thisVC.getReference().getBases()); altBases = ArrayUtils.addAll(altBases, thisVC.getAlternateAllele(0).getBases()); for( int locus = thisStart + refBases.length; locus < nextStart; locus++ ) { @@ -428,15 +427,11 @@ public class GenotypingEngine { refBases = ArrayUtils.add(refBases, refByte); altBases = ArrayUtils.add(altBases, refByte); } - if( nextVC.hasReferenceBaseForIndel() ) { - refBases = ArrayUtils.add(refBases, nextVC.getReferenceBaseForIndel()); - altBases = ArrayUtils.add(altBases, nextVC.getReferenceBaseForIndel()); - } refBases = ArrayUtils.addAll(refBases, nextVC.getReference().getBases()); altBases = ArrayUtils.addAll(altBases, nextVC.getAlternateAllele(0).getBases()); int iii = 0; - if( refBases.length == altBases.length && VCFAlleleClipper.needsPadding(thisVC) ) { // special case of insertion + deletion of same length creates an MNP --> trim padding bases off the allele + if( refBases.length == altBases.length ) { // special case of insertion + deletion of same length creates an MNP --> trim padding bases off the allele while( iii < refBases.length && refBases[iii] == altBases[iii] ) { iii++; } } final ArrayList mergedAlleles = new ArrayList(); @@ -530,10 +525,10 @@ public class GenotypingEngine { final int elementLength = ce.getLength(); switch( ce.getOperator() ) { case I: - final byte[] insertionBases = Arrays.copyOfRange( alignment, alignmentPos, alignmentPos + elementLength ); + final byte[] insertionBases = Arrays.copyOfRange( alignment, alignmentPos - 1, alignmentPos + elementLength ); // add padding base boolean allN = true; - for( final byte b : insertionBases ) { - if( b != (byte) 'N' ) { + for( int i = 1; i < insertionBases.length; i++ ) { // check all bases except for the padding base + if( insertionBases[i] != (byte) 'N' ) { allN = false; break; } @@ -541,14 +536,13 @@ public class GenotypingEngine { if( !allN ) { final ArrayList insertionAlleles = new ArrayList(); final int insertionStart = refLoc.getStart() + refPos - 1; + insertionAlleles.add( Allele.create(ref[refPos-1], true) ); if( haplotype != null && (haplotype.leftBreakPoint + alignmentStartHapwrtRef + refLoc.getStart() - 1 == insertionStart + elementLength + 1 || haplotype.rightBreakPoint + alignmentStartHapwrtRef + refLoc.getStart() - 1 == insertionStart + elementLength + 1) ) { - insertionAlleles.add( Allele.create(ref[refPos-1], true) ); insertionAlleles.add( SYMBOLIC_UNASSEMBLED_EVENT_ALLELE ); vcs.put(insertionStart, new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), insertionStart, insertionStart, insertionAlleles).make()); } else { - insertionAlleles.add( Allele.create(Allele.NULL_ALLELE_STRING, true) ); insertionAlleles.add( Allele.create(insertionBases, false) ); - vcs.put(insertionStart, new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), insertionStart, insertionStart, insertionAlleles).referenceBaseForIndel(ref[refPos-1]).make()); + vcs.put(insertionStart, new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), insertionStart, insertionStart, insertionAlleles).make()); } } @@ -558,7 +552,7 @@ public class GenotypingEngine { alignmentPos += elementLength; break; case D: - final byte[] deletionBases = Arrays.copyOfRange( ref, refPos, refPos + elementLength ); + final byte[] deletionBases = Arrays.copyOfRange( ref, refPos - 1, refPos + elementLength ); // add padding base final ArrayList deletionAlleles = new ArrayList(); final int deletionStart = refLoc.getStart() + refPos - 1; // BUGBUG: how often does this symbolic deletion allele case happen? @@ -569,8 +563,8 @@ public class GenotypingEngine { // vcs.put(deletionStart, new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), deletionStart, deletionStart, deletionAlleles).make()); //} else { deletionAlleles.add( Allele.create(deletionBases, true) ); - deletionAlleles.add( Allele.create(Allele.NULL_ALLELE_STRING, false) ); - vcs.put(deletionStart, new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), deletionStart, deletionStart + elementLength, deletionAlleles).referenceBaseForIndel(ref[refPos-1]).make()); + deletionAlleles.add( Allele.create(ref[refPos-1], false) ); + vcs.put(deletionStart, new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), deletionStart, deletionStart + elementLength, deletionAlleles).make()); //} refPos += elementLength; break; diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngineUnitTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngineUnitTest.java index 04bb3a753..4bcf5a0a0 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngineUnitTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/GenotypingEngineUnitTest.java @@ -262,8 +262,6 @@ public class GenotypingEngineUnitTest extends BaseTest { Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // SNP + ref + SNP = MNP with ref base gap thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("T","G").make(); @@ -274,11 +272,9 @@ public class GenotypingEngineUnitTest extends BaseTest { Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // insertion + SNP - thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("-","AAAAA").referenceBaseForIndel("T").make(); + thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("T","TAAAAA").make(); nextVC = new VariantContextBuilder().loc("2", 1705, 1705).alleles("C","G").make(); truthVC = new VariantContextBuilder().loc("2", 1703, 1705).alleles("TCC","TAAAAACG").source("merged").make(); mergedVC = GenotypingEngine.createMergedVariantContext(thisVC, nextVC, ref, refLoc); @@ -286,23 +282,19 @@ public class GenotypingEngineUnitTest extends BaseTest { Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // SNP + insertion thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("T","G").make(); - nextVC = new VariantContextBuilder().loc("2", 1705, 1705).alleles("-","AAAAA").referenceBaseForIndel("C").make(); + nextVC = new VariantContextBuilder().loc("2", 1705, 1705).alleles("C","CAAAAA").make(); truthVC = new VariantContextBuilder().loc("2", 1703, 1705).alleles("TCC","GCCAAAAA").source("merged").make(); mergedVC = GenotypingEngine.createMergedVariantContext(thisVC, nextVC, ref, refLoc); logger.warn(truthVC + " == " + mergedVC); Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // deletion + SNP - thisVC = new VariantContextBuilder().loc("2", 1703, 1704).alleles("C","-").referenceBaseForIndel("T").make(); + thisVC = new VariantContextBuilder().loc("2", 1703, 1704).alleles("TC","T").make(); nextVC = new VariantContextBuilder().loc("2", 1705, 1705).alleles("C","G").make(); truthVC = new VariantContextBuilder().loc("2", 1703, 1705).alleles("TCC","TG").source("merged").make(); mergedVC = GenotypingEngine.createMergedVariantContext(thisVC, nextVC, ref, refLoc); @@ -310,68 +302,56 @@ public class GenotypingEngineUnitTest extends BaseTest { Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // SNP + deletion thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("T","G").make(); - nextVC = new VariantContextBuilder().loc("2", 1705, 1706).alleles("G","-").referenceBaseForIndel("C").make(); + nextVC = new VariantContextBuilder().loc("2", 1705, 1706).alleles("CG","C").make(); truthVC = new VariantContextBuilder().loc("2", 1703, 1706).alleles("TCCG","GCC").source("merged").make(); mergedVC = GenotypingEngine.createMergedVariantContext(thisVC, nextVC, ref, refLoc); logger.warn(truthVC + " == " + mergedVC); Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // insertion + deletion = MNP - thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("-","A").referenceBaseForIndel("T").make(); - nextVC = new VariantContextBuilder().loc("2", 1705, 1706).alleles("G","-").referenceBaseForIndel("C").make(); + thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("T","TA").make(); + nextVC = new VariantContextBuilder().loc("2", 1705, 1706).alleles("CG","C").make(); truthVC = new VariantContextBuilder().loc("2", 1704, 1706).alleles("CCG","ACC").source("merged").make(); mergedVC = GenotypingEngine.createMergedVariantContext(thisVC, nextVC, ref, refLoc); logger.warn(truthVC + " == " + mergedVC); Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // insertion + deletion - thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("-","AAAAA").referenceBaseForIndel("T").make(); - nextVC = new VariantContextBuilder().loc("2", 1705, 1706).alleles("G","-").referenceBaseForIndel("C").make(); + thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("T","TAAAAA").make(); + nextVC = new VariantContextBuilder().loc("2", 1705, 1706).alleles("CG","C").make(); truthVC = new VariantContextBuilder().loc("2", 1703, 1706).alleles("TCCG","TAAAAACC").source("merged").make(); mergedVC = GenotypingEngine.createMergedVariantContext(thisVC, nextVC, ref, refLoc); logger.warn(truthVC + " == " + mergedVC); Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // insertion + insertion - thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("-","A").referenceBaseForIndel("T").make(); - nextVC = new VariantContextBuilder().loc("2", 1705, 1705).alleles("-","A").referenceBaseForIndel("C").make(); + thisVC = new VariantContextBuilder().loc("2", 1703, 1703).alleles("T","TA").make(); + nextVC = new VariantContextBuilder().loc("2", 1705, 1705).alleles("C","CA").make(); truthVC = new VariantContextBuilder().loc("2", 1703, 1705).alleles("TCC","TACCA").source("merged").make(); mergedVC = GenotypingEngine.createMergedVariantContext(thisVC, nextVC, ref, refLoc); logger.warn(truthVC + " == " + mergedVC); Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // deletion + deletion - thisVC = new VariantContextBuilder().loc("2", 1701, 1702).alleles("T","-").referenceBaseForIndel("A").make(); - nextVC = new VariantContextBuilder().loc("2", 1705, 1706).alleles("G","-").referenceBaseForIndel("C").make(); + thisVC = new VariantContextBuilder().loc("2", 1701, 1702).alleles("AT","A").make(); + nextVC = new VariantContextBuilder().loc("2", 1705, 1706).alleles("CG","C").make(); truthVC = new VariantContextBuilder().loc("2", 1701, 1706).alleles("ATTCCG","ATCC").source("merged").make(); mergedVC = GenotypingEngine.createMergedVariantContext(thisVC, nextVC, ref, refLoc); logger.warn(truthVC + " == " + mergedVC); Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); // complex + complex thisVC = new VariantContextBuilder().loc("2", 1703, 1704).alleles("TC","AAA").make(); @@ -382,8 +362,6 @@ public class GenotypingEngineUnitTest extends BaseTest { Assert.assertTrue(truthVC.hasSameAllelesAs(mergedVC)); Assert.assertEquals(truthVC.getStart(), mergedVC.getStart()); Assert.assertEquals(truthVC.getEnd(), mergedVC.getEnd()); - Assert.assertEquals(truthVC.hasReferenceBaseForIndel(), mergedVC.hasReferenceBaseForIndel()); - Assert.assertEquals(truthVC.getReferenceBaseForIndel(), mergedVC.getReferenceBaseForIndel()); } /** diff --git a/public/java/src/org/broadinstitute/sting/gatk/refdata/VariantContextAdaptors.java b/public/java/src/org/broadinstitute/sting/gatk/refdata/VariantContextAdaptors.java index fe069c2d9..dd1eea8a4 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/refdata/VariantContextAdaptors.java +++ b/public/java/src/org/broadinstitute/sting/gatk/refdata/VariantContextAdaptors.java @@ -163,43 +163,45 @@ public class VariantContextAdaptors { @Override public VariantContext convert(String name, Object input, ReferenceContext ref) { OldDbSNPFeature dbsnp = (OldDbSNPFeature)input; - if ( ! Allele.acceptableAlleleBases(dbsnp.getNCBIRefBase()) ) + + int index = dbsnp.getStart() - ref.getWindow().getStart() - 1; + if ( index < 0 ) + return null; // we weren't given enough reference context to create the VariantContext + + final byte refBaseForIndel = ref.getBases()[index]; + + Allele refAllele; + if ( dbsnp.getNCBIRefBase().equals("-") ) + refAllele = Allele.create(refBaseForIndel); + else if ( ! Allele.acceptableAlleleBases(dbsnp.getNCBIRefBase()) ) return null; - Allele refAllele = Allele.create(dbsnp.getNCBIRefBase(), true); + else + refAllele = Allele.create(refBaseForIndel + dbsnp.getNCBIRefBase(), true); - if ( isSNP(dbsnp) || isIndel(dbsnp) || isMNP(dbsnp) || dbsnp.getVariantType().contains("mixed") ) { - // add the reference allele - List alleles = new ArrayList(); - alleles.add(refAllele); - - // add all of the alt alleles - boolean sawNullAllele = refAllele.isNull(); - for ( String alt : getAlternateAlleleList(dbsnp) ) { - if ( ! Allele.acceptableAlleleBases(alt) ) { - //System.out.printf("Excluding dbsnp record %s%n", dbsnp); - return null; - } - Allele altAllele = Allele.create(alt, false); - alleles.add(altAllele); - if ( altAllele.isNull() ) - sawNullAllele = true; - } - - Map attributes = new HashMap(); - - int index = dbsnp.getStart() - ref.getWindow().getStart() - 1; - if ( index < 0 ) - return null; // we weren't given enough reference context to create the VariantContext - Byte refBaseForIndel = new Byte(ref.getBases()[index]); - - final VariantContextBuilder builder = new VariantContextBuilder(); - builder.source(name).id(dbsnp.getRsID()); - builder.loc(dbsnp.getChr(), dbsnp.getStart() - (sawNullAllele ? 1 : 0), dbsnp.getEnd() - (refAllele.isNull() ? 1 : 0)); - builder.alleles(alleles); - builder.referenceBaseForIndel(refBaseForIndel); - return builder.make(); - } else + boolean addPaddingBase; + if ( isSNP(dbsnp) || isMNP(dbsnp) ) + addPaddingBase = false; + else if ( isIndel(dbsnp) || dbsnp.getVariantType().contains("mixed") ) + addPaddingBase = true; + else return null; // can't handle anything else + + final List alleles = new ArrayList(); + alleles.add(refAllele); + + // add all of the alt alleles + for ( String alt : getAlternateAlleleList(dbsnp) ) { + if ( ! Allele.acceptableAlleleBases(alt) ) { + return null; + } + alleles.add(Allele.create((addPaddingBase ? refBaseForIndel : "") + alt, false)); + } + + final VariantContextBuilder builder = new VariantContextBuilder(); + builder.source(name).id(dbsnp.getRsID()); + builder.loc(dbsnp.getChr(), dbsnp.getStart() - (addPaddingBase ? 1 : 0), dbsnp.getEnd() - (addPaddingBase && refAllele.length() == 1 ? 1 : 0)); + builder.alleles(alleles); + return builder.make(); } } @@ -351,7 +353,7 @@ public class VariantContextAdaptors { long end = hapmap.getEnd(); if ( deletionLength > 0 ) end += deletionLength; - VariantContext vc = new VariantContextBuilder(name, hapmap.getChr(), hapmap.getStart(), end, alleles).id(hapmap.getName()).genotypes(genotypes).referenceBaseForIndel(refBaseForIndel).make(); + VariantContext vc = new VariantContextBuilder(name, hapmap.getChr(), hapmap.getStart(), end, alleles).id(hapmap.getName()).genotypes(genotypes).make(); return vc; } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java index 523aa81b1..261f6433b 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/annotator/DepthPerAlleleBySample.java @@ -42,10 +42,6 @@ import java.util.List; */ public class DepthPerAlleleBySample extends GenotypeAnnotation implements StandardAnnotation { - private static final String REF_ALLELE = "REF"; - - private static final String DEL = "DEL"; // constant, for speed: no need to create a key string for deletion allele every time - public void annotate(RefMetaDataTracker tracker, AnnotatorCompatible walker, ReferenceContext ref, AlignmentContext stratifiedContext, VariantContext vc, Genotype g, GenotypeBuilder gb) { if ( g == null || !g.isCalled() ) return; @@ -53,10 +49,10 @@ public class DepthPerAlleleBySample extends GenotypeAnnotation implements Standa if ( vc.isSNP() ) annotateSNP(stratifiedContext, vc, gb); else if ( vc.isIndel() ) - annotateIndel(stratifiedContext, vc, gb); + annotateIndel(stratifiedContext, ref.getBase(), vc, gb); } - private void annotateSNP(AlignmentContext stratifiedContext, VariantContext vc, GenotypeBuilder gb) { + private void annotateSNP(final AlignmentContext stratifiedContext, final VariantContext vc, final GenotypeBuilder gb) { HashMap alleleCounts = new HashMap(); for ( Allele allele : vc.getAlleles() ) @@ -77,62 +73,47 @@ public class DepthPerAlleleBySample extends GenotypeAnnotation implements Standa gb.AD(counts); } - private void annotateIndel(AlignmentContext stratifiedContext, VariantContext vc, GenotypeBuilder gb) { + private void annotateIndel(final AlignmentContext stratifiedContext, final byte refBase, final VariantContext vc, final GenotypeBuilder gb) { ReadBackedPileup pileup = stratifiedContext.getBasePileup(); if ( pileup == null ) return; - final HashMap alleleCounts = new HashMap(); - alleleCounts.put(REF_ALLELE, 0); + final HashMap alleleCounts = new HashMap(); final Allele refAllele = vc.getReference(); - for ( Allele allele : vc.getAlternateAlleles() ) { - - if ( allele.isNoCall() ) { - continue; // this does not look so good, should we die??? - } - - alleleCounts.put(getAlleleRepresentation(allele), 0); + for ( final Allele allele : vc.getAlleles() ) { + alleleCounts.put(allele, 0); } for ( PileupElement p : pileup ) { if ( p.isBeforeInsertion() ) { - final String b = p.getEventBases(); - if ( alleleCounts.containsKey(b) ) { - alleleCounts.put(b, alleleCounts.get(b)+1); + final Allele insertion = Allele.create(refBase + p.getEventBases(), false); + if ( alleleCounts.containsKey(insertion) ) { + alleleCounts.put(insertion, alleleCounts.get(insertion)+1); } } else if ( p.isBeforeDeletionStart() ) { - if ( p.getEventLength() == refAllele.length() ) { - // this is indeed the deletion allele recorded in VC - final String b = DEL; - if ( alleleCounts.containsKey(b) ) { - alleleCounts.put(b, alleleCounts.get(b)+1); - } + if ( p.getEventLength() == refAllele.length() + 1 ) { + // this is indeed the deletion allele recorded in VC + final Allele deletion = Allele.create(refBase); + if ( alleleCounts.containsKey(deletion) ) { + alleleCounts.put(deletion, alleleCounts.get(deletion)+1); } + } } else if ( p.getRead().getAlignmentEnd() > vc.getStart() ) { - alleleCounts.put(REF_ALLELE, alleleCounts.get(REF_ALLELE)+1); + alleleCounts.put(refAllele, alleleCounts.get(refAllele)+1); } } - int[] counts = new int[alleleCounts.size()]; - counts[0] = alleleCounts.get(REF_ALLELE); + final int[] counts = new int[alleleCounts.size()]; + counts[0] = alleleCounts.get(refAllele); for (int i = 0; i < vc.getAlternateAlleles().size(); i++) - counts[i+1] = alleleCounts.get( getAlleleRepresentation(vc.getAlternateAllele(i)) ); + counts[i+1] = alleleCounts.get( vc.getAlternateAllele(i) ); gb.AD(counts); } - private String getAlleleRepresentation(Allele allele) { - if ( allele.isNull() ) { // deletion wrt the ref - return DEL; - } else { // insertion, pass actual bases - return allele.getBaseString(); - } - - } - // public String getIndelBases() public List getKeyNames() { return Arrays.asList(VCFConstants.GENOTYPE_ALLELE_DEPTHS); } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/BeagleOutputToVCF.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/BeagleOutputToVCF.java index 627d561f6..c8abbfa5a 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/BeagleOutputToVCF.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/BeagleOutputToVCF.java @@ -247,8 +247,6 @@ public class BeagleOutputToVCF extends RodWalker { // Beagle always produces genotype strings based on the strings we input in the likelihood file. String refString = vc_input.getReference().getDisplayString(); - if (refString.length() == 0) // ref was null - refString = Allele.NULL_ALLELE_STRING; Allele bglAlleleA, bglAlleleB; diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/ProduceBeagleInput.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/ProduceBeagleInput.java index 14e92a066..470a1d477 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/ProduceBeagleInput.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/ProduceBeagleInput.java @@ -236,7 +236,7 @@ public class ProduceBeagleInput extends RodWalker { if ( markers != null ) markers.append(marker).append("\t").append(Integer.toString(markerCounter++)).append("\t"); for ( Allele allele : preferredVC.getAlleles() ) { String bglPrintString; - if (allele.isNoCall() || allele.isNull()) + if (allele.isNoCall()) bglPrintString = "-"; else bglPrintString = allele.getBaseString(); // get rid of * in case of reference allele diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/VariantsToBeagleUnphased.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/VariantsToBeagleUnphased.java index 6d83a1d2a..f338f0124 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/VariantsToBeagleUnphased.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/beagle/VariantsToBeagleUnphased.java @@ -146,7 +146,7 @@ public class VariantsToBeagleUnphased extends RodWalker { // write out the alleles at this site for ( Allele allele : vc.getAlleles() ) { - beagleOut.append(allele.isNoCall() || allele.isNull() ? "-" : allele.getBaseString()).append(" "); + beagleOut.append(allele.isNoCall() ? "-" : allele.getBaseString()).append(" "); } // write out sample level genotypes diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java index cef09a913..d2071a9fb 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/ConsensusAlleleCounter.java @@ -246,18 +246,19 @@ public class ConsensusAlleleCounter { // get ref bases of accurate deletion final int startIdxInReference = 1 + loc.getStart() - ref.getWindow().getStart(); stop = loc.getStart() + dLen; - final byte[] refBases = Arrays.copyOfRange(ref.getBases(), startIdxInReference, startIdxInReference + dLen); + final byte[] refBases = Arrays.copyOfRange(ref.getBases(), startIdxInReference - 1, startIdxInReference + dLen); // add reference padding if (Allele.acceptableAlleleBases(refBases, false)) { refAllele = Allele.create(refBases, true); - altAllele = Allele.create(Allele.NULL_ALLELE_STRING, false); + altAllele = Allele.create(ref.getBase(), false); } else continue; // don't go on with this allele if refBases are non-standard } else { // insertion case - if (Allele.acceptableAlleleBases(s, false)) { // don't allow N's in insertions - refAllele = Allele.create(Allele.NULL_ALLELE_STRING, true); - altAllele = Allele.create(s, false); + final String insertionBases = ref.getBase() + s; // add reference padding + if (Allele.acceptableAlleleBases(insertionBases, false)) { // don't allow N's in insertions + refAllele = Allele.create(ref.getBase(), true); + altAllele = Allele.create(insertionBases, false); stop = loc.getStart(); } else continue; // go on to next allele if consensus insertion has any non-standard base. @@ -267,7 +268,6 @@ public class ConsensusAlleleCounter { final VariantContextBuilder builder = new VariantContextBuilder().source(""); builder.loc(loc.getContig(), loc.getStart(), stop); builder.alleles(Arrays.asList(refAllele, altAllele)); - builder.referenceBaseForIndel(ref.getBase()); builder.noGenotypes(); if (doMultiAllelicCalls) { vcs.add(builder.make()); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java index 230d6c324..7eabe7a18 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsCalculationModel.java @@ -123,7 +123,7 @@ public class IndelGenotypeLikelihoodsCalculationModel extends GenotypeLikelihood final int endLoc = computeEndLocation(alleleList, loc,allelesArePadded); final int eventLength = getEventLength(alleleList); - final VariantContextBuilder builder = new VariantContextBuilder("UG_call", loc.getContig(), loc.getStart(), endLoc, alleleList).referenceBaseForIndel(ref.getBase()); + final VariantContextBuilder builder = new VariantContextBuilder("UG_call", loc.getContig(), loc.getStart(), endLoc, alleleList); // create the genotypes; no-call everyone for now GenotypesContext genotypes = GenotypesContext.create(); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java index 32564984a..d4c45e19d 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperEngine.java @@ -37,7 +37,6 @@ import org.broadinstitute.sting.gatk.walkers.annotator.VariantAnnotatorEngine; import org.broadinstitute.sting.utils.*; import org.broadinstitute.sting.utils.baq.BAQ; import org.broadinstitute.sting.utils.classloader.PluginManager; -import org.broadinstitute.sting.utils.codecs.vcf.VCFAlleleClipper; import org.broadinstitute.sting.utils.codecs.vcf.VCFConstants; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; import org.broadinstitute.sting.utils.exceptions.UserException; @@ -283,7 +282,7 @@ public class UnifiedGenotyperEngine { VariantContext vcInput = UnifiedGenotyperEngine.getVCFromAllelesRod(tracker, ref, rawContext.getLocation(), false, logger, UAC.alleles); if ( vcInput == null ) return null; - vc = new VariantContextBuilder("UG_call", ref.getLocus().getContig(), vcInput.getStart(), vcInput.getEnd(), vcInput.getAlleles()).referenceBaseForIndel(vcInput.getReferenceBaseForIndel()).make(); + vc = new VariantContextBuilder("UG_call", ref.getLocus().getContig(), vcInput.getStart(), vcInput.getEnd(), vcInput.getAlleles()).make(); } else { // deal with bad/non-standard reference bases if ( !Allele.acceptableAlleleBases(new byte[]{ref.getBase()}) ) @@ -408,11 +407,6 @@ public class UnifiedGenotyperEngine { builder.log10PError(phredScaledConfidence/-10.0); if ( ! passesCallThreshold(phredScaledConfidence) ) builder.filters(filter); - if ( limitedContext ) { - builder.referenceBaseForIndel(vc.getReferenceBaseForIndel()); - } else { - builder.referenceBaseForIndel(refContext.getBase()); - } // create the genotypes final GenotypesContext genotypes = afcm.get().subsetAlleles(vc, myAlleles, true,ploidy); @@ -491,10 +485,8 @@ public class UnifiedGenotyperEngine { builder.attributes(attributes); VariantContext vcCall = builder.make(); - // if we are subsetting alleles (either because there were too many or because some were not polymorphic) - // then we may need to trim the alleles (because the original VariantContext may have had to pad at the end). - if ( myAlleles.size() != vc.getAlleles().size() && !limitedContext ) // TODO - this function doesn't work with mixed records or records that started as mixed and then became non-mixed - vcCall = VCFAlleleClipper.reverseTrimAlleles(vcCall); + // TODO -- if we are subsetting alleles (either because there were too many or because some were not polymorphic) + // TODO -- then we may need to trim the alleles (because the original VariantContext may have had to pad at the end). if ( annotationEngine != null && !limitedContext ) { // Note: we want to use the *unfiltered* and *unBAQed* context for the annotations diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/indels/SomaticIndelDetector.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/indels/SomaticIndelDetector.java index 21db1412b..0c7e2ec5f 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/indels/SomaticIndelDetector.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/indels/SomaticIndelDetector.java @@ -1128,12 +1128,13 @@ public class SomaticIndelDetector extends ReadWalker { List alleles = new ArrayList(2); // actual observed (distinct!) alleles at the site List homref_alleles = null; // when needed, will contain two identical copies of ref allele - needed to generate hom-ref genotype + final byte referencePaddingBase = refBases[(int)start-1]; if ( call.getVariant() == null ) { - // we will need to cteate genotype with two (hom) ref alleles (below). + // we will need to create genotype with two (hom) ref alleles (below). // we can not use 'alleles' list here, since that list is supposed to contain // only *distinct* alleles observed at the site or VCFContext will frown upon us... - alleles.add( Allele.create(refBases[(int)start-1],true) ); + alleles.add( Allele.create(referencePaddingBase,true) ); homref_alleles = new ArrayList(2); homref_alleles.add( alleles.get(0)); homref_alleles.add( alleles.get(0)); @@ -1142,7 +1143,7 @@ public class SomaticIndelDetector extends ReadWalker { // (Genotype will tell us whether it is an actual call or not!) int event_length = call.getVariant().lengthOnRef(); if ( event_length < 0 ) event_length = 0; - fillAlleleList(alleles,call); + fillAlleleList(alleles,call,referencePaddingBase); stop += event_length; } @@ -1162,7 +1163,7 @@ public class SomaticIndelDetector extends ReadWalker { filters.add("NoCall"); } VariantContext vc = new VariantContextBuilder("IGv2_Indel_call", refName, start, stop, alleles) - .genotypes(genotypes).filters(filters).referenceBaseForIndel(refBases[(int)start-1]).make(); + .genotypes(genotypes).filters(filters).make(); vcf.add(vc); } @@ -1172,16 +1173,16 @@ public class SomaticIndelDetector extends ReadWalker { * @param l * @param call */ - private void fillAlleleList(List l, IndelPrecall call) { + private void fillAlleleList(List l, IndelPrecall call, byte referencePaddingBase) { int event_length = call.getVariant().lengthOnRef(); if ( event_length == 0 ) { // insertion - l.add( Allele.create(Allele.NULL_ALLELE_STRING,true) ); - l.add( Allele.create(call.getVariant().getBases(), false )); + l.add( Allele.create(referencePaddingBase,true) ); + l.add( Allele.create(referencePaddingBase + call.getVariant().getBases(), false )); } else { //deletion: - l.add( Allele.create(call.getVariant().getBases(), true )); - l.add( Allele.create(Allele.NULL_ALLELE_STRING,false) ); + l.add( Allele.create(referencePaddingBase + call.getVariant().getBases(), true )); + l.add( Allele.create(referencePaddingBase,false) ); } } @@ -1215,19 +1216,20 @@ public class SomaticIndelDetector extends ReadWalker { // } boolean homRefT = ( tCall.getVariant() == null ); boolean homRefN = ( nCall.getVariant() == null ); + final byte referencePaddingBase = refBases[(int)start-1]; if ( tCall.getVariant() == null && nCall.getVariant() == null) { // no indel at all ; create base-representation ref/ref alleles for genotype construction - alleles.add( Allele.create(refBases[(int)start-1],true) ); + alleles.add( Allele.create(referencePaddingBase,true) ); } else { // we got indel(s) int event_length = 0; if ( tCall.getVariant() != null ) { // indel in tumor event_length = tCall.getVariant().lengthOnRef(); - fillAlleleList(alleles, tCall); + fillAlleleList(alleles, tCall, referencePaddingBase); } else { event_length = nCall.getVariant().lengthOnRef(); - fillAlleleList(alleles, nCall); + fillAlleleList(alleles, nCall, referencePaddingBase); } if ( event_length > 0 ) stop += event_length; } @@ -1259,7 +1261,7 @@ public class SomaticIndelDetector extends ReadWalker { } VariantContext vc = new VariantContextBuilder("IGv2_Indel_call", refName, start, stop, alleles) - .genotypes(genotypes).filters(filters).attributes(attrs).referenceBaseForIndel(refBases[(int)start-1]).make(); + .genotypes(genotypes).filters(filters).attributes(attrs).make(); vcf.add(vc); } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/GenomeEvent.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/GenomeEvent.java index af6a52002..67ddc47ff 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/GenomeEvent.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/GenomeEvent.java @@ -26,7 +26,6 @@ package org.broadinstitute.sting.gatk.walkers.validation.validationsiteselector; import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.GenomeLocParser; -import org.broadinstitute.sting.utils.codecs.vcf.VCFConstants; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; import org.broadinstitute.sting.utils.variantcontext.Allele; import org.broadinstitute.sting.utils.variantcontext.VariantContext; @@ -40,14 +39,11 @@ public class GenomeEvent implements Comparable { final protected GenomeLoc loc; /** A set of the alleles segregating in this context */ final protected List alleles; - final protected Byte refBase; // final protected HashMap attributes; - public GenomeEvent(GenomeLocParser parser, final String contig, final int start, final int stop, final List alleles, HashMap attributes, - byte base) { + public GenomeEvent(GenomeLocParser parser, final String contig, final int start, final int stop, final List alleles, HashMap attributes) { this.loc = parser.createGenomeLoc(contig, start, stop); this.alleles = alleles; - this.refBase = base; // this.attributes = attributes; } @@ -68,7 +64,7 @@ public class GenomeEvent implements Comparable { public VariantContext createVariantContextFromEvent() { return new VariantContextBuilder("event", loc.getContig(), loc.getStart(), loc.getStop(), alleles) - .log10PError(0.0).referenceBaseForIndel(refBase).make(); + .log10PError(0.0).make(); } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/KeepAFSpectrumFrequencySelector.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/KeepAFSpectrumFrequencySelector.java index 4b68eed2e..7c1d63f02 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/KeepAFSpectrumFrequencySelector.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/KeepAFSpectrumFrequencySelector.java @@ -115,7 +115,7 @@ public class KeepAFSpectrumFrequencySelector extends FrequencyModeSelector { // create bare-bones event and log in corresponding bin // attributes contains AC,AF,AN pulled from original vc, and we keep them here and log in output file for bookkeeping purposes - GenomeEvent event = new GenomeEvent(parser, vc.getChr(), vc.getStart(), vc.getEnd(),vc.getAlleles(), attributes, vc.getReferenceBaseForIndel()); + GenomeEvent event = new GenomeEvent(parser, vc.getChr(), vc.getStart(), vc.getEnd(),vc.getAlleles(), attributes); binnedEventArray[binIndex].add(event); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/UniformSamplingFrequencySelector.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/UniformSamplingFrequencySelector.java index eda75d647..4019c5631 100644 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/UniformSamplingFrequencySelector.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/validation/validationsiteselector/UniformSamplingFrequencySelector.java @@ -65,7 +65,7 @@ public class UniformSamplingFrequencySelector extends FrequencyModeSelector { } // create bare-bones event and log in corresponding bin // attributes contains AC,AF,AN pulled from original vc, and we keep them here and log in output file for bookkeeping purposes - GenomeEvent event = new GenomeEvent(parser, vc.getChr(), vc.getStart(), vc.getEnd(),vc.getAlleles(), attributes, vc.getReferenceBaseForIndel()); + GenomeEvent event = new GenomeEvent(parser, vc.getChr(), vc.getStart(), vc.getEnd(),vc.getAlleles(), attributes); binnedEventArray.add(event); } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/evaluators/ThetaVariantEvaluator.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/evaluators/ThetaVariantEvaluator.java index 88bf3aef9..a509294ff 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/evaluators/ThetaVariantEvaluator.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/varianteval/evaluators/ThetaVariantEvaluator.java @@ -56,7 +56,7 @@ public class ThetaVariantEvaluator extends VariantEvaluator { //increment stats for pairwise mismatches for (Allele allele : genotype.getAlleles()) { - if (allele.isNonNull() && allele.isCalled()) { + if (allele.isCalled()) { String alleleString = allele.toString(); alleleCounts.putIfAbsent(alleleString, 0); alleleCounts.put(alleleString, alleleCounts.get(alleleString) + 1); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/LeftAlignVariants.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/LeftAlignVariants.java index c1755aa00..3f19e22d9 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/LeftAlignVariants.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/LeftAlignVariants.java @@ -226,6 +226,6 @@ public class LeftAlignVariants extends RodWalker { newGenotypes.add(new GenotypeBuilder(genotype).alleles(newAlleles).make()); } - return new VariantContextBuilder(vc).alleles(alleleMap.values()).genotypes(newGenotypes).referenceBaseForIndel(refBaseForIndel).make(); + return new VariantContextBuilder(vc).alleles(alleleMap.values()).genotypes(newGenotypes).make(); } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/LiftoverVariants.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/LiftoverVariants.java index 60d41abd5..094897edc 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/LiftoverVariants.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/LiftoverVariants.java @@ -116,7 +116,6 @@ public class LiftoverVariants extends RodWalker { if ( toInterval != null ) { // check whether the strand flips, and if so reverse complement everything - // TODO -- make this work for indels (difficult because the 'previous base' context needed will be changing based on indel type/size) if ( fromInterval.isPositiveStrand() != toInterval.isPositiveStrand() && vc.isPointEvent() ) { vc = VariantContextUtils.reverseComplement(vc); } @@ -129,11 +128,10 @@ public class LiftoverVariants extends RodWalker { .attribute("OriginalStart", fromInterval.getStart()).make(); } - VariantContext newVC = VCFAlleleClipper.createVariantContextWithPaddedAlleles(vc); - if ( originalVC.isSNP() && originalVC.isBiallelic() && VariantContextUtils.getSNPSubstitutionType(originalVC) != VariantContextUtils.getSNPSubstitutionType(newVC) ) { + if ( originalVC.isSNP() && originalVC.isBiallelic() && VariantContextUtils.getSNPSubstitutionType(originalVC) != VariantContextUtils.getSNPSubstitutionType(vc) ) { logger.warn(String.format("VCF at %s / %d => %s / %d is switching substitution type %s/%s to %s/%s", - originalVC.getChr(), originalVC.getStart(), newVC.getChr(), newVC.getStart(), - originalVC.getReference(), originalVC.getAlternateAllele(0), newVC.getReference(), newVC.getAlternateAllele(0))); + originalVC.getChr(), originalVC.getStart(), vc.getChr(), vc.getStart(), + originalVC.getReference(), originalVC.getAlternateAllele(0), vc.getReference(), vc.getAlternateAllele(0))); } writer.add(vc); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/ValidateVariants.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/ValidateVariants.java index 530258fe0..995e98931 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/ValidateVariants.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/ValidateVariants.java @@ -127,35 +127,16 @@ public class ValidateVariants extends RodWalker { return; // get the true reference allele - Allele reportedRefAllele = vc.getReference(); - Allele observedRefAllele = null; - // insertions - if ( vc.isSimpleInsertion() ) { - observedRefAllele = Allele.create(Allele.NULL_ALLELE_STRING); + final Allele reportedRefAllele = vc.getReference(); + final int refLength = reportedRefAllele.length(); + if ( refLength > 100 ) { + logger.info(String.format("Reference allele is too long (%d) at position %s:%d; skipping that record.", refLength, vc.getChr(), vc.getStart())); + return; } - // deletions - else if ( vc.isSimpleDeletion() || vc.isMNP() ) { - // we can't validate arbitrarily long deletions - if ( reportedRefAllele.length() > 100 ) { - logger.info(String.format("Reference allele is too long (%d) at position %s:%d; skipping that record.", reportedRefAllele.length(), vc.getChr(), vc.getStart())); - return; - } - // deletions are associated with the (position of) the last (preceding) non-deleted base; - // hence to get actually deleted bases we need offset = 1 - int offset = vc.isMNP() ? 0 : 1; - byte[] refBytes = ref.getBases(); - byte[] trueRef = new byte[reportedRefAllele.length()]; - for (int i = 0; i < reportedRefAllele.length(); i++) - trueRef[i] = refBytes[i+offset]; - observedRefAllele = Allele.create(trueRef, true); - } - // SNPs, etc. but not mixed types because they are too difficult - else if ( !vc.isMixed() ) { - byte[] refByte = new byte[1]; - refByte[0] = ref.getBase(); - observedRefAllele = Allele.create(refByte, true); - } + final byte[] observedRefBases = new byte[refLength]; + System.arraycopy(ref.getBases(), 0, observedRefBases, 0, refLength); + final Allele observedRefAllele = Allele.create(observedRefBases); // get the RS IDs Set rsIDs = null; @@ -168,10 +149,10 @@ public class ValidateVariants extends RodWalker { try { switch( type ) { case ALL: - vc.extraStrictValidation(observedRefAllele, ref.getBase(), rsIDs); + vc.extraStrictValidation(reportedRefAllele, observedRefAllele, rsIDs); break; case REF: - vc.validateReferenceBases(observedRefAllele, ref.getBase()); + vc.validateReferenceBases(reportedRefAllele, observedRefAllele); break; case IDS: vc.validateRSIDs(rsIDs); diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantsToTable.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantsToTable.java index 996ac75e7..4806b2ebc 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantsToTable.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantsToTable.java @@ -378,7 +378,7 @@ public class VariantsToTable extends RodWalker { getters.put("REF", new Getter() { public String get(VariantContext vc) { StringBuilder x = new StringBuilder(); - x.append(vc.getAlleleStringWithRefPadding(vc.getReference())); + x.append(vc.getReference()); return x.toString(); } }); @@ -390,7 +390,7 @@ public class VariantsToTable extends RodWalker { for ( int i = 0; i < n; i++ ) { if ( i != 0 ) x.append(","); - x.append(vc.getAlleleStringWithRefPadding(vc.getAlternateAllele(i))); + x.append(vc.getAlternateAllele(i)); } return x.toString(); } @@ -432,11 +432,8 @@ public class VariantsToTable extends RodWalker { private static Object splitAltAlleles(VariantContext vc) { final int numAltAlleles = vc.getAlternateAlleles().size(); if ( numAltAlleles == 1 ) - return vc.getAlleleStringWithRefPadding(vc.getAlternateAllele(0)); + return vc.getAlternateAllele(0); - final List alleles = new ArrayList(numAltAlleles); - for ( Allele allele : vc.getAlternateAlleles() ) - alleles.add(vc.getAlleleStringWithRefPadding(allele)); - return alleles; + return vc.getAlternateAlleles(); } } diff --git a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantsToVCF.java b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantsToVCF.java index e8c6794f2..cf568a62e 100755 --- a/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantsToVCF.java +++ b/public/java/src/org/broadinstitute/sting/gatk/walkers/variantutils/VariantsToVCF.java @@ -100,12 +100,6 @@ public class VariantsToVCF extends RodWalker { @Argument(fullName="sample", shortName="sample", doc="The sample name represented by the variant rod", required=false) protected String sampleName = null; - /** - * This argument is useful for fixing input VCFs with bad reference bases (the output will be a fixed version of the VCF). - */ - @Argument(fullName="fixRef", shortName="fixRef", doc="Fix common reference base in case there's an indel without padding", required=false) - protected boolean fixReferenceBase = false; - private Set allowedGenotypeFormatStrings = new HashSet(); private boolean wroteHeader = false; private Set samples; @@ -137,10 +131,6 @@ public class VariantsToVCF extends RodWalker { builder.genotypes(g); } - if ( fixReferenceBase ) { - builder.referenceBaseForIndel(ref.getBase()); - } - writeRecord(builder.make(), tracker, ref.getLocus()); } @@ -166,8 +156,8 @@ public class VariantsToVCF extends RodWalker { continue; Map alleleMap = new HashMap(2); - alleleMap.put(RawHapMapFeature.DELETION, Allele.create(Allele.NULL_ALLELE_STRING, dbsnpVC.isSimpleInsertion())); - alleleMap.put(RawHapMapFeature.INSERTION, Allele.create(((RawHapMapFeature)record).getAlleles()[1], !dbsnpVC.isSimpleInsertion())); + alleleMap.put(RawHapMapFeature.DELETION, Allele.create(ref.getBase(), dbsnpVC.isSimpleInsertion())); + alleleMap.put(RawHapMapFeature.INSERTION, Allele.create(ref.getBase() + ((RawHapMapFeature)record).getAlleles()[1], !dbsnpVC.isSimpleInsertion())); hapmap.setActualAlleles(alleleMap); // also, use the correct positioning for insertions diff --git a/public/java/src/org/broadinstitute/sting/utils/BaseUtils.java b/public/java/src/org/broadinstitute/sting/utils/BaseUtils.java index 393dd5735..0065f9258 100644 --- a/public/java/src/org/broadinstitute/sting/utils/BaseUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/BaseUtils.java @@ -431,6 +431,37 @@ public class BaseUtils { return new String(simpleComplement(bases.getBytes())); } + /** + * Returns the uppercased version of the bases + * + * @param bases the bases + * @return the upper cased version + */ + static public byte[] convertToUpperCase(final byte[] bases) { + for ( int i = 0; i < bases.length; i++ ) { + if ( (char)bases[i] >= 'a' ) + bases[i] = toUpperCaseBase(bases[i]); + } + return bases; + } + + static public byte toUpperCaseBase(final byte base) { + switch (base) { + case 'a': + return 'A'; + case 'c': + return 'C'; + case 'g': + return 'G'; + case 't': + return 'T'; + case 'n': + return 'N'; + default: + return base; + } + } + /** * Returns the index of the most common base in the basecounts array. To be used with * pileup.getBaseCounts. diff --git a/public/java/src/org/broadinstitute/sting/utils/codecs/bcf2/BCF2Codec.java b/public/java/src/org/broadinstitute/sting/utils/codecs/bcf2/BCF2Codec.java index 0b9654610..0f9cc34e7 100644 --- a/public/java/src/org/broadinstitute/sting/utils/codecs/bcf2/BCF2Codec.java +++ b/public/java/src/org/broadinstitute/sting/utils/codecs/bcf2/BCF2Codec.java @@ -305,27 +305,6 @@ public final class BCF2Codec implements FeatureCodec { builder.id(id); } - /** - * Annoying routine that deals with allele clipping from the BCF2 encoding to the standard - * GATK encoding. - * - * @param position - * @param ref - * @param unclippedAlleles - * @return - */ - @Requires({"position > 0", "ref != null && ref.length() > 0", "! unclippedAlleles.isEmpty()"}) - @Ensures("result.size() == unclippedAlleles.size()") - protected List clipAllelesIfNecessary(final int position, - final String ref, - final List unclippedAlleles) { - // the last argument of 1 allows us to safely ignore the end, because we are - // ultimately going to use the end in the record itself - final VCFAlleleClipper.ClippedAlleles clipped = VCFAlleleClipper.clipAlleles(position, ref, unclippedAlleles, 1); - if ( clipped.getError() != null ) error(clipped.getError()); - return clipped.getClippedAlleles(); - } - /** * Decode the alleles from this BCF2 file and put the results in builder * @param builder @@ -353,11 +332,9 @@ public final class BCF2Codec implements FeatureCodec { } assert ref != null; - alleles = clipAllelesIfNecessary(pos, ref, alleles); builder.alleles(alleles); assert ref.length() > 0; - builder.referenceBaseForIndel(ref.getBytes()[0]); return alleles; } diff --git a/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/AbstractVCFCodec.java b/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/AbstractVCFCodec.java index b3420514b..2b5695e3a 100755 --- a/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/AbstractVCFCodec.java +++ b/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/AbstractVCFCodec.java @@ -248,6 +248,7 @@ public abstract class AbstractVCFCodec extends AsciiFeatureCodec builder.id(parts[2]); final String ref = getCachedString(parts[3].toUpperCase()); + builder.stop(pos + ref.length() - 1); final String alts = getCachedString(parts[4].toUpperCase()); builder.log10PError(parseQual(parts[5])); @@ -257,8 +258,8 @@ public abstract class AbstractVCFCodec extends AsciiFeatureCodec builder.attributes(attrs); // get our alleles, filters, and setup an attribute map - final List rawAlleles = parseAlleles(ref, alts, lineNo); - final List alleles = updateBuilderAllelesAndStop(builder, ref, pos, rawAlleles, attrs); + final List alleles = parseAlleles(ref, alts, lineNo); + builder.alleles(alleles); // do we have genotyping data if (parts.length > NUM_STANDARD_FIELDS && includeGenotypes) { @@ -275,7 +276,6 @@ public abstract class AbstractVCFCodec extends AsciiFeatureCodec VariantContext vc = null; try { - builder.referenceBaseForIndel(ref.getBytes()[0]); vc = builder.make(); } catch (Exception e) { generateException(e.getMessage()); @@ -284,31 +284,6 @@ public abstract class AbstractVCFCodec extends AsciiFeatureCodec return vc; } - private final List updateBuilderAllelesAndStop(final VariantContextBuilder builder, - final String ref, - final int pos, - final List rawAlleles, - final Map attrs) { - int endForSymbolicAlleles = pos; // by default we use the pos - if ( attrs.containsKey(VCFConstants.END_KEY) ) { - // update stop with the end key if provided - try { - endForSymbolicAlleles = Integer.valueOf(attrs.get(VCFConstants.END_KEY).toString()); - } catch (Exception e) { - generateException("the END value in the INFO field is not valid"); - } - } - - // find out our current location, and clip the alleles down to their minimum length - final VCFAlleleClipper.ClippedAlleles clipped = VCFAlleleClipper.clipAlleles(pos, ref, rawAlleles, endForSymbolicAlleles); - if ( clipped.getError() != null ) - generateException(clipped.getError(), lineNo); - - builder.stop(clipped.getStop()); - builder.alleles(clipped.getClippedAlleles()); - return clipped.getClippedAlleles(); - } - /** * get the name of this codec * @return our set name diff --git a/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFAlleleClipper.java b/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFAlleleClipper.java deleted file mode 100644 index 40ba23d9d..000000000 --- a/public/java/src/org/broadinstitute/sting/utils/codecs/vcf/VCFAlleleClipper.java +++ /dev/null @@ -1,434 +0,0 @@ -/* - * Copyright (c) 2012, The Broad Institute - * - * Permission is hereby granted, free of charge, to any person - * obtaining a copy of this software and associated documentation - * files (the "Software"), to deal in the Software without - * restriction, including without limitation the rights to use, - * copy, modify, merge, publish, distribute, sublicense, and/or sell - * copies of the Software, and to permit persons to whom the - * Software is furnished to do so, subject to the following - * conditions: - * - * The above copyright notice and this permission notice shall be - * included in all copies or substantial portions of the Software. - * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES - * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND - * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT - * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, - * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING - * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR - * OTHER DEALINGS IN THE SOFTWARE. - */ - -package org.broadinstitute.sting.utils.codecs.vcf; - -import com.google.java.contract.Ensures; -import com.google.java.contract.Invariant; -import com.google.java.contract.Requires; -import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; -import org.broadinstitute.sting.utils.variantcontext.*; - -import java.util.*; - -/** - * All of the gross allele clipping and padding routines in one place - * - * Having attempted to understand / fix / document this code myself - * I can only conclude that this entire approach needs to be rethought. This - * code just doesn't work robustly with symbolic alleles, with multiple alleles, - * requires a special "reference base for indels" stored in the VariantContext - * whose correctness isn't enforced, and overall has strange special cases - * all over the place. - * - * The reason this code is so complex is due to symbolics and multi-alleleic - * variation, which frequently occur when combining variants from multiple - * VCF files. - * - * TODO rethink this class, make it clean, and make it easy to create, mix, and write out alleles - * TODO this code doesn't work with reverse clipped alleles (ATA / GTTA -> AT / GT) - * - * @author Mark DePristo - * @since 6/12 - */ -public final class VCFAlleleClipper { - private VCFAlleleClipper() { } - - /** - * Determine whether we should clip off the first base of all unclippped alleles or not - * - * Returns true if all of the alleles in unclippedAlleles share a common first base with - * ref0. Ref0 should be the first base of the reference allele UnclippedAlleles may - * contain the reference allele itself, or just the alternate alleles, it doesn't matter. - * - * The algorithm returns true if the first base should be clipped off, or false otherwise - * - * This algorithm works even in the presence of symbolic alleles, logically ignoring these - * values. It - * - * @param unclippedAlleles list of unclipped alleles to assay - * @param ref0 the first base of the reference allele - * @return true if we should clip the first base of unclippedAlleles - */ - @Requires("unclippedAlleles != null") - public static boolean shouldClipFirstBaseP(final List unclippedAlleles, - final byte ref0) { - boolean allSymbolicAlt = true; - - for ( final Allele a : unclippedAlleles ) { - if ( a.isSymbolic() ) { - continue; - } - - // already know we aren't symbolic, so we only need to decide if we have only seen a ref - if ( ! a.isReference() ) - allSymbolicAlt = false; - - if ( a.length() < 1 || (a.getBases()[0] != ref0) ) { - return false; - } - } - - // to reach here all alleles are consistent with clipping the first base matching ref0 - // but we don't clip if all ALT alleles are symbolic - return ! allSymbolicAlt; - } - - public static int computeReverseClipping(final List unclippedAlleles, - final byte[] ref, - final int forwardClipping, - final boolean allowFullClip) { - int clipping = 0; - boolean stillClipping = true; - - while ( stillClipping ) { - for ( final Allele a : unclippedAlleles ) { - if ( a.isSymbolic() ) - continue; - - // we need to ensure that we don't reverse clip out all of the bases from an allele because we then will have the wrong - // position set for the VariantContext (although it's okay to forward clip it all out, because the position will be fine). - if ( a.length() - clipping == 0 ) - return clipping - (allowFullClip ? 0 : 1); - - if ( a.length() - clipping <= forwardClipping || a.length() - forwardClipping == 0 ) { - stillClipping = false; - } - else if ( ref.length == clipping ) { - if ( allowFullClip ) - stillClipping = false; - else - return -1; - } - else if ( a.getBases()[a.length()-clipping-1] != ref[ref.length-clipping-1] ) { - stillClipping = false; - } - } - if ( stillClipping ) - clipping++; - } - - return clipping; - } - - /** - * Are the alleles describing a polymorphism substitution one base for another? - * - * @param alleles a list of alleles, must not be empty - * @return Return true if the length of any allele in alleles isn't 1 - */ - @Requires("!alleles.isEmpty()") - private static boolean isSingleNucleotideEvent(final List alleles) { - for ( final Allele a : alleles ) { - if ( a.length() != 1 ) - return false; - } - return true; - } - - /** - * clip the alleles, based on the reference, returning a ClippedAlleles object describing what happened - * - * The ClippedAlleles object contains the implied stop position of the alleles, given the provided start - * position, after clipping. It also contains the list of alleles, in the same order as the provided - * unclipped ones, that are the fully clipped version of the input alleles. If an error occurs - * during this option the getError() function returns a string describing the problem (for use in parsers). - * - * The basic operation are: - * - * single allele - * => stop == start and clipped == unclipped - * any number of single nucleotide events - * => stop == start and clipped == unclipped - * two alleles, second being symbolic - * => stop == start and clipped == unclipped - * Note in this case that the STOP should be computed by other means (from END in VCF, for example) - * Note that if there's more than two alleles and the second is a symbolic the code produces an error - * Any other case: - * The alleles are trimmed of any sequence shared at the end of the alleles. If N bases - * are common then the alleles will all be at least N bases shorter. - * The stop position returned is the start position + the length of the - * reverse trimmed only reference allele - 1. - * If the alleles all share a single common starting sequence (just one base is considered) - * then the alleles have this leading common base removed as well. - * - * TODO This code is gross and brittle and needs to be rethought from scratch - * - * @param start the unadjusted start position (pre-clipping) - * @param ref the reference string - * @param unclippedAlleles the list of unclipped alleles, including the reference allele - * @return the new reference end position of this event - */ - @Requires({"start > 0", "ref != null && ref.length() > 0", "!unclippedAlleles.isEmpty()"}) - @Ensures("result != null") - public static ClippedAlleles clipAlleles(final int start, - final String ref, - final List unclippedAlleles, - final int endForSymbolicAllele ) { - // no variation or single nucleotide events are by definition fully clipped - if ( unclippedAlleles.size() == 1 || isSingleNucleotideEvent(unclippedAlleles) ) - return new ClippedAlleles(start, unclippedAlleles, null); - - // we've got to sort out the clipping by looking at the alleles themselves - final byte firstRefBase = (byte) ref.charAt(0); - final boolean firstBaseIsClipped = shouldClipFirstBaseP(unclippedAlleles, firstRefBase); - final int forwardClipping = firstBaseIsClipped ? 1 : 0; - final int reverseClipping = computeReverseClipping(unclippedAlleles, ref.getBytes(), forwardClipping, false); - final boolean needsClipping = forwardClipping > 0 || reverseClipping > 0; - - if ( reverseClipping == -1 ) - return new ClippedAlleles("computeReverseClipping failed due to bad alleles"); - - boolean sawSymbolic = false; - List clippedAlleles; - if ( ! needsClipping ) { - // there's nothing to clip, so clippedAlleles are the original alleles - clippedAlleles = unclippedAlleles; - } else { - clippedAlleles = new ArrayList(unclippedAlleles.size()); - for ( final Allele a : unclippedAlleles ) { - if ( a.isSymbolic() ) { - sawSymbolic = true; - clippedAlleles.add(a); - } else { - final byte[] allele = Arrays.copyOfRange(a.getBases(), forwardClipping, a.getBases().length - reverseClipping); - if ( !Allele.acceptableAlleleBases(allele) ) - return new ClippedAlleles("Unparsable vcf record with bad allele [" + allele + "]"); - clippedAlleles.add(Allele.create(allele, a.isReference())); - } - } - } - - int stop = VariantContextUtils.computeEndFromAlleles(clippedAlleles, start, endForSymbolicAllele); - - // TODO - // TODO - // TODO COMPLETELY BROKEN CODE -- THE GATK CURRENTLY ENCODES THE STOP POSITION FOR CLIPPED ALLELES AS + 1 - // TODO ITS TRUE SIZE TO DIFFERENTIATE CLIPPED VS. UNCLIPPED ALLELES. NEEDS TO BE FIXED - // TODO - // TODO - if ( needsClipping && ! sawSymbolic && ! clippedAlleles.get(0).isNull() ) stop++; - // TODO - // TODO - // TODO COMPLETELY BROKEN CODE -- THE GATK CURRENTLY ENCODES THE STOP POSITION FOR CLIPPED ALLELES AS + 1 - // TODO ITS TRUE SIZE TO DIFFERENTIATE CLIPPED VS. UNCLIPPED ALLELES. NEEDS TO BE FIXED - // TODO - // TODO - - final Byte refBaseForIndel = firstBaseIsClipped ? firstRefBase : null; - return new ClippedAlleles(stop, clippedAlleles, refBaseForIndel); - } - - /** - * Returns true if the alleles in inputVC should have reference bases added for padding - * - * We need to pad a VC with a common base if the length of the reference allele is - * less than the length of the VariantContext. This happens because the position of - * e.g. an indel is always one before the actual event (as per VCF convention). - * - * @param inputVC the VC to evaluate, cannot be null - * @return true if - */ - public static boolean needsPadding(final VariantContext inputVC) { - // biallelic sites with only symbolic never need padding - if ( inputVC.isBiallelic() && inputVC.getAlternateAllele(0).isSymbolic() ) - return false; - - final int recordLength = inputVC.getEnd() - inputVC.getStart() + 1; - final int referenceLength = inputVC.getReference().length(); - - if ( referenceLength == recordLength ) - return false; - else if ( referenceLength == recordLength - 1 ) - return true; - else if ( !inputVC.hasSymbolicAlleles() ) - throw new IllegalArgumentException("Badly formed variant context at location " + String.valueOf(inputVC.getStart()) + - " in contig " + inputVC.getChr() + ". Reference length must be at most one base shorter than location size"); - else if ( inputVC.isMixed() && inputVC.hasSymbolicAlleles() ) - throw new IllegalArgumentException("GATK infrastructure limitation prevents needsPadding from working properly with VariantContexts containing a mixture of symbolic and concrete alleles at " + inputVC); - return false; - } - - public static Allele padAllele(final VariantContext vc, final Allele allele) { - assert needsPadding(vc); - - if ( allele.isSymbolic() ) - return allele; - else { - // get bases for current allele and create a new one with trimmed bases - final StringBuilder sb = new StringBuilder(); - sb.append((char)vc.getReferenceBaseForIndel().byteValue()); - sb.append(allele.getDisplayString()); - final String newBases = sb.toString(); - return Allele.create(newBases, allele.isReference()); - } - } - - public static VariantContext createVariantContextWithPaddedAlleles(VariantContext inputVC) { - final boolean padVC = needsPadding(inputVC); - - // nothing to do if we don't need to pad bases - if ( padVC ) { - if ( !inputVC.hasReferenceBaseForIndel() ) - throw new ReviewedStingException("Badly formed variant context at location " + inputVC.getChr() + ":" + inputVC.getStart() + "; no padded reference base is available."); - - final ArrayList alleles = new ArrayList(inputVC.getNAlleles()); - final Map unpaddedToPadded = inputVC.hasGenotypes() ? new HashMap(inputVC.getNAlleles()) : null; - - boolean paddedAtLeastOne = false; - for (final Allele a : inputVC.getAlleles()) { - final Allele padded = padAllele(inputVC, a); - paddedAtLeastOne = paddedAtLeastOne || padded != a; - alleles.add(padded); - if ( unpaddedToPadded != null ) unpaddedToPadded.put(a, padded); // conditional to avoid making unnecessary make - } - - if ( ! paddedAtLeastOne ) - throw new ReviewedStingException("VC was supposed to need padding but no allele was actually changed at location " + inputVC.getChr() + ":" + inputVC.getStart() + " with allele " + inputVC.getAlleles()); - - final VariantContextBuilder vcb = new VariantContextBuilder(inputVC); - vcb.alleles(alleles); - - // the position of the inputVC is one further, if it doesn't contain symbolic alleles - vcb.computeEndFromAlleles(alleles, inputVC.getStart(), inputVC.getEnd()); - - if ( inputVC.hasGenotypes() ) { - assert unpaddedToPadded != null; - - // now we can recreate new genotypes with trimmed alleles - final GenotypesContext genotypes = GenotypesContext.create(inputVC.getNSamples()); - for (final Genotype g : inputVC.getGenotypes() ) { - final List newGenotypeAlleles = new ArrayList(g.getAlleles().size()); - for (final Allele a : g.getAlleles()) { - newGenotypeAlleles.add( a.isCalled() ? unpaddedToPadded.get(a) : Allele.NO_CALL); - } - genotypes.add(new GenotypeBuilder(g).alleles(newGenotypeAlleles).make()); - } - vcb.genotypes(genotypes); - } - - return vcb.make(); - } - else - return inputVC; - - } - - public static VariantContext reverseTrimAlleles( final VariantContext inputVC ) { - // see if we need to trim common reference base from all alleles - - final int trimExtent = computeReverseClipping(inputVC.getAlleles(), inputVC.getReference().getDisplayString().getBytes(), 0, true); - if ( trimExtent <= 0 || inputVC.getAlleles().size() <= 1 ) - return inputVC; - - final List alleles = new ArrayList(); - final GenotypesContext genotypes = GenotypesContext.create(); - final Map originalToTrimmedAlleleMap = new HashMap(); - - for (final Allele a : inputVC.getAlleles()) { - if (a.isSymbolic()) { - alleles.add(a); - originalToTrimmedAlleleMap.put(a, a); - } else { - // get bases for current allele and create a new one with trimmed bases - final byte[] newBases = Arrays.copyOfRange(a.getBases(), 0, a.length()-trimExtent); - final Allele trimmedAllele = Allele.create(newBases, a.isReference()); - alleles.add(trimmedAllele); - originalToTrimmedAlleleMap.put(a, trimmedAllele); - } - } - - // now we can recreate new genotypes with trimmed alleles - for ( final Genotype genotype : inputVC.getGenotypes() ) { - final List originalAlleles = genotype.getAlleles(); - final List trimmedAlleles = new ArrayList(); - for ( final Allele a : originalAlleles ) { - if ( a.isCalled() ) - trimmedAlleles.add(originalToTrimmedAlleleMap.get(a)); - else - trimmedAlleles.add(Allele.NO_CALL); - } - genotypes.add(new GenotypeBuilder(genotype).alleles(trimmedAlleles).make()); - } - - return new VariantContextBuilder(inputVC).stop(inputVC.getStart() + alleles.get(0).length() + (inputVC.isMixed() ? -1 : 0)).alleles(alleles).genotypes(genotypes).make(); - } - - @Invariant("stop != -1 || error != null") // we're either an error or a meaningful result but not both - public static class ClippedAlleles { - private final int stop; - private final List clippedAlleles; - private final Byte refBaseForIndel; - private final String error; - - @Requires({"stop > 0", "clippedAlleles != null"}) - private ClippedAlleles(final int stop, final List clippedAlleles, final Byte refBaseForIndel) { - this.stop = stop; - this.clippedAlleles = clippedAlleles; - this.error = null; - this.refBaseForIndel = refBaseForIndel; - } - - @Requires("error != null") - private ClippedAlleles(final String error) { - this.stop = -1; - this.clippedAlleles = null; - this.refBaseForIndel = null; - this.error = error; - } - - /** - * Get an error if it occurred - * @return the error message, or null if no error occurred - */ - public String getError() { - return error; - } - - /** - * Get the stop position to use after the clipping as been applied, given the - * provided position to clipAlleles - * @return - */ - public int getStop() { - return stop; - } - - /** - * Get the clipped alleles themselves - * @return the clipped alleles in the order of the input unclipped alleles - */ - public List getClippedAlleles() { - return clippedAlleles; - } - - /** - * Returns the reference base we should use for indels, or null if none is appropriate - * @return - */ - public Byte getRefBaseForIndel() { - return refBaseForIndel; - } - } -} diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/Allele.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/Allele.java index 2e1770581..1947ef01e 100755 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/Allele.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/Allele.java @@ -1,9 +1,9 @@ package org.broadinstitute.sting.utils.variantcontext; -import java.util.ArrayList; +import org.broadinstitute.sting.utils.BaseUtils; + import java.util.Arrays; import java.util.Collection; -import java.util.List; /** * Immutable representation of an allele @@ -77,22 +77,19 @@ public class Allele implements Comparable { private static final byte[] EMPTY_ALLELE_BASES = new byte[0]; private boolean isRef = false; - private boolean isNull = false; private boolean isNoCall = false; private boolean isSymbolic = false; private byte[] bases = null; - public final static String NULL_ALLELE_STRING = "-"; public final static String NO_CALL_STRING = "."; /** A generic static NO_CALL allele for use */ // no public way to create an allele private Allele(byte[] bases, boolean isRef) { - // standardize our representation of null allele and bases + // null alleles are no longer allowed if ( wouldBeNullAllele(bases) ) { - bases = EMPTY_ALLELE_BASES; - isNull = true; + throw new IllegalArgumentException("Null alleles are not supported"); } else if ( wouldBeNoCallAllele(bases) ) { bases = EMPTY_ALLELE_BASES; isNoCall = true; @@ -101,8 +98,8 @@ public class Allele implements Comparable { isSymbolic = true; if ( isRef ) throw new IllegalArgumentException("Cannot tag a symbolic allele as the reference allele"); } -// else -// bases = new String(bases).toUpperCase().getBytes(); // todo -- slow performance + else + bases = BaseUtils.convertToUpperCase(bases); this.isRef = isRef; this.bases = bases; @@ -126,8 +123,6 @@ public class Allele implements Comparable { private final static Allele ALT_T = new Allele("T", false); private final static Allele REF_N = new Allele("N", true); private final static Allele ALT_N = new Allele("N", false); - private final static Allele REF_NULL = new Allele(NULL_ALLELE_STRING, true); - private final static Allele ALT_NULL = new Allele(NULL_ALLELE_STRING, false); public final static Allele NO_CALL = new Allele(NO_CALL_STRING, false); // --------------------------------------------------------------------------------------------------------- @@ -154,7 +149,6 @@ public class Allele implements Comparable { case '.': if ( isRef ) throw new IllegalArgumentException("Cannot tag a NoCall allele as the reference allele"); return NO_CALL; - case '-': return isRef ? REF_NULL : ALT_NULL; case 'A': case 'a' : return isRef ? REF_A : ALT_A; case 'C': case 'c' : return isRef ? REF_C : ALT_C; case 'G': case 'g' : return isRef ? REF_G : ALT_G; @@ -179,7 +173,7 @@ public class Allele implements Comparable { public static Allele extend(Allele left, byte[] right) { if (left.isSymbolic()) throw new IllegalArgumentException("Cannot extend a symbolic allele"); - byte[] bases = null; + byte[] bases; if ( left.length() == 0 ) bases = right; else { @@ -242,7 +236,10 @@ public class Allele implements Comparable { } public static boolean acceptableAlleleBases(byte[] bases, boolean allowNsAsAcceptable) { - if ( wouldBeNullAllele(bases) || wouldBeNoCallAllele(bases) || wouldBeSymbolicAllele(bases) ) + if ( wouldBeNullAllele(bases) ) + return false; + + if ( wouldBeNoCallAllele(bases) || wouldBeSymbolicAllele(bases) ) return true; for (byte base : bases ) { @@ -299,11 +296,6 @@ public class Allele implements Comparable { // // --------------------------------------------------------------------------------------------------------- - //Returns true if this is the null allele - public boolean isNull() { return isNull; } - // Returns true if this is not the null allele - public boolean isNonNull() { return ! isNull(); } - // Returns true if this is the NO_CALL allele public boolean isNoCall() { return isNoCall; } // Returns true if this is not the NO_CALL allele @@ -319,7 +311,7 @@ public class Allele implements Comparable { // Returns a nice string representation of this object public String toString() { - return (isNull() ? NULL_ALLELE_STRING : ( isNoCall() ? NO_CALL_STRING : getDisplayString() )) + (isReference() ? "*" : ""); + return ( isNoCall() ? NO_CALL_STRING : getDisplayString() ) + (isReference() ? "*" : ""); } /** @@ -384,27 +376,27 @@ public class Allele implements Comparable { * @return true if this and other are equal */ public boolean equals(Allele other, boolean ignoreRefState) { - return this == other || (isRef == other.isRef || ignoreRefState) && isNull == other.isNull && isNoCall == other.isNoCall && (bases == other.bases || Arrays.equals(bases, other.bases)); + return this == other || (isRef == other.isRef || ignoreRefState) && isNoCall == other.isNoCall && (bases == other.bases || Arrays.equals(bases, other.bases)); } /** * @param test bases to test against * - * @return true if this Alelle contains the same bases as test, regardless of its reference status; handles Null and NO_CALL alleles + * @return true if this Allele contains the same bases as test, regardless of its reference status; handles Null and NO_CALL alleles */ public boolean basesMatch(byte[] test) { return !isSymbolic && (bases == test || Arrays.equals(bases, test)); } /** * @param test bases to test against * - * @return true if this Alelle contains the same bases as test, regardless of its reference status; handles Null and NO_CALL alleles + * @return true if this Allele contains the same bases as test, regardless of its reference status; handles Null and NO_CALL alleles */ public boolean basesMatch(String test) { return basesMatch(test.toUpperCase().getBytes()); } /** * @param test allele to test against * - * @return true if this Alelle contains the same bases as test, regardless of its reference status; handles Null and NO_CALL alleles + * @return true if this Allele contains the same bases as test, regardless of its reference status; handles Null and NO_CALL alleles */ public boolean basesMatch(Allele test) { return basesMatch(test.getBases()); } @@ -421,10 +413,6 @@ public class Allele implements Comparable { // // --------------------------------------------------------------------------------------------------------- - public static Allele getMatchingAllele(Collection allAlleles, String alleleBases) { - return getMatchingAllele(allAlleles, alleleBases.getBytes()); - } - public static Allele getMatchingAllele(Collection allAlleles, byte[] alleleBases) { for ( Allele a : allAlleles ) { if ( a.basesMatch(alleleBases) ) { @@ -438,26 +426,6 @@ public class Allele implements Comparable { return null; // couldn't find anything } - public static List resolveAlleles(List possibleAlleles, List alleleStrings) { - List myAlleles = new ArrayList(alleleStrings.size()); - - for ( String alleleString : alleleStrings ) { - Allele allele = getMatchingAllele(possibleAlleles, alleleString); - - if ( allele == null ) { - if ( Allele.wouldBeNoCallAllele(alleleString.getBytes()) ) { - allele = create(alleleString); - } else { - throw new IllegalArgumentException("Allele " + alleleString + " not present in the list of alleles " + possibleAlleles); - } - } - - myAlleles.add(allele); - } - - return myAlleles; - } - public int compareTo(Allele other) { if ( isReference() && other.isNonReference() ) return -1; @@ -468,9 +436,6 @@ public class Allele implements Comparable { } public static boolean oneIsPrefixOfOther(Allele a1, Allele a2) { - if ( a1.isNull() || a2.isNull() ) - return true; - if ( a2.length() >= a1.length() ) return firstIsPrefixOfSecond(a1, a2); else diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContext.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContext.java index dcdd95d00..f298f1187 100755 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContext.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContext.java @@ -188,8 +188,6 @@ public class VariantContext implements Feature { // to enable tribble integratio @Deprecated // ID is no longer stored in the attributes map private final static String ID_KEY = "ID"; - private final Byte REFERENCE_BASE_FOR_INDEL; - public final static Set PASSES_FILTERS = Collections.unmodifiableSet(new LinkedHashSet()); /** The location of this VariantContext */ @@ -228,7 +226,6 @@ public class VariantContext implements Feature { // to enable tribble integratio // --------------------------------------------------------------------------------------------------------- public enum Validation { - REF_PADDING, ALLELES, GENOTYPES } @@ -250,7 +247,7 @@ public class VariantContext implements Feature { // to enable tribble integratio this(other.getSource(), other.getID(), other.getChr(), other.getStart(), other.getEnd(), other.getAlleles(), other.getGenotypes(), other.getLog10PError(), other.getFiltersMaybeNull(), - other.getAttributes(), other.REFERENCE_BASE_FOR_INDEL, + other.getAttributes(), other.fullyDecoded, NO_VALIDATION); } @@ -266,7 +263,6 @@ public class VariantContext implements Feature { // to enable tribble integratio * @param log10PError qual * @param filters filters: use null for unfiltered and empty set for passes filters * @param attributes attributes - * @param referenceBaseForIndel padded reference base * @param validationToPerform set of validation steps to take */ protected VariantContext(final String source, @@ -279,7 +275,6 @@ public class VariantContext implements Feature { // to enable tribble integratio final double log10PError, final Set filters, final Map attributes, - final Byte referenceBaseForIndel, final boolean fullyDecoded, final EnumSet validationToPerform ) { if ( contig == null ) { throw new IllegalArgumentException("Contig cannot be null"); } @@ -292,7 +287,6 @@ public class VariantContext implements Feature { // to enable tribble integratio this.ID = ID.equals(VCFConstants.EMPTY_ID_FIELD) ? VCFConstants.EMPTY_ID_FIELD : ID; this.commonInfo = new CommonInfo(source, log10PError, filters, attributes); - REFERENCE_BASE_FOR_INDEL = referenceBaseForIndel; // todo -- remove me when this check is no longer necessary if ( this.commonInfo.hasAttribute(ID_KEY) ) @@ -340,8 +334,9 @@ public class VariantContext implements Feature { // to enable tribble integratio * in this VC is returned as the set of alleles in the subContext, even if * some of those alleles aren't in the samples * - * @param sampleNames - * @return + * @param sampleNames the sample names + * @param rederiveAllelesFromGenotypes if true, returns the alleles to just those in use by the samples + * @return new VariantContext subsetting to just the given samples */ public VariantContext subContextFromSamples(Set sampleNames, final boolean rederiveAllelesFromGenotypes ) { if ( sampleNames.containsAll(getSampleNames()) ) { @@ -501,7 +496,7 @@ public class VariantContext implements Feature { // to enable tribble integratio */ public boolean isSimpleInsertion() { // can't just call !isSimpleDeletion() because of complex indels - return getType() == Type.INDEL && getReference().isNull() && isBiallelic(); + return getType() == Type.INDEL && isBiallelic() && getReference().length() < getAlternateAllele(0).length(); } /** @@ -509,7 +504,7 @@ public class VariantContext implements Feature { // to enable tribble integratio */ public boolean isSimpleDeletion() { // can't just call !isSimpleInsertion() because of complex indels - return getType() == Type.INDEL && getAlternateAllele(0).isNull() && isBiallelic(); + return getType() == Type.INDEL && isBiallelic() && getReference().length() > getAlternateAllele(0).length(); } /** @@ -553,22 +548,6 @@ public class VariantContext implements Feature { // to enable tribble integratio return ID; } - public boolean hasReferenceBaseForIndel() { - return REFERENCE_BASE_FOR_INDEL != null; - } - - // the indel base that gets stripped off for indels - public Byte getReferenceBaseForIndel() { - return REFERENCE_BASE_FOR_INDEL; - } - - public String getAlleleStringWithRefPadding(final Allele allele) { - if ( VCFAlleleClipper.needsPadding(this) ) - return VCFAlleleClipper.padAllele(this, allele).getDisplayString(); - else - return allele.getDisplayString(); - } - // --------------------------------------------------------------------------------------------------------- // @@ -808,8 +787,8 @@ public class VariantContext implements Feature { // to enable tribble integratio * Returns a map from sampleName -> Genotype for the genotype associated with sampleName. Returns a map * for consistency with the multi-get function. * - * @param sampleName - * @return + * @param sampleName the sample name + * @return mapping from sample name to genotype * @throws IllegalArgumentException if sampleName isn't bound to a genotype */ public GenotypesContext getGenotypes(String sampleName) { @@ -823,7 +802,7 @@ public class VariantContext implements Feature { // to enable tribble integratio * For testing convenience only * * @param sampleNames a unique list of sample names - * @return + * @return subsetting genotypes context * @throws IllegalArgumentException if sampleName isn't bound to a genotype */ protected GenotypesContext getGenotypes(Collection sampleNames) { @@ -1011,13 +990,13 @@ public class VariantContext implements Feature { // to enable tribble integratio /** * Run all extra-strict validation tests on a Variant Context object * - * @param reference the true reference allele - * @param paddedRefBase the reference base used for padding indels - * @param rsIDs the true dbSNP IDs + * @param reportedReference the reported reference allele + * @param observedReference the actual reference allele + * @param rsIDs the true dbSNP IDs */ - public void extraStrictValidation(Allele reference, Byte paddedRefBase, Set rsIDs) { + public void extraStrictValidation(final Allele reportedReference, final Allele observedReference, final Set rsIDs) { // validate the reference - validateReferenceBases(reference, paddedRefBase); + validateReferenceBases(reportedReference, observedReference); // validate the RS IDs validateRSIDs(rsIDs); @@ -1032,18 +1011,9 @@ public class VariantContext implements Feature { // to enable tribble integratio //checkReferenceTrack(); } - public void validateReferenceBases(Allele reference, Byte paddedRefBase) { - if ( reference == null ) - return; - - // don't validate if we're a complex event - if ( !isComplexIndel() && !reference.isNull() && !reference.basesMatch(getReference()) ) { - throw new TribbleException.InternalCodecException(String.format("the REF allele is incorrect for the record at position %s:%d, fasta says %s vs. VCF says %s", getChr(), getStart(), reference.getBaseString(), getReference().getBaseString())); - } - - // we also need to validate the padding base for simple indels - if ( hasReferenceBaseForIndel() && !getReferenceBaseForIndel().equals(paddedRefBase) ) { - throw new TribbleException.InternalCodecException(String.format("the padded REF base is incorrect for the record at position %s:%d, fasta says %s vs. VCF says %s", getChr(), getStart(), (char)paddedRefBase.byteValue(), (char)getReferenceBaseForIndel().byteValue())); + public void validateReferenceBases(final Allele reportedReference, final Allele observedReference) { + if ( reportedReference != null && !reportedReference.basesMatch(observedReference) ) { + throw new TribbleException.InternalCodecException(String.format("the REF allele is incorrect for the record at position %s:%d, fasta says %s vs. VCF says %s", getChr(), getStart(), observedReference.getBaseString(), reportedReference.getBaseString())); } } @@ -1135,7 +1105,6 @@ public class VariantContext implements Feature { // to enable tribble integratio for (final Validation val : validationToPerform ) { switch (val) { case ALLELES: validateAlleles(); break; - case REF_PADDING: validateReferencePadding(); break; case GENOTYPES: validateGenotypes(); break; default: throw new IllegalArgumentException("Unexpected validation mode " + val); } @@ -1164,20 +1133,11 @@ public class VariantContext implements Feature { // to enable tribble integratio } } - private void validateReferencePadding() { - if ( hasSymbolicAlleles() ) // symbolic alleles don't need padding... - return; - - boolean needsPadding = (getReference().length() == getEnd() - getStart()); // off by one because padded base was removed - - if ( needsPadding && !hasReferenceBaseForIndel() ) - throw new ReviewedStingException("Badly formed variant context at location " + getChr() + ":" + getStart() + "; no padded reference base was provided."); - } - private void validateAlleles() { - // check alleles - boolean alreadySeenRef = false, alreadySeenNull = false; - for ( Allele allele : alleles ) { + + boolean alreadySeenRef = false; + + for ( final Allele allele : alleles ) { // make sure there's only one reference allele if ( allele.isReference() ) { if ( alreadySeenRef ) throw new IllegalArgumentException("BUG: Received two reference tagged alleles in VariantContext " + alleles + " this=" + this); @@ -1187,24 +1147,14 @@ public class VariantContext implements Feature { // to enable tribble integratio if ( allele.isNoCall() ) { throw new IllegalArgumentException("BUG: Cannot add a no call allele to a variant context " + alleles + " this=" + this); } - - // make sure there's only one null allele - if ( allele.isNull() ) { - if ( alreadySeenNull ) throw new IllegalArgumentException("BUG: Received two null alleles in VariantContext " + alleles + " this=" + this); - alreadySeenNull = true; - } } // make sure there's one reference allele if ( ! alreadySeenRef ) throw new IllegalArgumentException("No reference allele found in VariantContext"); -// if ( getType() == Type.INDEL ) { -// if ( getReference().length() != (getLocation().size()-1) ) { - long length = (stop - start) + 1; - if ( ! hasSymbolicAlleles() - && ((getReference().isNull() && length != 1 ) - || (getReference().isNonNull() && (length - getReference().length() > 1)))) { + final long length = (stop - start) + 1; + if ( ! hasSymbolicAlleles() && length != getReference().length() ) { throw new IllegalStateException("BUG: GenomeLoc " + contig + ":" + start + "-" + stop + " has a size == " + length + " but the variation reference allele has length " + getReference().length() + " this = " + this); } } diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextBuilder.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextBuilder.java index f2375f6f9..d8ab4bd23 100644 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextBuilder.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextBuilder.java @@ -25,9 +25,6 @@ package org.broadinstitute.sting.utils.variantcontext; import com.google.java.contract.*; -import org.broad.tribble.Feature; -import org.broad.tribble.TribbleException; -import org.broad.tribble.util.ParsingUtils; import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.codecs.vcf.VCFConstants; import org.broadinstitute.sting.utils.exceptions.ReviewedStingException; @@ -74,7 +71,6 @@ public class VariantContextBuilder { private Set filters = null; private Map attributes = null; private boolean attributesCanBeModified = false; - private Byte referenceBaseForIndel = null; /** enum of what must be validated */ final private EnumSet toValidate = EnumSet.noneOf(VariantContext.Validation.class); @@ -117,7 +113,6 @@ public class VariantContextBuilder { this.genotypes = parent.genotypes; this.ID = parent.getID(); this.log10PError = parent.getLog10PError(); - this.referenceBaseForIndel = parent.getReferenceBaseForIndel(); this.source = parent.getSource(); this.start = parent.getStart(); this.stop = parent.getEnd(); @@ -132,7 +127,6 @@ public class VariantContextBuilder { this.genotypes = parent.genotypes; this.ID = parent.ID; this.log10PError = parent.log10PError; - this.referenceBaseForIndel = parent.referenceBaseForIndel; this.source = parent.source; this.start = parent.start; this.stop = parent.stop; @@ -362,21 +356,6 @@ public class VariantContextBuilder { return this; } - /** - * Tells us that the resulting VariantContext should use this byte for the reference base - * Null means no refBase is available - * @param referenceBaseForIndel - */ - public VariantContextBuilder referenceBaseForIndel(final Byte referenceBaseForIndel) { - this.referenceBaseForIndel = referenceBaseForIndel; - toValidate.add(VariantContext.Validation.REF_PADDING); - return this; - } - - public VariantContextBuilder referenceBaseForIndel(final String referenceBaseForIndel) { - return referenceBaseForIndel(referenceBaseForIndel.getBytes()[0]); - } - /** * Tells us that the resulting VariantContext should have source field set to source * @param source @@ -401,7 +380,6 @@ public class VariantContextBuilder { this.start = start; this.stop = stop; toValidate.add(VariantContext.Validation.ALLELES); - toValidate.add(VariantContext.Validation.REF_PADDING); return this; } @@ -416,7 +394,6 @@ public class VariantContextBuilder { this.start = loc.getStart(); this.stop = loc.getStop(); toValidate.add(VariantContext.Validation.ALLELES); - toValidate.add(VariantContext.Validation.REF_PADDING); return this; } @@ -440,7 +417,6 @@ public class VariantContextBuilder { public VariantContextBuilder start(final long start) { this.start = start; toValidate.add(VariantContext.Validation.ALLELES); - toValidate.add(VariantContext.Validation.REF_PADDING); return this; } @@ -517,6 +493,6 @@ public class VariantContextBuilder { public VariantContext make() { return new VariantContext(source, ID, contig, start, stop, alleles, genotypes, log10PError, filters, attributes, - referenceBaseForIndel, fullyDecoded, toValidate); + fullyDecoded, toValidate); } } diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java index d7e072980..e1a043e94 100755 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/VariantContextUtils.java @@ -64,9 +64,9 @@ public class VariantContextUtils { * Ensures that VC contains all of the samples in allSamples by adding missing samples to * the resulting VC with default diploid ./. genotypes * - * @param vc - * @param allSamples - * @return + * @param vc the VariantContext + * @param allSamples all of the samples needed + * @return a new VariantContext with missing samples added */ public static VariantContext addMissingSamples(final VariantContext vc, final Set allSamples) { // TODO -- what's the fastest way to do this calculation? @@ -376,9 +376,9 @@ public class VariantContextUtils { /** * @deprecated use variant context builder version instead - * @param vc - * @param keysToPreserve - * @return + * @param vc the variant context + * @param keysToPreserve the keys to preserve + * @return a pruned version of the original variant context */ @Deprecated public static VariantContext pruneVariantContext(final VariantContext vc, Collection keysToPreserve ) { @@ -486,14 +486,13 @@ public class VariantContextUtils { if ( genotypeMergeOptions == GenotypeMergeType.REQUIRE_UNIQUE ) verifyUniqueSampleNames(unsortedVCs); - final List prepaddedVCs = sortVariantContextsByPriority(unsortedVCs, priorityListOfVCs, genotypeMergeOptions); + final List preFilteredVCs = sortVariantContextsByPriority(unsortedVCs, priorityListOfVCs, genotypeMergeOptions); // Make sure all variant contexts are padded with reference base in case of indels if necessary final List VCs = new ArrayList(); - for (final VariantContext vc : prepaddedVCs) { - // also a reasonable place to remove filtered calls, if needed + for (final VariantContext vc : preFilteredVCs) { if ( ! filteredAreUncalled || vc.isNotFiltered() ) - VCs.add(VCFAlleleClipper.createVariantContextWithPaddedAlleles(vc)); + VCs.add(vc); } if ( VCs.size() == 0 ) // everything is filtered out and we're filteredAreUncalled return null; @@ -547,9 +546,6 @@ public class VariantContextUtils { filters.addAll(vc.getFilters()); - if ( referenceBaseForIndel == null ) - referenceBaseForIndel = vc.getReferenceBaseForIndel(); - // // add attributes // @@ -661,10 +657,9 @@ public class VariantContextUtils { builder.genotypes(genotypes); builder.log10PError(log10PError); builder.filters(filters).attributes(mergeInfoWithMaxAC ? attributesWithMaxAC : attributes); - builder.referenceBaseForIndel(referenceBaseForIndel); // Trim the padded bases of all alleles if necessary - final VariantContext merged = createVariantContextWithTrimmedAlleles(builder.make()); + final VariantContext merged = builder.make(); if ( printMessages && remapped ) System.out.printf("Remapped => %s%n", merged); return merged; } @@ -700,73 +695,6 @@ public class VariantContextUtils { return true; } - private static VariantContext createVariantContextWithTrimmedAlleles(VariantContext inputVC) { - // see if we need to trim common reference base from all alleles - boolean trimVC; - - // We need to trim common reference base from all alleles in all genotypes if a ref base is common to all alleles - Allele refAllele = inputVC.getReference(); - if (!inputVC.isVariant()) - trimVC = false; - else if (refAllele.isNull()) - trimVC = false; - else { - trimVC = VCFAlleleClipper.shouldClipFirstBaseP(inputVC.getAlternateAlleles(), (byte) inputVC.getReference().getDisplayString().charAt(0)); - } - - // nothing to do if we don't need to trim bases - if (trimVC) { - List alleles = new ArrayList(); - GenotypesContext genotypes = GenotypesContext.create(); - - Map originalToTrimmedAlleleMap = new HashMap(); - - for (final Allele a : inputVC.getAlleles()) { - if (a.isSymbolic()) { - alleles.add(a); - originalToTrimmedAlleleMap.put(a, a); - } else { - // get bases for current allele and create a new one with trimmed bases - byte[] newBases = Arrays.copyOfRange(a.getBases(), 1, a.length()); - Allele trimmedAllele = Allele.create(newBases, a.isReference()); - alleles.add(trimmedAllele); - originalToTrimmedAlleleMap.put(a, trimmedAllele); - } - } - - // detect case where we're trimming bases but resulting vc doesn't have any null allele. In that case, we keep original representation - // example: mixed records such as {TA*,TGA,TG} - boolean hasNullAlleles = false; - - for (final Allele a: originalToTrimmedAlleleMap.values()) { - if (a.isNull()) - hasNullAlleles = true; - } - - if (!hasNullAlleles) - return inputVC; - // now we can recreate new genotypes with trimmed alleles - for ( final Genotype genotype : inputVC.getGenotypes() ) { - - List originalAlleles = genotype.getAlleles(); - List trimmedAlleles = new ArrayList(); - for ( final Allele a : originalAlleles ) { - if ( a.isCalled() ) - trimmedAlleles.add(originalToTrimmedAlleleMap.get(a)); - else - trimmedAlleles.add(Allele.NO_CALL); - } - genotypes.add(new GenotypeBuilder(genotype).alleles(trimmedAlleles).make()); - - } - - final VariantContextBuilder builder = new VariantContextBuilder(inputVC); - return builder.alleles(alleles).genotypes(genotypes).referenceBaseForIndel(new Byte(inputVC.getReference().getBases()[0])).make(); - } - - return inputVC; - } - public static GenotypesContext stripPLs(GenotypesContext genotypes) { GenotypesContext newGs = GenotypesContext.create(genotypes.size()); @@ -979,7 +907,7 @@ public class VariantContextUtils { HashMap alleleMap = new HashMap(vc.getAlleles().size()); for ( Allele originalAllele : vc.getAlleles() ) { Allele newAllele; - if ( originalAllele.isNoCall() || originalAllele.isNull() ) + if ( originalAllele.isNoCall() ) newAllele = originalAllele; else newAllele = Allele.create(BaseUtils.simpleReverseComplement(originalAllele.getBases()), originalAllele.isReference()); diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/BCF2Writer.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/BCF2Writer.java index df2008e8e..b5da206ad 100644 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/BCF2Writer.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/BCF2Writer.java @@ -274,10 +274,7 @@ class BCF2Writer extends IndexingVariantContextWriter { } private void buildAlleles( VariantContext vc ) throws IOException { - final boolean needsPadding = VCFAlleleClipper.needsPadding(vc); for ( Allele allele : vc.getAlleles() ) { - if ( needsPadding ) - allele = VCFAlleleClipper.padAllele(vc, allele); final byte[] s = allele.getDisplayBases(); if ( s == null ) throw new ReviewedStingException("BUG: BCF2Writer encountered null padded allele" + allele); diff --git a/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/VCFWriter.java b/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/VCFWriter.java index 4548e026e..ea968e153 100755 --- a/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/VCFWriter.java +++ b/public/java/src/org/broadinstitute/sting/utils/variantcontext/writer/VCFWriter.java @@ -162,7 +162,6 @@ class VCFWriter extends IndexingVariantContextWriter { vc = new VariantContextBuilder(vc).noGenotypes().make(); try { - vc = VCFAlleleClipper.createVariantContextWithPaddedAlleles(vc); super.add(vc); Map alleleMap = buildAlleleMap(vc); diff --git a/public/java/test/org/broadinstitute/sting/utils/variantcontext/AlleleUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/variantcontext/AlleleUnitTest.java index ed9805d19..3bf020df7 100755 --- a/public/java/test/org/broadinstitute/sting/utils/variantcontext/AlleleUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/variantcontext/AlleleUnitTest.java @@ -37,8 +37,6 @@ import org.testng.annotations.Test; // public Allele(byte[] bases, boolean isRef) { // public Allele(boolean isRef) { // public Allele(String bases, boolean isRef) { -// public boolean isNullAllele() { return length() == 0; } -// public boolean isNonNullAllele() { return ! isNullAllele(); } // public boolean isReference() { return isRef; } // public boolean isNonReference() { return ! isReference(); } // public byte[] getBases() { return bases; } @@ -72,8 +70,6 @@ public class AlleleUnitTest { Assert.assertFalse(A.isReference()); Assert.assertTrue(A.basesMatch("A")); Assert.assertEquals(A.length(), 1); - Assert.assertTrue(A.isNonNull()); - Assert.assertFalse(A.isNull()); Assert.assertTrue(ARef.isReference()); Assert.assertFalse(ARef.isNonReference()); @@ -92,8 +88,8 @@ public class AlleleUnitTest { Assert.assertFalse(NoCall.isReference()); Assert.assertFalse(NoCall.basesMatch(".")); Assert.assertEquals(NoCall.length(), 0); - Assert.assertTrue(NoCall.isNonNull()); - Assert.assertFalse(NoCall.isNull()); + Assert.assertTrue(NoCall.isNoCall()); + Assert.assertFalse(NoCall.isCalled()); } @@ -111,8 +107,6 @@ public class AlleleUnitTest { Assert.assertFalse(del.basesMatch("-")); Assert.assertTrue(del.basesMatch("")); Assert.assertEquals(del.length(), 0); - Assert.assertFalse(del.isNonNull()); - Assert.assertTrue(del.isNull()); } diff --git a/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUnitTest.java index 1d290118f..11c75ed9a 100755 --- a/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUnitTest.java @@ -845,7 +845,6 @@ public class VariantContextUnitTest extends BaseTest { Assert.assertEquals(sub.getLog10PError(), vc.getLog10PError()); Assert.assertEquals(sub.getFilters(), vc.getFilters()); Assert.assertEquals(sub.getID(), vc.getID()); - Assert.assertEquals(sub.getReferenceBaseForIndel(), vc.getReferenceBaseForIndel()); Assert.assertEquals(sub.getAttributes(), vc.getAttributes()); Set expectedGenotypes = new HashSet();