use of org.broadinstitute.hellbender.tools.exome.ModeledSegment in project gatk-protected by broadinstitute.
the class PerformAlleleFractionSegmentation method doWork.
@Override
public Object doWork() {
final String sampleName = FilenameUtils.getBaseName(snpCountsFile.getAbsolutePath());
final AllelicPanelOfNormals allelicPoN = allelicPoNFile != null ? AllelicPanelOfNormals.read(allelicPoNFile) : AllelicPanelOfNormals.EMPTY_PON;
final AllelicCountCollection acc = new AllelicCountCollection(snpCountsFile);
final AlleleFractionSegmenter segmenter = new AlleleFractionSegmenter(initialNumStates, acc, allelicPoN);
final List<ModeledSegment> segments = segmenter.getModeledSegments();
SegmentUtils.writeModeledSegmentFile(outputSegmentsFile, segments, sampleName, true);
return "SUCCESS";
}
use of org.broadinstitute.hellbender.tools.exome.ModeledSegment in project gatk by broadinstitute.
the class CopyRatioSegmenterUnitTest method testSegmentation.
@Test
public void testSegmentation() {
final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));
final List<Double> trueWeights = Arrays.asList(0.2, 0.5, 0.3);
final List<Double> trueLog2CopyRatios = Arrays.asList(-2.0, 0.0, 1.4);
final double trueMemoryLength = 1e5;
final double trueStandardDeviation = 0.2;
final CopyRatioHMM trueModel = new CopyRatioHMM(trueLog2CopyRatios, trueWeights, trueMemoryLength, trueStandardDeviation);
final int chainLength = 10000;
final List<SimpleInterval> positions = randomPositions("chr1", chainLength, rng, trueMemoryLength / 4);
final List<Integer> trueStates = trueModel.generateHiddenStateChain(positions);
final List<Double> trueLog2CopyRatioSequence = trueStates.stream().map(n -> trueLog2CopyRatios.get(n)).collect(Collectors.toList());
final List<Double> data = trueLog2CopyRatioSequence.stream().map(cr -> generateData(trueStandardDeviation, cr, rng)).collect(Collectors.toList());
final List<Target> targets = positions.stream().map(Target::new).collect(Collectors.toList());
final ReadCountCollection rcc = new ReadCountCollection(targets, Arrays.asList("SAMPLE"), new Array2DRowRealMatrix(data.stream().mapToDouble(x -> x).toArray()));
final CopyRatioSegmenter segmenter = new CopyRatioSegmenter(10, rcc);
final List<ModeledSegment> segments = segmenter.getModeledSegments();
final double[] segmentCopyRatios = segments.stream().flatMap(s -> Collections.nCopies((int) s.getTargetCount(), s.getSegmentMeanInLog2CRSpace()).stream()).mapToDouble(x -> x).toArray();
final double averageCopyRatioError = IntStream.range(0, trueLog2CopyRatioSequence.size()).mapToDouble(n -> Math.abs(segmentCopyRatios[n] - trueLog2CopyRatioSequence.get(n))).average().getAsDouble();
Assert.assertEquals(averageCopyRatioError, 0, 0.025);
}
use of org.broadinstitute.hellbender.tools.exome.ModeledSegment in project gatk by broadinstitute.
the class PerformAlleleFractionSegmentationIntegrationTest method testCommandLine.
// checks that segmentation output is created -- only the unit test checks correctness of results
@Test
public void testCommandLine() {
final File snpFile = ALLELIC_COUNTS_FILE;
final File outputSegmentFile = createTempFile("segments", ".seg");
final int initialNumStates = 10;
final String[] arguments = { "-" + ExomeStandardArgumentDefinitions.TUMOR_ALLELIC_COUNTS_FILE_SHORT_NAME, snpFile.getAbsolutePath(), "-" + PerformAlleleFractionSegmentation.INITIAL_NUM_STATES_SHORT_NAME, Integer.toString(initialNumStates), "-" + ExomeStandardArgumentDefinitions.SEGMENT_FILE_SHORT_NAME, outputSegmentFile.getAbsolutePath() };
runCommandLine(arguments);
final List<ModeledSegment> segments = SegmentUtils.readModeledSegmentsFromSegmentFile(outputSegmentFile);
}
use of org.broadinstitute.hellbender.tools.exome.ModeledSegment in project gatk by broadinstitute.
the class PerformCopyRatioSegmentationIntegrationTest method testCommandLine.
// checks that segmentation output is created -- only the unit test checks correctness of results
@Test
public void testCommandLine() throws IOException {
final File tnCoverageFile = LOG2_TN_COVERAGE_FILE;
final File outputSegmentFile = createTempFile("segments", ".seg");
final int initialNumStates = 10;
final String[] arguments = { "-" + ExomeStandardArgumentDefinitions.TANGENT_NORMALIZED_COUNTS_FILE_SHORT_NAME, tnCoverageFile.getAbsolutePath(), "-" + PerformCopyRatioSegmentation.INITIAL_NUM_STATES_SHORT_NAME, Integer.toString(initialNumStates), "-" + ExomeStandardArgumentDefinitions.SEGMENT_FILE_SHORT_NAME, outputSegmentFile.getAbsolutePath() };
runCommandLine(arguments);
final List<ModeledSegment> segments = SegmentUtils.readModeledSegmentsFromSegmentFile(outputSegmentFile);
}
use of org.broadinstitute.hellbender.tools.exome.ModeledSegment in project gatk by broadinstitute.
the class AlleleFractionSegmenterUnitTest method testChromosomesOnDifferentSegments.
@Test
public void testChromosomesOnDifferentSegments() {
final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));
final double[] trueMinorAlleleFractions = new double[] { 0.12, 0.32, 0.5 };
final double trueMemoryLength = 1e5;
final AlleleFractionGlobalParameters trueParams = new AlleleFractionGlobalParameters(1.0, 0.01, 0.01);
// randomly set positions
final int chainLength = 100;
final List<SimpleInterval> positions = CopyRatioSegmenterUnitTest.randomPositions("chr1", chainLength, rng, trueMemoryLength / 4);
positions.addAll(CopyRatioSegmenterUnitTest.randomPositions("chr2", chainLength, rng, trueMemoryLength / 4));
positions.addAll(CopyRatioSegmenterUnitTest.randomPositions("chr3", chainLength, rng, trueMemoryLength / 4));
//fix everything to the same state 2
final int trueState = 2;
final List<Double> minorAlleleFractionSequence = Collections.nCopies(positions.size(), trueMinorAlleleFractions[trueState]);
final AllelicCountCollection counts = generateCounts(minorAlleleFractionSequence, positions, rng, trueParams);
final AlleleFractionSegmenter segmenter = new AlleleFractionSegmenter(10, counts, AllelicPanelOfNormals.EMPTY_PON);
final List<ModeledSegment> segments = segmenter.getModeledSegments();
//check that each chromosome has at least one segment
final int numDifferentContigsInSegments = (int) segments.stream().map(ModeledSegment::getContig).distinct().count();
Assert.assertEquals(numDifferentContigsInSegments, 3);
}
Aggregations