Search in sources :

Example 46 with Haplotype

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

the class AssemblyResultSet method calculateOriginalByTrimmedHaplotypes.

private Map<Haplotype, Haplotype> calculateOriginalByTrimmedHaplotypes(final AssemblyRegion trimmedAssemblyRegion) {
    if (debug) {
        logger.info("Trimming active region " + getRegionForGenotyping() + " with " + getHaplotypeCount() + " haplotypes");
    }
    final List<Haplotype> haplotypeList = getHaplotypeList();
    // trim down the haplotypes
    final Map<Haplotype, Haplotype> originalByTrimmedHaplotypes = trimDownHaplotypes(trimmedAssemblyRegion, haplotypeList);
    // create the final list of trimmed haplotypes
    final List<Haplotype> trimmedHaplotypes = new ArrayList<>(originalByTrimmedHaplotypes.keySet());
    // resort the trimmed haplotypes.
    Collections.sort(trimmedHaplotypes, Haplotype.SIZE_AND_BASE_ORDER);
    final Map<Haplotype, Haplotype> sortedOriginalByTrimmedHaplotypes = mapOriginalToTrimmed(originalByTrimmedHaplotypes, trimmedHaplotypes);
    if (debug) {
        logger.info("Trimmed region to " + trimmedAssemblyRegion.getSpan() + " size " + trimmedAssemblyRegion.getSpan().size() + " reduced number of haplotypes from " + haplotypeList.size() + " to only " + trimmedHaplotypes.size());
        for (final Haplotype remaining : trimmedHaplotypes) {
            logger.info("Remains: " + remaining + " cigar " + remaining.getCigar());
        }
    }
    return sortedOriginalByTrimmedHaplotypes;
}
Also used : Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

Example 47 with Haplotype

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

the class AssemblyBasedCallerUtils method assembleReads.

/**
     * High-level function that runs the assembler on the given region's reads,
     * returning a data structure with the resulting information needed
     * for further HC steps
     */
public static AssemblyResultSet assembleReads(final AssemblyRegion region, final List<VariantContext> givenAlleles, final AssemblyBasedCallerArgumentCollection argumentCollection, final SAMFileHeader header, final SampleList sampleList, final Logger logger, final ReferenceSequenceFile referenceReader, final ReadThreadingAssembler assemblyEngine) {
    finalizeRegion(region, argumentCollection.errorCorrectReads, argumentCollection.dontUseSoftClippedBases, (byte) (argumentCollection.minBaseQualityScore - 1), header, sampleList);
    if (argumentCollection.debug) {
        logger.info("Assembling " + region.getSpan() + " with " + region.size() + " reads:    (with overlap region = " + region.getExtendedSpan() + ")");
    }
    final byte[] fullReferenceWithPadding = region.getAssemblyRegionReference(referenceReader, REFERENCE_PADDING_FOR_ASSEMBLY);
    final SimpleInterval paddedReferenceLoc = getPaddedReferenceLoc(region, REFERENCE_PADDING_FOR_ASSEMBLY, referenceReader);
    final Haplotype referenceHaplotype = createReferenceHaplotype(region, paddedReferenceLoc, referenceReader);
    final ReadErrorCorrector readErrorCorrector = argumentCollection.errorCorrectReads ? new ReadErrorCorrector(argumentCollection.assemblerArgs.kmerLengthForReadErrorCorrection, HaplotypeCallerEngine.MIN_TAIL_QUALITY_WITH_ERROR_CORRECTION, argumentCollection.assemblerArgs.minObservationsForKmerToBeSolid, argumentCollection.debug, fullReferenceWithPadding) : null;
    try {
        final AssemblyResultSet assemblyResultSet = assemblyEngine.runLocalAssembly(region, referenceHaplotype, fullReferenceWithPadding, paddedReferenceLoc, givenAlleles, readErrorCorrector, header);
        assemblyResultSet.debugDump(logger);
        return assemblyResultSet;
    } catch (final Exception e) {
        // Capture any exception that might be thrown, and write out the assembly failure BAM if requested
        if (argumentCollection.captureAssemblyFailureBAM) {
            try (final SAMFileWriter writer = ReadUtils.createCommonSAMWriter(new File("assemblyFailure.bam"), null, header, false, false, false)) {
                for (final GATKRead read : region.getReads()) {
                    writer.addAlignment(read.convertToSAMRecord(header));
                }
            }
        }
        throw e;
    }
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) SAMFileWriter(htsjdk.samtools.SAMFileWriter) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile) CachingIndexedFastaSequenceFile(org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile) File(java.io.File) FileNotFoundException(java.io.FileNotFoundException) UserException(org.broadinstitute.hellbender.exceptions.UserException)

Example 48 with Haplotype

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

the class AssemblyBasedCallerUtils method realignReadsToTheirBestHaplotype.

/**
     * Returns a map with the original read as a key and the realigned read as the value.
     * <p>
     *     Missing keys or equivalent key and value pairs mean that the read was not realigned.
     * </p>
     * @return never {@code null}
     */
public static Map<GATKRead, GATKRead> realignReadsToTheirBestHaplotype(final ReadLikelihoods<Haplotype> originalReadLikelihoods, final Haplotype refHaplotype, final Locatable paddedReferenceLoc) {
    final Collection<ReadLikelihoods<Haplotype>.BestAllele<Haplotype>> bestAlleles = originalReadLikelihoods.bestAlleles();
    final Map<GATKRead, GATKRead> result = new HashMap<>(bestAlleles.size());
    for (final ReadLikelihoods<Haplotype>.BestAllele<Haplotype> bestAllele : bestAlleles) {
        final GATKRead originalRead = bestAllele.read;
        final Haplotype bestHaplotype = bestAllele.allele;
        final boolean isInformative = bestAllele.isInformative();
        final GATKRead realignedRead = AlignmentUtils.createReadAlignedToRef(originalRead, bestHaplotype, refHaplotype, paddedReferenceLoc.getStart(), isInformative);
        result.put(originalRead, realignedRead);
    }
    return result;
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype) ReadLikelihoods(org.broadinstitute.hellbender.utils.genotyper.ReadLikelihoods)

Example 49 with Haplotype

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

the class RandomLikelihoodCalculationEngine method computeReadLikelihoods.

@Override
public ReadLikelihoods<Haplotype> computeReadLikelihoods(final AssemblyResultSet assemblyResultSet, final SampleList samples, final Map<String, List<GATKRead>> reads) {
    Utils.nonNull(assemblyResultSet, "assemblyResultSet is null");
    Utils.nonNull(samples, "samples is null");
    Utils.nonNull(reads, "perSampleReadList is null");
    final AlleleList<Haplotype> haplotypes = new IndexedAlleleList<>(assemblyResultSet.getHaplotypeList());
    final ReadLikelihoods<Haplotype> result = new ReadLikelihoods<>(samples, haplotypes, reads);
    final Random rnd = Utils.getRandomGenerator();
    final int sampleCount = samples.numberOfSamples();
    final int alleleCount = haplotypes.numberOfAlleles();
    for (int i = 0; i < sampleCount; i++) {
        final List<GATKRead> sampleReads = result.sampleReads(i);
        final int readCount = sampleReads.size();
        final LikelihoodMatrix<Haplotype> sampleLikelihoods = result.sampleMatrix(i);
        for (int a = 0; a < alleleCount; a++) {
            for (int r = 0; r < readCount; r++) {
                sampleLikelihoods.set(a, r, -Math.abs(rnd.nextDouble()));
            }
        }
    }
    return result;
}
Also used : GATKRead(org.broadinstitute.hellbender.utils.read.GATKRead) Random(java.util.Random) Haplotype(org.broadinstitute.hellbender.utils.haplotype.Haplotype)

Example 50 with Haplotype

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

the class HaplotypeCallerEngine method referenceModelForNoVariation.

/**
     * Create an ref model result (ref model or no calls depending on mode) for an active region without any variation
     * (not is active, or assembled to just ref)
     *
     * @param region the region to return a no-variation result
     * @param needsToBeFinalized should the region be finalized before computing the ref model (should be false if already done)
     * @return a list of variant contexts (can be empty) to emit for this ref region
     */
private List<VariantContext> referenceModelForNoVariation(final AssemblyRegion region, final boolean needsToBeFinalized) {
    if (emitReferenceConfidence()) {
        //TODO - perhaps we can remove the last parameter of this method and the three lines bellow?
        if (needsToBeFinalized) {
            finalizeRegion(region);
        }
        filterNonPassingReads(region);
        final SimpleInterval paddedLoc = region.getExtendedSpan();
        final Haplotype refHaplotype = AssemblyBasedCallerUtils.createReferenceHaplotype(region, paddedLoc, referenceReader);
        final List<Haplotype> haplotypes = Collections.singletonList(refHaplotype);
        return referenceConfidenceModel.calculateRefConfidence(refHaplotype, haplotypes, paddedLoc, region, createDummyStratifiedReadMap(refHaplotype, samplesList, region), genotypingEngine.getPloidyModel(), Collections.emptyList());
    } else {
        return NO_CALLS;
    }
}
Also used : SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) 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