Java Code Examples for htsjdk.samtools.metrics.MetricsFile#addHistogram()

/**
 * Do the work after command line has been parsed. RuntimeException may be
 * thrown by this method, and are reported appropriately.
 *
 * @return program exit status.
 */
@Override
protected int doWork() {
	IOUtil.assertFileIsReadable(INPUT);
	IOUtil.assertFileIsWritable(OUTPUT);
	Map<String, ReadQualityMetrics> metricsMap = gatherMetrics(INPUT);
	MetricsFile<ReadQualityMetrics, Integer> outFile = new MetricsFile<>();
	outFile.addHistogram(metricsMap.get(this.GLOBAL).getHistogram());
       // Make sure the GLOBAL metrics is added first
       outFile.addMetric(metricsMap.remove(this.GLOBAL));
	for (ReadQualityMetrics metrics: metricsMap.values())
		outFile.addMetric(metrics);
	BufferedWriter w = IOUtil.openFileForBufferedWriting(OUTPUT);
	outFile.write(w);
	// close properly.
	try {
		w.close();
	} catch (IOException io) {
		throw new TranscriptomeException("Problem writing file", io);
	}
	return 0;
}
 

@Override
protected void finish() {
    // Generate a "Histogram" of mean quality and write it to the file
    final MetricsFile<?,Integer> metrics = getMetricsFile();
    metrics.addHistogram(q.getMeanQualityHistogram());
    if (!oq.isEmpty()) metrics.addHistogram(oq.getMeanQualityHistogram());
    metrics.write(OUTPUT);

    if (q.isEmpty() && oq.isEmpty()) {
        log.warn("No valid bases found in input file. No plot will be produced.");
    }
    else {
        // Now run R to generate a chart
        final int rResult = RExecutor.executeFromClasspath(
                "picard/analysis/meanQualityByCycle.R",
                OUTPUT.getAbsolutePath(),
                CHART_OUTPUT.getAbsolutePath(),
                INPUT.getName(),
                plotSubtitle);

        if (rResult != 0) {
            throw new PicardException("R script meanQualityByCycle.R failed with return code " + rResult);
        }
    }
}
 

private MetricsFile<?, Byte> makeMetrics(final Counts result) {
    // Built the Histograms out of the long[]s
    final Histogram<Byte> qHisto  = new Histogram<>("QUALITY", "COUNT_OF_Q");
    final Histogram<Byte> oqHisto = new Histogram<>("QUALITY", "COUNT_OF_OQ");

    for (int i=0; i< result.qCounts.length; ++i) {
        //Note: the checks are done to avoid work if there's nothing to add to the histogram
        if (result.qCounts[i]  > 0) qHisto.increment( (byte) i, (double) result.qCounts[i]);
        if (result.oqCounts[i] > 0) oqHisto.increment((byte) i, (double) result.oqCounts[i]);
    }

    final MetricsFile<?,Byte> metrics = getMetricsFile();
    metrics.addHistogram(qHisto);
    if (!oqHisto.isEmpty()) {
        metrics.addHistogram(oqHisto);
    }
    return metrics;
}
 

private void writeSummary (final Histogram<Integer> h) {

		MetricsFile<TrimMetric, Integer> mf = new MetricsFile<TrimMetric, Integer>();
		mf.addHistogram(h);
		TrimMetric tm=new TrimMetric(h);
		mf.addMetric(tm);
		mf.write(this.OUTPUT_SUMMARY);
	}
 

private void writeSummary(final Histogram<Integer> h) {

		MetricsFile<TrimMetric, Integer> mf = new MetricsFile<>();
		mf.addHistogram(h);
		TrimMetric tm = new TrimMetric(h);
		mf.addMetric(tm);
		mf.write(this.OUTPUT_SUMMARY);
	}
 

/**
 * Writes the metrics given by the libraryIdGenerator to the outputFile.
 *
 * @param libraryIdGenerator A {@link LibraryIdGenerator} object that contains the map from library to {@link DuplicationMetrics} for
 *                           that library
 * @param metricsFile        An empty {@link MetricsFile} object that will be filled, with "finalized" metrics and written out.
 *                           It needs to be generated from a non-static context so that various commandline information is
 *                           added to the header when {@link CommandLineProgram#getMetricsFile()} is called.
 * @param outputFile         The file to write the metrics to
 */
static public void finalizeAndWriteMetrics(final LibraryIdGenerator libraryIdGenerator, final MetricsFile<DuplicationMetrics, Double> metricsFile, final File outputFile) {
    final Map<String, DuplicationMetrics> metricsByLibrary = libraryIdGenerator.getMetricsByLibraryMap();
    final Histogram<Short> opticalDuplicatesByLibraryId = libraryIdGenerator.getOpticalDuplicatesByLibraryIdMap();
    final Map<String, Short> libraryIds = libraryIdGenerator.getLibraryIdsMap();

    // Write out the metrics
    for (final Map.Entry<String, DuplicationMetrics> entry : metricsByLibrary.entrySet()) {
        final String libraryName = entry.getKey();
        final DuplicationMetrics metrics = entry.getValue();

        metrics.READ_PAIRS_EXAMINED = metrics.READ_PAIRS_EXAMINED / 2;
        metrics.READ_PAIR_DUPLICATES = metrics.READ_PAIR_DUPLICATES / 2;

        // Add the optical dupes to the metrics
        final Short libraryId = libraryIds.get(libraryName);
        if (libraryId != null) {
            final Histogram.Bin<Short> bin = opticalDuplicatesByLibraryId.get(libraryId);
            if (bin != null) {
                metrics.READ_PAIR_OPTICAL_DUPLICATES = (long) bin.getValue();
            }
        }
        metrics.calculateDerivedFields();
        metricsFile.addMetric(metrics);
    }

    if (metricsByLibrary.size() == 1) {
        metricsFile.setHistogram(metricsByLibrary.values().iterator().next().calculateRoiHistogram());
    }

    // Add set size histograms - the set size counts are printed on adjacent columns to the ROI metric.
    metricsFile.addHistogram(libraryIdGenerator.getDuplicateCountHist());
    metricsFile.addHistogram(libraryIdGenerator.getOpticalDuplicateCountHist());
    metricsFile.addHistogram(libraryIdGenerator.getNonOpticalDuplicateCountHist());

    metricsFile.write(outputFile);
}
 

/**
 * Adds collected metrics and depth histogram to file
 * @param file MetricsFile for result of collector's work
 * @param dupeFilter         counting filter for duplicate reads
 * @param adapterFilter      counting filter for adapter reads
 * @param mapqFilter         counting filter for mapping quality
 * @param pairFilter         counting filter for reads without a mapped mate pair
 */
public void addToMetricsFile(final MetricsFile<WgsMetrics, Integer> file,
        final boolean includeBQHistogram,
        final CountingFilter dupeFilter,
        final CountingFilter adapterFilter,
        final CountingFilter mapqFilter,
        final CountingPairedFilter pairFilter) {
    final WgsMetrics metrics = getMetrics(dupeFilter, adapterFilter, mapqFilter, pairFilter);

    // add them to the file
    file.addMetric(metrics);
    file.addHistogram(getHighQualityDepthHistogram());
    if (includeBQHistogram) addBaseQHistogram(file);
}
 

@Override
public void addMetricsToFile(final MetricsFile<RnaSeqMetrics, Integer> file) {
    // Compute metrics based on coverage of top 1000 genes
    final Histogram<Integer> normalizedCovByPos = computeCoverageMetrics();
    file.addMetric(metrics);
    file.addHistogram(normalizedCovByPos);
}
 

@Override
public void addToMetricsFile(final MetricsFile<WgsMetrics, Integer> file,
                             final boolean includeBQHistogram,
                             final CountingFilter dupeFilter,
                             final CountingFilter adapterFilter,
                             final CountingFilter mapqFilter,
                             final CountingPairedFilter pairFilter) {
    highQualityDepthHistogram = getDepthHistogram();
    highQualityDepthHistogramNonZero = getDepthHistogramNonZero();

    // calculate metrics the same way as in CollectWgsMetrics
    final WgsMetricsWithNonZeroCoverage metrics = (WgsMetricsWithNonZeroCoverage) getMetrics(dupeFilter, adapterFilter, mapqFilter, pairFilter);
    metrics.CATEGORY = WgsMetricsWithNonZeroCoverage.Category.WHOLE_GENOME;

    // set count of the coverage-zero bin to 0 and re-calculate metrics
    // note we don't need to update the base quality histogram; there are no bases in the depth = 0 bin
    highQualityDepthHistogramArray[0] = 0;
    unfilteredDepthHistogramArray[0] = 0;

    final WgsMetricsWithNonZeroCoverage metricsNonZero = (WgsMetricsWithNonZeroCoverage) getMetrics(dupeFilter, adapterFilter, mapqFilter, pairFilter);
    metricsNonZero.CATEGORY = WgsMetricsWithNonZeroCoverage.Category.NON_ZERO_REGIONS;

    file.addMetric(metrics);
    file.addMetric(metricsNonZero);
    file.addHistogram(highQualityDepthHistogram);
    file.addHistogram(highQualityDepthHistogramNonZero);

    if (includeBQHistogram) {
        addBaseQHistogram(file);
    }
}
 

@Override
protected void finish() {
    // Built the Histograms out of the long[]s
    final Histogram<Byte> qHisto  = new Histogram<Byte>("QUALITY", "COUNT_OF_Q");
    final Histogram<Byte> oqHisto = new Histogram<Byte>("QUALITY", "COUNT_OF_OQ");

    for (int i=0; i< qCounts.length; ++i) {
        if (qCounts[i]  > 0) qHisto.increment( (byte) i, (double) qCounts[i]);
        if (oqCounts[i] > 0) oqHisto.increment((byte) i, (double) oqCounts[i]);
    }

    final MetricsFile<?,Byte> metrics = getMetricsFile();
    metrics.addHistogram(qHisto);
    if (!oqHisto.isEmpty()) metrics.addHistogram(oqHisto);
    metrics.write(OUTPUT);
    if (qHisto.isEmpty() && oqHisto.isEmpty()) {
        log.warn("No valid bases found in input file. No plot will be produced.");
    }
    else {
        // Now run R to generate a chart
        final int rResult = RExecutor.executeFromClasspath(
                "picard/analysis/qualityScoreDistribution.R",
                OUTPUT.getAbsolutePath(),
                CHART_OUTPUT.getAbsolutePath(),
                INPUT.getName(),
                this.plotSubtitle);

        if (rResult != 0) {
            throw new PicardException("R script qualityScoreDistribution.R failed with return code " + rResult);
        }
    }
}
 

private MetricsFile<?,Integer> finish(final HistogramGenerator q, final HistogramGenerator oq) {
    // Generate a "Histogram" of mean quality and write it to the file
    final MetricsFile<?,Integer> metrics = getMetricsFile();
    metrics.addHistogram(q.getMeanQualityHistogram());
    if (!oq.isEmpty()) {
        metrics.addHistogram(oq.getMeanQualityHistogram());
    }
    return metrics;
}
 

private void writeSummary(final BeadSynthesisErrorsSummaryMetric summary, final File out) {
	MetricsFile<BeadSynthesisErrorsSummaryMetric, Integer> outFile = new MetricsFile<>();
	outFile.addMetric(summary);
	outFile.addHistogram(summary.getHistogram());
	outFile.write(out);
}
 

protected void addBaseQHistogram(final MetricsFile<WgsMetrics, Integer> file) {
    file.addHistogram(getUnfilteredBaseQHistogram());
}
 

public void addMetricsToFile(final MetricsFile<InsertSizeMetrics,Integer> file) {
    // get the number of inserts, and the maximum and minimum keys across, across all orientations
    for (final Histogram<Integer> h : this.histograms.values()) {
        totalInserts += h.getCount();
    }
    if (0 == totalInserts) return; // nothing to store

    for(final Map.Entry<SamPairUtil.PairOrientation, Histogram<Integer>> entry : histograms.entrySet()) {
        final SamPairUtil.PairOrientation pairOrientation = entry.getKey();
        final Histogram<Integer> histogram = entry.getValue();
        final double total = histogram.getCount();

        // Only include a category if it has a sufficient percentage of the data in it
        if( total >= totalInserts * minimumPct ) {
            final InsertSizeMetrics metrics = new InsertSizeMetrics();
            metrics.SAMPLE             = this.sample;
            metrics.LIBRARY            = this.library;
            metrics.READ_GROUP         = this.readGroup;
            metrics.PAIR_ORIENTATION   = pairOrientation;
            if (!histogram.isEmpty()) {
                metrics.READ_PAIRS = (long) total;
                metrics.MAX_INSERT_SIZE = (int) histogram.getMax();
                metrics.MIN_INSERT_SIZE = (int) histogram.getMin();
                metrics.MEDIAN_INSERT_SIZE = histogram.getMedian();
                metrics.MODE_INSERT_SIZE = histogram.getMode();
                metrics.MEDIAN_ABSOLUTE_DEVIATION = histogram.getMedianAbsoluteDeviation();

                final double median = histogram.getMedian();
                double covered = 0;
                double low = median;
                double high = median;

                while (low >= histogram.getMin()-1 || high <= histogram.getMax()+1) {
                    final Histogram.Bin<Integer> lowBin = histogram.get((int) low);
                    if (lowBin != null) covered += lowBin.getValue();

                    if (low != high) {
                        final Histogram.Bin<Integer> highBin = histogram.get((int) high);
                        if (highBin != null) covered += highBin.getValue();
                    }

                    final double percentCovered = covered / total;
                    final int distance = (int) (high - low) + 1;
                    if (percentCovered >= 0.1 && metrics.WIDTH_OF_10_PERCENT == 0) metrics.WIDTH_OF_10_PERCENT = distance;
                    if (percentCovered >= 0.2 && metrics.WIDTH_OF_20_PERCENT == 0) metrics.WIDTH_OF_20_PERCENT = distance;
                    if (percentCovered >= 0.3 && metrics.WIDTH_OF_30_PERCENT == 0) metrics.WIDTH_OF_30_PERCENT = distance;
                    if (percentCovered >= 0.4 && metrics.WIDTH_OF_40_PERCENT == 0) metrics.WIDTH_OF_40_PERCENT = distance;
                    if (percentCovered >= 0.5 && metrics.WIDTH_OF_50_PERCENT == 0) metrics.WIDTH_OF_50_PERCENT = distance;
                    if (percentCovered >= 0.6 && metrics.WIDTH_OF_60_PERCENT == 0) metrics.WIDTH_OF_60_PERCENT = distance;
                    if (percentCovered >= 0.7 && metrics.WIDTH_OF_70_PERCENT == 0) metrics.WIDTH_OF_70_PERCENT = distance;
                    if (percentCovered >= 0.8 && metrics.WIDTH_OF_80_PERCENT == 0) metrics.WIDTH_OF_80_PERCENT = distance;
                    if (percentCovered >= 0.9 && metrics.WIDTH_OF_90_PERCENT == 0) metrics.WIDTH_OF_90_PERCENT = distance;
                    if (percentCovered >= 0.95 && metrics.WIDTH_OF_95_PERCENT == 0) metrics.WIDTH_OF_95_PERCENT = distance;
                    if (percentCovered >= 0.99 && metrics.WIDTH_OF_99_PERCENT == 0) metrics.WIDTH_OF_99_PERCENT = distance;

                    --low;
                    ++high;
                }
            }

            // Trim the Histogram down to get rid of outliers that would make the chart useless.
            final Histogram<Integer> trimmedHistogram = histogram; // alias it
            trimmedHistogram.trimByWidth(getWidthToTrimTo(metrics));

            if (!trimmedHistogram.isEmpty()) {
                metrics.MEAN_INSERT_SIZE = trimmedHistogram.getMean();
                metrics.STANDARD_DEVIATION = trimmedHistogram.getStandardDeviation();
            }

            file.addHistogram(trimmedHistogram);
            file.addMetric(metrics);
        }
    }
}
 

@Override
public void addMetricsToFile(final MetricsFile<METRIC_TYPE, Integer> hsMetricsComparableMetricsFile) {
    hsMetricsComparableMetricsFile.addMetric(convertMetric(this.metrics));
    hsMetricsComparableMetricsFile.addHistogram(highQualityDepthHistogram);
    hsMetricsComparableMetricsFile.addHistogram(unfilteredBaseQHistogram);
}