Search in sources :

Example 76 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project mixcr by milaboratory.

the class FullSeqAssemblerTest method testRandom1.

@Test
public void testRandom1() throws Exception {
    CloneFraction[] clones = { new CloneFraction(750, masterSeq1WT), // V: S346:G->T
    new CloneFraction(1000, masterSeq1VSub1), // J: D55:A
    new CloneFraction(1000, masterSeq1VDel1JDel1), // J: D62:C
    new CloneFraction(500, masterSeq1VDel1JDelVSub2) };
    Well19937c rand = new Well19937c();
    rand.setSeed(12345);
    RandomDataGenerator rdg = new RandomDataGenerator(rand);
    List<SequenceRead> readsOrig = new ArrayList<>();
    int readLength = 100;
    int id = -1;
    for (CloneFraction clone : clones) {
        for (int i = 0; i < clone.count; i++) {
            // Left read with CDR3
            ++id;
            readsOrig.add(new PairedRead(new SingleReadImpl(id, new NSequenceWithQuality(clone.seq.getRangeFromCDR3Begin(-rand.nextInt(readLength - clone.seq.cdr3Part), readLength)), "R1_" + id), new SingleReadImpl(id, new NSequenceWithQuality(clone.seq.getRangeFromCDR3End(rdg.nextInt(-clone.seq.cdr3Part / 2, clone.seq.jPart), readLength).getReverseComplement()), "R2_" + id)));
            ++id;
            readsOrig.add(new PairedRead(new SingleReadImpl(id, new NSequenceWithQuality(clone.seq.getRangeFromCDR3Begin(rdg.nextInt(-clone.seq.vPart, clone.seq.cdr3Part / 2 - readLength), readLength)), "R1_" + id), new SingleReadImpl(id, new NSequenceWithQuality(clone.seq.getRangeFromCDR3Begin(-rand.nextInt(readLength - clone.seq.cdr3Part), readLength)).getReverseComplement(), "R2_" + id)));
        }
    }
    // readsOrig = Arrays.asList(setReadId(0, readsOrig.get(12)), setReadId(1, readsOrig.get(13)));
    int[] perm = rdg.nextPermutation(readsOrig.size(), readsOrig.size());
    List<SequenceRead> reads = new ArrayList<>();
    for (int i = 0; i < readsOrig.size(); i++) reads.add(readsOrig.get(perm[i]));
    RunMiXCR.RunMiXCRAnalysis params = new RunMiXCR.RunMiXCRAnalysis(new SequenceReaderCloseable<SequenceRead>() {

        int counter = 0;

        @Override
        public void close() {
        }

        @Override
        public long getNumberOfReads() {
            return counter;
        }

        @Override
        public synchronized SequenceRead take() {
            if (counter == reads.size())
                return null;
            return reads.get(counter++);
        }
    }, true);
    params.alignerParameters = VDJCParametersPresets.getByName("rna-seq");
    params.alignerParameters.setSaveOriginalReads(true);
    params.alignerParameters.setVAlignmentParameters(params.alignerParameters.getVAlignerParameters().setGeneFeatureToAlign(GeneFeature.VTranscriptWithP));
    RunMiXCR.AlignResult align = RunMiXCR.align(params);
    // // TODO exception for translation
    // for (VDJCAlignments al : align.alignments) {
    // for (int i = 0; i < al.numberOfTargets(); i++) {
    // System.out.println(VDJCAlignmentsFormatter.getTargetAsMultiAlignment(al, i));
    // System.out.println();
    // }
    // System.out.println();
    // System.out.println(" ================================================ ");
    // System.out.println();
    // }
    RunMiXCR.AssembleResult assemble = RunMiXCR.assemble(align);
    Assert.assertEquals(1, assemble.cloneSet.size());
    CloneFactory cloneFactory = new CloneFactory(align.parameters.cloneAssemblerParameters.getCloneFactoryParameters(), align.parameters.cloneAssemblerParameters.getAssemblingFeatures(), align.usedGenes, align.parameters.alignerParameters.getFeaturesToAlignMap());
    FullSeqAssembler agg = new FullSeqAssembler(cloneFactory, DEFAULT_PARAMETERS, assemble.cloneSet.get(0), align.parameters.alignerParameters);
    FullSeqAssembler.RawVariantsData prep = agg.calculateRawData(() -> CUtils.asOutputPort(align.alignments.stream().filter(a -> a.getFeature(GeneFeature.CDR3) != null).collect(Collectors.toList())));
    List<Clone> clns = new ArrayList<>(new CloneSet(Arrays.asList(agg.callVariants(prep))).getClones());
    clns.sort(Comparator.comparingDouble(Clone::getCount).reversed());
    System.out.println("# Clones: " + clns.size());
    id = 0;
    for (Clone clone : clns) {
        clone = clone.setId(id++);
        System.out.println(clone.numberOfTargets());
        System.out.println(clone.getCount());
        System.out.println(clone.getFraction());
        System.out.println(clone.getBestHit(GeneType.Variable).getAlignment(0).getAbsoluteMutations());
        System.out.println(clone.getBestHit(GeneType.Joining).getAlignment(0).getAbsoluteMutations());
        System.out.println();
    // ActionExportClonesPretty.outputCompact(System.out, clone);
    }
}
Also used : java.util(java.util) SequenceRead(com.milaboratory.core.io.sequence.SequenceRead) Well44497b(org.apache.commons.math3.random.Well44497b) SequenceQuality(com.milaboratory.core.sequence.SequenceQuality) Clone(com.milaboratory.mixcr.basictypes.Clone) GeneFeature(io.repseq.core.GeneFeature) CloneFactory(com.milaboratory.mixcr.assembler.CloneFactory) Main(com.milaboratory.mixcr.cli.Main) StreamSupport(java.util.stream.StreamSupport) PairedRead(com.milaboratory.core.io.sequence.PairedRead) RunMiXCR(com.milaboratory.mixcr.util.RunMiXCR) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) VDJCAlignmentsFormatter(com.milaboratory.mixcr.basictypes.VDJCAlignmentsFormatter) CUtils(cc.redberry.pipe.CUtils) Test(org.junit.Test) Collectors(java.util.stream.Collectors) TIntHashSet(gnu.trove.set.hash.TIntHashSet) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) Well19937c(org.apache.commons.math3.random.Well19937c) CloneSet(com.milaboratory.mixcr.basictypes.CloneSet) SingleReadImpl(com.milaboratory.core.io.sequence.SingleReadImpl) GeneType(io.repseq.core.GeneType) NSequenceWithQuality(com.milaboratory.core.sequence.NSequenceWithQuality) ActionExportClonesPretty(com.milaboratory.mixcr.cli.ActionExportClonesPretty) VDJCParametersPresets(com.milaboratory.mixcr.vdjaligners.VDJCParametersPresets) Assert(org.junit.Assert) SequenceReaderCloseable(com.milaboratory.core.io.sequence.SequenceReaderCloseable) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) SingleReadImpl(com.milaboratory.core.io.sequence.SingleReadImpl) PairedRead(com.milaboratory.core.io.sequence.PairedRead) Well19937c(org.apache.commons.math3.random.Well19937c) CloneSet(com.milaboratory.mixcr.basictypes.CloneSet) NSequenceWithQuality(com.milaboratory.core.sequence.NSequenceWithQuality) Clone(com.milaboratory.mixcr.basictypes.Clone) RandomDataGenerator(org.apache.commons.math3.random.RandomDataGenerator) SequenceRead(com.milaboratory.core.io.sequence.SequenceRead) RunMiXCR(com.milaboratory.mixcr.util.RunMiXCR) CloneFactory(com.milaboratory.mixcr.assembler.CloneFactory) Test(org.junit.Test)

Example 77 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project accumulo by apache.

the class RolllingStatsTest method testZipf.

@Test
public void testZipf() {
    ZipfDistribution zd = new ZipfDistribution(new Well19937c(42), 1000, 2);
    testDistribrution(() -> zd.sample() * 100);
}
Also used : Well19937c(org.apache.commons.math3.random.Well19937c) ZipfDistribution(org.apache.commons.math3.distribution.ZipfDistribution) Test(org.junit.Test)

Example 78 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project ignite by apache.

the class BootstrappedDatasetBuilder method build.

/**
 * {@inheritDoc}
 */
@Override
public BootstrappedDatasetPartition build(LearningEnvironment env, Iterator<UpstreamEntry<K, V>> upstreamData, long upstreamDataSize, EmptyContext ctx) {
    BootstrappedVector[] dataset = new BootstrappedVector[Math.toIntExact(upstreamDataSize)];
    int cntr = 0;
    PoissonDistribution poissonDistribution = new PoissonDistribution(new Well19937c(env.randomNumbersGenerator().nextLong()), subsampleSize, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    while (upstreamData.hasNext()) {
        UpstreamEntry<K, V> nextRow = upstreamData.next();
        LabeledVector<Double> vecAndLb = preprocessor.apply(nextRow.getKey(), nextRow.getValue());
        Vector features = vecAndLb.features();
        Double lb = vecAndLb.label();
        int[] repetitionCounters = new int[samplesCnt];
        Arrays.setAll(repetitionCounters, i -> poissonDistribution.sample());
        dataset[cntr++] = new BootstrappedVector(features, lb, repetitionCounters);
    }
    return new BootstrappedDatasetPartition(dataset);
}
Also used : PoissonDistribution(org.apache.commons.math3.distribution.PoissonDistribution) Well19937c(org.apache.commons.math3.random.Well19937c) LabeledVector(org.apache.ignite.ml.structures.LabeledVector) Vector(org.apache.ignite.ml.math.primitives.vector.Vector)

Example 79 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project repseqio by repseqio.

the class ExportCloneSequencesAction method go.

@Override
public void go(ActionHelper helper) throws Exception {
    Chains chains = params.getChains();
    GeneFeature geneFeature = params.getGeneFeature();
    RandomGenerator random = new Well19937c(1232434);
    try (GRepertoireReader input = new GRepertoireReader(createBufferedReader(params.getInput()));
        FastaWriter<NucleotideSequence> output = createSingleFastaWriter(params.getOutput())) {
        List<DescriptionExtractor> extractors = params.getExtractors(input.getLibrary());
        long i = 0;
        for (GClone clone : CUtils.it(input)) {
            int f = params.factor == null ? 1 : randomizedRound(clone.abundance * params.factor, random);
            for (int j = 0; j < f; j++) for (Map.Entry<String, GGene> e : clone.genes.entrySet()) if (chains.contains(e.getKey())) {
                StringBuilder descriptionLine = new StringBuilder("GClone");
                for (DescriptionExtractor extractor : extractors) descriptionLine.append("|").append(extractor.extract(clone, e.getValue(), e.getKey()));
                output.write(new FastaRecord<>(i++, descriptionLine.toString(), e.getValue().getFeature(geneFeature)));
            }
        }
    }
}
Also used : GeneFeature(io.repseq.core.GeneFeature) Chains(io.repseq.core.Chains) Well19937c(org.apache.commons.math3.random.Well19937c) RandomGenerator(org.apache.commons.math3.random.RandomGenerator) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence)

Example 80 with Well19937c

use of org.apache.commons.math3.random.Well19937c in project repseqio by repseqio.

the class BasicGCloneModelTest method test1.

@Test
public void test1() throws Exception {
    GGeneModel geneModel = new BasicGGeneModel(new IndependentVDJCGenesModel(b("TRBV12-2*00", 1.0).put("TRBV12-3*00", 0.0).get(), b("TRBD1*00", 1.0).put("TRBD2*00", 0.0).get(), b("TRBJ1-2*00", 1.0).put("TRBJ1-3*00", 0.0).get(), b("TRBC1*00", 1.0).put("TRBC2*00", 0.0).get()), new IndependentVDJTrimmingModel(new CommonCategoricalGeneTrimmingModel(b(-1, 1.0).get()), new CommonCategoricalDGeneTrimmingModel(b("-2|-3", 1.0).get()), new CommonCategoricalGeneTrimmingModel(b(-4, 1.0).get())), new FixedInsertModel(new NucleotideSequence("ATTA")), new FixedInsertModel(new NucleotideSequence("GACA")));
    BasicGCloneModel model = new BasicGCloneModel(new VDJCLibraryId("default", 9606), new FixedRealModel(1.0), b("TRB", geneModel).get());
    TestUtil.assertJson(model);
    GCloneGenerator gen = model.create(new Well19937c(123), VDJCLibraryRegistry.getDefault());
    GClone clone = gen.sample();
    TestUtil.assertJson(clone);
    GGene trb = clone.genes.get("TRB");
    assertEquals(trb.vdjcGenes.v.getName(), "TRBV12-2*00");
    assertEquals(trb.vdjcGenes.d.getName(), "TRBD1*00");
    assertEquals(trb.vdjcGenes.j.getName(), "TRBJ1-2*00");
    assertEquals(trb.vdjcGenes.c.getName(), "TRBC1*00");
    assertEquals(new VDJTrimming(-1, -4, -2, -3), trb.vdjTrimming);
    assertEquals(new NucleotideSequence("ATTA"), trb.vInsert);
    assertEquals(new NucleotideSequence("GACA"), trb.djInsert);
}
Also used : VDJCLibraryId(io.repseq.core.VDJCLibraryId) Well19937c(org.apache.commons.math3.random.Well19937c) VDJTrimming(io.repseq.gen.VDJTrimming) GGene(io.repseq.gen.GGene) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) GClone(io.repseq.gen.GClone) Test(org.junit.Test)

Aggregations

Well19937c (org.apache.commons.math3.random.Well19937c)81 RandomDataGenerator (org.apache.commons.math3.random.RandomDataGenerator)42 ArrayList (java.util.ArrayList)31 RandomGenerator (org.apache.commons.math3.random.RandomGenerator)28 Test (org.junit.Test)23 FakeGradientFunction (gdsc.smlm.function.FakeGradientFunction)17 DoubleEquality (gdsc.core.utils.DoubleEquality)7 Gaussian2DFunction (gdsc.smlm.function.gaussian.Gaussian2DFunction)6 DenseMatrix64F (org.ejml.data.DenseMatrix64F)6 TimingService (gdsc.core.test.TimingService)4 Statistics (gdsc.core.utils.Statistics)4 StoredDataStatistics (gdsc.core.utils.StoredDataStatistics)4 Gradient1Function (gdsc.smlm.function.Gradient1Function)4 ValueProcedure (gdsc.smlm.function.ValueProcedure)4 ErfGaussian2DFunction (gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction)4 TooManyEvaluationsException (org.apache.commons.math3.exception.TooManyEvaluationsException)4 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)3 PseudoRandomGenerator (gdsc.core.utils.PseudoRandomGenerator)3 PrecomputedGradient1Function (gdsc.smlm.function.PrecomputedGradient1Function)3 SingleFreeCircularGaussian2DFunction (gdsc.smlm.function.gaussian.SingleFreeCircularGaussian2DFunction)3