use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class AssemblyRegionTestDataSet method assemblyResultSet.
public AssemblyResultSet assemblyResultSet() {
if (assemblyResultSet == null) {
final ReadThreadingGraph rtg = new ReadThreadingGraph(kmerSize);
rtg.addSequence("anonymous", getReference().getBytes(), true);
for (final String haplotype : haplotypesStrings()) {
rtg.addSequence("anonymous", haplotype.getBytes(), false);
}
rtg.buildGraphIfNecessary();
if (rtg.hasCycles())
throw new RuntimeException("there is cycles in the reference with kmer size " + kmerSize + ". Don't use this size for the benchmark or change the reference");
List<Haplotype> haplotypeList = haplotypeList();
assemblyResultSet = new AssemblyResultSet();
final AssemblyResult ar = new AssemblyResult((haplotypeList.size() > 1 ? AssemblyResult.Status.ASSEMBLED_SOME_VARIATION : AssemblyResult.Status.JUST_ASSEMBLED_REFERENCE), rtg.toSequenceGraph(), rtg);
for (final Haplotype h : haplotypeList) assemblyResultSet.add(h, ar);
}
return assemblyResultSet;
}
use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class AssemblyRegionTestDataSet method expandAllHaplotypeCombinations.
private List<Haplotype> expandAllHaplotypeCombinations(final String civarString, final String reference) {
final Civar civar = Civar.fromCharSequence(civarString);
final List<Civar> unrolledCivars = civar.optionalizeAll().unroll();
List<Haplotype> result = new ArrayList<>(unrolledCivars.size());
for (final Civar c : unrolledCivars) {
final String baseString = c.applyTo(reference);
final Haplotype haplotype = new Haplotype(baseString.getBytes(), baseString.equals(reference));
haplotype.setGenomeLocation(genomeLocParser.createGenomeLoc("1", 1, reference.length()));
try {
haplotype.setCigar(c.toCigar(reference.length()));
} catch (final RuntimeException ex) {
c.applyTo(reference);
c.toCigar(reference.length());
throw new RuntimeException("" + c + " " + ex.getMessage(), ex);
}
result.add(haplotype);
}
return result;
}
use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class PairHMMUnitTest method testLikelihoodsFromHaplotypes.
@Test(dataProvider = "JustHMMProvider")
public void testLikelihoodsFromHaplotypes(final PairHMM hmm) {
final int readSize = 10;
final int refSize = 20;
final byte[] readBases = Utils.dupBytes((byte) 'A', readSize);
final byte[] refBases = Utils.dupBytes((byte) 'A', refSize);
final byte baseQual = 20;
final byte insQual = 100;
final byte gcp = 100;
final Haplotype refH = new Haplotype(refBases, true);
final byte[] readQuals = Utils.dupBytes(baseQual, readBases.length);
final List<GATKRead> reads = Arrays.asList(ArtificialReadUtils.createArtificialRead(readBases, readQuals, readBases.length + "M"));
final Map<GATKRead, byte[]> gpcs = buildGapContinuationPenalties(reads, gcp);
hmm.computeLog10Likelihoods(matrix(Arrays.asList(refH)), Collections.emptyList(), gpcs);
Assert.assertEquals(hmm.getLogLikelihoodArray(), null);
hmm.computeLog10Likelihoods(matrix(Arrays.asList(refH)), reads, gpcs);
final double expected = getExpectedMatchingLogLikelihood(readBases, refBases, baseQual, insQual);
final double[] la = hmm.getLogLikelihoodArray();
Assert.assertEquals(la.length, 1);
Assert.assertEquals(la[0], expected, 1e-3, "Likelihoods should sum to just the error prob of the read " + String.format("readSize=%d refSize=%d", readSize, refSize));
}
use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class VectorPairHMMUnitTest method testLikelihoodsFromHaplotypes.
@Test(dataProvider = "JustHMMProvider")
public void testLikelihoodsFromHaplotypes(final PairHMM hmm, Boolean loaded) {
// skip if not loaded
if (!loaded.booleanValue()) {
throw new SkipException("AVX PairHMM is not supported on this system or the library is not available");
}
BasicInputParser parser = null;
try {
parser = new BasicInputParser(true, new FileInputStream(pairHMMTestData));
} catch (FileNotFoundException e) {
Assert.fail("PairHMM test data not found : " + pairHMMTestData);
}
while (parser.hasNext()) {
String[] tokens = parser.next();
final Haplotype hap = new Haplotype(tokens[0].getBytes(), true);
final byte[] bases = tokens[1].getBytes();
final byte[] baseQuals = normalize(tokens[2].getBytes(), 6);
final byte[] insertionQuals = normalize(tokens[3].getBytes());
final byte[] deletionQuals = normalize(tokens[4].getBytes());
final byte[] gcp = normalize(tokens[5].getBytes());
final double expectedResult = Double.parseDouble(tokens[6]);
final int readLength = bases.length;
final GATKRead read = ArtificialReadUtils.createArtificialRead(bases, baseQuals, readLength + "M");
ReadUtils.setInsertionBaseQualities(read, insertionQuals);
ReadUtils.setDeletionBaseQualities(read, deletionQuals);
final Map<GATKRead, byte[]> gpcs = new LinkedHashMap<>(readLength);
gpcs.put(read, gcp);
hmm.initialize(Arrays.asList(hap), null, 0, 0);
hmm.computeLog10Likelihoods(matrix(Arrays.asList(hap)), Arrays.asList(read), gpcs);
final double[] la = hmm.getLogLikelihoodArray();
Assert.assertEquals(la[0], expectedResult, 1e-5, "Likelihood not in expected range.");
}
hmm.close();
}
use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.
the class ReferenceConfidenceModel method createReferenceHaplotype.
/**
* Create a reference haplotype for an active region
*
* @param activeRegion the active region
* @param refBases the ref bases
* @param paddedReferenceLoc the location spanning of the refBases -- can be longer than activeRegion.getLocation()
* @return a reference haplotype
*/
public static Haplotype createReferenceHaplotype(final AssemblyRegion activeRegion, final byte[] refBases, final SimpleInterval paddedReferenceLoc) {
Utils.nonNull(activeRegion, "null region");
Utils.nonNull(refBases, "null refBases");
Utils.nonNull(paddedReferenceLoc, "null paddedReferenceLoc");
final int alignmentStart = activeRegion.getExtendedSpan().getStart() - paddedReferenceLoc.getStart();
if (alignmentStart < 0) {
throw new IllegalStateException("Bad alignment start in createReferenceHaplotype " + alignmentStart);
}
final Haplotype refHaplotype = new Haplotype(refBases, true);
refHaplotype.setAlignmentStartHapwrtRef(alignmentStart);
final Cigar c = new Cigar();
c.add(new CigarElement(refHaplotype.getBases().length, CigarOperator.M));
refHaplotype.setCigar(c);
return refHaplotype;
}
Aggregations