read.setBaseQualities(Utils.dupBytes((byte)30, readLength));
reads.add(read);
}
// create a pileup with all reads having offset 0
final ReadBackedPileup pileup = new ReadBackedPileupImpl(myLocation, reads, 0);
Allele base_A = Allele.create(BaseUtils.Base.A.base);
Allele base_C = Allele.create(BaseUtils.Base.C.base);
Allele base_T = Allele.create(BaseUtils.Base.T.base);
PerReadAlleleLikelihoodMap perReadAlleleLikelihoodMap = new PerReadAlleleLikelihoodMap();
int idx = 0;
for ( final PileupElement e : pileup ) {
for ( final Allele allele : Arrays.asList(base_A,base_C,base_T) ) {
Double likelihood;
if ( idx % 2 == 0 )
likelihood = allele == base_A ? -0.04 : -3.0;
else
likelihood = allele == base_C ? -0.04 : -3.0;
perReadAlleleLikelihoodMap.add(e,allele,likelihood);
}
idx++;
}
Assert.assertEquals(perReadAlleleLikelihoodMap.size(),pileup.depthOfCoverage());
Assert.assertEquals(perReadAlleleLikelihoodMap.getAlleleStratifiedReadMap().keySet().size(),3);
Map<Allele,List<GATKSAMRecord>> halfAhalfC = perReadAlleleLikelihoodMap.getAlleleStratifiedReadMap();
Assert.assertEquals(halfAhalfC.get(base_A).size(),pileup.depthOfCoverage()/2);
Assert.assertEquals(halfAhalfC.get(base_C).size(),pileup.depthOfCoverage()/2);
Assert.assertEquals(halfAhalfC.get(base_T).size(),0);
// make sure the likelihoods are retrievable
idx = 0;
for ( final PileupElement e : pileup ) {
Assert.assertTrue(perReadAlleleLikelihoodMap.containsPileupElement(e));
Map<Allele,Double> likelihoods = perReadAlleleLikelihoodMap.getLikelihoodsAssociatedWithPileupElement(e);
for ( final Allele allele : Arrays.asList(base_A,base_C,base_T) ) {
Double expLik;
if ( idx % 2 == 0 )
expLik = allele == base_A ? -0.04 : -3.0;
else
expLik = allele == base_C ? -0.04 : -3.0;
Assert.assertEquals(likelihoods.get(allele),expLik);
}
idx++;
}
// and test downsampling for good measure
final List<GATKSAMRecord> excessReads = new LinkedList<GATKSAMRecord>();
int prevSize = perReadAlleleLikelihoodMap.size();
for ( int i = 0; i < 10 ; i++ ) {
final GATKSAMRecord read = ArtificialSAMUtils.createArtificialRead(header, "myExcessRead" + i, 0, 1, readLength);
final byte[] bases = Utils.dupBytes((byte)'A', readLength);
bases[0] = (byte)(i % 2 == 0 ? 'A' : 'C'); // every other base is a C
// set the read's bases and quals
read.setReadBases(bases);
read.setBaseQualities(Utils.dupBytes((byte)30, readLength));
for ( final Allele allele : Arrays.asList(base_A,base_C,base_T) ) {
perReadAlleleLikelihoodMap.add(read,allele,allele==base_A ? -0.04 : -3.0);
}
Assert.assertEquals(perReadAlleleLikelihoodMap.size(),1+prevSize);
prevSize = perReadAlleleLikelihoodMap.size();
}
Assert.assertEquals(perReadAlleleLikelihoodMap.size(),pileup.depthOfCoverage()+10);
Assert.assertEquals(perReadAlleleLikelihoodMap.getAlleleStratifiedReadMap().get(base_A).size(),60);
perReadAlleleLikelihoodMap.performPerAlleleDownsampling(0.1);
Assert.assertEquals(perReadAlleleLikelihoodMap.size(),(int) (0.9*(pileup.depthOfCoverage()+10)));
Map<Allele,List<GATKSAMRecord>> downsampledStrat = perReadAlleleLikelihoodMap.getAlleleStratifiedReadMap();
Assert.assertEquals(downsampledStrat.get(base_A).size(),(int) (pileup.depthOfCoverage()/2) - 1);
Assert.assertEquals(downsampledStrat.get(base_C).size(),(int) (pileup.depthOfCoverage()/2));
Assert.assertEquals(downsampledStrat.get(base_T).size(),0);
}