Java Code Examples for htsjdk.samtools.metrics.MetricsFile

The following examples show how to use htsjdk.samtools.metrics.MetricsFile. These examples are extracted from open source projects. 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 Project: Drop-seq   Source File: GatherReadQualityMetrics.java    License: 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 Project: picard   Source File: ExtractIlluminaBarcodesTest.java    License: MIT License 6 votes vote down vote up
/**
 * Testing the quality thresholding. Looking at a single barcode (ACAGTG) with a min quality of 25 and no mismatches
 */
@Test(dataProvider = "qualityBarcodeData")
public void testQualityBarcodes(final int quality,
                                final int maxMismatches, final int perfectMatches, final int oneMismatch,
                                final String testName) throws Exception {
    final File metricsFile = File.createTempFile("qual.", ".metrics");
    metricsFile.deleteOnExit();

    final String[] args = new String[]{
            "BASECALLS_DIR=" + qual.getPath(),
            "LANE=" + 1,
            "READ_STRUCTURE=25T8B25T",
            "METRICS_FILE=" + metricsFile.getPath(),
            "MINIMUM_BASE_QUALITY=" + quality,
            "MAX_MISMATCHES=" + maxMismatches,
            "BARCODE=CAATAGTC"
    };

    Assert.assertEquals(runPicardCommandLine(args), 0);
    final MetricsFile<ExtractIlluminaBarcodes.BarcodeMetric, Integer> result = new MetricsFile<>();
    result.read(new FileReader(metricsFile));
    Assert.assertEquals(result.getMetrics().get(0).PERFECT_MATCHES, perfectMatches, "Got wrong number of perfect matches for test: '" + testName + "'");
    Assert.assertEquals(result.getMetrics().get(0).ONE_MISMATCH_MATCHES, oneMismatch, "Got wrong number of one-mismatch matches for test: '" + testName + "'");
}
 
Example 3
Source Project: picard   Source File: TheoreticalSensitivityTest.java    License: MIT License 6 votes vote down vote up
@Test(dataProvider = "TheoreticalSensitivityConstantDepthDataProvider")
public void testSensitivityAtConstantDepth(final double expected, final File metricsFile, final double alleleFraction, final int depth, final int sampleSize, final double tolerance) throws Exception {
    // This tests Theoretical Sensitivity assuming a uniform depth with a distribution of base quality scores.
    // Because this only tests sensitivity at a constant depth, we use this for testing the code at high depths.
    final MetricsFile<?, Integer> metrics = new MetricsFile<>();
    try (final FileReader metricsFileReader = new FileReader(metricsFile)) {
        metrics.read(metricsFileReader);
    }

    final List<Histogram<Integer>> histograms = metrics.getAllHistograms();
    final Histogram<Integer> qualityHistogram = histograms.get(1);

    // We ensure that even using different random seeds we converge to roughly the same value.
    for (int i = 0; i < 3; i++) {
        double result = TheoreticalSensitivity.sensitivityAtConstantDepth(depth, qualityHistogram, 3, sampleSize, alleleFraction, i);
        Assert.assertEquals(result, expected, tolerance);
    }
}
 
Example 4
Source Project: picard   Source File: CollectIlluminaBasecallingMetricsTest.java    License: MIT License 6 votes vote down vote up
@Test
public void testNovaseqIndexedRun() throws Exception {
    final MetricsFile<IlluminaBasecallingMetrics, Integer> metricsFile = runIt(1, "151T8B8B151T",
            new File("testdata/picard/illumina/151T8B8B151T_cbcl/Data/Intensities/BaseCalls"), null, true);
    final IlluminaBasecallingMetrics laneMetric = metricsFile.getMetrics().get(0);
    Assert.assertEquals(laneMetric.LANE, "1");
    Assert.assertEquals(laneMetric.MOLECULAR_BARCODE_SEQUENCE_1, "CACCTAGTACTCGAGT");
    Assert.assertEquals(laneMetric.MOLECULAR_BARCODE_NAME, "SA_CACCTAGTACTCGAGT");
    Assert.assertEquals(laneMetric.MEAN_CLUSTERS_PER_TILE, 1.0);
    Assert.assertEquals(laneMetric.SD_CLUSTERS_PER_TILE, 0.0);
    Assert.assertEquals(laneMetric.MEAN_PF_CLUSTERS_PER_TILE,1.0);
    Assert.assertEquals(laneMetric.SD_PF_CLUSTERS_PER_TILE, 0.0);
    Assert.assertEquals(laneMetric.MEAN_PCT_PF_CLUSTERS_PER_TILE, 100.0);
    Assert.assertEquals(laneMetric.SD_PCT_PF_CLUSTERS_PER_TILE, 0.0);
    Assert.assertEquals(laneMetric.TOTAL_BASES, 302);
    Assert.assertEquals(laneMetric.TOTAL_READS, 2);
    Assert.assertEquals(laneMetric.PF_BASES, 302);
    Assert.assertEquals(laneMetric.PF_READS, 2);


    Assert.assertEquals(metricsFile.getMetrics().size(),6);
}
 
Example 5
Source Project: picard   Source File: GenotypeConcordance.java    License: MIT License 6 votes vote down vote up
/**
* Outputs the detailed statistics tables for SNP and Indel match categories.
**/
public static void outputDetailMetricsFile(final VariantContext.Type variantType, final MetricsFile<GenotypeConcordanceDetailMetrics,?> genotypeConcordanceDetailMetricsFile,
                                           final GenotypeConcordanceCounts counter, final String truthSampleName, final String callSampleName,
                                           final boolean missingSitesHomRef, final boolean outputAllRows) {
    final GenotypeConcordanceSchemeFactory schemeFactory = new GenotypeConcordanceSchemeFactory();
    final GenotypeConcordanceScheme scheme = schemeFactory.getScheme(missingSitesHomRef);
    scheme.validateScheme();
    for (final TruthState truthState : TruthState.values()) {
        for (final CallState callState : CallState.values()) {
            final long count = counter.getCount(truthState, callState);
            final String contingencyValues = scheme.getContingencyStateString(truthState, callState);
            if (count > 0 || outputAllRows) {
                final GenotypeConcordanceDetailMetrics detailMetrics = new GenotypeConcordanceDetailMetrics();
                detailMetrics.VARIANT_TYPE = variantType;
                detailMetrics.TRUTH_SAMPLE = truthSampleName;
                detailMetrics.CALL_SAMPLE = callSampleName;
                detailMetrics.TRUTH_STATE = truthState;
                detailMetrics.CALL_STATE = callState;
                detailMetrics.COUNT = count;
                detailMetrics.CONTINGENCY_VALUES = contingencyValues;
                genotypeConcordanceDetailMetricsFile.addMetric(detailMetrics);
            }
        }
    }
}
 
Example 6
Source Project: picard   Source File: ClusterCrosscheckMetrics.java    License: MIT License 6 votes vote down vote up
@Override
protected int doWork() {
    IOUtil.assertFileIsReadable(INPUT);
    if(OUTPUT != null) IOUtil.assertFileIsWritable(OUTPUT);

    final MetricsFile<CrosscheckMetric, ?> metricsFile = getMetricsFile();

    try {
        metricsFile.read(new FileReader(INPUT));
    } catch (FileNotFoundException e) {
        e.printStackTrace();
        return 1;
    }

    clusterMetrics(metricsFile.getMetrics()).write(OUTPUT);

    return 0;
}
 
Example 7
Source Project: picard   Source File: CollectSamErrorMetricsTest.java    License: MIT License 6 votes vote down vote up
@Test(dataProvider = "OneCovariateErrorMetricsDataProvider")
public void testOneCovariateErrorMetrics(final String errorSubscript, final File samFile, final int priorQ, BaseErrorMetric expectedMetric) {

    ErrorMetric.setPriorError(QualityUtil.getErrorProbabilityFromPhredScore(priorQ));
    expectedMetric.calculateDerivedFields();

    // Note that soft clipped bases are not counted
    List<BaseErrorMetric> metrics = MetricsFile.readBeans(new File(errorMetrics.get(samFile).getAbsolutePath() + errorSubscript));

    BaseErrorMetric metric = metrics
            .stream()
            .filter(m -> m.COVARIATE.equals(expectedMetric.COVARIATE))
            .findAny()
            .orElseThrow(() -> new AssertionError("didn't find metric with COVARIATE==" + expectedMetric.COVARIATE));

    Assert.assertEquals(metric, expectedMetric);
}
 
Example 8
Source Project: picard   Source File: TheoreticalSensitivityTest.java    License: MIT License 6 votes vote down vote up
@Test(dataProvider = "TheoreticalSensitivityDataProvider")
public void testSensitivity(final double expected, final File metricsFile, final double alleleFraction, final int sampleSize) throws Exception {
    // This tests Theoretical Sensitivity using distributions on both base quality scores
    // and the depth histogram.

    // We use a pretty forgiving tolerance here because for these tests
    // we are not using large enough sample sizes to converge.
    final double tolerance = 0.02;

    final MetricsFile<?, Integer> metrics = new MetricsFile<>();
    try (final FileReader metricsFileReader = new FileReader(metricsFile)) {
        metrics.read(metricsFileReader);
    }

    final List<Histogram<Integer>> histograms = metrics.getAllHistograms();
    final Histogram<Integer> depthHistogram = histograms.get(0);
    final Histogram<Integer> qualityHistogram = histograms.get(1);

    final double result = TheoreticalSensitivity.theoreticalSensitivity(depthHistogram, qualityHistogram, sampleSize, 3, alleleFraction);

    Assert.assertEquals(result, expected, tolerance);
}
 
Example 9
Source Project: picard   Source File: CollectIlluminaBasecallingMetricsTest.java    License: MIT License 6 votes vote down vote up
@Test
public void testNonIndexedRunLane1() throws Exception {
    final MetricsFile<IlluminaBasecallingMetrics, Integer> metricsFile = runIt(1, "125T125T",
            new File(rootTestDir, "125T125T/Data/Intensities/BaseCalls"), null, false);
    final IlluminaBasecallingMetrics laneMetric = metricsFile.getMetrics().get(0);

    Assert.assertEquals(laneMetric.LANE, "1");
    Assert.assertEquals(laneMetric.MOLECULAR_BARCODE_SEQUENCE_1, null);
    Assert.assertEquals(laneMetric.MOLECULAR_BARCODE_NAME, null);
    Assert.assertEquals(laneMetric.MEAN_CLUSTERS_PER_TILE, 2000.0);
    Assert.assertEquals(laneMetric.SD_CLUSTERS_PER_TILE, 0.0);
    Assert.assertEquals(laneMetric.MEAN_PF_CLUSTERS_PER_TILE,1863.0);
    Assert.assertEquals(laneMetric.SD_PF_CLUSTERS_PER_TILE, 0.0);
    Assert.assertEquals(laneMetric.MEAN_PCT_PF_CLUSTERS_PER_TILE, 93.15);
    Assert.assertEquals(laneMetric.SD_PCT_PF_CLUSTERS_PER_TILE, 0.0);
    Assert.assertEquals(laneMetric.TOTAL_BASES, 500000);
    Assert.assertEquals(laneMetric.TOTAL_READS, 4000);
    Assert.assertEquals(laneMetric.PF_BASES, 465750);
    Assert.assertEquals(laneMetric.PF_READS, 3726);


    Assert.assertEquals(metricsFile.getMetrics().size(),1);
}
 
Example 10
Source Project: picard   Source File: CollectBaseDistributionByCycle.java    License: MIT License 6 votes vote down vote up
@Override
protected void finish() {
    final MetricsFile<BaseDistributionByCycleMetrics, ?> metrics = getMetricsFile();
    hist.addToMetricsFile(metrics);
    metrics.write(OUTPUT);
    if (hist.isEmpty()) {
        log.warn("No valid bases found in input file. No plot will be produced.");
    } else {
        final int rResult = RExecutor.executeFromClasspath("picard/analysis/baseDistributionByCycle.R",
                OUTPUT.getAbsolutePath(),
                CHART_OUTPUT.getAbsolutePath(),
                INPUT.getName(),
                plotSubtitle);
        if (rResult != 0) {
            throw new PicardException("R script nucleotideDistributionByCycle.R failed with return code " + rResult);
        }
    }
}
 
Example 11
Source Project: picard   Source File: CollectAlignmentSummaryMetricsTest.java    License: MIT License 6 votes vote down vote up
@Test
public void testZeroLengthReads() throws IOException {
    final File input = new File(TEST_DATA_DIR, "summary_alignment_stats_test2.sam");
    final File outfile = File.createTempFile("alignmentMetrics", ".txt");
    outfile.deleteOnExit();
    final String[] args = new String[]{
            "INPUT=" + input.getAbsolutePath(),
            "OUTPUT=" + outfile.getAbsolutePath(),
            "COLLECT_ALIGNMENT_INFORMATION=false"
    };
    Assert.assertEquals(runPicardCommandLine(args), 0);

    final MetricsFile<AlignmentSummaryMetrics, Comparable<?>> output = new MetricsFile<>();
    try (FileReader reader = new FileReader(outfile)) {
        output.read(reader);
    }
    for (final AlignmentSummaryMetrics metrics : output.getMetrics()) {
        // test that it doesn't blow up
    }
}
 
Example 12
Source Project: picard   Source File: AlignmentSummaryMetricsCollector.java    License: MIT License 6 votes vote down vote up
@Override
public void addMetricsToFile(final MetricsFile<AlignmentSummaryMetrics, Comparable<?>> file) {
    if (firstOfPairCollector.getMetrics().TOTAL_READS > 0) {
        // override how bad cycle is determined for paired reads, it should be
        // the sum of first and second reads
        pairCollector.getMetrics().BAD_CYCLES = firstOfPairCollector.getMetrics().BAD_CYCLES +
                secondOfPairCollector.getMetrics().BAD_CYCLES;

        file.addMetric(firstOfPairCollector.getMetrics());
        file.addMetric(secondOfPairCollector.getMetrics());
        file.addMetric(pairCollector.getMetrics());
    }

    //if there are no reads in any category then we will returned an unpaired alignment summary metric with all zero values
    if (unpairedCollector.getMetrics().TOTAL_READS > 0 || firstOfPairCollector.getMetrics().TOTAL_READS == 0) {
        file.addMetric(unpairedCollector.getMetrics());
    }
}
 
Example 13
Source Project: picard   Source File: CollectTargetedMetricsTest.java    License: MIT License 6 votes vote down vote up
@Test(dataProvider = "targetedIntervalDataProvider")
public void runCollectTargetedMetricsTest(final File input, final File outfile, final File perTargetOutfile, final String referenceFile,
                            final String targetIntervals, final int sampleSize) throws IOException {

    final String[] args = new String[] {
            "TARGET_INTERVALS=" + targetIntervals,
            "INPUT=" + input.getAbsolutePath(),
            "OUTPUT=" + outfile.getAbsolutePath(),
            "REFERENCE_SEQUENCE=" + referenceFile,
            "PER_TARGET_COVERAGE=" + perTargetOutfile.getAbsolutePath(),
            "LEVEL=ALL_READS",
            "AMPLICON_INTERVALS=" + targetIntervals,
            "SAMPLE_SIZE=" + sampleSize
    };

    Assert.assertEquals(runPicardCommandLine(args), 0);

    final MetricsFile<TargetedPcrMetrics, Comparable<?>> output = new MetricsFile<>();
    output.read(new FileReader(outfile));

    for (final TargetedPcrMetrics metrics : output.getMetrics()) {
        Assert.assertEquals(metrics.TOTAL_READS, numReads * 2);
        Assert.assertEquals(metrics.HET_SNP_SENSITIVITY, .997972, .02);
    }
}
 
Example 14
Source Project: picard   Source File: TheoreticalSensitivityTest.java    License: MIT License 6 votes vote down vote up
@Test(dataProvider = "sumOfGaussiansDataProvider")
public void testDrawSumOfQScores(final File metricsFile, final int altDepth, final double tolerance) throws Exception {
    final MetricsFile<TheoreticalSensitivityMetrics, Integer> metrics = new MetricsFile<>();
    try (final FileReader metricsFileReader = new FileReader(metricsFile)) {
        metrics.read(metricsFileReader);
    }

    final List<Histogram<Integer>> histograms = metrics.getAllHistograms();

    final Histogram<Integer> qualityHistogram = histograms.get(1);
    final TheoreticalSensitivity.RouletteWheel qualityRW = new TheoreticalSensitivity.RouletteWheel(TheoreticalSensitivity.trimDistribution(TheoreticalSensitivity.normalizeHistogram(qualityHistogram)));

    final Random randomNumberGenerator = new Random(51);

    // Calculate mean and deviation of quality score distribution to enable Gaussian sampling below
    final double averageQuality = qualityHistogram.getMean();
    final double standardDeviationQuality = qualityHistogram.getStandardDeviation();

    for (int k = 0; k < 1; k++) {
        int sumOfQualitiesFull = IntStream.range(0, altDepth).map(n -> qualityRW.draw()).sum();
        int sumOfQualities = TheoreticalSensitivity.drawSumOfQScores(altDepth, averageQuality, standardDeviationQuality, randomNumberGenerator.nextGaussian());

        Assert.assertEquals(sumOfQualitiesFull, sumOfQualities, sumOfQualitiesFull * tolerance);
    }
}
 
Example 15
Source Project: picard   Source File: CollectSamErrorMetricsTest.java    License: MIT License 6 votes vote down vote up
@Test(dataProvider = "oneCovariateIndelErrorMetricsDataProvider")
public void testOneCovariateIndelErrorMetrics(final String errorSubscript, final File samFile, final int priorQ, BaseErrorMetric expectedMetric) {

    final File errorByAll = new File(errorMetrics.get(samFile).getAbsolutePath() + errorSubscript);
    errorByAll.deleteOnExit();

    ErrorMetric.setPriorError(QualityUtil.getErrorProbabilityFromPhredScore(priorQ));
    expectedMetric.calculateDerivedFields();

    // Note that soft clipped bases are not counted
    List<BaseErrorMetric> metrics = MetricsFile.readBeans(errorByAll);

    BaseErrorMetric metric = metrics
            .stream()
            .filter(m -> m.COVARIATE.equals(expectedMetric.COVARIATE))
            .findAny()
            .orElseThrow(() -> new AssertionError("didn't find metric with COVARIATE==" + expectedMetric.COVARIATE));

    Assert.assertEquals(metric, expectedMetric);
}
 
Example 16
protected void saveResults(final MetricsFile<?, Integer> metrics, final SAMFileHeader readsHeader, final String inputFileName) {
    MetricsUtils.saveMetrics(metrics, out);

    if (metrics.getAllHistograms().isEmpty()) {
        logger.warn("No valid bases found in input file.");
    } else if (chartOutput != null) {
        // Now run R to generate a chart

        // If we're working with a single library, assign that library's name
        // as a suffix to the plot title
        final List<SAMReadGroupRecord> readGroups = readsHeader.getReadGroups();

        /*
         * A subtitle for the plot, usually corresponding to a library.
         */
        String plotSubtitle = "";
        if (readGroups.size() == 1) {
            plotSubtitle = StringUtil.asEmptyIfNull(readGroups.get(0).getLibrary());
        }
        final RScriptExecutor executor = new RScriptExecutor();
        executor.addScript(getBaseDistributionByCycleRScriptResource());
        executor.addArgs(out, chartOutput.getAbsolutePath(), inputFileName, plotSubtitle);
        executor.exec();
    }
}
 
Example 17
Source Project: Drop-seq   Source File: SpermSeqMarkDuplicates.java    License: MIT License 5 votes vote down vote up
@Override
protected int doWork() {
	for (File i: INPUT)
		IOUtil.assertFileIsReadable(i);
	IOUtil.assertFileIsWritable(OUTPUT);
	IOUtil.assertFileIsWritable(OUTPUT_STATS);

	// for reporting, what is the aggregate output name
	globalMetrics.CELL_BARCODE=AGGREGATE_NAME;

	List<String> cellBarcodes = getCellBarcodes();

       // no need to maintain input sort because records must be re-sorted to do dupe marking
	SamHeaderAndIterator headerAndIterator = SamFileMergeUtil.mergeInputs(INPUT, false, samReaderFactory);
	SAMFileHeader h= headerAndIterator.header;
       h.setSortOrder(SORT_ORDER);
	SAMFileWriter writer= new SAMFileWriterFactory().makeSAMOrBAMWriter(h, false, OUTPUT);

	Map<String, PCRDuplicateMetrics> metricsMap=null;

	metricsMap=processReadsByPosition(cellBarcodes, headerAndIterator, writer);

	log.info("Processing done, writing final BAM - this requires resorting the output BAM into genomic order");
	writer.close();

	MetricsFile<PCRDuplicateMetrics, String> file = new MetricsFile<>();
	this.globalMetrics.calculateStats();
	file.addMetric(this.globalMetrics);
	for (String key: cellBarcodes) {
		PCRDuplicateMetrics cell = metricsMap.get(key);
		if (cell!=null) {
			cell.calculateStats();
			file.addMetric(cell);
		}

	}
	file.write(this.OUTPUT_STATS);
	return 0;
}
 
Example 18
Source Project: Drop-seq   Source File: TagReadWithGeneFunction.java    License: MIT License 5 votes vote down vote up
@Override
protected int doWork() {
	IOUtil.assertFileIsReadable(this.INPUT);
	IOUtil.assertFileIsReadable(this.ANNOTATIONS_FILE);
	if (this.SUMMARY!=null) IOUtil.assertFileIsWritable(this.SUMMARY);
	IOUtil.assertFileIsWritable(this.OUTPUT);

	SamReader inputSam = SamReaderFactory.makeDefault().open(INPUT);

	SAMFileHeader header = inputSam.getFileHeader();
	SamHeaderUtil.addPgRecord(header, this);
	SAMSequenceDictionary bamDict = header.getSequenceDictionary();

       final OverlapDetector<Gene> geneOverlapDetector = GeneAnnotationReader.loadAnnotationsFile(ANNOTATIONS_FILE, bamDict);
       SAMFileWriter writer= new SAMFileWriterFactory().makeSAMOrBAMWriter(header, true, OUTPUT);

       for (SAMRecord r: inputSam) {
       	pl.record(r);

       	if (!r.getReadUnmappedFlag())
			// r=	setGeneExons(r, geneOverlapDetector, this.ALLOW_MULTI_GENE_READS);
       		r= setAnnotations(r, geneOverlapDetector, this.ALLOW_MULTI_GENE_READS);
       	writer.addAlignment(r);
       }

	CloserUtil.close(inputSam);
	writer.close();
	// if (this.USE_STRAND_INFO) log.info(this.metrics.toString());
	if (SUMMARY==null) return 0;

	//process summary
	MetricsFile<ReadTaggingMetric, Integer> outFile = new MetricsFile<ReadTaggingMetric, Integer>();
	outFile.addMetric(this.metrics);
	outFile.write(this.SUMMARY);
	return 0;

}
 
Example 19
Source Project: Drop-seq   Source File: DigitalExpression.java    License: MIT License 5 votes vote down vote up
public static void writeSummary(final Collection<DESummary> summaryCollection, final MetricsFile<DESummary, Integer> summaryMetricsFile, final File outFile) {
    List<DESummary> sc = new ArrayList<>(summaryCollection);
    Collections.sort(sc, DigitalExpression.TRANSCRIPT_ORDER_DESCENDING);
    for (DESummary z: sc)
        summaryMetricsFile.addMetric(z);
    summaryMetricsFile.write(outFile);
}
 
Example 20
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 21
Source Project: Drop-seq   Source File: FilterBam.java    License: MIT License 5 votes vote down vote up
/**
 * Write the summary output file of reads accepted/rejected.
 * @param summaryOutputFile
 * @param metrics
 */
private void writeSummary (File summaryOutputFile, FilteredReadsMetric metrics) {
	if (summaryOutputFile!=null) {
		MetricsFile<FilteredReadsMetric, Integer> outSummary = getMetricsFile();
		outSummary.addMetric(metrics);
		outSummary.write(summaryOutputFile);		
	}
}
 
Example 22
Source Project: picard   Source File: GenotypeConcordanceTest.java    License: MIT License 5 votes vote down vote up
private void assertMetricsFileEqual(final File actualMetricsFile, final File expectedMetricsFile) throws FileNotFoundException {
    // Actual metrics file
    final MetricsFile<GenotypeConcordanceSummaryMetrics, Comparable<?>> actual = new MetricsFile<GenotypeConcordanceSummaryMetrics, Comparable<?>>();
    actual.read(new FileReader(actualMetricsFile));

    // Expected metrics file
    final MetricsFile<GenotypeConcordanceSummaryMetrics, Comparable<?>> expected = new MetricsFile<GenotypeConcordanceSummaryMetrics, Comparable<?>>();
    expected.read(new FileReader(expectedMetricsFile));

    // Note - cannot use .equals as it calls .areHeadersEqual and they are not since the timestamp (at a minimum is different)
    Assert.assertTrue(expected.areMetricsEqual(actual));
    Assert.assertTrue(expected.areHistogramsEqual(actual));
}
 
Example 23
Source Project: Drop-seq   Source File: TestUtils.java    License: MIT License 5 votes vote down vote up
public static boolean testMetricsFilesEqual(final File expected, final File actual) throws FileNotFoundException {
	final MetricsFile<?, ?> metricsA = new MetricsFile();
	final MetricsFile<?, ?> metricsB = new MetricsFile();
	metricsA.read(new FileReader(expected));
	metricsB.read(new FileReader(actual));
	return metricsA.areMetricsEqual(metricsB) && metricsA.areHistogramsEqual(metricsB);
}
 
Example 24
@Override
public void onTraversalStart() {

    // Gets around issue 2274 in gatk public
    if (transitions.size() == 0) {
        transitions.add(DEFAULT_ARTIFACT_MODE);
    }

    // Sort the input artifacts argument
    transitions.sort(null);

    final MetricsFile<SequencingArtifactMetrics.PreAdapterDetailMetrics, Comparable<?>> mf = new MetricsFile<>();

    try {
        mf.read(new FileReader(preAdapterMetricsFile));
    } catch (final FileNotFoundException fnfe) {
        throw new UserException("Could not find file: " + preAdapterMetricsFile.getAbsolutePath());
    }

    // Collect all of the transitions that were specified in the parameters.
    relevantTransitions.addAll(transitions.stream().map(s -> convertParameterToTransition(s)).collect(Collectors.toSet()));

    // Get the PreAdapterQ score from the picard metrics tool, which gives an indication of how badly infested the bam file is.
    transitionToPreAdapterScoreMap = PreAdapterOrientationScorer.scoreOrientationBiasMetricsOverContext(mf.getMetrics());
    logger.info("preAdapter scores:");
    transitionToPreAdapterScoreMap.keySet().stream().forEach(k -> logger.info(k + ": " + transitionToPreAdapterScoreMap.get(k)));

    setupVCFWriter();
}
 
Example 25
/**
 * Initialize the collector with input arguments;
 */
@Override
public void initialize(
        final QualityYieldMetricsArgumentCollection inputArgs,
        final SAMFileHeader samHeader,
        final List<Header> defaultHeaders) {
    metricsFile = new MetricsFile<QualityYieldMetrics, Integer>();
    defaultHeaders.stream().forEach(h -> metricsFile.addHeader(h));
    this.args = inputArgs;
}
 
Example 26
public static List<Histogram<Integer>> sumHistogramsFromFiles(final List<MetricsFile<?,Integer>> metricsFiles, final boolean ref){
    Utils.nonNull(metricsFiles, "files may not be null");
    if (metricsFiles.isEmpty()) {
        return Collections.emptyList();
    }

    final List<Histogram<Integer>> histogramList = metricsFiles.get(0).getAllHistograms();
    if (ref){
        Utils.validate(histogramList.size() == F1R2FilterConstants.NUM_KMERS,
                "The list of ref histograms need to include all kmers as enforced by CollectF1R2Counts");
        Utils.validate(histogramList.stream().allMatch(h -> F1R2FilterConstants.ALL_KMERS.contains(h.getValueLabel())),
                "a histogram contains an unsupported, non-kmer header");
    } else {
        Utils.validate(histogramList.size() == F1R2FilterConstants.NUM_KMERS * F1R2FilterConstants.numAltHistogramsPerContext,
                "The list of alt histograms missing some (kmer, alt allele, f1r2) triple");
    }

    for (int i = 1; i < metricsFiles.size(); i++){
        final List<Histogram<Integer>> ithHistograms = metricsFiles.get(i).getAllHistograms();
        for (final Histogram<Integer> jthHistogram : ithHistograms){
            final String refContext = jthHistogram.getValueLabel();
            final Optional<Histogram<Integer>> hist = histogramList.stream().filter(h -> h.getValueLabel().equals(refContext)).findAny();
            Utils.validate(hist.isPresent(),"Missing histogram header for: " + refContext);

            hist.get().addHistogram(jthHistogram);
        }
    }
    return histogramList;
}
 
Example 27
Source Project: picard   Source File: CollectIlluminaLaneMetrics.java    License: MIT License 5 votes vote down vote up
public static File writePhasingMetrics(final Map<Integer, ? extends Collection<Tile>> laneTiles, final File outputDirectory,
                                       final String outputPrefix, final MetricsFile<MetricBase, Comparable<?>> phasingMetricsFile,
                                       final String fileExtension, final boolean isNovaSeq) {
    laneTiles.forEach((key, value) -> IlluminaPhasingMetrics.getPhasingMetricsForTiles(key.longValue(),
            value, !isNovaSeq).forEach(phasingMetricsFile::addMetric));

    return writeMetrics(phasingMetricsFile, outputDirectory, outputPrefix, IlluminaPhasingMetrics.getExtension() + fileExtension);
}
 
Example 28
Source Project: picard   Source File: CollectIlluminaLaneMetrics.java    License: MIT License 5 votes vote down vote up
private static File writeMetrics(final MetricsFile<MetricBase, Comparable<?>> metricsFile, final File outputDirectory,
                                 final String outputPrefix, final String outputExtension) {
    final File outputFile = new File(outputDirectory, String.format("%s.%s", outputPrefix, outputExtension));
    LOG.info(String.format("Writing %s lane metrics to %s ...", metricsFile.getMetrics().size(), outputFile));
    metricsFile.write(outputFile);
    return outputFile;
}
 
Example 29
Source Project: picard   Source File: CollectHsMetricsTest.java    License: MIT License 5 votes vote down vote up
/** Read back the first metrics record in an hs metrics file. */
private HsMetrics readMetrics(final File f) {
    try {
        final MetricsFile<HsMetrics, Comparable<?>> mFile = new MetricsFile<HsMetrics, Comparable<?>>();
        mFile.read(new FileReader(f));
        return mFile.getMetrics().get(0);
    }
    catch (IOException ioe) {
         throw new RuntimeIOException(ioe);
    }
}
 
Example 30
Source Project: picard   Source File: CollectSamErrorMetrics.java    License: MIT License 5 votes vote down vote up
/**
 * Given a "full" BaseCalculator, get all of its metrics and put them in the appropriate Metrics file.
 *
 * @param locusAggregator a BaseCalculator that has been "loaded up" with bases.
 */
private void writeMetricsFileForAggregator(final BaseErrorAggregation locusAggregator) {
    final MetricsFile<ErrorMetric, Integer> file = getMetricsFile();

    ErrorMetric.setPriorError(QualityUtil.getErrorProbabilityFromPhredScore(PRIOR_Q));

    for (final ErrorMetric metric : locusAggregator.getMetrics()) {
        metric.calculateDerivedFields();
        file.addMetric(metric);
    }

    file.write(new File(OUTPUT + "." + locusAggregator.getSuffix()));
}