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

The following examples show how to use htsjdk.samtools.metrics.MetricsFile#addHistogram() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: GatherReadQualityMetrics.java    From Drop-seq with MIT License 6 votes vote down vote up
/**
 * 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;
}
 
Example 2
Source File: MeanQualityByCycle.java    From picard with MIT License 6 votes vote down vote up
@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);
        }
    }
}
 
Example 3
Source File: QualityScoreDistributionSpark.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
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;
}
 
Example 4
Source File: TrimStartingSequence.java    From Drop-seq with MIT License 5 votes vote down vote up
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);
	}
 
Example 5
Source File: PolyATrimmer.java    From Drop-seq with MIT License 5 votes vote down vote up
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);
	}
 
Example 6
Source File: AbstractMarkDuplicatesCommandLineProgram.java    From picard with MIT License 5 votes vote down vote up
/**
 * 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);
}
 
Example 7
Source File: AbstractWgsMetricsCollector.java    From picard with MIT License 5 votes vote down vote up
/**
 * 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);
}
 
Example 8
Source File: RnaSeqMetricsCollector.java    From picard with MIT License 5 votes vote down vote up
@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);
}
 
Example 9
Source File: CollectWgsMetricsWithNonZeroCoverage.java    From picard with MIT License 5 votes vote down vote up
@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);
    }
}
 
Example 10
Source File: QualityScoreDistribution.java    From picard with MIT License 5 votes vote down vote up
@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);
        }
    }
}
 
Example 11
Source File: MeanQualityByCycleSpark.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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;
}
 
Example 12
Source File: DetectBeadSynthesisErrors.java    From Drop-seq with MIT License 4 votes vote down vote up
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);
}
 
Example 13
Source File: AbstractWgsMetricsCollector.java    From picard with MIT License 4 votes vote down vote up
protected void addBaseQHistogram(final MetricsFile<WgsMetrics, Integer> file) {
    file.addHistogram(getUnfilteredBaseQHistogram());
}
 
Example 14
Source File: InsertSizeMetricsCollector.java    From picard with MIT License 4 votes vote down vote up
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);
        }
    }
}
 
Example 15
Source File: TargetMetricsCollector.java    From picard with MIT License 4 votes vote down vote up
@Override
public void addMetricsToFile(final MetricsFile<METRIC_TYPE, Integer> hsMetricsComparableMetricsFile) {
    hsMetricsComparableMetricsFile.addMetric(convertMetric(this.metrics));
    hsMetricsComparableMetricsFile.addHistogram(highQualityDepthHistogram);
    hsMetricsComparableMetricsFile.addHistogram(unfilteredBaseQHistogram);
}