Search in sources :

Example 1 with Target

use of com.milaboratory.core.Target 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);
}
Also used : VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) VDJCAlignmentsWriter(com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter) VDJCAligner(com.milaboratory.mixcr.vdjaligners.VDJCAligner) ParallelProcessor(cc.redberry.pipe.blocks.ParallelProcessor) Target(com.milaboratory.core.Target) Indexer(cc.redberry.pipe.util.Indexer) CanReportProgress(com.milaboratory.util.CanReportProgress) SingleFastqWriter(com.milaboratory.core.io.sequence.fastq.SingleFastqWriter) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) PairedEndReadsLayout(com.milaboratory.core.PairedEndReadsLayout) VDJCAlignmentResult(com.milaboratory.mixcr.vdjaligners.VDJCAlignmentResult) PairedFastqWriter(com.milaboratory.core.io.sequence.fastq.PairedFastqWriter) Chunk(cc.redberry.pipe.util.Chunk) ProcessException(com.milaboratory.cli.ProcessException) CountLimitingOutputPort(cc.redberry.pipe.util.CountLimitingOutputPort) SequenceWriter(com.milaboratory.core.io.sequence.SequenceWriter) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) SequenceRead(com.milaboratory.core.io.sequence.SequenceRead)

Example 2 with Target

use of com.milaboratory.core.Target 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 3 with Target

use of com.milaboratory.core.Target in project mixcr by milaboratory.

the class VDJCAlignerS method processStrict.

private VDJCAlignmentResult<SingleRead> processStrict(SingleRead input) {
    Target[] targets = parameters.getReadsLayout().createTargets(input);
    // Algorithm below relies on this fact
    assert targets.length <= 2;
    boolean anyIsFull = false, anyHasJ = false, anyHasV = false;
    KVJResultsForSingle[] results = new KVJResultsForSingle[targets.length];
    for (int i = 0; i < results.length; i++) {
        results[i] = align(targets[i]);
        results[i].performChainFilteringIfNeeded();
        anyIsFull |= results[i].hasKVAndJHits();
        anyHasJ |= results[i].hasKJHits();
        anyHasV |= results[i].hasKVHits();
    }
    if (!anyIsFull) {
        if (!anyHasJ && !anyHasV)
            onFailedAlignment(input, VDJCAlignmentFailCause.NoHits);
        else if (!anyHasJ)
            onFailedAlignment(input, VDJCAlignmentFailCause.NoJHits);
        else
            onFailedAlignment(input, VDJCAlignmentFailCause.NoVHits);
        return new VDJCAlignmentResult<>(input);
    }
    KVJResultsForSingle topResult = null;
    // Calculating best result
    for (KVJResultsForSingle result : results) if (result.hasKVAndJHits())
        topResult = topResult != null ? null : result;
    if (topResult == null) {
        // Finalizing alignment for both results to determine who is the best
        for (KVJResultsForSingle result : results) {
            result.alignDC();
            if (topResult == null || topResult.sumScore() < result.sumScore())
                topResult = result;
        }
    } else
        // Align C and D genes only for best result
        topResult.alignDC();
    // Checking minimal sum score
    if (topResult.sumScore() < parameters.getMinSumScore()) {
        onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
        return new VDJCAlignmentResult<>(input);
    }
    // Filtering hits basing on minSumScore
    topResult.calculateHits(parameters.getMinSumScore(), parameters.getMaxHits());
    if (topResult.hasVAndJHits()) {
        VDJCAlignments alignment = topResult.toVDJCAlignments(input.getId(), input);
        onSuccessfulAlignment(input, alignment);
        return new VDJCAlignmentResult<>(input, alignment);
    } else {
        onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
        return new VDJCAlignmentResult<>(input);
    }
}
Also used : Target(com.milaboratory.core.Target) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments)

Example 4 with Target

use of com.milaboratory.core.Target in project mixcr by milaboratory.

the class VDJCAlignerS method processPartial.

private VDJCAlignmentResult<SingleRead> processPartial(SingleRead input) {
    Target[] targets = parameters.getReadsLayout().createTargets(input);
    KVJResultsForSingle[] results = new KVJResultsForSingle[targets.length];
    for (int i = 0; i < results.length; i++) {
        results[i] = align(targets[i]);
        results[i].performChainFilteringIfNeeded();
    }
    KVJResultsForSingle topResult = null;
    for (KVJResultsForSingle result : results) {
        if (!result.hasNoKVNorKJHits())
            result.alignDC();
        if (topResult == null || topResult.sumScore() < result.sumScore())
            topResult = result;
    }
    if (topResult.hasNoKVNorKJHits()) {
        onFailedAlignment(input, VDJCAlignmentFailCause.NoHits);
        return new VDJCAlignmentResult<>(input);
    }
    // Checking minimal sum score
    if (topResult.sumScore() < parameters.getMinSumScore()) {
        onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
        return new VDJCAlignmentResult<>(input);
    }
    // Filtering hits basing on minSumScore
    topResult.calculateHits(parameters.getMinSumScore(), parameters.getMaxHits());
    VDJCAlignments alignment = topResult.toVDJCAlignments(input.getId(), input);
    // Final check
    if (!parameters.getAllowNoCDR3PartAlignments()) {
        // CDR3 Begin / End
        if (!alignment.getPartitionedTarget(0).getPartitioning().isAvailable(reqPointL) && !alignment.getPartitionedTarget(0).getPartitioning().isAvailable(reqPointR)) {
            onFailedAlignment(input, VDJCAlignmentFailCause.NoCDR3Parts);
            return new VDJCAlignmentResult<>(input);
        }
    }
    // Read successfully aligned
    onSuccessfulAlignment(input, alignment);
    return new VDJCAlignmentResult<>(input, alignment);
}
Also used : Target(com.milaboratory.core.Target) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments)

Aggregations

Target (com.milaboratory.core.Target)4 VDJCAlignments (com.milaboratory.mixcr.basictypes.VDJCAlignments)4 SequenceRead (com.milaboratory.core.io.sequence.SequenceRead)2 ParallelProcessor (cc.redberry.pipe.blocks.ParallelProcessor)1 Chunk (cc.redberry.pipe.util.Chunk)1 CountLimitingOutputPort (cc.redberry.pipe.util.CountLimitingOutputPort)1 Indexer (cc.redberry.pipe.util.Indexer)1 ProcessException (com.milaboratory.cli.ProcessException)1 PairedEndReadsLayout (com.milaboratory.core.PairedEndReadsLayout)1 PairedRead (com.milaboratory.core.io.sequence.PairedRead)1 SequenceWriter (com.milaboratory.core.io.sequence.SequenceWriter)1 SingleRead (com.milaboratory.core.io.sequence.SingleRead)1 SingleReadImpl (com.milaboratory.core.io.sequence.SingleReadImpl)1 PairedFastqWriter (com.milaboratory.core.io.sequence.fastq.PairedFastqWriter)1 SingleFastqWriter (com.milaboratory.core.io.sequence.fastq.SingleFastqWriter)1 NSequenceWithQuality (com.milaboratory.core.sequence.NSequenceWithQuality)1 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)1 SequenceHistory (com.milaboratory.mixcr.basictypes.SequenceHistory)1 VDJCAlignmentsWriter (com.milaboratory.mixcr.basictypes.VDJCAlignmentsWriter)1 AlignedTarget (com.milaboratory.mixcr.partialassembler.AlignedTarget)1