Search in sources :

Example 1 with SingleReadImpl

use of com.milaboratory.core.io.sequence.SingleReadImpl in project mixcr by milaboratory.

the class VDJCAlignerPVFirstTest method test2.

@Test
@Ignore
public void test2() throws Exception {
    PairedRead read1 = new PairedRead(new SingleReadImpl(0, new NSequenceWithQuality(new NucleotideSequence("GCTGTGTATTACTGTGCAAGAGGGCCCCAAGAAAATAGTGGTTATTACTACGGGTTTGACTACTGGGGCCAGGGA"), SequenceQuality.GOOD_QUALITY_VALUE), "206"), new SingleReadImpl(0, new NSequenceWithQuality(new NucleotideSequence("GGCGCCAGGGGGAAGACCGATGGGCCCTTGGTGGAGGCTGAGGAGACGGTGACCAGGGTTCCCTGGCCCCAGTAG"), SequenceQuality.GOOD_QUALITY_VALUE), "206"));
    PairedRead read2 = new PairedRead(new SingleReadImpl(1, new NSequenceWithQuality(new NucleotideSequence("GCTGTGTATTACTGTGCAAGAGGGCCCCAAGAAAATAGTGGTTATTACTACGGGTTTGACTACTGGGGCCAGGGA"), SequenceQuality.GOOD_QUALITY_VALUE), "11621"), new SingleReadImpl(1, new NSequenceWithQuality(new NucleotideSequence("GGCGCCAGGGGGAAGACCGATGGGCCCTTGGTGGAGGCTGAGGAGACGGTGACCAGGGTTCCCTGGCCCCAGTAG"), SequenceQuality.GOOD_QUALITY_VALUE), "11621"));
    RunMiXCR.RunMiXCRAnalysis params = new RunMiXCR.RunMiXCRAnalysis(read1);
    RunMiXCR.AlignResult align = RunMiXCR.align(params);
    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();
    }
}
Also used : SingleReadImpl(com.milaboratory.core.io.sequence.SingleReadImpl) NSequenceWithQuality(com.milaboratory.core.sequence.NSequenceWithQuality) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) PairedRead(com.milaboratory.core.io.sequence.PairedRead) RunMiXCR(com.milaboratory.mixcr.util.RunMiXCR) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) Ignore(org.junit.Ignore) Test(org.junit.Test)

Example 2 with SingleReadImpl

use of com.milaboratory.core.io.sequence.SingleReadImpl in project mixcr by milaboratory.

the class VDJCAlignerSTest method test2.

@Test
@Ignore
public void test2() throws Exception {
    // @
    // GCTGTGTATTACTGTGCAAGAGGGCCCCAAGAAAATAGTGGTTATTACTACGGGTTTGACTACTGGGGCCAGGGAACCCTGGTCACCGTCTCCTCAGCCTCCACCAAGGGCCCATCGGTCTTCCCCCTGGCGCC
    // +
    // CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
    VDJCAlignerParameters parameters = VDJCParametersPresets.getByName("default");
    VDJCAlignerS aligner = new VDJCAlignerS(parameters);
    for (VDJCGene gene : VDJCLibraryRegistry.getDefault().getLibrary("default", "hs").getGenes(Chains.IGH)) if (parameters.containsRequiredFeature(gene))
        aligner.addGene(gene);
    SingleReadImpl read = new SingleReadImpl(0, new NSequenceWithQuality(new NucleotideSequence("GCTGTGTATTACTGTGCAAGAGGGCCCCAAGAAAATAGTGGTTATTACTACGGGTTTGACTACTGGGGCCAGGGAACCCTGGTCACCGTCTCCTCAGCCTCCACCAAGGGCCCATCGGTCTTCCCCCTGGCGCC"), SequenceQuality.GOOD_QUALITY_VALUE), "");
    RandomUtil.getThreadLocalRandom().setSeed(29);
    VDJCAlignmentResult<SingleRead> result = aligner.process0(read);
}
Also used : SingleReadImpl(com.milaboratory.core.io.sequence.SingleReadImpl) NSequenceWithQuality(com.milaboratory.core.sequence.NSequenceWithQuality) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) VDJCGene(io.repseq.core.VDJCGene) SingleRead(com.milaboratory.core.io.sequence.SingleRead) Ignore(org.junit.Ignore) Test(org.junit.Test)

Example 3 with SingleReadImpl

use of com.milaboratory.core.io.sequence.SingleReadImpl in project mixcr by milaboratory.

the class VDJCAlignerPVFirst method process0.

@Override
protected VDJCAlignmentResult<PairedRead> process0(final PairedRead input) {
    Target[] targets = getTargets(input);
    // Creates helper classes for each PTarget
    PAlignmentHelper[] helpers = createInitialHelpers(targets);
    VDJCAlignmentResult<PairedRead> result = parameters.getAllowPartialAlignments() ? processPartial(input, helpers) : processStrict(input, helpers);
    // if sAligner == null (which means --no-merge option), no merge will be performed
    if (result.alignment != null && sAligner != null) {
        final VDJCAlignments alignment = result.alignment;
        final TargetMerger.TargetMergingResult mergeResult = alignmentsMerger.merge(new AlignedTarget(alignment, 0), new AlignedTarget(alignment, 1), false);
        if (mergeResult.failedDueInconsistentAlignments()) {
            GeneType geneType = mergeResult.getFailedMergedGeneType();
            int removeId = alignment.getBestHit(geneType).getAlignment(0).getScore() > alignment.getBestHit(geneType).getAlignment(1).getScore() ? 1 : 0;
            if (listener != null)
                listener.onTopHitSequenceConflict(input, alignment, geneType);
            return new VDJCAlignmentResult<>(input, alignment.removeBestHitAlignment(geneType, removeId));
        } else if (mergeResult.isSuccessful()) {
            assert mergeResult.isUsingAlignments();
            NSequenceWithQuality alignedTarget = mergeResult.getResult().getTarget();
            SingleRead sRead = new SingleReadImpl(input.getId(), alignedTarget, "");
            VDJCAlignmentResult<SingleRead> sResult = sAligner.process0(sRead);
            if (sResult.alignment == null)
                return result;
            VDJCAlignments sAlignment = sResult.alignment.setHistory(new SequenceHistory[] { new SequenceHistory.Merge(SequenceHistory.OverlapType.AlignmentOverlap, result.alignment.getHistory(0), result.alignment.getHistory(1), mergeResult.getOffset(), mergeResult.getMismatched()) }, new SequenceRead[] { input });
            if (listener != null)
                listener.onSuccessfulAlignmentOverlap(input, sAlignment);
            return new VDJCAlignmentResult<>(input, sAlignment);
        }
    }
    return result;
}
Also used : SingleReadImpl(com.milaboratory.core.io.sequence.SingleReadImpl) TargetMerger(com.milaboratory.mixcr.partialassembler.TargetMerger) SequenceHistory(com.milaboratory.mixcr.basictypes.SequenceHistory) PairedRead(com.milaboratory.core.io.sequence.PairedRead) AlignedTarget(com.milaboratory.mixcr.partialassembler.AlignedTarget) Target(com.milaboratory.core.Target) AlignedTarget(com.milaboratory.mixcr.partialassembler.AlignedTarget) NSequenceWithQuality(com.milaboratory.core.sequence.NSequenceWithQuality) SequenceRead(com.milaboratory.core.io.sequence.SequenceRead) SingleRead(com.milaboratory.core.io.sequence.SingleRead) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments)

Example 4 with SingleReadImpl

use of com.milaboratory.core.io.sequence.SingleReadImpl 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 5 with SingleReadImpl

use of com.milaboratory.core.io.sequence.SingleReadImpl in project mixcr by milaboratory.

the class FullSeqAssemblerTest method test1.

@Test
public void test1() throws Exception {
    int len = 140;
    PairedRead read1 = new PairedRead(new SingleReadImpl(0, new NSequenceWithQuality(masterSeq1WT.getRangeFromCDR3Begin(-20, len)), "R1"), new SingleReadImpl(0, new NSequenceWithQuality(masterSeq1WT.getRangeFromCDR3Begin(-200, len).getReverseComplement()), "R2"));
    PairedRead read2 = new PairedRead(new SingleReadImpl(1, new NSequenceWithQuality(masterSeq1WT.getRangeFromCDR3Begin(-30, len)), "R1"), new SingleReadImpl(1, new NSequenceWithQuality(masterSeq1WT.getRangeFromCDR3Begin(-150, len).getReverseComplement()), "R2"));
    RunMiXCR.RunMiXCRAnalysis params = new RunMiXCR.RunMiXCRAnalysis(read1, read2);
    // [-200, -60]  [-20, 120]
    // [-150, 110]
    // 
    // [-200, -150], [110, 120] = 60
    // [-60, -20] = 40
    params.alignerParameters = VDJCParametersPresets.getByName("rna-seq");
    params.alignerParameters.setSaveOriginalReads(true);
    params.cloneAssemblerParameters.updateFrom(params.alignerParameters);
    RunMiXCR.AlignResult align = RunMiXCR.align(params);
    // for (VDJCAlignments al : align.alignments) {
    // for (int i = 0; i < al.numberOfTargets(); i++)
    // System.out.println(VDJCAlignmentsFormatter.getTargetAsMultiAlignment(al, i));
    // System.out.println();
    // }
    RunMiXCR.AssembleResult assemble = RunMiXCR.assemble(align);
    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);
    PointSequence[] r2s = agg.toPointSequences(align.alignments.get(1));
    TIntHashSet p2 = new TIntHashSet(Arrays.stream(r2s).mapToInt(s -> s.point).toArray());
    Assert.assertEquals(261 - masterSeq1WT.cdr3Part, p2.size());
    PointSequence[] r1s = agg.toPointSequences(align.alignments.get(0));
    TIntHashSet p1 = new TIntHashSet(Arrays.stream(r1s).mapToInt(s -> s.point).toArray());
    Assert.assertEquals(281 - masterSeq1WT.cdr3Part, p1.size());
    FullSeqAssembler.RawVariantsData prep = agg.calculateRawData(() -> CUtils.asOutputPort(align.alignments));
    long uniq1 = StreamSupport.stream(CUtils.it(prep.createPort()).spliterator(), false).mapToInt(l -> l[0]).filter(c -> c == 0xFFFFFFFF).count();
    long uniq2 = StreamSupport.stream(CUtils.it(prep.createPort()).spliterator(), false).mapToInt(l -> l[1]).filter(c -> c == 0xFFFFFFFF).count();
    Assert.assertEquals(40, uniq1);
    Assert.assertEquals(60, uniq2);
    for (Clone clone : new CloneSet(Arrays.asList(agg.callVariants(prep))).getClones()) {
        ActionExportClonesPretty.outputCompact(System.out, clone);
        System.out.println();
        System.out.println(" ================================================ ");
        System.out.println();
    }
}
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) TIntHashSet(gnu.trove.set.hash.TIntHashSet) CloneSet(com.milaboratory.mixcr.basictypes.CloneSet) NSequenceWithQuality(com.milaboratory.core.sequence.NSequenceWithQuality) RunMiXCR(com.milaboratory.mixcr.util.RunMiXCR) CloneFactory(com.milaboratory.mixcr.assembler.CloneFactory) Clone(com.milaboratory.mixcr.basictypes.Clone) Test(org.junit.Test)

Aggregations

SingleReadImpl (com.milaboratory.core.io.sequence.SingleReadImpl)7 NSequenceWithQuality (com.milaboratory.core.sequence.NSequenceWithQuality)6 VDJCAlignments (com.milaboratory.mixcr.basictypes.VDJCAlignments)6 PairedRead (com.milaboratory.core.io.sequence.PairedRead)5 SequenceRead (com.milaboratory.core.io.sequence.SequenceRead)5 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)5 Test (org.junit.Test)5 RunMiXCR (com.milaboratory.mixcr.util.RunMiXCR)4 CUtils (cc.redberry.pipe.CUtils)3 SequenceReaderCloseable (com.milaboratory.core.io.sequence.SequenceReaderCloseable)3 SequenceQuality (com.milaboratory.core.sequence.SequenceQuality)3 CloneFactory (com.milaboratory.mixcr.assembler.CloneFactory)3 Clone (com.milaboratory.mixcr.basictypes.Clone)3 CloneSet (com.milaboratory.mixcr.basictypes.CloneSet)3 VDJCAlignmentsFormatter (com.milaboratory.mixcr.basictypes.VDJCAlignmentsFormatter)3 ActionExportClonesPretty (com.milaboratory.mixcr.cli.ActionExportClonesPretty)3 Main (com.milaboratory.mixcr.cli.Main)3 VDJCParametersPresets (com.milaboratory.mixcr.vdjaligners.VDJCParametersPresets)3 TIntHashSet (gnu.trove.set.hash.TIntHashSet)3 GeneFeature (io.repseq.core.GeneFeature)3