Search in sources :

Example 16 with Haplotype

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;
}
Also used : ReadThreadingGraph(org.broadinstitute.hellbender.tools.walkers.haplotypecaller.readthreading.ReadThreadingGraph) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

Example 17 with Haplotype

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;
}
Also used : ArrayList(java.util.ArrayList) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

Example 18 with Haplotype

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));
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 19 with Haplotype

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();
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) BasicInputParser(org.broadinstitute.hellbender.utils.text.parsers.BasicInputParser) FileNotFoundException(java.io.FileNotFoundException) FileInputStream(java.io.FileInputStream) SkipException(org.testng.SkipException) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 20 with Haplotype

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;
}
Also used : Cigar(htsjdk.samtools.Cigar) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) CigarElement(htsjdk.samtools.CigarElement)

Aggregations

Haplotype (org.broadinstitute.hellbender.utils.haplotype.Haplotype)88 Test (org.testng.annotations.Test)31 GATKRead (org.broadinstitute.hellbender.utils.read.GATKRead)28 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)28 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)12 DataProvider (org.testng.annotations.DataProvider)10 VariantContext (htsjdk.variant.variantcontext.VariantContext)9 KBestHaplotype (org.broadinstitute.hellbender.tools.walkers.haplotypecaller.graphs.KBestHaplotype)8 ArrayList (java.util.ArrayList)7 Cigar (htsjdk.samtools.Cigar)6 AssemblyRegion (org.broadinstitute.hellbender.engine.AssemblyRegion)6 HomogeneousPloidyModel (org.broadinstitute.hellbender.tools.walkers.genotyper.HomogeneousPloidyModel)6 IndependentSampleGenotypesModel (org.broadinstitute.hellbender.tools.walkers.genotyper.IndependentSampleGenotypesModel)6 PloidyModel (org.broadinstitute.hellbender.tools.walkers.genotyper.PloidyModel)6 ReadThreadingGraph (org.broadinstitute.hellbender.tools.walkers.haplotypecaller.readthreading.ReadThreadingGraph)6 SampleList (org.broadinstitute.hellbender.utils.genotyper.SampleList)6 EventMap (org.broadinstitute.hellbender.utils.haplotype.EventMap)6 CigarElement (htsjdk.samtools.CigarElement)4 Allele (htsjdk.variant.variantcontext.Allele)4 File (java.io.File)4