use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk by broadinstitute.
the class JointAFCRSegmenterUnitTest method testSegmentation.
@Test
public void testSegmentation() {
final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));
// probability that a datum is a het i.e. #hets / (#hets + #targets)
final double hetProportion = 0.25;
final List<Double> trueWeights = Arrays.asList(0.2, 0.5, 0.3);
final double[] trueMinorAlleleFractions = new double[] { 0.12, 0.32, 0.5 };
final double[] trueLog2CopyRatios = new double[] { -2.0, 0.0, 1.7 };
final List<AFCRHiddenState> trueJointStates = IntStream.range(0, trueLog2CopyRatios.length).mapToObj(n -> new AFCRHiddenState(trueMinorAlleleFractions[n], trueLog2CopyRatios[n])).collect(Collectors.toList());
final double trueMemoryLength = 1e5;
final double trueCauchyWidth = 0.2;
final int initialNumCRStates = 20;
final int initialNumAFStates = 20;
final AlleleFractionGlobalParameters trueAFParams = new AlleleFractionGlobalParameters(1.0, 0.01, 0.01);
final JointAFCRHMM trueJointModel = new JointAFCRHMM(trueJointStates, trueWeights, trueMemoryLength, trueAFParams, AllelicPanelOfNormals.EMPTY_PON, trueCauchyWidth);
// generate joint truth
final int chainLength = 10000;
final List<SimpleInterval> positions = CopyRatioSegmenterUnitTest.randomPositions("chr1", chainLength, rng, trueMemoryLength / 4);
final List<Integer> trueHiddenStates = trueJointModel.generateHiddenStateChain(positions);
final List<AFCRHiddenState> trueAFCRSequence = trueHiddenStates.stream().map(trueJointModel::getHiddenStateValue).collect(Collectors.toList());
final double[] trueCopyRatioSequence = trueAFCRSequence.stream().mapToDouble(AFCRHiddenState::getLog2CopyRatio).toArray();
final double[] trueAlleleFractionSequence = trueAFCRSequence.stream().mapToDouble(AFCRHiddenState::getMinorAlleleFraction).toArray();
// generate separate af and cr data
final GammaDistribution biasGenerator = AlleleFractionSegmenterUnitTest.getGammaDistribution(trueAFParams, rng);
final double outlierProbability = trueAFParams.getOutlierProbability();
final AllelicCountCollection afData = new AllelicCountCollection();
final List<Double> crData = new ArrayList<>();
final List<Target> crTargets = new ArrayList<>();
for (int n = 0; n < positions.size(); n++) {
final SimpleInterval position = positions.get(n);
final AFCRHiddenState jointState = trueAFCRSequence.get(n);
final double minorFraction = jointState.getMinorAlleleFraction();
final double log2CopyRatio = jointState.getLog2CopyRatio();
if (rng.nextDouble() < hetProportion) {
// het datum
afData.add(AlleleFractionSegmenterUnitTest.generateAllelicCount(minorFraction, position, rng, biasGenerator, outlierProbability));
} else {
//target datum
crTargets.add(new Target(position));
crData.add(CopyRatioSegmenterUnitTest.generateData(trueCauchyWidth, log2CopyRatio, rng));
}
}
final ReadCountCollection rcc = new ReadCountCollection(crTargets, Arrays.asList("SAMPLE"), new Array2DRowRealMatrix(crData.stream().mapToDouble(x -> x).toArray()));
final JointAFCRSegmenter segmenter = JointAFCRSegmenter.createJointSegmenter(initialNumCRStates, rcc, initialNumAFStates, afData, AllelicPanelOfNormals.EMPTY_PON);
final TargetCollection<SimpleInterval> tc = new HashedListTargetCollection<>(positions);
final List<Pair<SimpleInterval, AFCRHiddenState>> segmentation = segmenter.findSegments();
final List<ACNVModeledSegment> jointSegments = segmentation.stream().map(pair -> {
final SimpleInterval position = pair.getLeft();
final AFCRHiddenState jointState = pair.getRight();
final PosteriorSummary crSummary = PerformJointSegmentation.errorlessPosterior(jointState.getLog2CopyRatio());
final PosteriorSummary afSummary = PerformJointSegmentation.errorlessPosterior(jointState.getMinorAlleleFraction());
return new ACNVModeledSegment(position, crSummary, afSummary);
}).collect(Collectors.toList());
final double[] segmentCopyRatios = jointSegments.stream().flatMap(s -> Collections.nCopies(tc.targetCount(s.getInterval()), s.getSegmentMeanPosteriorSummary().getCenter()).stream()).mapToDouble(x -> x).toArray();
final double[] segmentMinorFractions = jointSegments.stream().flatMap(s -> Collections.nCopies(tc.targetCount(s.getInterval()), s.getMinorAlleleFractionPosteriorSummary().getCenter()).stream()).mapToDouble(x -> x).toArray();
final double averageMinorFractionError = Arrays.stream(MathArrays.ebeSubtract(trueAlleleFractionSequence, segmentMinorFractions)).map(Math::abs).average().getAsDouble();
final double averageCopyRatioError = Arrays.stream(MathArrays.ebeSubtract(trueCopyRatioSequence, segmentCopyRatios)).map(Math::abs).average().getAsDouble();
Assert.assertEquals(averageMinorFractionError, 0, 0.04);
Assert.assertEquals(averageCopyRatioError, 0, 0.04);
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk-protected by broadinstitute.
the class PreAdapterOrientationScorer method countOrientationBiasMetricsOverContext.
/**
* Gets PreAdapterQ collapsed over contexts.
*
* rows {PRO, CON} x cols {ref, alt}
* @param metrics metrics usually read from a picard preadapter detail file.
* @return mapping to score all orientation bias artifact modes to a PreAdapterQ score. This score can be used as a bam-file level score for level of artifacts.
*/
@VisibleForTesting
static Map<Transition, RealMatrix> countOrientationBiasMetricsOverContext(final List<SequencingArtifactMetrics.PreAdapterDetailMetrics> metrics) {
Utils.nonNull(metrics, "Input metrics cannot be null");
// Artifact mode to a matrix
final Map<Transition, RealMatrix> result = new HashMap<>();
// Collapse over context
for (SequencingArtifactMetrics.PreAdapterDetailMetrics metric : metrics) {
final Transition key = Transition.transitionOf(metric.REF_BASE, metric.ALT_BASE);
result.putIfAbsent(key, new Array2DRowRealMatrix(2, 2));
final RealMatrix preAdapterCountMatrix = result.get(key);
preAdapterCountMatrix.addToEntry(PreAdapterOrientationScorer.PRO, PreAdapterOrientationScorer.ALT, metric.PRO_ALT_BASES);
preAdapterCountMatrix.addToEntry(PreAdapterOrientationScorer.CON, PreAdapterOrientationScorer.ALT, metric.CON_ALT_BASES);
preAdapterCountMatrix.addToEntry(PreAdapterOrientationScorer.PRO, PreAdapterOrientationScorer.REF, metric.PRO_REF_BASES);
preAdapterCountMatrix.addToEntry(PreAdapterOrientationScorer.CON, PreAdapterOrientationScorer.REF, metric.CON_REF_BASES);
}
return result;
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk-protected by broadinstitute.
the class MatrixSummaryUtilsUnitTest method testGetRowAndColumnMedians.
@Test(dataProvider = "rowAndColumnMedianTestData")
public void testGetRowAndColumnMedians(double[][] testdata, double[] gtRowMedian, double[] gtColumnMedian) {
final RealMatrix tmp = new Array2DRowRealMatrix(testdata);
PoNTestUtils.assertEqualsDoubleArrays(MatrixSummaryUtils.getRowMedians(tmp), gtRowMedian);
PoNTestUtils.assertEqualsDoubleArrays(MatrixSummaryUtils.getColumnMedians(tmp), gtColumnMedian);
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk-protected by broadinstitute.
the class HDF5LibraryUnitTest method testCreateLargeMatrix.
@Test
public void testCreateLargeMatrix() {
// Creates a large PoN of junk values and simply tests that these can be written and read.
// Make a big, fake set of read counts.
final int numRows = 2500000;
final int numCols = 10;
final double mean = 3e-7;
final double sigma = 1e-9;
final RealMatrix bigCounts = createMatrixOfGaussianValues(numRows, numCols, mean, sigma);
final File tempOutputHD5 = IOUtils.createTempFile("big-ol-", ".hd5");
final HDF5File hdf5File = new HDF5File(tempOutputHD5, HDF5File.OpenMode.CREATE);
final String hdf5Path = "/test/m";
hdf5File.makeDoubleMatrix(hdf5Path, bigCounts.getData());
hdf5File.close();
final HDF5File hdf5FileForReading = new HDF5File(tempOutputHD5, HDF5File.OpenMode.READ_ONLY);
final double[][] result = hdf5FileForReading.readDoubleMatrix(hdf5Path);
final RealMatrix resultAsRealMatrix = new Array2DRowRealMatrix(result);
Assert.assertTrue(resultAsRealMatrix.getRowDimension() == numRows);
Assert.assertTrue(resultAsRealMatrix.getColumnDimension() == numCols);
final RealMatrix readMatrix = new Array2DRowRealMatrix(result);
PoNTestUtils.assertEqualsMatrix(readMatrix, bigCounts, false);
}
use of org.apache.commons.math3.linear.Array2DRowRealMatrix in project gatk-protected by broadinstitute.
the class GATKProtectedMathUtilsTest method testRowInf.
@Test
public void testRowInf() {
final double[][] array = { { 1, 2, Double.POSITIVE_INFINITY }, { 5, 5, 5 }, { 7, 8, 9 }, { -15, 2, 12 } };
final double[] guess = GATKProtectedMathUtils.rowSums(new Array2DRowRealMatrix(array));
final double[] gt = { Double.POSITIVE_INFINITY, 15, 24, -1 };
Assert.assertEquals(guess.length, 4);
Assert.assertEquals(guess[1], gt[1]);
Assert.assertEquals(guess[2], gt[2]);
Assert.assertEquals(guess[3], gt[3]);
Assert.assertTrue(Double.isInfinite(guess[0]));
}
Aggregations