Merge branch 'master' of ssh://nickel.broadinstitute.org/humgen/gsa-scr1/gsa-engineering/git/unstable

This commit is contained in:
Ryan Poplin 2012-02-03 17:14:20 -05:00
commit f94d547e97
8 changed files with 486 additions and 151 deletions

View File

@ -0,0 +1,22 @@
package org.broadinstitute.sting.gatk.walkers.diagnostics.targets;
/**
* Short one line description of the walker.
*
* @author Mauricio Carneiro
* @since 2/1/12
*/
public enum CallableStatus {
/** the reference base was an N, which is not considered callable the GATK */
REF_N,
/** the base satisfied the min. depth for calling but had less than maxDepth to avoid having EXCESSIVE_COVERAGE */
CALLABLE,
/** absolutely no reads were seen at this locus, regardless of the filtering parameters */
NO_COVERAGE,
/** there were less than min. depth bases at the locus, after applying filters */
LOW_COVERAGE,
/** more than -maxDepth read at the locus, indicating some sort of mapping problem */
EXCESSIVE_COVERAGE,
/** more than --maxFractionOfReadsWithLowMAPQ at the locus, indicating a poor mapping quality of the reads */
POOR_QUALITY
}

View File

@ -0,0 +1,172 @@
package org.broadinstitute.sting.gatk.walkers.diagnostics.targets;
import org.broad.tribble.Feature;
import org.broadinstitute.sting.commandline.Argument;
import org.broadinstitute.sting.commandline.Input;
import org.broadinstitute.sting.commandline.IntervalBinding;
import org.broadinstitute.sting.commandline.Output;
import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
import org.broadinstitute.sting.gatk.walkers.By;
import org.broadinstitute.sting.gatk.walkers.DataSource;
import org.broadinstitute.sting.gatk.walkers.LocusWalker;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocComparator;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.exceptions.UserException;
import java.io.PrintStream;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import java.util.TreeSet;
/**
* Short one line description of the walker.
*
* <p>
* [Long description of the walker]
* </p>
*
*
* <h2>Input</h2>
* <p>
* [Description of the Input]
* </p>
*
* <h2>Output</h2>
* <p>
* [Description of the Output]
* </p>
*
* <h2>Examples</h2>
* <pre>
* java
* -jar GenomeAnalysisTK.jar
* -T [walker name]
* </pre>
*
* @author Mauricio Carneiro
* @since 2/1/12
*/
@By(value = DataSource.READS)
public class DiagnoseTargets extends LocusWalker<Long, Long> {
@Input(fullName = "interval_track", shortName = "int", doc = "", required = true)
private IntervalBinding<Feature> intervalTrack = null;
@Output
private PrintStream out = System.out;
@Argument(fullName = "expand_interval", shortName = "exp", doc = "", required = false)
private int expandInterval = 50;
@Argument(fullName = "minimum_base_quality", shortName = "mbq", doc = "", required = false)
private int minimumBaseQuality = 20;
@Argument(fullName = "minimum_mapping_quality", shortName = "mmq", doc = "", required = false)
private int minimumMappingQuality = 20;
@Argument(fullName = "minimum_coverage", shortName = "mincov", doc = "", required = false)
private int minimumCoverage = 5;
@Argument(fullName = "maximum_coverage", shortName = "maxcov", doc = "", required = false)
private int maximumCoverage = 700;
private TreeSet<GenomeLoc> intervalList = null; // The list of intervals of interest (plus expanded intervals if user wants them)
private HashMap<GenomeLoc, IntervalStatistics> intervalMap = null; // interval => statistics
private Iterator<GenomeLoc> intervalListIterator; // An iterator to go over all the intervals provided as we traverse the genome
private GenomeLoc currentInterval = null; // The "current" interval loaded and being filled with statistics
private IntervalStatistics currentIntervalStatistics = null; // The "current" interval loaded and being filled with statistics
private GenomeLocParser parser; // just an object to allow us to create genome locs (for the expanded intervals)
@Override
public void initialize() {
super.initialize();
if (intervalTrack == null)
throw new UserException("This tool currently only works if you provide an interval track");
parser = new GenomeLocParser(getToolkit().getMasterSequenceDictionary()); // Important to initialize the parser before creating the intervals below
List<GenomeLoc> originalList = intervalTrack.getIntervals(getToolkit()); // The original list of targets provided by the user that will be expanded or not depending on the options provided
intervalList = new TreeSet<GenomeLoc>(new GenomeLocComparator());
intervalMap = new HashMap<GenomeLoc, IntervalStatistics>(originalList.size() * 2);
for (GenomeLoc interval : originalList)
addAndExpandIntervalToLists(interval);
intervalListIterator = intervalList.iterator();
}
@Override
public Long map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
GenomeLoc refLocus = ref.getLocus();
while (currentInterval == null || currentInterval.isBefore(refLocus)) {
if (!intervalListIterator.hasNext())
return 0L;
currentInterval = intervalListIterator.next();
currentIntervalStatistics = intervalMap.get(currentInterval);
}
if (currentInterval.isPast(refLocus))
return 0L;
byte[] mappingQualities = context.getBasePileup().getMappingQuals();
byte[] baseQualities = context.getBasePileup().getQuals();
int coverage = context.getBasePileup().getBaseAndMappingFilteredPileup(minimumBaseQuality, minimumMappingQuality).depthOfCoverage();
int rawCoverage = context.size();
IntervalStatisticLocus locusData = new IntervalStatisticLocus(mappingQualities, baseQualities, coverage, rawCoverage);
currentIntervalStatistics.addLocus(refLocus, locusData);
return 1L;
}
@Override
public Long reduceInit() {
return 0L;
}
@Override
public Long reduce(Long value, Long sum) {
return sum + value;
}
@Override
public void onTraversalDone(Long result) {
super.onTraversalDone(result);
out.println("Interval\tCallStatus\tCOV\tAVG");
for (GenomeLoc interval : intervalList) {
IntervalStatistics stats = intervalMap.get(interval);
out.println(String.format("%s\t%s\t%d\t%f", interval, stats.callableStatus(), stats.totalCoverage(), stats.averageCoverage()));
}
}
private GenomeLoc createIntervalBefore(GenomeLoc interval) {
int start = Math.max(interval.getStart() - expandInterval, 0);
int stop = Math.max(interval.getStart() - 1, 0);
return parser.createGenomeLoc(interval.getContig(), interval.getContigIndex(), start, stop);
}
private GenomeLoc createIntervalAfter(GenomeLoc interval) {
int contigLimit = getToolkit().getSAMFileHeader().getSequenceDictionary().getSequence(interval.getContigIndex()).getSequenceLength();
int start = Math.min(interval.getStop() + 1, contigLimit);
int stop = Math.min(interval.getStop() + expandInterval, contigLimit);
return parser.createGenomeLoc(interval.getContig(), interval.getContigIndex(), start, stop);
}
private void addAndExpandIntervalToLists(GenomeLoc interval) {
if (expandInterval > 0) {
GenomeLoc before = createIntervalBefore(interval);
GenomeLoc after = createIntervalAfter(interval);
intervalList.add(before);
intervalList.add(after);
intervalMap.put(before, new IntervalStatistics(before, minimumCoverage, maximumCoverage, minimumMappingQuality, minimumBaseQuality));
intervalMap.put(after, new IntervalStatistics(after, minimumCoverage, maximumCoverage, minimumMappingQuality, minimumBaseQuality));
}
intervalList.add(interval);
intervalMap.put(interval, new IntervalStatistics(interval, minimumCoverage, maximumCoverage, minimumMappingQuality, minimumBaseQuality));
}
}

View File

@ -0,0 +1,34 @@
package org.broadinstitute.sting.gatk.walkers.diagnostics.targets;
/**
* The definition of a locus for the DiagnoseTargets walker statistics calculation
*
* @author Mauricio Carneiro
* @since 2/3/12
*/
class IntervalStatisticLocus {
private final byte[] mappingQuality;
private final byte[] baseQuality;
private final int coverage;
private final int rawCoverage;
public IntervalStatisticLocus(byte[] mappingQuality, byte[] baseQuality, int coverage, int rawCoverage) {
this.mappingQuality = mappingQuality;
this.baseQuality = baseQuality;
this.coverage = coverage;
this.rawCoverage = rawCoverage;
}
public IntervalStatisticLocus() {
this(new byte[1], new byte[1], 0, 0);
}
public int getCoverage() {
return coverage;
}
public int getRawCoverage() {
return rawCoverage;
}
}

View File

@ -0,0 +1,122 @@
package org.broadinstitute.sting.gatk.walkers.diagnostics.targets;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.exceptions.ReviewedStingException;
import java.util.ArrayList;
import java.util.HashMap;
/**
* Short one line description of the walker.
*
* @author Mauricio Carneiro
* @since 2/1/12
*/
class IntervalStatistics {
private final GenomeLoc interval;
private final ArrayList<IntervalStatisticLocus> loci;
private final int minimumCoverageThreshold;
private final int maximumCoverageThreshold;
private final int minimumMappingQuality;
private final int minimumBaseQuality;
private int preComputedTotalCoverage = -1; // avoids re-calculating the total sum (-1 means we haven't pre-computed it yet)
private IntervalStatistics(GenomeLoc interval, ArrayList<IntervalStatisticLocus> loci, int minimumCoverageThreshold, int maximumCoverageThreshold, int minimumMappingQuality, int minimumBaseQuality) {
this.interval = interval;
this.loci = loci;
this.minimumCoverageThreshold = minimumCoverageThreshold;
this.maximumCoverageThreshold = maximumCoverageThreshold;
this.minimumMappingQuality = minimumMappingQuality;
this.minimumBaseQuality = minimumBaseQuality;
}
public IntervalStatistics(GenomeLoc interval, int minimumCoverageThreshold, int maximumCoverageThreshold, int minimumMappingQuality, int minimumBaseQuality) {
this(interval, new ArrayList<IntervalStatisticLocus>(interval.size()), minimumCoverageThreshold, maximumCoverageThreshold, minimumMappingQuality, minimumBaseQuality);
// Initialize every loci (this way we don't have to worry about non-existent loci in the object
for (int i = 0; i < interval.size(); i++)
this.loci.add(i, new IntervalStatisticLocus());
}
public long totalCoverage() {
if (preComputedTotalCoverage < 0)
calculateTotalCoverage();
return preComputedTotalCoverage;
}
public double averageCoverage() {
if (preComputedTotalCoverage < 0)
calculateTotalCoverage();
return (double) preComputedTotalCoverage / loci.size();
}
/**
* Calculates the callable status of the entire interval
*
* @return the callable status of the entire interval
*/
public CallableStatus callableStatus() {
long max = -1;
CallableStatus maxCallableStatus = null;
HashMap<CallableStatus, Integer> statusCounts = new HashMap<CallableStatus, Integer>(CallableStatus.values().length);
// initialize the statusCounts with all callable states
for (CallableStatus key : CallableStatus.values())
statusCounts.put(key, 0);
// calculate the callable status for each locus
for (int i = 0; i < loci.size(); i++) {
CallableStatus status = callableStatus(i);
int count = statusCounts.get(status) + 1;
statusCounts.put(status, count);
if (count > max) {
max = count;
maxCallableStatus = status;
}
}
return maxCallableStatus;
}
public void addLocus(GenomeLoc locus, IntervalStatisticLocus locusData) {
if (!interval.containsP(locus))
throw new ReviewedStingException(String.format("Locus %s is not part of the Interval", locus));
int locusIndex = locus.getStart() - interval.getStart();
loci.add(locusIndex, locusData);
}
/**
* returns the callable status of this locus without taking the reference base into account.
*
* @param locusIndex location in the genome to inquire (only one locus)
* @return the callable status of a locus
*/
private CallableStatus callableStatus(int locusIndex) {
if (loci.get(locusIndex).getCoverage() > maximumCoverageThreshold)
return CallableStatus.EXCESSIVE_COVERAGE;
if (loci.get(locusIndex).getCoverage() >= minimumCoverageThreshold)
return CallableStatus.CALLABLE;
if (loci.get(locusIndex).getRawCoverage() >= minimumCoverageThreshold)
return CallableStatus.POOR_QUALITY;
if (loci.get(locusIndex).getRawCoverage() > 0)
return CallableStatus.LOW_COVERAGE;
return CallableStatus.NO_COVERAGE;
}
private void calculateTotalCoverage() {
preComputedTotalCoverage = 0;
for (IntervalStatisticLocus locus : loci)
preComputedTotalCoverage += locus.getCoverage();
}
}

View File

@ -436,7 +436,7 @@ public class GenomeLoc implements Comparable<GenomeLoc>, Serializable, HasGenome
* never be < 1.
*/
@Ensures("result > 0")
public long size() {
public int size() {
return stop - start + 1;
}

View File

@ -49,7 +49,6 @@ public class MathUtils {
* high precision
*/
/**
* Private constructor. No instantiating this class!
*/
@ -60,48 +59,48 @@ public class MathUtils {
// under/overflow checking, so this shouldn't be used in the general case (but is fine
// if one is already make those checks before calling in to the rounding).
public static int fastRound(double d) {
return (d > 0) ? (int)(d + 0.5d) : (int)(d - 0.5d);
return (d > 0) ? (int) (d + 0.5d) : (int) (d - 0.5d);
}
public static double approximateLog10SumLog10(final double[] vals) {
return approximateLog10SumLog10(vals, vals.length);
return approximateLog10SumLog10(vals, vals.length);
}
public static double approximateLog10SumLog10(final double[] vals, final int endIndex) {
final int maxElementIndex = MathUtils.maxElementIndex(vals, endIndex);
double approxSum = vals[maxElementIndex];
if ( approxSum == Double.NEGATIVE_INFINITY )
final int maxElementIndex = MathUtils.maxElementIndex(vals, endIndex);
double approxSum = vals[maxElementIndex];
if (approxSum == Double.NEGATIVE_INFINITY)
return approxSum;
for ( int i = 0; i < endIndex; i++ ) {
if ( i == maxElementIndex || vals[i] == Double.NEGATIVE_INFINITY )
continue;
for (int i = 0; i < endIndex; i++) {
if (i == maxElementIndex || vals[i] == Double.NEGATIVE_INFINITY)
continue;
final double diff = approxSum - vals[i];
if ( diff < MathUtils.MAX_JACOBIAN_TOLERANCE ) {
// See notes from the 2-inout implementation below
final int ind = fastRound(diff / MathUtils.JACOBIAN_LOG_TABLE_STEP); // hard rounding
approxSum += MathUtils.jacobianLogTable[ind];
}
}
final double diff = approxSum - vals[i];
if (diff < MathUtils.MAX_JACOBIAN_TOLERANCE) {
// See notes from the 2-inout implementation below
final int ind = fastRound(diff / MathUtils.JACOBIAN_LOG_TABLE_STEP); // hard rounding
approxSum += MathUtils.jacobianLogTable[ind];
}
}
return approxSum;
}
public static double approximateLog10SumLog10(double small, double big) {
// make sure small is really the smaller value
if ( small > big ) {
if (small > big) {
final double t = big;
big = small;
small = t;
}
if ( small == Double.NEGATIVE_INFINITY || big == Double.NEGATIVE_INFINITY )
if (small == Double.NEGATIVE_INFINITY || big == Double.NEGATIVE_INFINITY)
return big;
final double diff = big - small;
if ( diff >= MathUtils.MAX_JACOBIAN_TOLERANCE )
final double diff = big - small;
if (diff >= MathUtils.MAX_JACOBIAN_TOLERANCE)
return big;
// OK, so |y-x| < tol: we use the following identity then:
@ -138,10 +137,15 @@ public class MathUtils {
return size;
}
public static double average(Collection<Integer> x) {
return (double) sum(x) / x.size();
}
public static double average(Collection<Number> numbers, boolean ignoreNan) {
if (ignoreNan) {
return sum(numbers, true) / nonNanSize(numbers);
} else {
}
else {
return sum(numbers, false) / nonNanSize(numbers);
}
}
@ -172,10 +176,17 @@ public class MathUtils {
public static double sum(double[] values) {
double s = 0.0;
for (double v : values) s += v;
for (double v : values)
s += v;
return s;
}
public static long sum(int[] x) {
long total = 0;
for (int v : x)
total += v;
return total;
}
/**
* Calculates the log10 cumulative sum of an array with log10 probabilities
@ -218,21 +229,23 @@ public class MathUtils {
public static double sumDoubles(List<Double> values) {
double s = 0.0;
for (double v : values) s += v;
for (double v : values)
s += v;
return s;
}
public static int sumIntegers(List<Integer> values) {
int s = 0;
for (int v : values) s += v;
for (int v : values)
s += v;
return s;
}
public static double sumLog10(double[] log10values) {
return Math.pow(10.0, log10sumLog10(log10values));
// double s = 0.0;
// for ( double v : log10values) s += Math.pow(10.0, v);
// return s;
// double s = 0.0;
// for ( double v : log10values) s += Math.pow(10.0, v);
// return s;
}
public static double log10sumLog10(double[] log10values) {
@ -445,7 +458,6 @@ public class MathUtils {
return Math.sqrt(rms);
}
/**
* calculate the Root Mean Square of an array of integers
*
@ -506,7 +518,6 @@ public class MathUtils {
return result;
}
/**
* normalizes the log10-based array. ASSUMES THAT ALL ARRAY ENTRIES ARE <= 0 (<= 1 IN REAL-SPACE).
*
@ -543,7 +554,8 @@ public class MathUtils {
sum += normalized[i];
for (int i = 0; i < array.length; i++) {
double x = normalized[i] / sum;
if (takeLog10OfOutput) x = Math.log10(x);
if (takeLog10OfOutput)
x = Math.log10(x);
normalized[i] = x;
}
@ -565,7 +577,8 @@ public class MathUtils {
sum += normalized[i];
for (int i = 0; i < array.size(); i++) {
double x = normalized[i] / sum;
if (takeLog10OfOutput) x = Math.log10(x);
if (takeLog10OfOutput)
x = Math.log10(x);
normalized[i] = x;
}
@ -587,11 +600,12 @@ public class MathUtils {
}
public static int maxElementIndex(final double[] array) {
return maxElementIndex(array, array.length);
return maxElementIndex(array, array.length);
}
public static int maxElementIndex(final double[] array, final int endIndex) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
if (array == null)
throw new IllegalArgumentException("Array cannot be null!");
int maxI = -1;
for (int i = 0; i < endIndex; i++) {
@ -603,11 +617,12 @@ public class MathUtils {
}
public static int maxElementIndex(final int[] array) {
return maxElementIndex(array, array.length);
return maxElementIndex(array, array.length);
}
public static int maxElementIndex(final int[] array, int endIndex) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
if (array == null)
throw new IllegalArgumentException("Array cannot be null!");
int maxI = -1;
for (int i = 0; i < endIndex; i++) {
@ -635,7 +650,8 @@ public class MathUtils {
}
public static int minElementIndex(double[] array) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
if (array == null)
throw new IllegalArgumentException("Array cannot be null!");
int minI = -1;
for (int i = 0; i < array.length; i++) {
@ -647,7 +663,8 @@ public class MathUtils {
}
public static int minElementIndex(byte[] array) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
if (array == null)
throw new IllegalArgumentException("Array cannot be null!");
int minI = -1;
for (int i = 0; i < array.length; i++) {
@ -659,7 +676,8 @@ public class MathUtils {
}
public static int minElementIndex(int[] array) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
if (array == null)
throw new IllegalArgumentException("Array cannot be null!");
int minI = -1;
for (int i = 0; i < array.length; i++) {
@ -671,20 +689,26 @@ public class MathUtils {
}
public static int arrayMaxInt(List<Integer> array) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
if (array.size() == 0) throw new IllegalArgumentException("Array size cannot be 0!");
if (array == null)
throw new IllegalArgumentException("Array cannot be null!");
if (array.size() == 0)
throw new IllegalArgumentException("Array size cannot be 0!");
int m = array.get(0);
for (int e : array) m = Math.max(m, e);
for (int e : array)
m = Math.max(m, e);
return m;
}
public static double arrayMaxDouble(List<Double> array) {
if (array == null) throw new IllegalArgumentException("Array cannot be null!");
if (array.size() == 0) throw new IllegalArgumentException("Array size cannot be 0!");
if (array == null)
throw new IllegalArgumentException("Array cannot be null!");
if (array.size() == 0)
throw new IllegalArgumentException("Array size cannot be 0!");
double m = array.get(0);
for (double e : array) m = Math.max(m, e);
for (double e : array)
m = Math.max(m, e);
return m;
}
@ -722,6 +746,13 @@ public class MathUtils {
return average(vals, vals.size());
}
public static double average(int[] x) {
int sum = 0;
for (int v : x)
sum += v;
return (double) sum / x.length;
}
public static byte average(byte[] vals) {
int sum = 0;
for (byte v : vals) {
@ -798,7 +829,6 @@ public class MathUtils {
return permutation;
}
public static int[] permuteArray(int[] array, Integer[] permutation) {
int[] output = new int[array.length];
for (int i = 0; i < output.length; i++) {
@ -839,7 +869,6 @@ public class MathUtils {
return output;
}
/**
* Draw N random elements from list.
*/
@ -905,7 +934,8 @@ public class MathUtils {
public static <T> int countOccurrences(T x, List<T> l) {
int count = 0;
for (T y : l) {
if (x.equals(y)) count++;
if (x.equals(y))
count++;
}
return count;
@ -1013,9 +1043,11 @@ public class MathUtils {
for (Comparable y : list) {
if (x.compareTo(y) > 0) {
lessThanX.add(y);
} else if (x.compareTo(y) < 0) {
}
else if (x.compareTo(y) < 0) {
greaterThanX.add(y);
} else
}
else
equalToX.add(y);
}
@ -1028,7 +1060,6 @@ public class MathUtils {
}
public static Object getMedian(List<Comparable> list) {
return orderStatisticSearch((int) Math.ceil(list.size() / 2), list);
}
@ -1058,10 +1089,12 @@ public class MathUtils {
if (quality < qk) {
lessThanQReads.add(read);
lessThanQOffsets.add(offset);
} else if (quality > qk) {
}
else if (quality > qk) {
greaterThanQReads.add(read);
greaterThanQOffsets.add(offset);
} else {
}
else {
equalToQReads.add(reads.get(iter));
}
}
@ -1079,6 +1112,13 @@ public class MathUtils {
return getQScoreOrderStatistic(reads, offsets, (int) Math.floor(reads.size() / 2.));
}
public static long sum(Collection<Integer> x) {
long sum = 0;
for (int v : x)
sum += v;
return sum;
}
/**
* A utility class that computes on the fly average and standard deviation for a stream of numbers.
* The number of observations does not have to be known in advance, and can be also very big (so that
@ -1184,8 +1224,7 @@ public class MathUtils {
log10Cache[k] = Math.log10(k);
for (int k = 0; k < JACOBIAN_LOG_TABLE_SIZE; k++) {
jacobianLogTable[k] = Math.log10(1.0 + Math.pow(10.0, -((double) k)
* JACOBIAN_LOG_TABLE_STEP));
jacobianLogTable[k] = Math.log10(1.0 + Math.pow(10.0, -((double) k) * JACOBIAN_LOG_TABLE_STEP));
}
}
@ -1232,7 +1271,8 @@ public class MathUtils {
else if (diff >= 0) {
int ind = (int) (diff * INV_JACOBIAN_LOG_TABLE_STEP + 0.5);
return x + jacobianLogTable[ind];
} else {
}
else {
int ind = (int) (-diff * INV_JACOBIAN_LOG_TABLE_STEP + 0.5);
return y + jacobianLogTable[ind];
}
@ -1273,71 +1313,7 @@ public class MathUtils {
/**
* Constants to simplify the log gamma function calculation.
*/
private static final double
zero = 0.0,
one = 1.0,
half = .5,
a0 = 7.72156649015328655494e-02,
a1 = 3.22467033424113591611e-01,
a2 = 6.73523010531292681824e-02,
a3 = 2.05808084325167332806e-02,
a4 = 7.38555086081402883957e-03,
a5 = 2.89051383673415629091e-03,
a6 = 1.19270763183362067845e-03,
a7 = 5.10069792153511336608e-04,
a8 = 2.20862790713908385557e-04,
a9 = 1.08011567247583939954e-04,
a10 = 2.52144565451257326939e-05,
a11 = 4.48640949618915160150e-05,
tc = 1.46163214496836224576e+00,
tf = -1.21486290535849611461e-01,
tt = -3.63867699703950536541e-18,
t0 = 4.83836122723810047042e-01,
t1 = -1.47587722994593911752e-01,
t2 = 6.46249402391333854778e-02,
t3 = -3.27885410759859649565e-02,
t4 = 1.79706750811820387126e-02,
t5 = -1.03142241298341437450e-02,
t6 = 6.10053870246291332635e-03,
t7 = -3.68452016781138256760e-03,
t8 = 2.25964780900612472250e-03,
t9 = -1.40346469989232843813e-03,
t10 = 8.81081882437654011382e-04,
t11 = -5.38595305356740546715e-04,
t12 = 3.15632070903625950361e-04,
t13 = -3.12754168375120860518e-04,
t14 = 3.35529192635519073543e-04,
u0 = -7.72156649015328655494e-02,
u1 = 6.32827064025093366517e-01,
u2 = 1.45492250137234768737e+00,
u3 = 9.77717527963372745603e-01,
u4 = 2.28963728064692451092e-01,
u5 = 1.33810918536787660377e-02,
v1 = 2.45597793713041134822e+00,
v2 = 2.12848976379893395361e+00,
v3 = 7.69285150456672783825e-01,
v4 = 1.04222645593369134254e-01,
v5 = 3.21709242282423911810e-03,
s0 = -7.72156649015328655494e-02,
s1 = 2.14982415960608852501e-01,
s2 = 3.25778796408930981787e-01,
s3 = 1.46350472652464452805e-01,
s4 = 2.66422703033638609560e-02,
s5 = 1.84028451407337715652e-03,
s6 = 3.19475326584100867617e-05,
r1 = 1.39200533467621045958e+00,
r2 = 7.21935547567138069525e-01,
r3 = 1.71933865632803078993e-01,
r4 = 1.86459191715652901344e-02,
r5 = 7.77942496381893596434e-04,
r6 = 7.32668430744625636189e-06,
w0 = 4.18938533204672725052e-01,
w1 = 8.33333333333329678849e-02,
w2 = -2.77777777728775536470e-03,
w3 = 7.93650558643019558500e-04,
w4 = -5.95187557450339963135e-04,
w5 = 8.36339918996282139126e-04,
w6 = -1.63092934096575273989e-03;
private static final double zero = 0.0, one = 1.0, half = .5, a0 = 7.72156649015328655494e-02, a1 = 3.22467033424113591611e-01, a2 = 6.73523010531292681824e-02, a3 = 2.05808084325167332806e-02, a4 = 7.38555086081402883957e-03, a5 = 2.89051383673415629091e-03, a6 = 1.19270763183362067845e-03, a7 = 5.10069792153511336608e-04, a8 = 2.20862790713908385557e-04, a9 = 1.08011567247583939954e-04, a10 = 2.52144565451257326939e-05, a11 = 4.48640949618915160150e-05, tc = 1.46163214496836224576e+00, tf = -1.21486290535849611461e-01, tt = -3.63867699703950536541e-18, t0 = 4.83836122723810047042e-01, t1 = -1.47587722994593911752e-01, t2 = 6.46249402391333854778e-02, t3 = -3.27885410759859649565e-02, t4 = 1.79706750811820387126e-02, t5 = -1.03142241298341437450e-02, t6 = 6.10053870246291332635e-03, t7 = -3.68452016781138256760e-03, t8 = 2.25964780900612472250e-03, t9 = -1.40346469989232843813e-03, t10 = 8.81081882437654011382e-04, t11 = -5.38595305356740546715e-04, t12 = 3.15632070903625950361e-04, t13 = -3.12754168375120860518e-04, t14 = 3.35529192635519073543e-04, u0 = -7.72156649015328655494e-02, u1 = 6.32827064025093366517e-01, u2 = 1.45492250137234768737e+00, u3 = 9.77717527963372745603e-01, u4 = 2.28963728064692451092e-01, u5 = 1.33810918536787660377e-02, v1 = 2.45597793713041134822e+00, v2 = 2.12848976379893395361e+00, v3 = 7.69285150456672783825e-01, v4 = 1.04222645593369134254e-01, v5 = 3.21709242282423911810e-03, s0 = -7.72156649015328655494e-02, s1 = 2.14982415960608852501e-01, s2 = 3.25778796408930981787e-01, s3 = 1.46350472652464452805e-01, s4 = 2.66422703033638609560e-02, s5 = 1.84028451407337715652e-03, s6 = 3.19475326584100867617e-05, r1 = 1.39200533467621045958e+00, r2 = 7.21935547567138069525e-01, r3 = 1.71933865632803078993e-01, r4 = 1.86459191715652901344e-02, r5 = 7.77942496381893596434e-04, r6 = 7.32668430744625636189e-06, w0 = 4.18938533204672725052e-01, w1 = 8.33333333333329678849e-02, w2 = -2.77777777728775536470e-03, w3 = 7.93650558643019558500e-04, w4 = -5.95187557450339963135e-04, w5 = 8.36339918996282139126e-04, w6 = -1.63092934096575273989e-03;
/**
* Efficient rounding functions to simplify the log gamma function calculation
@ -1368,14 +1344,17 @@ public class MathUtils {
/* purge off +-inf, NaN, +-0, and negative arguments */
int ix = hx & 0x7fffffff;
if (ix >= 0x7ff00000) return Double.POSITIVE_INFINITY;
if ((ix | lx) == 0 || hx < 0) return Double.NaN;
if (ix >= 0x7ff00000)
return Double.POSITIVE_INFINITY;
if ((ix | lx) == 0 || hx < 0)
return Double.NaN;
if (ix < 0x3b900000) { /* |x|<2**-70, return -log(|x|) */
return -Math.log(x);
}
/* purge off 1 and 2 */
if ((((ix - 0x3ff00000) | lx) == 0) || (((ix - 0x40000000) | lx) == 0)) r = 0;
if ((((ix - 0x3ff00000) | lx) == 0) || (((ix - 0x40000000) | lx) == 0))
r = 0;
/* for x < 2.0 */
else if (ix < 0x40000000) {
if (ix <= 0x3feccccc) { /* lgamma(x) = lgamma(x+1)-log(x) */
@ -1383,22 +1362,27 @@ public class MathUtils {
if (ix >= 0x3FE76944) {
y = one - x;
i = 0;
} else if (ix >= 0x3FCDA661) {
}
else if (ix >= 0x3FCDA661) {
y = x - (tc - one);
i = 1;
} else {
}
else {
y = x;
i = 2;
}
} else {
}
else {
r = zero;
if (ix >= 0x3FFBB4C3) {
y = 2.0 - x;
i = 0;
} /* [1.7316,2] */ else if (ix >= 0x3FF3B4C4) {
} /* [1.7316,2] */
else if (ix >= 0x3FF3B4C4) {
y = x - tc;
i = 1;
} /* [1.23,1.73] */ else {
} /* [1.23,1.73] */
else {
y = x - one;
i = 2;
}
@ -1426,7 +1410,8 @@ public class MathUtils {
p2 = one + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * v5))));
r += (-0.5 * y + p1 / p2);
}
} else if (ix < 0x40200000) { /* x < 8.0 */
}
else if (ix < 0x40200000) { /* x < 8.0 */
i = (int) x;
t = zero;
y = x - (double) i;
@ -1449,13 +1434,15 @@ public class MathUtils {
break;
}
/* 8.0 <= x < 2**58 */
} else if (ix < 0x43900000) {
}
else if (ix < 0x43900000) {
t = Math.log(x);
z = one / x;
y = z * z;
w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * w6)))));
r = (x - half) * (t - one) + w;
} else
}
else
/* 2**58 <= x <= inf */
r = x * (Math.log(x) - one);
return r;
@ -1490,7 +1477,6 @@ public class MathUtils {
return log10BinomialCoefficient(n, k) + log10p * k + log10OneMinusP * (n - k);
}
/**
* Calculates the log10 of the multinomial coefficient. Designed to prevent
* overflows even with very large numbers.
@ -1534,7 +1520,6 @@ public class MathUtils {
return log10Gamma(x + 1);
}
/**
* Adds two arrays together and returns a new array with the sum.
*
@ -1572,17 +1557,18 @@ public class MathUtils {
/**
* Vector operations
*
* @param v1 first numerical array
* @param v2 second numerical array
* @return a new array with the elements added
* @return a new array with the elements added
*/
public static <E extends Number> Double[] vectorSum(E v1[], E v2[]) {
if (v1.length != v2.length)
throw new UserException("BUG: vectors v1, v2 of different size in vectorSum()");
Double[] result = new Double[v1.length];
for (int k=0; k < v1.length; k++)
result[k] = v1[k].doubleValue()+v2[k].doubleValue();
for (int k = 0; k < v1.length; k++)
result[k] = v1[k].doubleValue() + v2[k].doubleValue();
return result;
}
@ -1590,19 +1576,19 @@ public class MathUtils {
public static <E extends Number> Double[] scalarTimesVector(E a, E[] v1) {
Double result[] = new Double[v1.length];
for (int k=0; k < v1.length; k++)
result[k] = a.doubleValue()*v1[k].doubleValue();
for (int k = 0; k < v1.length; k++)
result[k] = a.doubleValue() * v1[k].doubleValue();
return result;
}
public static <E extends Number> Double dotProduct(E[] v1, E[] v2) {
public static <E extends Number> Double dotProduct(E[] v1, E[] v2) {
if (v1.length != v2.length)
throw new UserException("BUG: vectors v1, v2 of different size in vectorSum()");
Double result = 0.0;
for (int k=0; k < v1.length; k++)
result += v1[k].doubleValue() *v2[k].doubleValue();
for (int k = 0; k < v1.length; k++)
result += v1[k].doubleValue() * v2[k].doubleValue();
return result;
@ -1610,7 +1596,7 @@ public class MathUtils {
public static double[] vectorLog10(double v1[]) {
double result[] = new double[v1.length];
for (int k=0; k < v1.length; k++)
for (int k = 0; k < v1.length; k++)
result[k] = Math.log10(v1[k]);
return result;
@ -1620,7 +1606,7 @@ public class MathUtils {
// todo - silly overloading, just because Java can't unbox/box arrays of primitive types, and we can't do generics with primitive types!
public static Double[] vectorLog10(Double v1[]) {
Double result[] = new Double[v1.length];
for (int k=0; k < v1.length; k++)
for (int k = 0; k < v1.length; k++)
result[k] = Math.log10(v1[k]);
return result;

View File

@ -314,10 +314,10 @@ public class IntervalUtils {
* @param reference The reference for the intervals.
* @return A map of contig names with their sizes.
*/
public static Map<String, Long> getContigSizes(File reference) {
public static Map<String, Integer> getContigSizes(File reference) {
ReferenceDataSource referenceSource = new ReferenceDataSource(reference);
List<GenomeLoc> locs = GenomeLocSortedSet.createSetFromSequenceDictionary(referenceSource.getReference().getSequenceDictionary()).toList();
Map<String, Long> lengths = new LinkedHashMap<String, Long>();
Map<String, Integer> lengths = new LinkedHashMap<String, Integer>();
for (GenomeLoc loc: locs)
lengths.put(loc.getContig(), loc.size());
return lengths;

View File

@ -8,13 +8,12 @@ import org.broadinstitute.sting.BaseTest;
import org.broadinstitute.sting.commandline.IntervalBinding;
import org.broadinstitute.sting.gatk.GenomeAnalysisEngine;
import org.broadinstitute.sting.gatk.datasources.reference.ReferenceDataSource;
import org.broadinstitute.sting.utils.GenomeLocSortedSet;
import org.testng.Assert;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.GenomeLoc;
import org.broadinstitute.sting.utils.GenomeLocParser;
import org.broadinstitute.sting.utils.GenomeLocSortedSet;
import org.broadinstitute.sting.utils.exceptions.UserException;
import org.broadinstitute.sting.utils.fasta.CachingIndexedFastaSequenceFile;
import org.testng.Assert;
import org.testng.annotations.BeforeClass;
import org.testng.annotations.DataProvider;
import org.testng.annotations.Test;
@ -341,7 +340,7 @@ public class IntervalUtilsUnitTest extends BaseTest {
@Test
public void testGetContigLengths() {
Map<String, Long> lengths = IntervalUtils.getContigSizes(new File(BaseTest.hg18Reference));
Map<String, Integer> lengths = IntervalUtils.getContigSizes(new File(BaseTest.hg18Reference));
Assert.assertEquals((long)lengths.get("chr1"), 247249719);
Assert.assertEquals((long)lengths.get("chr2"), 242951149);
Assert.assertEquals((long)lengths.get("chr3"), 199501827);