Search in sources :

Example 51 with GATKException

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");
    }
}
Also used : GATKException(org.broadinstitute.hellbender.exceptions.GATKException)

Example 52 with GATKException

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;
    }
}
Also used : SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) GATKException(org.broadinstitute.hellbender.exceptions.GATKException)

Example 53 with GATKException

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;
}
Also used : FileNotFoundException(java.io.FileNotFoundException) ReferenceSequenceFile(htsjdk.samtools.reference.ReferenceSequenceFile) GATKException(org.broadinstitute.hellbender.exceptions.GATKException) PrintWriter(java.io.PrintWriter)

Example 54 with GATKException

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;
}
Also used : SAMException(htsjdk.samtools.SAMException) SAMSequenceRecord(htsjdk.samtools.SAMSequenceRecord) ReferenceSequence(htsjdk.samtools.reference.ReferenceSequence) GATKException(org.broadinstitute.hellbender.exceptions.GATKException)

Example 55 with GATKException

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);
}
Also used : Cigar(htsjdk.samtools.Cigar) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VariantContext(htsjdk.variant.variantcontext.VariantContext) GATKException(org.broadinstitute.hellbender.exceptions.GATKException) CigarElement(htsjdk.samtools.CigarElement)

Aggregations

GATKException (org.broadinstitute.hellbender.exceptions.GATKException)96 IOException (java.io.IOException)19 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)13 CigarElement (htsjdk.samtools.CigarElement)12 ArrayList (java.util.ArrayList)10 UserException (org.broadinstitute.hellbender.exceptions.UserException)10 SAMSequenceRecord (htsjdk.samtools.SAMSequenceRecord)8 Cigar (htsjdk.samtools.Cigar)7 File (java.io.File)6 SAMFileHeader (htsjdk.samtools.SAMFileHeader)5 OutputStream (java.io.OutputStream)5 GATKRead (org.broadinstitute.hellbender.utils.read.GATKRead)5 VisibleForTesting (com.google.common.annotations.VisibleForTesting)4 Utils (org.broadinstitute.hellbender.utils.Utils)4 Tuple2 (scala.Tuple2)4 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)3 BufferedOutputStream (java.io.BufferedOutputStream)3 InputStream (java.io.InputStream)3 BigInteger (java.math.BigInteger)3 java.util (java.util)3