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;
}
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;
}
}
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;
}
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;
}
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;
}
}
Aggregations