Examples with VDJCAlignmentsWriter com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter used on opensource projects

Example 1 with VDJCAlignmentsWriter

use of com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter in project mixcr by milaboratory.

the class ActionAlign method go.

@Override
@SuppressWarnings("unchecked")
public void go(ActionHelper helper) throws Exception {
    // FIXME remove in 2.2
    if (actionParameters.printNonFunctionalWarnings())
        System.out.println("WARNING: -wf / --non-functional-warnings option is deprecated, will be removed in 2.2 " + "release. Use --verbose instead.");
    // Saving initial timestamp
    long beginTimestamp = System.currentTimeMillis();
    // Getting aligner parameters
    VDJCAlignerParameters alignerParameters = actionParameters.getAlignerParameters();
    // FIXME remove in 2.3
    if (actionParameters.getSaveOriginalReads()) {
        System.out.println("WARNING: -g / --save-reads option is deprecated, will be removed in 2.3 " + "release. Use -OsaveOriginalReads=true.");
        alignerParameters.setSaveOriginalReads(true);
    }
    // FIXME remove in 2.3
    if (actionParameters.getSaveReadDescription()) {
        System.out.println("WARNING: -a / --save-description option is deprecated, will be removed in 2.3 " + "release. Use -OsaveOriginalReads=true.");
        alignerParameters.setSaveOriginalReads(true);
    }
    if (!actionParameters.overrides.isEmpty()) {
        // Perform parameters overriding
        alignerParameters = JsonOverrider.override(alignerParameters, VDJCAlignerParameters.class, actionParameters.overrides);
        if (alignerParameters == null)
            throw new ProcessException("Failed to override some parameter.");
    }
    // FIXME remove in 2.2
    if (actionParameters.allowDifferentVJLoci != null && actionParameters.allowDifferentVJLoci) {
        System.out.println("Warning: usage of --diff-loci is deprecated. Use -OallowChimeras=true instead.");
        alignerParameters.setAllowChimeras(true);
    }
    // Creating aligner
    VDJCAligner aligner = VDJCAligner.createAligner(alignerParameters, actionParameters.isInputPaired(), !actionParameters.getNoMerge());
    // Detect if automatic featureToAlign correction is required
    int totalV = 0, totalVErrors = 0, hasVRegion = 0;
    GeneFeature correctingFeature = alignerParameters.getVAlignerParameters().getGeneFeatureToAlign().hasReversedRegions() ? GeneFeature.VRegionWithP : GeneFeature.VRegion;
    VDJCLibrary library = VDJCLibraryRegistry.getDefault().getLibrary(actionParameters.library, actionParameters.species);
    System.out.println("Reference library: " + library.getLibraryId());
    // Printing library level warnings, if specified for the library
    if (!library.getWarnings().isEmpty()) {
        System.out.println("Library warnings:");
        for (String l : library.getWarnings()) System.out.println(l);
    }
    // Printing citation notice, if specified for the library
    if (!library.getCitations().isEmpty()) {
        System.out.println("Please cite:");
        for (String l : library.getCitations()) System.out.println(l);
    }
    for (VDJCGene gene : library.getGenes(actionParameters.getChains())) {
        if (gene.getGeneType() == GeneType.Variable) {
            totalV++;
            if (!alignerParameters.containsRequiredFeature(gene)) {
                totalVErrors++;
                if (gene.getPartitioning().isAvailable(correctingFeature))
                    hasVRegion++;
            }
        }
    }
    // Performing V featureToAlign correction if needed
    if (totalVErrors > totalV * 0.9 && hasVRegion > totalVErrors * 0.8) {
        System.out.println("WARNING: forcing -OvParameters.geneFeatureToAlign=" + GeneFeature.encode(correctingFeature) + " since current gene feature (" + GeneFeature.encode(alignerParameters.getVAlignerParameters().getGeneFeatureToAlign()) + ") is absent in " + Util.PERCENT_FORMAT.format(100.0 * totalVErrors / totalV) + "% of V genes.");
        alignerParameters.getVAlignerParameters().setGeneFeatureToAlign(correctingFeature);
    }
    int numberOfExcludedNFGenes = 0;
    int numberOfExcludedFGenes = 0;
    for (VDJCGene gene : library.getGenes(actionParameters.getChains())) {
        NucleotideSequence featureSequence = alignerParameters.extractFeatureToAlign(gene);
        // exclusionReason is null ==> gene is not excluded
        String exclusionReason = null;
        if (featureSequence == null)
            exclusionReason = "absent " + GeneFeature.encode(alignerParameters.getFeatureToAlign(gene.getGeneType()));
        else if (featureSequence.containsWildcards())
            exclusionReason = "wildcard symbols in " + GeneFeature.encode(alignerParameters.getFeatureToAlign(gene.getGeneType()));
        if (exclusionReason == null)
            // If there are no reasons to exclude the gene, adding it to aligner
            aligner.addGene(gene);
        else {
            if (gene.isFunctional()) {
                ++numberOfExcludedFGenes;
                if (actionParameters.verbose())
                    System.out.println("WARNING: Functional gene " + gene.getName() + " excluded due to " + exclusionReason);
            } else
                ++numberOfExcludedNFGenes;
        }
    }
    if (actionParameters.printWarnings() && numberOfExcludedFGenes > 0)
        System.out.println("WARNING: " + numberOfExcludedFGenes + " functional genes were excluded, re-run " + "with --verbose option to see the list of excluded genes and exclusion reason.");
    if (actionParameters.verbose() && numberOfExcludedNFGenes > 0)
        System.out.println("WARNING: " + numberOfExcludedNFGenes + " non-functional genes excluded.");
    if (aligner.getVGenesToAlign().isEmpty())
        throw new ProcessException("No V genes to align. Aborting execution. See warnings for more info " + "(turn on verbose warnings by adding --verbose option).");
    if (aligner.getJGenesToAlign().isEmpty())
        throw new ProcessException("No J genes to align. Aborting execution. See warnings for more info " + "(turn on verbose warnings by adding --verbose option).");
    AlignerReport report = new AlignerReport();
    report.setStartMillis(beginTimestamp);
    report.setInputFiles(actionParameters.getInputsForReport());
    report.setOutputFiles(actionParameters.getOutputsForReport());
    report.setCommandLine(helper.getCommandLineArguments());
    // Attaching report to aligner
    aligner.setEventsListener(report);
    try (SequenceReaderCloseable<? extends SequenceRead> reader = actionParameters.createReader();
        VDJCAlignmentsWriter writer = actionParameters.getOutputName().equals(".") ? null : new VDJCAlignmentsWriter(actionParameters.getOutputName());
        SequenceWriter notAlignedWriter = actionParameters.failedReadsR1 == null ? null : (actionParameters.isInputPaired() ? new PairedFastqWriter(actionParameters.failedReadsR1, actionParameters.failedReadsR2) : new SingleFastqWriter(actionParameters.failedReadsR1))) {
        if (writer != null)
            writer.header(aligner);
        OutputPort<? extends SequenceRead> sReads = reader;
        CanReportProgress progress = (CanReportProgress) reader;
        if (actionParameters.limit != 0) {
            sReads = new CountLimitingOutputPort<>(sReads, actionParameters.limit);
            progress = SmartProgressReporter.extractProgress((CountLimitingOutputPort<?>) sReads);
        }
        final boolean writeAllResults = actionParameters.getWriteAllResults();
        EnumMap<GeneType, VDJCHit[]> emptyHits = new EnumMap<>(GeneType.class);
        for (GeneType gt : GeneType.values()) if (alignerParameters.getGeneAlignerParameters(gt) != null)
            emptyHits.put(gt, new VDJCHit[0]);
        final PairedEndReadsLayout readsLayout = alignerParameters.getReadsLayout();
        SmartProgressReporter.startProgressReport("Alignment", progress);
        OutputPort<Chunk<? extends SequenceRead>> mainInputReads = CUtils.buffered((OutputPort) chunked(sReads, 64), 16);
        OutputPort<VDJCAlignmentResult> alignments = unchunked(new ParallelProcessor(mainInputReads, chunked(aligner), actionParameters.threads));
        for (VDJCAlignmentResult result : CUtils.it(new OrderedOutputPort<>(alignments, new Indexer<VDJCAlignmentResult>() {

            @Override
            public long getIndex(VDJCAlignmentResult o) {
                return o.read.getId();
            }
        }))) {
            VDJCAlignments alignment = result.alignment;
            SequenceRead read = result.read;
            if (alignment == null) {
                if (writeAllResults) // Creating empty alignment object if alignment for current read failed
                {
                    Target target = readsLayout.createTargets(read)[0];
                    alignment = new VDJCAlignments(emptyHits, target.targets, SequenceHistory.RawSequence.of(read.getId(), target), alignerParameters.isSaveOriginalReads() ? new SequenceRead[] { read } : null);
                } else {
                    if (notAlignedWriter != null)
                        notAlignedWriter.write(result.read);
                    continue;
                }
            }
            if (alignment.isChimera())
                report.onChimera();
            if (writer != null)
                writer.write(alignment);
        }
        if (writer != null)
            writer.setNumberOfProcessedReads(reader.getNumberOfReads());
    }
    report.setFinishMillis(System.currentTimeMillis());
    // 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);
}

Example 2 with VDJCAlignmentsWriter

use of com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter 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);
    }
}

Example 3 with VDJCAlignmentsWriter

use of com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter in project mixcr by milaboratory.

the class ActionMergeAlignments method go.

@Override
public void go(ActionHelper helper) throws Exception {
    try (MultiReader reader = new MultiReader(parameters.getInputFileNames());
        VDJCAlignmentsWriter writer = new VDJCAlignmentsWriter(parameters.getOutputFileName())) {
        reader.initNextReader();
        SmartProgressReporter.startProgressReport("Merging", reader);
        writer.header(reader.currentInnerReader.getParameters(), reader.currentInnerReader.getUsedGenes());
        for (VDJCAlignments record : CUtils.it(reader)) writer.write(record);
        writer.setNumberOfProcessedReads(reader.readIdOffset.get());
    }
}

Example 4 with VDJCAlignmentsWriter

use of com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter in project mixcr by milaboratory.

the class CloneAssemblerRunnerTest method runFullPipeline.

private static CloneSet runFullPipeline(String... fastqFiles) throws IOException, InterruptedException {
    // building alignments
    VDJCAlignerParameters alignerParameters = VDJCParametersPresets.getByName("default");
    VDJCAligner aligner = fastqFiles.length == 1 ? new VDJCAlignerS(alignerParameters) : new VDJCAlignerWithMerge(alignerParameters);
    for (VDJCGene gene : VDJCLibraryRegistry.getDefault().getLibrary("mi", "hs").getGenes(Chains.IGH)) if (alignerParameters.containsRequiredFeature(gene))
        aligner.addGene(gene);
    SequenceReader reader;
    if (fastqFiles.length == 1)
        reader = new SingleFastqReader(CloneAssemblerRunnerTest.class.getClassLoader().getResourceAsStream(fastqFiles[0]), true);
    else
        reader = new PairedFastqReader(CloneAssemblerRunnerTest.class.getClassLoader().getResourceAsStream(fastqFiles[0]), CloneAssemblerRunnerTest.class.getClassLoader().getResourceAsStream(fastqFiles[1]), true);
    // write alignments to byte array
    ByteArrayOutputStream alignmentsSerialized = new ByteArrayOutputStream();
    try (VDJCAlignmentsWriter writer = new VDJCAlignmentsWriter(alignmentsSerialized)) {
        writer.header(aligner);
        for (Object read : CUtils.it(reader)) {
            VDJCAlignmentResult result = (VDJCAlignmentResult) aligner.process((SequenceRead) read);
            if (result.alignment != null)
                writer.write(result.alignment);
        }
    }
    AlignmentsProvider alignmentsProvider = AlignmentsProvider.Util.createProvider(alignmentsSerialized.toByteArray(), VDJCLibraryRegistry.getDefault());
    LinearGapAlignmentScoring<NucleotideSequence> scoring = new LinearGapAlignmentScoring<>(NucleotideSequence.ALPHABET, 5, -9, -12);
    CloneFactoryParameters factoryParameters = new CloneFactoryParameters(new VJCClonalAlignerParameters(0.8f, scoring, 5), new VJCClonalAlignerParameters(0.8f, scoring, 5), null, new DAlignerParameters(GeneFeature.DRegion, 0.85f, 30.0f, 3, scoring));
    CloneAssemblerParameters assemblerParameters = new CloneAssemblerParameters(new GeneFeature[] { GeneFeature.CDR3 }, 12, QualityAggregationType.Average, new CloneClusteringParameters(2, 1, TreeSearchParameters.ONE_MISMATCH, new RelativeConcentrationFilter(1.0E-6)), factoryParameters, true, true, false, 0.4, true, (byte) 20, .8, "2 of 6", (byte) 15);
    System.out.println(GlobalObjectMappers.toOneLine(assemblerParameters));
    CloneAssemblerRunner assemblerRunner = new CloneAssemblerRunner(alignmentsProvider, new CloneAssembler(assemblerParameters, true, aligner.getUsedGenes(), alignerParameters), 2);
    SmartProgressReporter.startProgressReport(assemblerRunner);
    assemblerRunner.run();
    CloneSet cloneSet = assemblerRunner.getCloneSet(null);
    ByteArrayOutputStream bos = new ByteArrayOutputStream();
    CloneSetIO.write(cloneSet, bos);
    CloneSet cloneSetDeserialized = CloneSetIO.readClns(new ByteArrayInputStream(bos.toByteArray()));
    assertCSEquals(cloneSet, cloneSetDeserialized);
    OutputPortCloseable<ReadToCloneMapping> rrr = assemblerRunner.assembler.getAssembledReadsPort();
    ReadToCloneMapping take;
    while ((take = rrr.take()) != null) System.out.println(take);
    return cloneSet;
}

Example 5 with VDJCAlignmentsWriter

use of com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter in project mixcr by milaboratory.

the class ActionFilterAlignments method go.

@Override
public void go(ActionHelper helper) throws Exception {
    try (VDJCAlignmentsReader reader = parameters.getInput();
        VDJCAlignmentsWriter writer = parameters.getOutput()) {
        CanReportProgress progress = reader;
        OutputPort<VDJCAlignments> sReads = reader;
        if (parameters.limit != 0) {
            sReads = new CountLimitingOutputPort<>(sReads, parameters.limit);
            progress = SmartProgressReporter.extractProgress((CountLimitingOutputPort<?>) sReads);
        }
        writer.header(reader.getParameters(), reader.getUsedGenes());
        SmartProgressReporter.startProgressReport("Filtering", progress);
        int total = 0, passed = 0;
        final AlignmentsFilter filter = parameters.getFilter();
        for (VDJCAlignments al : CUtils.it(CUtils.buffered(sReads, 2048))) {
            ++total;
            if (filter.accept(al)) {
                writer.write(al);
                ++passed;
            }
        }
        writer.setNumberOfProcessedReads(reader.getNumberOfReads());
        System.out.printf("%s alignments analysed\n", total);
        System.out.printf("%s alignments written (%.1f%%)\n", passed, 100.0 * passed / total);
    }
}