use of org.broadinstitute.hellbender.exceptions.GATKException in project gatk by broadinstitute.
the class Tranche method tranchesString.
/**
* Returns an appropriately formatted string representing the raw tranches file on disk.
*
* @param tranches
* @return
*/
public static String tranchesString(final List<Tranche> tranches) {
try (final ByteArrayOutputStream bytes = new ByteArrayOutputStream();
final PrintStream stream = new PrintStream(bytes)) {
if (tranches.size() > 1)
Collections.sort(tranches, new TrancheTruthSensitivityComparator());
stream.println("# Variant quality score tranches file");
stream.println("# Version number " + CURRENT_VERSION);
stream.println("targetTruthSensitivity,numKnown,numNovel,knownTiTv,novelTiTv,minVQSLod,filterName,model,accessibleTruthSites,callsAtTruthSites,truthSensitivity");
Tranche prev = null;
for (Tranche t : tranches) {
stream.printf("%.2f,%d,%d,%.4f,%.4f,%.4f,VQSRTranche%s%.2fto%.2f,%s,%d,%d,%.4f%n", t.targetTruthSensitivity, t.numKnown, t.numNovel, t.knownTiTv, t.novelTiTv, t.minVQSLod, t.model.toString(), (prev == null ? 0.0 : prev.targetTruthSensitivity), t.targetTruthSensitivity, t.model.toString(), t.accessibleTruthSites, t.callsAtTruthSites, t.getTruthSensitivity());
prev = t;
}
return bytes.toString();
} catch (IOException e) {
throw new GATKException("IOException while converting tranche to a string");
}
}
use of org.broadinstitute.hellbender.exceptions.GATKException in project gatk by broadinstitute.
the class MRUCachingSAMSequenceDictionary method updateCache.
/**
* The key algorithm. Given a new record, update the last used record, contig
* name, and index.
*
* @param contig the contig we want to look up. If null, index is used instead
* @param index the contig index we want to look up. Only used if contig is null
* @throws GATKException if index isn't present in the dictionary
* @return the SAMSequenceRecord for contig / index
*/
private SAMSequenceRecord updateCache(final String contig, int index) {
SAMSequenceRecord rec = contig == null ? dict.getSequence(index) : dict.getSequence(contig);
if (rec == null) {
throw new GATKException("BUG: requested unknown contig=" + contig + " index=" + index);
} else {
lastSSR = rec;
lastContig = rec.getSequenceName();
lastIndex = rec.getSequenceIndex();
return rec;
}
}
use of org.broadinstitute.hellbender.exceptions.GATKException in project gatk by broadinstitute.
the class ValidateSamFile method doWork.
@Override
protected Object doWork() {
IOUtil.assertFileIsReadable(INPUT);
ReferenceSequenceFile reference = null;
if (REFERENCE_SEQUENCE != null) {
IOUtil.assertFileIsReadable(REFERENCE_SEQUENCE);
reference = ReferenceSequenceFileFactory.getReferenceSequenceFile(REFERENCE_SEQUENCE);
}
final PrintWriter out;
if (OUTPUT != null) {
IOUtil.assertFileIsWritable(OUTPUT);
try {
out = new PrintWriter(OUTPUT);
} catch (FileNotFoundException e) {
// we already asserted this so we should not get here
throw new GATKException("Unexpected exception", e);
}
} else {
out = new PrintWriter(System.out);
}
final SamReaderFactory factory = SamReaderFactory.makeDefault().referenceSequence(REFERENCE_SEQUENCE).validationStringency(ValidationStringency.SILENT).enable(SamReaderFactory.Option.VALIDATE_CRC_CHECKSUMS);
final SamReader samReader = factory.open(INPUT);
if (samReader.type() != SamReader.Type.BAM_TYPE)
VALIDATE_INDEX = false;
factory.setOption(SamReaderFactory.Option.CACHE_FILE_BASED_INDEXES, VALIDATE_INDEX);
factory.reapplyOptions(samReader);
final SamFileValidator validator = new SamFileValidator(out, MAX_OPEN_TEMP_FILES);
validator.setErrorsToIgnore(IGNORE);
if (IGNORE_WARNINGS) {
validator.setIgnoreWarnings(IGNORE_WARNINGS);
}
if (MODE == Mode.SUMMARY) {
validator.setVerbose(false, 0);
} else {
validator.setVerbose(true, MAX_OUTPUT);
}
if (IS_BISULFITE_SEQUENCED) {
validator.setBisulfiteSequenced(IS_BISULFITE_SEQUENCED);
}
if (VALIDATE_INDEX) {
validator.setIndexValidationStringency(BamIndexValidator.IndexValidationStringency.EXHAUSTIVE);
}
if (IOUtil.isRegularPath(INPUT)) {
// Do not check termination if reading from a stream
validator.validateBamFileTermination(INPUT);
}
boolean isValid = false;
switch(MODE) {
case SUMMARY:
isValid = validator.validateSamFileSummary(samReader, reference);
break;
case VERBOSE:
isValid = validator.validateSamFileVerbose(samReader, reference);
break;
}
out.flush();
return isValid;
}
use of org.broadinstitute.hellbender.exceptions.GATKException in project gatk by broadinstitute.
the class CachingIndexedFastaSequenceFile method getSubsequenceAt.
/**
* Gets the subsequence of the contig in the range [start,stop]
*
* Uses the sequence cache if possible, or updates the cache to handle the request. If the range
* is larger than the cache itself, just loads the sequence directly, not changing the cache at all
*
* @param contig Contig whose subsequence to retrieve.
* @param start inclusive, 1-based start of region.
* @param stop inclusive, 1-based stop of region.
* @return The partial reference sequence associated with this range. If preserveCase is false, then
* all of the bases in the ReferenceSequence returned by this method will be upper cased.
*/
@Override
public ReferenceSequence getSubsequenceAt(final String contig, long start, final long stop) {
final ReferenceSequence result;
if ((stop - start) >= cacheSize) {
cacheMisses++;
result = super.getSubsequenceAt(contig, start, stop);
if (!preserveCase)
StringUtil.toUpperCase(result.getBases());
if (!preserveIUPAC)
BaseUtils.convertIUPACtoN(result.getBases(), true, start < 1);
} else {
// todo -- potential optimization is to check if contig.name == contig, as this in general will be true
SAMSequenceRecord contigInfo = super.getSequenceDictionary().getSequence(contig);
if (contigInfo == null) {
throw new UserException.MissingContigInSequenceDictionary(contig, super.getSequenceDictionary());
}
if (stop > contigInfo.getSequenceLength())
throw new SAMException("Query asks for data past end of contig. Query contig " + contig + " start:" + start + " stop:" + stop + " contigLength:" + contigInfo.getSequenceLength());
if (start < cache.start || stop > cache.stop || cache.seq == null || cache.seq.getContigIndex() != contigInfo.getSequenceIndex()) {
cacheMisses++;
cache.start = Math.max(start - cacheMissBackup, 0);
cache.stop = Math.min(start + cacheSize + cacheMissBackup, contigInfo.getSequenceLength());
cache.seq = super.getSubsequenceAt(contig, cache.start, cache.stop);
// convert all of the bases in the sequence to upper case if we aren't preserving cases
if (!preserveCase)
StringUtil.toUpperCase(cache.seq.getBases());
if (!preserveIUPAC)
BaseUtils.convertIUPACtoN(cache.seq.getBases(), true, cache.start == 0);
} else {
cacheHits++;
}
// at this point we determine where in the cache we want to extract the requested subsequence
final int cacheOffsetStart = (int) (start - cache.start);
final int cacheOffsetStop = (int) (stop - start + cacheOffsetStart + 1);
try {
result = new ReferenceSequence(cache.seq.getName(), cache.seq.getContigIndex(), Arrays.copyOfRange(cache.seq.getBases(), cacheOffsetStart, cacheOffsetStop));
} catch (ArrayIndexOutOfBoundsException e) {
throw new GATKException(String.format("BUG: bad array indexing. Cache start %d and end %d, request start %d end %d, offset start %d and end %d, base size %d", cache.start, cache.stop, start, stop, cacheOffsetStart, cacheOffsetStop, cache.seq.getBases().length), e);
}
}
// for debugging -- print out our efficiency if requested
if (PRINT_EFFICIENCY && (getCacheHits() + getCacheMisses()) % PRINT_FREQUENCY == 0)
printEfficiency(Level.INFO);
return result;
}
use of org.broadinstitute.hellbender.exceptions.GATKException in project gatk by broadinstitute.
the class EventMap method processCigarForInitialEvents.
protected void processCigarForInitialEvents() {
final Cigar cigar = haplotype.getCigar();
final byte[] alignment = haplotype.getBases();
int refPos = haplotype.getAlignmentStartHapwrtRef();
if (refPos < 0) {
return;
}
// Protection against SW failures
final List<VariantContext> proposedEvents = new ArrayList<>();
int alignmentPos = 0;
for (int cigarIndex = 0; cigarIndex < cigar.numCigarElements(); cigarIndex++) {
final CigarElement ce = cigar.getCigarElement(cigarIndex);
final int elementLength = ce.getLength();
switch(ce.getOperator()) {
case I:
{
if (refPos > 0) {
// protect against trying to create insertions/deletions at the beginning of a contig
final List<Allele> insertionAlleles = new ArrayList<>();
final int insertionStart = refLoc.getStart() + refPos - 1;
final byte refByte = ref[refPos - 1];
if (BaseUtils.isRegularBase(refByte)) {
insertionAlleles.add(Allele.create(refByte, true));
}
if (cigarIndex == 0 || cigarIndex == cigar.numCigarElements() - 1) {
// if the insertion isn't completely resolved in the haplotype, skip it
// note this used to emit SYMBOLIC_UNASSEMBLED_EVENT_ALLELE but that seems dangerous
} else {
byte[] insertionBases = {};
// add the padding base
insertionBases = ArrayUtils.add(insertionBases, ref[refPos - 1]);
insertionBases = ArrayUtils.addAll(insertionBases, Arrays.copyOfRange(alignment, alignmentPos, alignmentPos + elementLength));
if (BaseUtils.isAllRegularBases(insertionBases)) {
insertionAlleles.add(Allele.create(insertionBases, false));
}
}
if (insertionAlleles.size() == 2) {
// found a proper ref and alt allele
proposedEvents.add(new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), insertionStart, insertionStart, insertionAlleles).make());
}
}
alignmentPos += elementLength;
break;
}
case S:
{
alignmentPos += elementLength;
break;
}
case D:
{
if (refPos > 0) {
// protect against trying to create insertions/deletions at the beginning of a contig
// add padding base
final byte[] deletionBases = Arrays.copyOfRange(ref, refPos - 1, refPos + elementLength);
final List<Allele> deletionAlleles = new ArrayList<>();
final int deletionStart = refLoc.getStart() + refPos - 1;
final byte refByte = ref[refPos - 1];
if (BaseUtils.isRegularBase(refByte) && BaseUtils.isAllRegularBases(deletionBases)) {
deletionAlleles.add(Allele.create(deletionBases, true));
deletionAlleles.add(Allele.create(refByte, false));
proposedEvents.add(new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), deletionStart, deletionStart + elementLength, deletionAlleles).make());
}
}
refPos += elementLength;
break;
}
case M:
case EQ:
case X:
{
for (int iii = 0; iii < elementLength; iii++) {
final byte refByte = ref[refPos];
final byte altByte = alignment[alignmentPos];
if (refByte != altByte) {
// SNP!
if (BaseUtils.isRegularBase(refByte) && BaseUtils.isRegularBase(altByte)) {
final List<Allele> snpAlleles = new ArrayList<>();
snpAlleles.add(Allele.create(refByte, true));
snpAlleles.add(Allele.create(altByte, false));
proposedEvents.add(new VariantContextBuilder(sourceNameToAdd, refLoc.getContig(), refLoc.getStart() + refPos, refLoc.getStart() + refPos, snpAlleles).make());
}
}
refPos++;
alignmentPos++;
}
break;
}
case N:
case H:
case P:
default:
throw new GATKException("Unsupported cigar operator created during SW alignment: " + ce.getOperator());
}
}
for (final VariantContext proposedEvent : proposedEvents) addVC(proposedEvent, true);
}
Aggregations