org.apache.commons.math3.linear.RealMatrix#walkInOptimizedOrder

Source File: PCATangentNormalizationUtils.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

6 votes

/**
 * Calculate the beta-hats that best fit case read counts given the panel of normals.
 *
 * @param normalsPseudoinverse the log-normalized or reduced-panel pseudoinverse from a panel of normals
 * @param input a {@code TxS} matrix where {@code T} is the number of targets and {@code S} the number of count groups (e.g. case samples).
 * @return never {@code null} an {@code NxS} matrix, where N is the number of samples in
 *  the panel and S the original name of count groups.
 */
@VisibleForTesting
public static RealMatrix calculateBetaHats(final RealMatrix normalsPseudoinverse,
                                           final RealMatrix input,
                                           final double epsilon) {
    Utils.nonNull(normalsPseudoinverse, "Normals inverse matrix cannot be null.");
    Utils.nonNull(input, "Input counts cannot be null.");
    Utils.validateArg(epsilon > 0, String.format("Invalid epsilon value, must be > 0: %f", epsilon));
    final double targetThreshold = (Math.log(epsilon) / Math.log(2)) + 1;

    // copy case samples in order to mask targets in-place and mask (set to zero) targets with coverage below threshold
    final RealMatrix maskedInput = input.copy();
    maskedInput.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(final int row, final int column, final double value) {
            return value > targetThreshold ? value : 0;
        }
    });

    return normalsPseudoinverse.multiply(maskedInput);
}

Source File: HDF5PCACoveragePoNCreationUtils.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

6 votes

/**
 * Final pre-panel normalization that consists of dividing all counts by the median of
 * its column and log it with base 2.
 * <p>
 *     The normalization occurs in-place.
 * </p>
 *
 * @param readCounts the input counts to normalize.
 */
@VisibleForTesting
static void normalizeAndLogReadCounts(final ReadCountCollection readCounts, final Logger logger) {
    final RealMatrix counts = readCounts.counts();
    final Median medianCalculator = new Median();

    final double[] medians = IntStream.range(0, counts.getColumnDimension()).mapToDouble(col -> medianCalculator.evaluate(counts.getColumn(col))).toArray();
    counts.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(final int row, final int column, final double value) {
            return Math.log(Math.max(EPSILON, value / medians[column])) * INV_LN_2;
        }
    });

    logger.info("Counts normalized by the column median and log2'd.");
}

Source File: HDF5PCACoveragePoNCreationUtils.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

6 votes

/**
 * Calculates the median of column medians and subtract it from all counts.
 * @param readCounts the input counts to center.
 * @return the median of medians that has been subtracted from all counts.
 */
@VisibleForTesting
static double subtractMedianOfMedians(final ReadCountCollection readCounts, final Logger logger) {
    final RealMatrix counts = readCounts.counts();
    final Median medianCalculator = new Median();
    final double[] columnMedians = MatrixSummaryUtils.getColumnMedians(counts);

    final double medianOfMedians = medianCalculator.evaluate(columnMedians);
    counts.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(final int row, final int column, final double value) {
            return value - medianOfMedians;
        }
    });
    logger.info(String.format("Counts centered around the median of medians %.2f", medianOfMedians));
    return medianOfMedians;
}

Source File: GCBiasSimulatedData.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

6 votes

public static Pair<ReadCountCollection, double[]> simulatedData(final int numTargets, final int numSamples) {
    final List<Target> phonyTargets = SimulatedTargets.phonyTargets(numTargets);
    final List<String> phonySamples = SimulatedSamples.phonySamples(numSamples);

    final Random random = new Random(13);
    final double[] gcContentByTarget = IntStream.range(0, numTargets)
            .mapToDouble(n -> 0.5 + 0.2*random.nextGaussian())
            .map(x -> Math.min(x,0.95)).map(x -> Math.max(x,0.05)).toArray();
    final double[] gcBiasByTarget = Arrays.stream(gcContentByTarget).map(QUADRATIC_GC_BIAS_CURVE::apply).toArray();

    // model mainly GC bias with a small random amount of non-GC bias
    // thus noise after GC correction should be nearly zero
    final RealMatrix counts = new Array2DRowRealMatrix(numTargets, numSamples);
    counts.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(final int target, final int column, final double value) {
            return gcBiasByTarget[target]*(1.0 + 0.01*random.nextDouble());
        }
    });
    final ReadCountCollection rcc = new ReadCountCollection(phonyTargets, phonySamples, counts);
    return new ImmutablePair<>(rcc, gcContentByTarget);
}

Source File: WeightedIntDiGraph.java From pacaya with Apache License 2.0

6 votes

public static WeightedIntDiGraph fromMatrix(RealMatrix m) {
    WeightedIntDiGraph g = new WeightedIntDiGraph();
    m.walkInOptimizedOrder(new RealMatrixPreservingVisitor() {
        
        @Override
        public void visit(int row, int column, double value) {
            if (value != 0.0) {
                g.addEdge(row, column, value);
            }
        }
        
        @Override
        public void start(int rows, int columns, int startRow, int endRow, int startColumn, int endColumn) {
            // do nothing
        }
        
        @Override
        public double end() {
            // do nothing
            return 0;
        }
    });
    return g;
}

Source File: GCBiasCorrectorUnitTest.java From gatk with BSD 3-Clause "New" or "Revised" License

6 votes

private static Pair<RealMatrix, double[]> simulateData(final int numSamples,
                                                       final int numIntervals) {
    final Random random = new Random(RANDOM_SEED);
    final double[] intervalGCContent = IntStream.range(0, numIntervals)
            .mapToDouble(n -> 0.5 + 0.2 * random.nextGaussian())
            .map(x -> Math.min(x, 0.95)).map(x -> Math.max(x, 0.05))
            .toArray();
    final double[] intervalGCBias = Arrays.stream(intervalGCContent)
            .map(QUADRATIC_GC_BIAS_CURVE::apply)
            .toArray();

    //model GC bias along with a small random amount of uniform non-GC-bias noise;
    //remaining noise after GC-bias correction should be only arise from the latter
    final RealMatrix readCounts = new Array2DRowRealMatrix(numSamples, numIntervals);
    readCounts.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(final int sample, final int interval, final double value) {
            return (int) (MEAN_READ_DEPTH * intervalGCBias[interval] * (1. + NON_GC_BIAS_NOISE_LEVEL * random.nextDouble()));
        }
    });
    return new ImmutablePair<>(readCounts, intervalGCContent);
}

Source File: SomaticGenotypingEngine.java From gatk with BSD 3-Clause "New" or "Revised" License

5 votes

public static <EVIDENCE> RealMatrix getAsRealMatrix(final LikelihoodMatrix<EVIDENCE, Allele> matrix) {
    final RealMatrix result = new Array2DRowRealMatrix(matrix.numberOfAlleles(), matrix.evidenceCount());
    result.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(int row, int column, double value) {
            return matrix.get(row, column);
        }
    });
    return result;
}

Source File: ReadCountCollectionUtils.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

5 votes

/**
 * Impute zero counts to the median of non-zero values in the enclosing target row.
 *
 * <p>The imputation is done in-place, thus the input matrix is well be modified as a result of this call.</p>
 *
 * @param readCounts the input and output read-count matrix.
 */
public static void imputeZeroCountsAsTargetMedians(final ReadCountCollection readCounts, final Logger logger) {

    final RealMatrix counts = readCounts.counts();
    final int targetCount = counts.getRowDimension();

    final Median medianCalculator = new Median();
    int totalCounts = counts.getColumnDimension() * counts.getRowDimension();

    // Get the number of zeroes contained in the counts.
    final long totalZeroCounts = IntStream.range(0, targetCount)
            .mapToLong(t -> DoubleStream.of(counts.getRow(t))
                    .filter(c -> c == 0.0).count()).sum();

    // Get the median of each row, not including any entries that are zero.
    final double[] medians = IntStream.range(0, targetCount)
            .mapToDouble(t -> medianCalculator.evaluate(
                    DoubleStream.of(counts.getRow(t))
                            .filter(c -> c != 0.0)
                            .toArray()
            )).toArray();

    // Change any zeros in the counts to the median for the row.
    counts.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(final int row, final int column, final double value) {
            return value != 0 ? value : medians[row];
        }
    });

    if (totalZeroCounts > 0) {
        final double totalZeroCountsPercentage = 100.0 * (totalZeroCounts / totalCounts);
        logger.info(String.format("Some 0.0 counts (%d out of %d, %.2f%%) were imputed to their enclosing target " +
                "non-zero median", totalZeroCounts, totalZeroCounts, totalZeroCountsPercentage));
    } else {
        logger.info("No count is 0.0 thus no count needed to be imputed.");
    }
}

Source File: SVDDenoisingUtils.java From gatk with BSD 3-Clause "New" or "Revised" License

5 votes

private static void transformToFractionalCoverage(final RealMatrix matrix) {
    logger.info("Transforming read counts to fractional coverage...");
    final double[] sampleSums = IntStream.range(0, matrix.getRowDimension())
            .mapToDouble(r -> MathUtils.sum(matrix.getRow(r))).toArray();
    matrix.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(int sampleIndex, int intervalIndex, double value) {
            return value / sampleSums[sampleIndex];
        }
    });
}

Source File: SVDDenoisingUtils.java From gatk with BSD 3-Clause "New" or "Revised" License

5 votes

/**
 * Preprocess (i.e., transform to fractional coverage and correct GC bias)
 * and standardize read counts for samples when no panel of normals is available.
 * The original {@code readCounts} has dimensions 1 x intervals and is not modified.
 * If {@code intervalGCContent} is null, GC-bias correction will not be performed
 */
public static RealMatrix preprocessAndStandardizeSample(final double[] readCounts,
                                                        final double[] intervalGCContent) {
    Utils.nonNull(readCounts);
    Utils.validateArg(intervalGCContent == null || readCounts.length == intervalGCContent.length,
            "Number of intervals for read counts must match those for GC-content annotations.");

    RealMatrix result = new Array2DRowRealMatrix(new double[][]{readCounts});

    //preprocess (transform to fractional coverage, correct GC bias) copy in place
    logger.info("Preprocessing read counts...");
    transformToFractionalCoverage(result);
    performOptionalGCBiasCorrection(result, intervalGCContent);
    logger.info("Sample read counts preprocessed.");

    //standardize copy in place
    logger.info("Standardizing read counts...");
    divideBySampleMedianAndTransformToLog2(result);
    logger.info("Subtracting sample median...");
    final double[] sampleLog2Medians = MatrixSummaryUtils.getRowMedians(result);
    result.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() {
        @Override
        public double visit(int sampleIndex, int intervalIndex, double value) {
            return value - sampleLog2Medians[sampleIndex];
        }
    });
    logger.info("Sample read counts standardized.");

    return result;
}