finders.add(haplotypeFinder);
final Iterator<KBestHaplotype> bestHaplotypes = haplotypeFinder.iterator(numBestHaplotypesPerGraph);
while (bestHaplotypes.hasNext()) {
final KBestHaplotype kBestHaplotype = bestHaplotypes.next();
final Haplotype h = kBestHaplotype.haplotype();
if( !returnHaplotypes.contains(h) ) {
final Cigar cigar = CigarUtils.calculateCigar(refHaplotype.getBases(),h.getBases());
if ( cigar == null ) {
failedCigars++; // couldn't produce a meaningful alignment of haplotype to reference, fail quietly
continue;
} else if( cigar.isEmpty() ) {
throw new IllegalStateException("Smith-Waterman alignment failure. Cigar = " + cigar + " with reference length " + cigar.getReferenceLength() +
" but expecting reference length of " + refHaplotype.getCigar().getReferenceLength());
} else if ( pathIsTooDivergentFromReference(cigar) || cigar.getReferenceLength() < MIN_HAPLOTYPE_REFERENCE_LENGTH ) {
// N cigar elements means that a bubble was too divergent from the reference so skip over this path
continue;
} else if( cigar.getReferenceLength() != refHaplotype.getCigar().getReferenceLength() ) { // SW failure
throw new IllegalStateException("Smith-Waterman alignment failure. Cigar = " + cigar + " with reference length "
+ cigar.getReferenceLength() + " but expecting reference length of " + refHaplotype.getCigar().getReferenceLength()
+ " ref = " + refHaplotype + " path " + new String(h.getBases()));
}
h.setCigar(cigar);
h.setAlignmentStartHapwrtRef(activeRegionStart);
h.setGenomeLocation(activeRegionWindow);
returnHaplotypes.add(h);
assemblyResultSet.add(h, assemblyResultByGraph.get(graph));
if ( debug )
logger.info("Adding haplotype " + h.getCigar() + " from graph with kmer " + graph.getKmerSize());
}
}
}
// Make sure that the ref haplotype is amongst the return haplotypes and calculate its score as
// the first returned by any finder.
if (!returnHaplotypes.contains(refHaplotype)) {
double refScore = Double.NaN;
for (final KBestHaplotypeFinder finder : finders) {
final double candidate = finder.score(refHaplotype);
if (Double.isNaN(candidate)) continue;
refScore = candidate;
break;
}
refHaplotype.setScore(refScore);
returnHaplotypes.add(refHaplotype);
}
if (failedCigars != 0)
logger.debug(String.format("failed to align some haplotypes (%d) back to the reference (loc=%s); these will be ignored.",failedCigars,refLoc.toString()));
if( debug ) {
if( returnHaplotypes.size() > 1 ) {
logger.info("Found " + returnHaplotypes.size() + " candidate haplotypes of " + returnHaplotypes.size() + " possible combinations to evaluate every read against.");
} else {
logger.info("Found only the reference haplotype in the assembly graph.");
}
for( final Haplotype h : returnHaplotypes ) {
logger.info( h.toString() );
logger.info( "> Cigar = " + h.getCigar() + " : " + h.getCigar().getReferenceLength() + " score " + h.getScore() + " ref " + h.isReference());
}
}
return new ArrayList<>(returnHaplotypes);