diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/PoolGenotypeLikelihoodsUnitTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/PoolGenotypeLikelihoodsUnitTest.java index 9aab24998..74abb6b11 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/PoolGenotypeLikelihoodsUnitTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/PoolGenotypeLikelihoodsUnitTest.java @@ -27,7 +27,6 @@ package org.broadinstitute.sting.gatk.walkers.genotyper; import net.sf.samtools.SAMUtils; import org.apache.log4j.Logger; import org.broadinstitute.sting.gatk.contexts.AlignmentContext; -import org.broadinstitute.sting.gatk.contexts.ReferenceContext; import org.broadinstitute.sting.gatk.walkers.Walker; import org.broadinstitute.sting.utils.BaseUtils; import org.broadinstitute.sting.utils.MathUtils; @@ -290,15 +289,17 @@ public class PoolGenotypeLikelihoodsUnitTest { } - @Test + // TODO -- Guillermo, this test cannot work because the ArtificialReadPileupTestProvider returns a position of chr1:5, which is less than + // TODO -- HAPLOTYPE_SIZE in IndelGenotypeLikelihoodsCalculationModel.getHaplotypeMapFromAlleles() so the HaplotypeMap is not populated. + @Test (enabled = false) public void testIndelErrorModel() { final ArtificialReadPileupTestProvider refPileupTestProvider = new ArtificialReadPileupTestProvider(1,"ref"); final byte refByte = refPileupTestProvider.getRefByte(); - final String altBases = refByte + "TCA"; + final String altBases = (char)refByte + "TCA"; final String refSampleName = refPileupTestProvider.getSampleNames().get(0); final List trueAlleles = new ArrayList(); trueAlleles.add(Allele.create(refByte, true)); - trueAlleles.add(Allele.create(refByte + "TC", false)); + trueAlleles.add(Allele.create((char)refByte + "TC", false)); final String fw = new String(refPileupTestProvider.getReferenceContext().getForwardBases()); final VariantContext refInsertionVC = new VariantContextBuilder("test","chr1",refPileupTestProvider.getReferenceContext().getLocus().getStart(), @@ -392,9 +393,6 @@ public class PoolGenotypeLikelihoodsUnitTest { final byte refByte = readPileupTestProvider.getRefByte(); final byte altByte = refByte == (byte)'T'? (byte) 'C': (byte)'T'; - final int refIdx = BaseUtils.simpleBaseToBaseIndex(refByte); - final int altIdx = BaseUtils.simpleBaseToBaseIndex(altByte); - final List allAlleles = new ArrayList(); // this contains only ref Allele up to now final Set laneIDs = new TreeSet(); laneIDs.add(GenotypeLikelihoodsCalculationModel.DUMMY_LANE); @@ -411,11 +409,17 @@ public class PoolGenotypeLikelihoodsUnitTest { for (String laneID : laneIDs) noisyErrorModels.put(laneID, Q30ErrorModel); + final int refIdx = 0; + int altIdx = 2; + + // ref allele must be first + allAlleles.add(Allele.create(refByte, true)); for (byte b: BaseUtils.BASES) { - if (refByte == b) - allAlleles.add(Allele.create(b,true)); - else + if (refByte != b) { + if (b == altByte) + altIdx = allAlleles.size(); allAlleles.add(Allele.create(b, false)); + } } PrintStream out = null; diff --git a/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/ArtificialReadPileupTestProvider.java b/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/ArtificialReadPileupTestProvider.java index a911718c1..17149220a 100644 --- a/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/ArtificialReadPileupTestProvider.java +++ b/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/ArtificialReadPileupTestProvider.java @@ -38,9 +38,6 @@ import org.broadinstitute.sting.utils.pileup.ReadBackedPileup; import org.broadinstitute.sting.utils.pileup.ReadBackedPileupImpl; import org.broadinstitute.sting.utils.sam.ArtificialSAMUtils; import org.broadinstitute.sting.utils.sam.GATKSAMRecord; -import org.broadinstitute.sting.utils.variantcontext.Allele; -import org.broadinstitute.sting.utils.variantcontext.VariantContext; -import org.broadinstitute.sting.utils.variantcontext.VariantContextBuilder; import java.util.*; @@ -103,39 +100,27 @@ public class ArtificialReadPileupTestProvider { boolean addBaseErrors, int phredScaledBaseErrorRate) { // RefMetaDataTracker tracker = new RefMetaDataTracker(null,referenceContext); - ArrayList vcAlleles = new ArrayList(); - String refBase = refBases.substring(offset,offset+1); // referenceContext.getBase()? - Allele refAllele, altAllele; + String refAllele, altAllele; if (eventLength == 0) { // SNP case - refAllele = Allele.create(refBase,true); - altAllele = Allele.create(altBases.substring(0,1), false); + refAllele = new String(new byte[]{referenceContext.getBase()}); + altAllele = new String(altBases.substring(0,1)); } else if (eventLength>0){ // insertion - refAllele = Allele.create(refBase,true); - altAllele = Allele.create(refBase + altBases.substring(0,eventLength), false); + refAllele = ""; + altAllele = altBases.substring(0,eventLength); } else { // deletion - refAllele = Allele.create(refBases.substring(offset,offset+Math.abs(eventLength)),true); - altAllele = Allele.create(refBase, false); + refAllele = refBases.substring(offset,offset+Math.abs(eventLength)); + altAllele = ""; } - int stop = loc.getStart(); - vcAlleles.add(refAllele); - vcAlleles.add(altAllele); - - final VariantContextBuilder builder = new VariantContextBuilder().source(""); - builder.loc(loc.getContig(), loc.getStart(), stop); - builder.alleles(vcAlleles); - builder.noGenotypes(); - - final VariantContext vc = builder.make(); Map contexts = new HashMap(); for (String sample: sampleNames) { - AlignmentContext context = new AlignmentContext(loc, generateRBPForVariant(loc,vc, altBases, numReadsPerAllele, sample, addBaseErrors, phredScaledBaseErrorRate)); + AlignmentContext context = new AlignmentContext(loc, generateRBPForVariant(loc, refAllele, altAllele, altBases, numReadsPerAllele, sample, addBaseErrors, phredScaledBaseErrorRate)); contexts.put(sample,context); } @@ -149,73 +134,71 @@ public class ArtificialReadPileupTestProvider { rg.setSample(name); return rg; } - private ReadBackedPileup generateRBPForVariant( GenomeLoc loc, VariantContext vc, String altBases, + + private ReadBackedPileup generateRBPForVariant( GenomeLoc loc, String refAllele, String altAllele, String altBases, int[] numReadsPerAllele, String sample, boolean addErrors, int phredScaledErrorRate) { List pileupElements = new ArrayList(); - int readStart = contigStart; int offset = (contigStop-contigStart+1)/2; - int refAlleleLength = 0; - int readCounter = 0; - int alleleCounter = 0; - for (Allele allele: vc.getAlleles()) { - if (allele.isReference()) - refAlleleLength = allele.getBases().length; - - int alleleLength = allele.getBases().length; - - for ( int d = 0; d < numReadsPerAllele[alleleCounter]; d++ ) { - byte[] readBases = trueHaplotype(allele, offset, refAlleleLength); - if (addErrors) - addBaseErrors(readBases, phredScaledErrorRate); - - byte[] readQuals = new byte[readBases.length]; - Arrays.fill(readQuals, (byte)phredScaledErrorRate); - - GATKSAMRecord read = new GATKSAMRecord(header); - read.setBaseQualities(readQuals); - read.setReadBases(readBases); - read.setReadName(artificialReadName+readCounter++); - - boolean isBeforeDeletion = false, isBeforeInsertion = false; - if (allele.isReference()) - read.setCigarString(readBases.length + "M"); - else { - isBeforeDeletion = alleleLengthrefAlleleLength; - if (isBeforeDeletion || isBeforeInsertion) - read.setCigarString(offset+"M"+ alleleLength + (isBeforeDeletion?"D":"I") + - (readBases.length-offset)+"M"); - else // SNP case - read.setCigarString(readBases.length+"M"); - } - - int eventLength = (isBeforeDeletion?refAlleleLength:(isBeforeInsertion?alleleLength:0)); - read.setReadPairedFlag(false); - read.setAlignmentStart(readStart); - read.setMappingQuality(artificialMappingQuality); - read.setReferenceName(loc.getContig()); - read.setReadNegativeStrandFlag(false); - read.setAttribute("RG", sampleRG(sample).getReadGroupId()); - - - pileupElements.add(new PileupElement(read,offset,false,isBeforeDeletion, false, isBeforeInsertion,false,false,altBases.substring(0,alleleLength),eventLength)); - } - alleleCounter++; - } + int refAlleleLength = refAllele.length(); + pileupElements.addAll(createPileupElements(refAllele, loc, numReadsPerAllele[0], sample, contigStart, offset, altBases, addErrors, phredScaledErrorRate, refAlleleLength, true)); + pileupElements.addAll(createPileupElements(altAllele, loc, numReadsPerAllele[1], sample, contigStart, offset, altBases, addErrors, phredScaledErrorRate, refAlleleLength, false)); return new ReadBackedPileupImpl(loc,pileupElements); } - private byte[] trueHaplotype(Allele allele, int offset, int refAlleleLength) { + private List createPileupElements(String allele, GenomeLoc loc, int numReadsPerAllele, String sample, int readStart, int offset, String altBases, boolean addErrors, int phredScaledErrorRate, int refAlleleLength, boolean isReference) { + + int alleleLength = allele.length(); + List pileupElements = new ArrayList(); + + int readCounter = 0; + for ( int d = 0; d < numReadsPerAllele; d++ ) { + byte[] readBases = trueHaplotype(allele, offset, refAlleleLength); + if (addErrors) + addBaseErrors(readBases, phredScaledErrorRate); + + byte[] readQuals = new byte[readBases.length]; + Arrays.fill(readQuals, (byte)phredScaledErrorRate); + + GATKSAMRecord read = new GATKSAMRecord(header); + read.setBaseQualities(readQuals); + read.setReadBases(readBases); + read.setReadName(artificialReadName+readCounter++); + + boolean isBeforeDeletion = false, isBeforeInsertion = false; + if (isReference) + read.setCigarString(readBases.length + "M"); + else { + isBeforeDeletion = alleleLengthrefAlleleLength; + if (isBeforeDeletion || isBeforeInsertion) + read.setCigarString(offset+"M"+ alleleLength + (isBeforeDeletion?"D":"I") + + (readBases.length-offset)+"M"); + else // SNP case + read.setCigarString(readBases.length+"M"); + } + + int eventLength = (isBeforeDeletion?refAlleleLength:(isBeforeInsertion?alleleLength:0)); + read.setReadPairedFlag(false); + read.setAlignmentStart(readStart); + read.setMappingQuality(artificialMappingQuality); + read.setReferenceName(loc.getContig()); + read.setReadNegativeStrandFlag(false); + read.setAttribute("RG", sampleRG(sample).getReadGroupId()); + + + pileupElements.add(new PileupElement(read,offset,false,isBeforeDeletion, false, isBeforeInsertion,false,false,altBases.substring(0,alleleLength),eventLength)); + } + + return pileupElements; + } + + private byte[] trueHaplotype(String allele, int offset, int refAlleleLength) { // create haplotype based on a particular allele - String prefix = refBases.substring(offset); - String alleleBases = new String(allele.getBases()); + String prefix = refBases.substring(0, offset); String postfix = refBases.substring(offset+refAlleleLength,refBases.length()); - return (prefix+alleleBases+postfix).getBytes(); - - - + return (prefix+allele+postfix).getBytes(); } private void addBaseErrors(final byte[] readBases, final int phredScaledErrorRate) { diff --git a/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsUnitTest.java b/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsUnitTest.java index c7ef51d0c..dfd4bc525 100644 --- a/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/IndelGenotypeLikelihoodsUnitTest.java @@ -70,7 +70,7 @@ public class IndelGenotypeLikelihoodsUnitTest extends BaseTest { List alleles = getConsensusAlleles(eventLength,true,10,0.1, altBases); Assert.assertEquals(alleles.size(),2); - Assert.assertEquals(alleles.get(1).getBaseString(), altBases.substring(0,eventLength)); + Assert.assertEquals(alleles.get(1).getBaseString().substring(1), altBases.substring(0,eventLength)); @@ -79,7 +79,7 @@ public class IndelGenotypeLikelihoodsUnitTest extends BaseTest { eventLength = 3; alleles = getConsensusAlleles(eventLength,false,10,0.1, altBases); Assert.assertEquals(alleles.size(),2); - Assert.assertEquals(alleles.get(0).getBaseString(), refBases.substring(pileupProvider.offset,pileupProvider.offset+eventLength)); + Assert.assertEquals(alleles.get(0).getBaseString().substring(1), refBases.substring(pileupProvider.offset,pileupProvider.offset+eventLength)); // same with min Reads = 11 alleles = getConsensusAlleles(eventLength,false,11,0.1, altBases); @@ -97,7 +97,7 @@ public class IndelGenotypeLikelihoodsUnitTest extends BaseTest { alleles = getConsensusAlleles(eventLength,true,10,0.1, altBases); Assert.assertEquals(alleles.size(),2); - Assert.assertEquals(alleles.get(1).getBaseString(), altBases.substring(0,eventLength)); + Assert.assertEquals(alleles.get(1).getBaseString().substring(1), altBases.substring(0,eventLength)); altBases = "CCTCNTGAGA"; eventLength = 5; diff --git a/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUtilsUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUtilsUnitTest.java index 8c86a54de..8ba11db02 100644 --- a/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUtilsUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/variantcontext/VariantContextUtilsUnitTest.java @@ -39,7 +39,7 @@ import java.io.FileNotFoundException; import java.util.*; public class VariantContextUtilsUnitTest extends BaseTest { - Allele Aref, T, C, delRef, Cref, ATC, ATCATC; + Allele Aref, T, C, Cref, ATC, ATCATC; private GenomeLocParser genomeLocParser; @BeforeSuite @@ -56,7 +56,6 @@ public class VariantContextUtilsUnitTest extends BaseTest { // alleles Aref = Allele.create("A", true); Cref = Allele.create("C", true); - delRef = Allele.create("-", true); T = Allele.create("T"); C = Allele.create("C"); ATC = Allele.create("ATC"); @@ -156,28 +155,23 @@ public class VariantContextUtilsUnitTest extends BaseTest { Arrays.asList(Aref, C), Arrays.asList(Aref, T, C)); // in order of appearence - // The following is actually a pathological case - there's no way on a vcf to represent a null allele that's non-variant. - // The code converts this (correctly) to a single-base non-variant vc with whatever base was there as a reference. - new MergeAllelesTest(Arrays.asList(delRef), - Arrays.asList(Cref)); + new MergeAllelesTest(Arrays.asList(Aref), + Arrays.asList(Aref, ATC), + Arrays.asList(Aref, ATC)); - new MergeAllelesTest(Arrays.asList(delRef), - Arrays.asList(delRef, ATC), - Arrays.asList(delRef, ATC)); - - new MergeAllelesTest(Arrays.asList(delRef), - Arrays.asList(delRef, ATC, ATCATC), - Arrays.asList(delRef, ATC, ATCATC)); + new MergeAllelesTest(Arrays.asList(Aref), + Arrays.asList(Aref, ATC, ATCATC), + Arrays.asList(Aref, ATC, ATCATC)); // alleles in the order we see them - new MergeAllelesTest(Arrays.asList(delRef, ATCATC), - Arrays.asList(delRef, ATC, ATCATC), - Arrays.asList(delRef, ATCATC, ATC)); + new MergeAllelesTest(Arrays.asList(Aref, ATCATC), + Arrays.asList(Aref, ATC, ATCATC), + Arrays.asList(Aref, ATCATC, ATC)); // same - new MergeAllelesTest(Arrays.asList(delRef, ATC), - Arrays.asList(delRef, ATCATC), - Arrays.asList(delRef, ATC, ATCATC)); + new MergeAllelesTest(Arrays.asList(Aref, ATC), + Arrays.asList(Aref, ATCATC), + Arrays.asList(Aref, ATC, ATCATC)); return MergeAllelesTest.getTests(MergeAllelesTest.class); }