Search in sources :

Example 11 with Haplotype

use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.

the class VectorLoglessPairHMM method computeLog10Likelihoods.

/**
     * {@inheritDoc}
     */
@Override
public void computeLog10Likelihoods(final LikelihoodMatrix<Haplotype> logLikelihoods, final List<GATKRead> processedReads, final Map<GATKRead, byte[]> gcp) {
    if (processedReads.isEmpty()) {
        return;
    }
    if (doProfiling) {
        startTime = System.nanoTime();
    }
    int readListSize = processedReads.size();
    int numHaplotypes = logLikelihoods.numberOfAlleles();
    ReadDataHolder[] readDataArray = new ReadDataHolder[readListSize];
    int idx = 0;
    for (GATKRead read : processedReads) {
        readDataArray[idx] = new ReadDataHolder();
        readDataArray[idx].readBases = read.getBases();
        readDataArray[idx].readQuals = read.getBaseQualities();
        readDataArray[idx].insertionGOP = ReadUtils.getBaseInsertionQualities(read);
        readDataArray[idx].deletionGOP = ReadUtils.getBaseDeletionQualities(read);
        readDataArray[idx].overallGCP = gcp.get(read);
        ++idx;
    }
    //to store results
    mLogLikelihoodArray = new double[readListSize * numHaplotypes];
    if (doProfiling) {
        threadLocalSetupTimeDiff = (System.nanoTime() - startTime);
    }
    //for(reads)
    //   for(haplotypes)
    //       compute_full_prob()
    pairHmm.computeLikelihoods(readDataArray, mHaplotypeDataArray, mLogLikelihoodArray);
    int readIdx = 0;
    for (int r = 0; r < readListSize; r++) {
        int hapIdx = 0;
        for (final Haplotype haplotype : logLikelihoods.alleles()) {
            //Since the order of haplotypes in the List<Haplotype> and alleleHaplotypeMap is different,
            //get idx of current haplotype in the list and use this idx to get the right likelihoodValue
            final int idxInsideHaplotypeList = haplotypeToHaplotypeListIdxMap.get(haplotype);
            logLikelihoods.set(hapIdx, r, mLogLikelihoodArray[readIdx + idxInsideHaplotypeList]);
            ++hapIdx;
        }
        readIdx += numHaplotypes;
    }
    if (doProfiling) {
        threadLocalPairHMMComputeTimeDiff = (System.nanoTime() - startTime);
        pairHMMComputeTime += threadLocalPairHMMComputeTimeDiff;
        pairHMMSetupTime += threadLocalSetupTimeDiff;
    }
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) ReadDataHolder(org.broadinstitute.gatk.nativebindings.pairhmm.ReadDataHolder)

Example 12 with Haplotype

use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.

the class ReadThreadingAssemblerUnitTest method testAssemblyWithVariant.

private void testAssemblyWithVariant(final ReadThreadingAssembler assembler, final byte[] refBases, final SimpleInterval loc, final int nReadsToUse, final VariantContext site) {
    final String preRef = new String(refBases).substring(0, site.getStart());
    final String postRef = new String(refBases).substring(site.getEnd() + 1, refBases.length);
    final byte[] altBases = (preRef + site.getAlternateAllele(0).getBaseString() + postRef).getBytes();
    //        logger.warn("ref " + new String(refBases));
    //        logger.warn("alt " + new String(altBases));
    final List<GATKRead> reads = new LinkedList<>();
    for (int i = 0; i < nReadsToUse; i++) {
        final byte[] bases = altBases.clone();
        final byte[] quals = Utils.dupBytes((byte) 30, altBases.length);
        final String cigar = altBases.length + "M";
        final GATKRead read = ArtificialReadUtils.createArtificialRead(header, loc.getContig(), loc.getContig(), loc.getStart(), bases, quals, cigar);
        reads.add(read);
    }
    final Haplotype refHaplotype = new Haplotype(refBases, true);
    final Haplotype altHaplotype = new Haplotype(altBases, false);
    final List<Haplotype> haplotypes = assemble(assembler, refBases, loc, reads);
    Assert.assertEquals(haplotypes, Arrays.asList(refHaplotype, altHaplotype));
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) KBestHaplotype(org.broadinstitute.hellbender.tools.walkers.haplotypecaller.graphs.KBestHaplotype)

Example 13 with Haplotype

use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.

the class ReadThreadingAssemblerUnitTest method testAssembleRef.

@Test(dataProvider = "AssembleIntervalsData")
public void testAssembleRef(final ReadThreadingAssembler assembler, final SimpleInterval loc, final int nReadsToUse) {
    final byte[] refBases = seq.getSubsequenceAt(loc.getContig(), loc.getStart(), loc.getEnd()).getBases();
    final List<GATKRead> reads = new LinkedList<>();
    for (int i = 0; i < nReadsToUse; i++) {
        final byte[] bases = refBases.clone();
        final byte[] quals = Utils.dupBytes((byte) 30, refBases.length);
        final String cigar = refBases.length + "M";
        final GATKRead read = ArtificialReadUtils.createArtificialRead(header, loc.getContig(), loc.getContig(), loc.getStart(), bases, quals, cigar);
        reads.add(read);
    }
    // TODO -- generalize to all assemblers
    final Haplotype refHaplotype = new Haplotype(refBases, true);
    final List<Haplotype> haplotypes = assemble(assembler, refBases, loc, reads);
    Assert.assertEquals(haplotypes, Collections.singletonList(refHaplotype));
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) KBestHaplotype(org.broadinstitute.hellbender.tools.walkers.haplotypecaller.graphs.KBestHaplotype) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 14 with Haplotype

use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.

the class ReadThreadingAssemblerUnitTest method testSimpleAssembly.

@Test(dataProvider = "SimpleAssemblyTestData")
public void testSimpleAssembly(final String name, final ReadThreadingAssembler assembler, final SimpleInterval loc, final String ref, final String alt) {
    final byte[] refBases = ref.getBytes();
    final byte[] altBases = alt.getBytes();
    final List<GATKRead> reads = new LinkedList<>();
    for (int i = 0; i < 20; i++) {
        final byte[] bases = altBases.clone();
        final byte[] quals = Utils.dupBytes((byte) 30, altBases.length);
        final String cigar = altBases.length + "M";
        final GATKRead read = ArtificialReadUtils.createArtificialRead(header, loc.getContig(), loc.getContig(), loc.getStart(), bases, quals, cigar);
        reads.add(read);
    }
    final Haplotype refHaplotype = new Haplotype(refBases, true);
    final Haplotype altHaplotype = new Haplotype(altBases, false);
    final List<Haplotype> haplotypes = assemble(assembler, refBases, loc, reads);
    Assert.assertTrue(haplotypes.size() > 0, "Failed to find ref haplotype");
    Assert.assertEquals(haplotypes.get(0), refHaplotype);
    Assert.assertEquals(haplotypes.size(), 2, "Failed to find single alt haplotype");
    Assert.assertEquals(haplotypes.get(1), altHaplotype);
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) KBestHaplotype(org.broadinstitute.hellbender.tools.walkers.haplotypecaller.graphs.KBestHaplotype) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 15 with Haplotype

use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk by broadinstitute.

the class AssemblyRegionTestDataSet method cigarToHaplotype.

private Haplotype cigarToHaplotype(final String reference, final String cigar, final int offset, final boolean global) {
    final String sequence = applyCigar(reference, cigar, offset, global);
    final Haplotype haplotype = new Haplotype(sequence.getBytes(), reference.equals(sequence));
    haplotype.setGenomeLocation(genomeLocParser.createGenomeLoc("chr1", 1, reference.length()));
    haplotype.setCigar(Civar.fromCharSequence(cigar).toCigar(reference.length()));
    return haplotype;
}
Also used : Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

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