Examples with PairOrientation htsjdk.samtools.SamPairUtil.PairOrientation used on opensource projects

Example 1 with PairOrientation

use of htsjdk.samtools.SamPairUtil.PairOrientation in project gatk by broadinstitute.

the class CollectJumpingLibraryMetrics method doWork.

/**
     * Calculates the detailed statistics about the jumping library and then generates the results.
     */
@Override
protected Object doWork() {
    for (File f : INPUT) {
        IOUtil.assertFileIsReadable(f);
    }
    IOUtil.assertFileIsWritable(OUTPUT);
    Histogram<Integer> innieHistogram = new Histogram<>();
    Histogram<Integer> outieHistogram = new Histogram<>();
    int fragments = 0;
    int innies = 0;
    int outies = 0;
    int innieDupes = 0;
    int outieDupes = 0;
    int crossChromPairs = 0;
    int superSized = 0;
    int tandemPairs = 0;
    double chimeraSizeMinimum = Math.max(getOutieMode(), (double) CHIMERA_KB_MIN);
    for (File f : INPUT) {
        SamReader reader = SamReaderFactory.makeDefault().open(f);
        if (reader.getFileHeader().getSortOrder() != SAMFileHeader.SortOrder.coordinate) {
            throw new UserException("SAM file must " + f.getName() + " must be sorted in coordinate order");
        }
        for (SAMRecord sam : reader) {
            // We're getting all our info from the first of each pair.
            if (!sam.getFirstOfPairFlag()) {
                continue;
            }
            // Ignore unmapped read pairs
            if (sam.getReadUnmappedFlag()) {
                if (!sam.getMateUnmappedFlag()) {
                    fragments++;
                    continue;
                }
                // If both ends are unmapped and we've hit unaligned reads we're done
                if (sam.getReferenceIndex() == SAMRecord.NO_ALIGNMENT_REFERENCE_INDEX) {
                    break;
                }
                continue;
            }
            if (sam.getMateUnmappedFlag()) {
                fragments++;
                continue;
            }
            // Ignore low-quality reads.  If we don't have the mate mapping quality, assume it's OK
            if ((sam.getAttribute(SAMTag.MQ.name()) != null && sam.getIntegerAttribute(SAMTag.MQ.name()) < MINIMUM_MAPPING_QUALITY) || sam.getMappingQuality() < MINIMUM_MAPPING_QUALITY) {
                continue;
            }
            final int absInsertSize = Math.abs(sam.getInferredInsertSize());
            if (absInsertSize > chimeraSizeMinimum) {
                superSized++;
            } else if (sam.getMateNegativeStrandFlag() == sam.getReadNegativeStrandFlag()) {
                tandemPairs++;
            } else if (!sam.getMateReferenceIndex().equals(sam.getReferenceIndex())) {
                crossChromPairs++;
            } else {
                final PairOrientation pairOrientation = SamPairUtil.getPairOrientation(sam);
                if (pairOrientation == PairOrientation.RF) {
                    outieHistogram.increment(absInsertSize);
                    outies++;
                    if (sam.getDuplicateReadFlag()) {
                        outieDupes++;
                    }
                } else if (pairOrientation == PairOrientation.FR) {
                    innieHistogram.increment(absInsertSize);
                    innies++;
                    if (sam.getDuplicateReadFlag()) {
                        innieDupes++;
                    }
                } else {
                    throw new IllegalStateException("This should never happen");
                }
            }
        }
        CloserUtil.close(reader);
    }
    MetricsFile<JumpingLibraryMetrics, Integer> metricsFile = getMetricsFile();
    JumpingLibraryMetrics metrics = new JumpingLibraryMetrics();
    metrics.JUMP_PAIRS = outies;
    metrics.JUMP_DUPLICATE_PAIRS = outieDupes;
    metrics.JUMP_DUPLICATE_PCT = outies != 0 ? outieDupes / (double) outies : 0;
    metrics.JUMP_LIBRARY_SIZE = (outies > 0 && outieDupes > 0) ? DuplicationMetrics.estimateLibrarySize(outies, outies - outieDupes) : 0;
    outieHistogram.trimByTailLimit(TAIL_LIMIT);
    metrics.JUMP_MEAN_INSERT_SIZE = outieHistogram.getMean();
    metrics.JUMP_STDEV_INSERT_SIZE = outieHistogram.getStandardDeviation();
    metrics.NONJUMP_PAIRS = innies;
    metrics.NONJUMP_DUPLICATE_PAIRS = innieDupes;
    metrics.NONJUMP_DUPLICATE_PCT = innies != 0 ? innieDupes / (double) innies : 0;
    metrics.NONJUMP_LIBRARY_SIZE = (innies > 0 && innieDupes > 0) ? DuplicationMetrics.estimateLibrarySize(innies, innies - innieDupes) : 0;
    innieHistogram.trimByTailLimit(TAIL_LIMIT);
    metrics.NONJUMP_MEAN_INSERT_SIZE = innieHistogram.getMean();
    metrics.NONJUMP_STDEV_INSERT_SIZE = innieHistogram.getStandardDeviation();
    metrics.CHIMERIC_PAIRS = crossChromPairs + superSized + tandemPairs;
    metrics.FRAGMENTS = fragments;
    double totalPairs = outies + innies + metrics.CHIMERIC_PAIRS;
    metrics.PCT_JUMPS = totalPairs != 0 ? outies / totalPairs : 0;
    metrics.PCT_NONJUMPS = totalPairs != 0 ? innies / totalPairs : 0;
    metrics.PCT_CHIMERAS = totalPairs != 0 ? metrics.CHIMERIC_PAIRS / totalPairs : 0;
    metricsFile.addMetric(metrics);
    metricsFile.write(OUTPUT);
    return null;
}