/*** * Immutable method that divides the current Histogram by an input Histogram and generates a new one * Throws an exception if the bins don't match up exactly * @param divisorHistogram * @return * @throws IllegalArgumentException if the keySet of this histogram is not equal to the keySet of the given divisorHistogram */ public Histogram<K> divideByHistogram(final Histogram<K> divisorHistogram) { final Histogram<K> output = new Histogram<K>(); if (!this.keySet().equals(divisorHistogram.keySet())) throw new IllegalArgumentException("Attempting to divide Histograms with non-identical bins"); for (final K key : this.keySet()){ final Bin<K> dividend = this.get(key); final Bin<K> divisor = divisorHistogram.get(key); output.increment(key, dividend.getValue()/divisor.getValue()); } return output; }
public static InsertSizeDistribution create(Histogram<Integer> insertSizeHistogram) { if (insertSizeHistogram == null) return null; int[] insertSize = new int[insertSizeHistogram.size()]; double[] count = new double[insertSizeHistogram.size()]; double total = insertSizeHistogram.getSumOfValues(); int i = 0; Set<Integer> keys = insertSizeHistogram.keySet(); for (Integer key : keys) { insertSize[i] = key; count[i] = insertSizeHistogram.get(key).getValue(); i++; } return new InsertSizeDistribution(insertSize, count, (long)total); } public InsertSizeDistribution(int[] singletons, double[] readCounts) {
public void finish() { cytoConversionRate = nCytoTotal == 0 ? 0 : nCytoConverted / (double)nCytoTotal; nCpgSeen = (int)cpgTotal.getSumOfValues(); nCpgConverted = (int)cpgConverted.getSumOfValues(); cpgConversionRate = nCpgSeen == 0 ? 0 : nCpgConverted / (double)nCpgSeen; coverageMean = cpgTotal.getMeanBinSize(); coverageMedian = (int)cpgTotal.getMedianBinSize(); }
/** * Calculates the mean bin size */ public double getMeanBinSize() { return (getSumOfValues() / size()); }
totalInserts += h.getCount(); final SamPairUtil.PairOrientation pairOrientation = entry.getKey(); final Histogram<Integer> histogram = entry.getValue(); final double total = histogram.getCount(); 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(); final Histogram.Bin<Integer> highBin = histogram.get((int) high); if (highBin != null) covered += highBin.getValue(); trimmedHistogram.trimByWidth(getWidthToTrimTo(metrics));
/** * Outputs validation summary report to out. * * @param samReader records to validate * @param reference if null, NM tag validation is skipped * @return boolean true if there are no validation errors, otherwise false */ public boolean validateSamFileSummary(final SamReader samReader, final ReferenceSequenceFile reference) { init(reference, samReader.getFileHeader()); validateSamFile(samReader, out); boolean result = errorsByType.isEmpty(); if (errorsByType.getCount() > 0) { // Convert to a histogram with String IDs so that WARNING: or ERROR: can be prepended to the error type. final Histogram<String> errorsAndWarningsByType = new Histogram<>("Error Type", "Count"); for (final Histogram.Bin<Type> bin : errorsByType.values()) { errorsAndWarningsByType.increment(bin.getId().getHistogramString(), bin.getValue()); } final MetricsFile<ValidationMetrics, String> metricsFile = new MetricsFile<>(); errorsByType.setBinLabel("Error Type"); errorsByType.setValueLabel("Count"); metricsFile.setHistogram(errorsAndWarningsByType); metricsFile.write(out); } cleanup(); return result; }
for (final Histogram.Bin<Integer> bin : highQualityDepthHistogram.values()) { final int depth = Math.min((int) bin.getIdValue(), coverageCap); depthHistogramArray[depth] += bin.getValue(); GENOME_TERRITORY = (long) highQualityDepthHistogram.getSumOfValues(); MEAN_COVERAGE = highQualityDepthHistogram.getMean(); SD_COVERAGE = highQualityDepthHistogram.getStandardDeviation(); MEDIAN_COVERAGE = highQualityDepthHistogram.getMedian(); MAD_COVERAGE = highQualityDepthHistogram.getMedianAbsoluteDeviation();
/** Gets the median absolute deviation of the distribution. */ public double getMedianAbsoluteDeviation() { final double median = getMedian(); final Histogram<Double> deviations = new Histogram<>(); for (final Bin<K> bin : values()) { final double dev = abs(bin.getIdValue() - median); deviations.increment(dev, bin.getValue()); } return deviations.getMedian(); }
Histogram<Integer> innieHistogram = new Histogram<Integer>(); Histogram<Integer> outieHistogram = new Histogram<Integer>(); final PairOrientation pairOrientation = SamPairUtil.getPairOrientation(sam); if (pairOrientation == PairOrientation.RF) { outieHistogram.increment(absInsertSize); outies++; if (sam.getDuplicateReadFlag()) { innieHistogram.increment(absInsertSize); innies++; if (sam.getDuplicateReadFlag()) { 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;
final List<Histogram<HKEY>> nonEmptyHistograms = new ArrayList<Histogram<HKEY>>(); for (final Histogram<HKEY> histo : this.histograms) { if (!histo.isEmpty()) nonEmptyHistograms.add(histo); final java.util.Set<HKEY> keys = new TreeSet<HKEY>(nonEmptyHistograms.get(0).comparator()); for (final Histogram<HKEY> histo : nonEmptyHistograms) { if (histo != null) keys.addAll(histo.keySet()); out.append(HISTO_HEADER + nonEmptyHistograms.get(0).keySet().iterator().next().getClass().getName()); out.newLine(); out.append(StringUtil.assertCharactersNotInString(nonEmptyHistograms.get(0).getBinLabel(), '\t', '\n')); for (final Histogram<HKEY> histo : nonEmptyHistograms) { out.append(SEPARATOR); out.append(StringUtil.assertCharactersNotInString(histo.getValueLabel(), '\t', '\n')); final Histogram.Bin<HKEY> bin = histo.get(key); final double value = (bin == null ? 0 : bin.getValue());
@Test public void testMad() { final int[] is = {4,4,4,4,5,5,5,5,6,6,6,6,6,6,7,7,7,7,7,7,7,8,8,8,8,8,8,8,8}; final Histogram<Integer> histo = new Histogram<>(); for (final int i : is) histo.increment(i); Assert.assertEquals(7d, histo.getMedian()); Assert.assertEquals(1d, histo.getMedianAbsoluteDeviation()); Assert.assertTrue(abs(histo.estimateSdViaMad() - 1.4826) < 0.0001); }
@Test(dataProvider = "histogramData") public void testHistogramFunctions(final int[] values, final double mean, final double stdev, final Integer trimByWidth) { final Histogram<Integer> histo = new Histogram<>(); for (int value : values) { histo.increment(value); } if (trimByWidth != null) histo.trimByWidth(trimByWidth); final double m = histo.getMean(); final double sd = histo.getStandardDeviation(); Assert.assertEquals(round(mean), round(m), "Means are not equal"); Assert.assertEquals(round(stdev), round(sd), "Stdevs are not equal"); }
public static double[] normalizeHistogram(final Histogram<Integer> histogram) { if (histogram == null) throw new PicardException("Histogram is null and cannot be normalized"); final double histogramSumOfValues = histogram.getSumOfValues(); final double[] normalizedHistogram = new double[histogram.size()]; for (int i = 0; i < histogram.size(); i++) { if (histogram.get(i) != null) { normalizedHistogram[i] = histogram.get(i).getValue() / histogramSumOfValues; } } return normalizedHistogram; }
public void onComplete() { //summarize read data if (metrics.TOTAL_READS > 0) { metrics.PCT_PF_READS = (double) metrics.PF_READS / (double) metrics.TOTAL_READS; metrics.PCT_ADAPTER = this.adapterReads / (double) metrics.PF_READS; metrics.MEAN_READ_LENGTH = readLengthHistogram.getMean(); //Calculate BAD_CYCLES metrics.BAD_CYCLES = 0; for (final Histogram.Bin<Integer> cycleBin : badCycleHistogram.values()) { final double badCyclePercentage = cycleBin.getValue() / metrics.TOTAL_READS; if (badCyclePercentage >= 0.8) { metrics.BAD_CYCLES++; } } if(doRefMetrics) { if (metrics.PF_READS > 0) metrics.PCT_PF_READS_ALIGNED = (double) metrics.PF_READS_ALIGNED / (double) metrics.PF_READS; if (metrics.PF_READS_ALIGNED > 0) metrics.PCT_READS_ALIGNED_IN_PAIRS = (double) metrics.READS_ALIGNED_IN_PAIRS / (double) metrics.PF_READS_ALIGNED; if (metrics.PF_READS_ALIGNED > 0) metrics.PCT_PF_READS_IMPROPER_PAIRS = (double) metrics.PF_READS_IMPROPER_PAIRS / (double) metrics.PF_READS_ALIGNED; if (metrics.PF_READS_ALIGNED > 0) metrics.STRAND_BALANCE = numPositiveStrand / (double) metrics.PF_READS_ALIGNED; if (this.chimerasDenominator > 0) metrics.PCT_CHIMERAS = this.chimeras / (double) this.chimerasDenominator; if (nonBisulfiteAlignedBases > 0) metrics.PF_MISMATCH_RATE = mismatchHistogram.getSum() / (double) nonBisulfiteAlignedBases; metrics.PF_HQ_MEDIAN_MISMATCHES = hqMismatchHistogram.getMedian(); if (hqNonBisulfiteAlignedBases > 0) metrics.PF_HQ_ERROR_RATE = hqMismatchHistogram.getSum() / (double) hqNonBisulfiteAlignedBases; if (metrics.PF_ALIGNED_BASES > 0) metrics.PF_INDEL_RATE = this.indels / (double) metrics.PF_ALIGNED_BASES; } } }