Search in sources :

Example 11 with VDJCAlignerParameters

use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.

the class ActionAssemble method go.

@Override
public void go(ActionHelper helper) throws Exception {
    // Saving initial timestamp
    long beginTimestamp = System.currentTimeMillis();
    // Checking consistency between actionParameters.doWriteClnA() value and file extension
    if ((actionParameters.getOutputFileName().toLowerCase().endsWith(".clna") && !actionParameters.doWriteClnA()) || (actionParameters.getOutputFileName().toLowerCase().endsWith(".clns") && actionParameters.doWriteClnA()))
        System.out.println("WARNING: Unexpected file extension, use .clns extension for clones-only (normal) output and\n" + ".clna if -a / --write-alignments options specified.");
    // Extracting V/D/J/C gene list from input vdjca file
    final List<VDJCGene> genes;
    final VDJCAlignerParameters alignerParameters;
    try (VDJCAlignmentsReader reader = new VDJCAlignmentsReader(actionParameters.getInputFileName(), VDJCLibraryRegistry.getDefault())) {
        genes = reader.getUsedGenes();
        // Saving aligner parameters to correct assembler parameters
        alignerParameters = reader.getParameters();
    }
    AlignmentsProvider alignmentsProvider = AlignmentsProvider.Util.createProvider(actionParameters.getInputFileName(), VDJCLibraryRegistry.getDefault());
    CloneAssemblerParameters assemblerParameters = actionParameters.getCloneAssemblerParameters();
    // set aligner parameters
    assemblerParameters.updateFrom(alignerParameters);
    // Overriding JSON parameters
    if (!actionParameters.overrides.isEmpty()) {
        assemblerParameters = JsonOverrider.override(assemblerParameters, CloneAssemblerParameters.class, actionParameters.overrides);
        if (assemblerParameters == null) {
            System.err.println("Failed to override some parameter.");
            System.exit(1);
        }
    }
    // Performing assembly
    try (CloneAssembler assembler = new CloneAssembler(assemblerParameters, actionParameters.doWriteClnA() || actionParameters.events != null, genes, alignerParameters.getFeaturesToAlignMap())) {
        // Creating event listener to collect run statistics
        CloneAssemblerReport report = new CloneAssemblerReport();
        report.setStartMillis(beginTimestamp);
        report.setInputFiles(new String[] { actionParameters.getInputFileName() });
        report.setOutputFiles(new String[] { actionParameters.getOutputFileName() });
        report.setCommandLine(helper.getCommandLineArguments());
        assembler.setListener(report);
        // Running assembler
        CloneAssemblerRunner assemblerRunner = new CloneAssemblerRunner(alignmentsProvider, assembler, actionParameters.threads);
        SmartProgressReporter.startProgressReport(assemblerRunner);
        assemblerRunner.run();
        // Getting results
        final CloneSet cloneSet = assemblerRunner.getCloneSet(alignerParameters);
        // Passing final cloneset to assemble last pieces of statistics for report
        report.onClonesetFinished(cloneSet);
        // Writing results
        if (actionParameters.doWriteClnA())
            try (ClnAWriter writer = new ClnAWriter(actionParameters.getOutputFileName())) {
                // writer will supply current stage and completion percent to the progress reporter
                SmartProgressReporter.startProgressReport(writer);
                // Writing clone block
                writer.writeClones(cloneSet);
                // Pre-soring alignments
                try (AlignmentsMappingMerger merged = new AlignmentsMappingMerger(alignmentsProvider.create(), assembler.getAssembledReadsPort())) {
                    writer.sortAlignments(merged, assembler.getAlignmentsCount());
                }
                writer.writeAlignmentsAndIndex();
            }
        else
            try (CloneSetIO.CloneSetWriter writer = new CloneSetIO.CloneSetWriter(cloneSet, actionParameters.getOutputFileName())) {
                SmartProgressReporter.startProgressReport(writer);
                writer.write();
            }
        // Writing report
        report.setFinishMillis(System.currentTimeMillis());
        assert cloneSet.getClones().size() == report.getCloneCount();
        report.setTotalReads(alignmentsProvider.getTotalNumberOfReads());
        // Writing report to stout
        System.out.println("============= Report ==============");
        Util.writeReportToStdout(report);
        if (actionParameters.report != null)
            Util.writeReport(actionParameters.report, report);
        if (actionParameters.jsonReport != null)
            Util.writeJsonReport(actionParameters.jsonReport, report);
        // Writing raw events (not documented feature)
        if (actionParameters.events != null)
            try (PipeWriter<ReadToCloneMapping> writer = new PipeWriter<>(actionParameters.events)) {
                CUtils.drain(assembler.getAssembledReadsPort(), writer);
            }
    }
}
Also used : VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) VDJCAlignmentsReader(com.milaboratory.mixcr.basictypes.VDJCAlignmentsReader) CloneSetIO(com.milaboratory.mixcr.basictypes.CloneSetIO) PipeWriter(com.milaboratory.primitivio.PipeWriter) CloneSet(com.milaboratory.mixcr.basictypes.CloneSet) ClnAWriter(com.milaboratory.mixcr.basictypes.ClnAWriter) VDJCGene(io.repseq.core.VDJCGene)

Example 12 with VDJCAlignerParameters

use of com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters in project mixcr by milaboratory.

the class PartialAlignmentsAssemblerTest method createTestData.

public static InputTestData createTestData(long seed) throws Exception {
    EnumMap<GeneType, String> geneNames = new EnumMap<GeneType, String>(GeneType.class) {

        {
            put(Variable, "TRBV20-1*00");
            put(Diversity, "TRBD2*00");
            put(Joining, "TRBJ2-6*00");
            put(Constant, "TRBC2*00");
        }
    };
    // config
    RandomGenerator rnd = RandomUtil.getThreadLocalRandom();
    rnd.setSeed(seed);
    final VDJCAlignerParameters defaultFeatures = VDJCParametersPresets.getByName("default");
    defaultFeatures.getVAlignerParameters().setGeneFeatureToAlign(VRegion);
    defaultFeatures.getDAlignerParameters().setGeneFeatureToAlign(DRegion);
    defaultFeatures.getJAlignerParameters().setGeneFeatureToAlign(JRegion);
    // used alleles
    EnumMap<GeneType, VDJCGene> genes = new EnumMap<>(GeneType.class);
    // germline parts of sequences
    EnumMap<GeneType, NucleotideSequence> germlineRegions = gtMap();
    // left, right cut of germline
    EnumMap<GeneType, int[]> germlineCuts = gtMap();
    // begin, end positions in assembled sequence
    EnumMap<GeneType, int[]> refPositions = gtMap();
    // single assembled sequence
    SequenceBuilder<NucleotideSequence> referenceBuilder = NucleotideSequence.ALPHABET.createBuilder();
    NucleotideSequence VDJunction = TestUtil.randomSequence(NucleotideSequence.ALPHABET, 3, 10);
    NucleotideSequence DJJunction = TestUtil.randomSequence(NucleotideSequence.ALPHABET, 3, 10);
    for (GeneType gt : GeneType.VDJC_REFERENCE) {
        VDJCGene gene = VDJCLibraryRegistry.getDefault().getLibrary("default", "hs").get(geneNames.get(gt));
        NucleotideSequence seq = gene.getFeature(defaultFeatures.getFeatureToAlign(gt));
        int[] cuts = null;
        switch(gt) {
            case Variable:
                cuts = new int[] { 0, rnd.nextInt(gene.getFeature(GermlineVCDR3Part).size() - 5) };
                break;
            case Diversity:
                cuts = new int[] { rnd.nextInt(seq.size() / 3), rnd.nextInt(seq.size() / 3) };
                break;
            case Joining:
                cuts = new int[] { rnd.nextInt(gene.getFeature(GermlineJCDR3Part).size() - 5), 0 };
                break;
            case Constant:
                cuts = new int[] { 0, rnd.nextInt(seq.size() / 2) };
                break;
        }
        NucleotideSequence gSeq = seq.getRange(cuts[0], seq.size() - cuts[1]);
        int[] positions = new int[2];
        positions[0] = referenceBuilder.size();
        referenceBuilder.append(gSeq);
        positions[1] = referenceBuilder.size();
        if (gt == Variable)
            referenceBuilder.append(VDJunction);
        if (gt == Diversity)
            referenceBuilder.append(DJJunction);
        genes.put(gt, gene);
        germlineCuts.put(gt, cuts);
        germlineRegions.put(gt, gSeq);
        refPositions.put(gt, positions);
    }
    NucleotideSequence VJJunction = NucleotideSequence.ALPHABET.createBuilder().append(VDJunction).append(germlineRegions.get(Diversity)).append(DJJunction).createAndDestroy();
    NucleotideSequence reference = referenceBuilder.createAndDestroy();
    return new InputTestData(genes, germlineRegions, germlineCuts, refPositions, VDJunction, DJJunction, reference, VJJunction);
}
Also used : VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) RandomGenerator(org.apache.commons.math3.random.RandomGenerator) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) VDJCGene(io.repseq.core.VDJCGene) GeneType(io.repseq.core.GeneType) EnumMap(java.util.EnumMap)

Aggregations

VDJCAlignerParameters (com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters)12 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)6 VDJCGene (io.repseq.core.VDJCGene)6 VDJCAlignments (com.milaboratory.mixcr.basictypes.VDJCAlignments)5 Test (org.junit.Test)5 ArrayList (java.util.ArrayList)4 ParallelProcessor (cc.redberry.pipe.blocks.ParallelProcessor)3 TargetBuilder (com.milaboratory.mixcr.tests.TargetBuilder)3 Well44497b (org.apache.commons.math3.random.Well44497b)3 Chunk (cc.redberry.pipe.util.Chunk)2 Indexer (cc.redberry.pipe.util.Indexer)2 ProcessException (com.milaboratory.cli.ProcessException)2 CloneAssemblerParameters (com.milaboratory.mixcr.assembler.CloneAssemblerParameters)2 PartialAlignmentsAssemblerAligner (com.milaboratory.mixcr.partialassembler.PartialAlignmentsAssemblerAligner)2 VDJCMultiRead (com.milaboratory.mixcr.partialassembler.VDJCMultiRead)2 PrimitivI (com.milaboratory.primitivio.PrimitivI)2 CUtils (cc.redberry.pipe.CUtils)1 OutputPort (cc.redberry.pipe.OutputPort)1 CountLimitingOutputPort (cc.redberry.pipe.util.CountLimitingOutputPort)1 DynamicParameter (com.beust.jcommander.DynamicParameter)1