Search in sources :

Example 1 with VDJCGene

use of io.repseq.core.VDJCGene in project mixcr by milaboratory.

the class ActionAssembleContigs method go.

@Override
public void go(ActionHelper helper) throws Exception {
    // TODO FIX!!!!!!!!!!!!!
    if (parameters.threads > 1)
        throw new ParameterException("Multithreaded processing is not supported yet.");
    long beginTimestamp = System.currentTimeMillis();
    FullSeqAssemblerParameters p = FullSeqAssemblerParameters.getByName("default");
    if (!parameters.overrides.isEmpty()) {
        // Perform parameters overriding
        p = JsonOverrider.override(p, FullSeqAssemblerParameters.class, parameters.overrides);
        if (p == null)
            throw new ProcessException("Failed to override some parameter.");
    }
    final FullSeqAssemblerReport report = new FullSeqAssemblerReport();
    FullSeqAssemblerParameters assemblerParameters = p;
    int totalClonesCount = 0;
    List<VDJCGene> genes;
    VDJCAlignerParameters alignerParameters;
    CloneAssemblerParameters cloneAssemblerParameters;
    try (ClnAReader reader = new ClnAReader(parameters.getInputFileName(), VDJCLibraryRegistry.getDefault());
        PrimitivO tmpOut = new PrimitivO(new BufferedOutputStream(new FileOutputStream(parameters.getOutputFileName())))) {
        final CloneFactory cloneFactory = new CloneFactory(reader.getAssemblerParameters().getCloneFactoryParameters(), reader.getAssemblingFeatures(), reader.getGenes(), reader.getAlignerParameters().getFeaturesToAlignMap());
        alignerParameters = reader.getAlignerParameters();
        cloneAssemblerParameters = reader.getAssemblerParameters();
        genes = reader.getGenes();
        IOUtil.registerGeneReferences(tmpOut, genes, alignerParameters);
        ClnAReader.CloneAlignmentsPort cloneAlignmentsPort = reader.clonesAndAlignments();
        SmartProgressReporter.startProgressReport("Assembling", cloneAlignmentsPort);
        OutputPort<Clone[]> parallelProcessor = new ParallelProcessor<>(cloneAlignmentsPort, cloneAlignments -> {
            FullSeqAssembler fullSeqAssembler = new FullSeqAssembler(cloneFactory, assemblerParameters, cloneAlignments.clone, alignerParameters);
            fullSeqAssembler.setReport(report);
            FullSeqAssembler.RawVariantsData rawVariantsData = fullSeqAssembler.calculateRawData(() -> {
                try {
                    return cloneAlignments.alignments();
                } catch (IOException e) {
                    throw new RuntimeException(e);
                }
            });
            return fullSeqAssembler.callVariants(rawVariantsData);
        }, parameters.threads);
        for (Clone[] clones : CUtils.it(parallelProcessor)) {
            totalClonesCount += clones.length;
            for (Clone cl : clones) tmpOut.writeObject(cl);
        }
        assert report.getInitialCloneCount() == reader.numberOfClones();
    }
    assert report.getFinalCloneCount() == totalClonesCount;
    assert report.getFinalCloneCount() >= report.getInitialCloneCount();
    Clone[] clones = new Clone[totalClonesCount];
    try (PrimitivI tmpIn = new PrimitivI(new BufferedInputStream(new FileInputStream(parameters.getOutputFileName())))) {
        IOUtil.registerGeneReferences(tmpIn, genes, alignerParameters);
        int i = 0;
        for (Clone clone : CUtils.it(new PipeDataInputReader<>(Clone.class, tmpIn, totalClonesCount))) clones[i++] = clone;
    }
    Arrays.sort(clones, Comparator.comparingDouble(c -> -c.getCount()));
    for (int i = 0; i < clones.length; i++) clones[i] = clones[i].setId(i);
    CloneSet cloneSet = new CloneSet(Arrays.asList(clones), genes, alignerParameters.getFeaturesToAlignMap(), alignerParameters, cloneAssemblerParameters);
    try (CloneSetIO.CloneSetWriter writer = new CloneSetIO.CloneSetWriter(cloneSet, parameters.getOutputFileName())) {
        SmartProgressReporter.startProgressReport(writer);
        writer.write();
    }
    ReportWrapper reportWrapper = new ReportWrapper(command(), report);
    reportWrapper.setStartMillis(beginTimestamp);
    reportWrapper.setInputFiles(parameters.getInputFileName());
    reportWrapper.setOutputFiles(parameters.getOutputFileName());
    reportWrapper.setCommandLine(helper.getCommandLineArguments());
    reportWrapper.setFinishMillis(System.currentTimeMillis());
    // Writing report to stout
    System.out.println("============= Report ==============");
    Util.writeReportToStdout(report);
    if (parameters.report != null)
        Util.writeReport(parameters.report, reportWrapper);
    if (parameters.jsonReport != null)
        Util.writeJsonReport(parameters.jsonReport, reportWrapper);
}
Also used : Parameters(com.beust.jcommander.Parameters) java.util(java.util) ParameterException(com.beust.jcommander.ParameterException) Action(com.milaboratory.cli.Action) Parameter(com.beust.jcommander.Parameter) CloneAssemblerParameters(com.milaboratory.mixcr.assembler.CloneAssemblerParameters) com.milaboratory.mixcr.basictypes(com.milaboratory.mixcr.basictypes) FullSeqAssemblerParameters(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssemblerParameters) ParallelProcessor(cc.redberry.pipe.blocks.ParallelProcessor) VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) PositiveInteger(com.beust.jcommander.validators.PositiveInteger) CloneFactory(com.milaboratory.mixcr.assembler.CloneFactory) FullSeqAssembler(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssembler) ActionHelper(com.milaboratory.cli.ActionHelper) PrimitivO(com.milaboratory.primitivio.PrimitivO) CUtils(cc.redberry.pipe.CUtils) FullSeqAssemblerReport(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssemblerReport) PrimitivI(com.milaboratory.primitivio.PrimitivI) OutputPort(cc.redberry.pipe.OutputPort) PipeDataInputReader(com.milaboratory.primitivio.PipeDataInputReader) java.io(java.io) VDJCGene(io.repseq.core.VDJCGene) SmartProgressReporter(com.milaboratory.util.SmartProgressReporter) DynamicParameter(com.beust.jcommander.DynamicParameter) VDJCLibraryRegistry(io.repseq.core.VDJCLibraryRegistry) ProcessException(com.milaboratory.cli.ProcessException) ActionParametersWithOutput(com.milaboratory.cli.ActionParametersWithOutput) VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) FullSeqAssemblerParameters(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssemblerParameters) CloneAssemblerParameters(com.milaboratory.mixcr.assembler.CloneAssemblerParameters) ParallelProcessor(cc.redberry.pipe.blocks.ParallelProcessor) ParameterException(com.beust.jcommander.ParameterException) PrimitivO(com.milaboratory.primitivio.PrimitivO) FullSeqAssembler(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssembler) FullSeqAssemblerReport(com.milaboratory.mixcr.assembler.fullseq.FullSeqAssemblerReport) ProcessException(com.milaboratory.cli.ProcessException) VDJCGene(io.repseq.core.VDJCGene) CloneFactory(com.milaboratory.mixcr.assembler.CloneFactory) PrimitivI(com.milaboratory.primitivio.PrimitivI)

Example 2 with VDJCGene

use of io.repseq.core.VDJCGene in project mixcr by milaboratory.

the class IOUtil method registerGeneReferences.

public static void registerGeneReferences(PrimitivI input, List<VDJCGene> genes, HasFeatureToAlign featuresToAlign) {
    // Putting genes references and feature sequences to be serialized/deserialized as references
    for (VDJCGene gene : genes) {
        input.putKnownReference(gene);
        // Also put sequences of certain gene features of genes as known references if required
        if (featuresToAlign != null) {
            GeneFeature featureToAlign = featuresToAlign.getFeatureToAlign(gene.getGeneType());
            if (featureToAlign == null)
                continue;
            NucleotideSequence featureSequence = gene.getFeature(featureToAlign);
            if (featureSequence == null)
                continue;
            input.putKnownReference(featureSequence);
        }
    }
}
Also used : GeneFeature(io.repseq.core.GeneFeature) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) VDJCGene(io.repseq.core.VDJCGene)

Example 3 with VDJCGene

use of io.repseq.core.VDJCGene in project mixcr by milaboratory.

the class IOUtil method readGeneReferences.

public static List<VDJCGene> readGeneReferences(PrimitivI input, VDJCLibraryRegistry registry) {
    // Reading gene ids
    int count = input.readInt();
    List<VDJCGene> genes = new ArrayList<>(count);
    for (int i = 0; i < count; i++) {
        VDJCGeneId id = input.readObject(VDJCGeneId.class);
        VDJCGene gene = registry.getGene(id);
        if (gene == null)
            throw new RuntimeException("Gene not found: " + id);
        genes.add(gene);
    }
    return genes;
}
Also used : ArrayList(java.util.ArrayList) VDJCGene(io.repseq.core.VDJCGene) VDJCGeneId(io.repseq.core.VDJCGeneId)

Example 4 with VDJCGene

use of io.repseq.core.VDJCGene 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 5 with VDJCGene

use of io.repseq.core.VDJCGene in project mixcr by milaboratory.

the class VDJCAlignerSTest method testSerialization1.

@Test
public void testSerialization1() throws Exception {
    VDJCAlignerParameters parameters = VDJCParametersPresets.getByName("default");
    // LociLibrary ll = LociLibraryManager.getDefault().getLibrary("mi");
    ByteArrayOutputStream bos = new ByteArrayOutputStream();
    List<VDJCAlignments> alignemntsList = new ArrayList<>();
    int header;
    try (SingleFastqReader reader = new SingleFastqReader(VDJCAlignerSTest.class.getClassLoader().getResourceAsStream("sequences/sample_IGH_R1.fastq"), true)) {
        VDJCAlignerS aligner = new VDJCAlignerS(parameters);
        for (VDJCGene gene : VDJCLibraryRegistry.getDefault().getLibrary("default", "hs").getGenes(Chains.IGH)) if (parameters.containsRequiredFeature(gene))
            aligner.addGene(gene);
        try (VDJCAlignmentsWriter writer = new VDJCAlignmentsWriter(bos)) {
            writer.header(aligner);
            header = bos.size();
            for (SingleRead read : CUtils.it(reader)) {
                VDJCAlignmentResult<SingleRead> result = aligner.process(read);
                if (result.alignment != null) {
                    writer.write(result.alignment);
                    alignemntsList.add(result.alignment);
                }
            }
        }
    }
    Assert.assertTrue(alignemntsList.size() > 10);
    System.out.println("Bytes per alignment: " + (bos.size() - header) / alignemntsList.size());
    try (VDJCAlignmentsReader reader = new VDJCAlignmentsReader(new ByteArrayInputStream(bos.toByteArray()))) {
        int i = 0;
        for (VDJCAlignments alignments : CUtils.it(reader)) Assert.assertEquals(alignemntsList.get(i++), alignments);
    }
}
Also used : SingleFastqReader(com.milaboratory.core.io.sequence.fastq.SingleFastqReader) VDJCAlignmentsReader(com.milaboratory.mixcr.basictypes.VDJCAlignmentsReader) ArrayList(java.util.ArrayList) VDJCAlignmentsWriter(com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter) ByteArrayOutputStream(java.io.ByteArrayOutputStream) ByteArrayInputStream(java.io.ByteArrayInputStream) VDJCGene(io.repseq.core.VDJCGene) SingleRead(com.milaboratory.core.io.sequence.SingleRead) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) Test(org.junit.Test)

Aggregations

VDJCGene (io.repseq.core.VDJCGene)27 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)15 ArrayList (java.util.ArrayList)8 Test (org.junit.Test)8 VDJCAlignments (com.milaboratory.mixcr.basictypes.VDJCAlignments)7 VDJCHit (com.milaboratory.mixcr.basictypes.VDJCHit)7 VDJCAlignerParameters (com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters)6 GeneFeature (io.repseq.core.GeneFeature)6 VDJCLibrary (io.repseq.core.VDJCLibrary)6 GeneType (io.repseq.core.GeneType)5 VDJCLibraryRegistry (io.repseq.core.VDJCLibraryRegistry)5 SingleRead (com.milaboratory.core.io.sequence.SingleRead)3 NSequenceWithQuality (com.milaboratory.core.sequence.NSequenceWithQuality)3 ReferencePoint (io.repseq.core.ReferencePoint)3 VDJCGeneId (io.repseq.core.VDJCGeneId)3 ByteArrayOutputStream (java.io.ByteArrayOutputStream)3 Pattern (java.util.regex.Pattern)3 ParallelProcessor (cc.redberry.pipe.blocks.ParallelProcessor)2 AlignmentHit (com.milaboratory.core.alignment.batch.AlignmentHit)2 PairedRead (com.milaboratory.core.io.sequence.PairedRead)2