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