Examples with NucleotideSequence com.milaboratory.core.sequence.NucleotideSequence used on opensource projects

Example 1 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence 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 NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class ActionExportAlignmentsPretty method outputVerbose.

public void outputVerbose(PrintStream output, final VDJCAlignments alignments) {
    output.println(">>> Read ids: " + Arrays.toString(alignments.getReadIds()).replace("[", "").replace("]", ""));
    output.println();
    output.println(">>> Target sequences (input sequences):");
    output.println();
    for (int i = 0; i < alignments.numberOfTargets(); i++) {
        output.println("Sequence" + i + ":");
        final VDJCPartitionedSequence partitionedTarget = alignments.getPartitionedTarget(i);
        printGeneFeatures(new Filter<GeneFeature>() {

            @Override
            public boolean accept(GeneFeature object) {
                return partitionedTarget.getPartitioning().isAvailable(object);
            }
        }, output, "Contains features: ");
        output.println();
        output.print(new NSequenceWithQualityPrintHelper(alignments.getTarget(i), LINE_OFFSET, LINE_LENGTH));
    }
    if (alignments.numberOfTargets() > 1) {
        // Printing a set of available gene features for a full read
        output.println(">>> Gene features that can be extracted from this paired-read: ");
        printGeneFeatures(new Filter<GeneFeature>() {

            @Override
            public boolean accept(GeneFeature object) {
                return alignments.getFeature(object) != null;
            }
        }, output, "");
    }
    output.println();
    for (GeneType geneType : GeneType.values()) {
        output.println(">>> Alignments with " + geneType.getLetter() + " gene:");
        output.println();
        boolean exists = false;
        VDJCHit[] hits = alignments.getHits(geneType);
        if (hits.length > 0) {
            hits = actionParameters.isOnlyTop() ? new VDJCHit[] { hits[0] } : hits;
            for (VDJCHit hit : hits) {
                exists = true;
                output.println(hit.getGene().getName() + " (total score = " + hit.getScore() + ")");
                for (int i = 0; i < alignments.numberOfTargets(); i++) {
                    Alignment<NucleotideSequence> alignment = hit.getAlignment(i);
                    if (alignment == null)
                        continue;
                    output.println("Alignment of Sequence" + i + " (score = " + (alignment == null ? "NaN" : alignment.getScore()) + "):");
                    if (alignment != null) {
                        for (AlignmentHelper subHelper : alignment.getAlignmentHelper().split(LINE_LENGTH, LINE_OFFSET)) {
                            output.println(subHelper.toStringWithSeq2Quality(alignments.getTarget(i).getQuality()));
                            output.println();
                        }
                        if (actionParameters.printGeneSequence()) {
                            output.println("Gene sequence:");
                            output.println(alignment.getSequence1());
                            output.println();
                        }
                    }
                }
            }
        }
        if (!exists) {
            output.println("No hits.");
            output.println();
        }
    }
    char[] ll = new char[94];
    Arrays.fill(ll, '=');
    output.println(ll);
    output.println();
}

Example 3 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class ActionClonesDiff method getKey.

private CKey getKey(Clone clone) {
    final NucleotideSequence[] clonalSequence = new NucleotideSequence[clone.numberOfTargets()];
    for (int i = 0; i < clonalSequence.length; i++) clonalSequence[i] = clone.getTarget(i).getSequence();
    final VDJCGeneId v = params.useV() ? getBestGene(clone, GeneType.Variable) : null;
    final VDJCGeneId j = params.useJ() ? getBestGene(clone, GeneType.Joining) : null;
    final VDJCGeneId c = params.useC() ? getBestGene(clone, GeneType.Constant) : null;
    return new CKey(clonalSequence, v, j, c);
}

Example 4 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class VDJCAlignmentsFormatter method drawAASequence.

public static void drawAASequence(MultiAlignmentHelper helper, SequencePartitioning partitioning, NucleotideSequence target) {
    List<RangeTranslationParameters> trParams = partitioning.getTranslationParameters(target.size());
    char[] line = new char[helper.size()];
    Arrays.fill(line, ' ');
    for (RangeTranslationParameters trParam : trParams) {
        NucleotideSequence mainSequence = target.getRange(trParam.range);
        NucleotideSequence leftover = trParam.codonLeftoverRange == null ? null : target.getRange(trParam.codonLeftoverRange);
        NucleotideSequence bigSeq = leftover == null ? mainSequence : trParam.leftIncompleteCodonRange() != null ? leftover.concatenate(mainSequence) : mainSequence.concatenate(leftover);
        AminoAcidSequence aa = AminoAcidSequence.translate(bigSeq, trParam.translationParameters);
        int aaPosition = 0;
        int ntPosition = trParam.range.getFrom() + AminoAcidSequence.convertAAPositionToNt(aaPosition, mainSequence.size(), trParam.translationParameters);
        if (aa.codeAt(aaPosition) == AminoAcidAlphabet.INCOMPLETE_CODON) {
            line[helper.subjectToAlignmentPosition(ntPosition)] = // '_'
            AminoAcidSequence.ALPHABET.codeToSymbol(aa.codeAt(aaPosition));
            ++aaPosition;
        }
        while (aaPosition < aa.size() && aa.codeAt(aaPosition) != AminoAcidAlphabet.INCOMPLETE_CODON) {
            ntPosition = trParam.range.getFrom() + AminoAcidSequence.convertAAPositionToNt(aaPosition, mainSequence.size(), trParam.translationParameters);
            line[helper.subjectToAlignmentPosition(ntPosition + 1)] = AminoAcidSequence.ALPHABET.codeToSymbol(aa.codeAt(aaPosition));
            ++aaPosition;
        }
        if (aaPosition < aa.size() && aa.codeAt(aaPosition) == AminoAcidAlphabet.INCOMPLETE_CODON) {
            ntPosition = trParam.range.getFrom() + AminoAcidSequence.convertAAPositionToNt(aaPosition, mainSequence.size(), trParam.translationParameters);
            line[helper.subjectToAlignmentPosition(ntPosition)] = AminoAcidSequence.ALPHABET.codeToSymbol(aa.codeAt(aaPosition));
        }
    }
    helper.addAnnotationString("", new String(line));
}

Example 5 with NucleotideSequence

use of com.milaboratory.core.sequence.NucleotideSequence in project mixcr by milaboratory.

the class VDJCHit method getIdentity.

public float getIdentity() {
    float identity = 0;
    int tSize = 0;
    for (Alignment<NucleotideSequence> alignment : alignments) {
        if (alignment == null)
            continue;
        AlignmentHelper h = alignment.getAlignmentHelper();
        identity += h.identity() * h.size();
        tSize += h.size();
    }
    return identity / tSize;
}