Search in sources :

Example 41 with Haplotype

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

the class ReadThreadingAssembler method createGraph.

/**
     * Creates the sequence graph for the given kmerSize
     *
     * @param reads            reads to use
     * @param refHaplotype     reference haplotype
     * @param kmerSize         kmer size
     * @param activeAlleleHaplotypes the GGA haplotypes to inject into the graph
     * @param allowLowComplexityGraphs if true, do not check for low-complexity graphs
     * @param allowNonUniqueKmersInRef if true, do not fail if the reference has non-unique kmers
     * @return sequence graph or null if one could not be created (e.g. because it contains cycles or too many paths or is low complexity)
     */
private AssemblyResult createGraph(final Iterable<GATKRead> reads, final Haplotype refHaplotype, final int kmerSize, final Iterable<Haplotype> activeAlleleHaplotypes, final boolean allowLowComplexityGraphs, final boolean allowNonUniqueKmersInRef, final SAMFileHeader header) {
    if (refHaplotype.length() < kmerSize) {
        // happens in cases where the assembled region is just too small
        return new AssemblyResult(AssemblyResult.Status.FAILED, null, null);
    }
    if (!allowNonUniqueKmersInRef && !ReadThreadingGraph.determineNonUniqueKmers(new ReadThreadingGraph.SequenceForKmers("ref", refHaplotype.getBases(), 0, refHaplotype.getBases().length, 1, true), kmerSize).isEmpty()) {
        if (debug) {
            logger.info("Not using kmer size of " + kmerSize + " in read threading assembler because reference contains non-unique kmers");
        }
        return null;
    }
    final ReadThreadingGraph rtgraph = new ReadThreadingGraph(kmerSize, debugGraphTransformations, MIN_BASE_QUALITY_TO_USE_IN_ASSEMBLY, numPruningSamples);
    rtgraph.setThreadingStartOnlyAtExistingVertex(!recoverDanglingBranches);
    // add the reference sequence to the graph
    rtgraph.addSequence("ref", refHaplotype.getBases(), true);
    // add the artificial GGA haplotypes to the graph
    int hapCount = 0;
    for (final Haplotype h : activeAlleleHaplotypes) {
        rtgraph.addSequence("activeAllele" + hapCount++, h.getBases(), GGA_MODE_ARTIFICIAL_COUNTS, false);
    }
    // Next pull kmers out of every read and throw them on the graph
    for (final GATKRead read : reads) {
        rtgraph.addRead(read, header);
    }
    // actually build the read threading graph
    rtgraph.buildGraphIfNecessary();
    // sanity check: make sure there are no cycles in the graph
    if (rtgraph.hasCycles()) {
        if (debug) {
            logger.info("Not using kmer size of " + kmerSize + " in read threading assembler because it contains a cycle");
        }
        return null;
    }
    // sanity check: make sure the graph had enough complexity with the given kmer
    if (!allowLowComplexityGraphs && rtgraph.isLowComplexity()) {
        if (debug) {
            logger.info("Not using kmer size of " + kmerSize + " in read threading assembler because it does not produce a graph with enough complexity");
        }
        return null;
    }
    return getAssemblyResult(refHaplotype, kmerSize, rtgraph);
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) AssemblyResult(org.broadinstitute.hellbender.tools.walkers.haplotypecaller.AssemblyResult) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

Example 42 with Haplotype

use of org.broadinstitute.hellbender.utils.haplotype.Haplotype in project gatk-protected 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)

Example 43 with Haplotype

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

the class KBestHaplotype method haplotype.

/**
     * Returns the solution haplotype.
     *
     * @return never {@code null}.
     */
public final Haplotype haplotype() {
    if (haplotype != null) {
        return haplotype;
    }
    haplotype = new Haplotype(bases(), isReference());
    if (score() > 0) {
        throw new IllegalStateException("score cannot be greater than 0: " + score());
    }
    haplotype.setScore(score());
    return haplotype;
}
Also used : Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

Example 44 with Haplotype

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

the class AssemblyResultSet method trimTo.

/**
     * Trims an assembly result set down based on a new set of trimmed haplotypes.
     *
     * @param trimmedAssemblyRegion the trimmed down active region.
     *
     * @throws NullPointerException if any argument in {@code null} or
     *      if there are {@code null} entries in {@code originalByTrimmedHaplotypes} for trimmed haplotype keys.
     * @throws IllegalArgumentException if there is no reference haplotype amongst the trimmed ones.
     *
     * @return never {@code null}, a new trimmed assembly result set.
     */
public AssemblyResultSet trimTo(final AssemblyRegion trimmedAssemblyRegion) {
    final Map<Haplotype, Haplotype> originalByTrimmedHaplotypes = calculateOriginalByTrimmedHaplotypes(trimmedAssemblyRegion);
    if (refHaplotype == null) {
        throw new IllegalStateException("refHaplotype is null");
    }
    Utils.nonNull(trimmedAssemblyRegion);
    final AssemblyResultSet result = new AssemblyResultSet();
    for (final Haplotype trimmed : originalByTrimmedHaplotypes.keySet()) {
        final Haplotype original = originalByTrimmedHaplotypes.get(trimmed);
        if (original == null) {
            throw new IllegalStateException("all trimmed haplotypes must have an original one");
        }
        final AssemblyResult as = assemblyResultByHaplotype.get(original);
        if (as == null) {
            result.add(trimmed);
        } else {
            result.add(trimmed, as);
        }
    }
    result.setRegionForGenotyping(trimmedAssemblyRegion);
    result.setFullReferenceWithPadding(fullReferenceWithPadding);
    result.setPaddedReferenceLoc(paddedReferenceLoc);
    if (result.refHaplotype == null) {
        throw new IllegalStateException("missing reference haplotype in the trimmed set");
    }
    result.wasTrimmed = true;
    return result;
}
Also used : Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

Example 45 with Haplotype

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

the class AssemblyResultSet method trimDownHaplotypes.

private Map<Haplotype, Haplotype> trimDownHaplotypes(final AssemblyRegion trimmedAssemblyRegion, final List<Haplotype> haplotypeList) {
    final Map<Haplotype, Haplotype> originalByTrimmedHaplotypes = new HashMap<>();
    for (final Haplotype h : haplotypeList) {
        final Haplotype trimmed = h.trim(trimmedAssemblyRegion.getExtendedSpan());
        if (trimmed != null) {
            if (originalByTrimmedHaplotypes.containsKey(trimmed)) {
                if (trimmed.isReference()) {
                    originalByTrimmedHaplotypes.remove(trimmed);
                    originalByTrimmedHaplotypes.put(trimmed, h);
                }
            } else {
                originalByTrimmedHaplotypes.put(trimmed, h);
            }
        } else if (h.isReference()) {
            throw new IllegalStateException("trimming eliminates the reference haplotype");
        } else if (debug) {
            logger.info("Throwing out haplotype " + h + " with cigar " + h.getCigar() + " because it starts with or ends with an insertion or deletion when trimmed to " + trimmedAssemblyRegion.getExtendedSpan());
        }
    }
    return originalByTrimmedHaplotypes;
}
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