Examples with SequenceRead com.milaboratory.core.io.sequence.SequenceRead used on opensource projects

Example 1 with SequenceRead

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

use of com.milaboratory.core.io.sequence.SequenceRead 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;
}

Example 3 with SequenceRead

use of com.milaboratory.core.io.sequence.SequenceRead in project mixcr by milaboratory.

the class FieldExtractors method getFields.

public static synchronized Field[] getFields() {
    if (descriptors == null) {
        List<Field> descriptorsList = new ArrayList<>();
        // Number of targets
        descriptorsList.add(new PL_O("-targets", "Export number of targets", "Number of targets", "numberOfTargets") {

            @Override
            protected String extract(VDJCObject object) {
                return Integer.toString(object.numberOfTargets());
            }
        });
        // Best hits
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new PL_O("-" + Character.toLowerCase(l) + "Hit", "Export best " + l + " hit", "Best " + l + " hit", "best" + l + "Hit") {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit bestHit = object.getBestHit(type);
                    if (bestHit == null)
                        return NULL;
                    return bestHit.getGene().getName();
                }
            });
        }
        // Best gene
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new PL_O("-" + Character.toLowerCase(l) + "Gene", "Export best " + l + " hit gene name (e.g. TRBV12-3 for TRBV12-3*00)", "Best " + l + " gene", "best" + l + "Gene") {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit bestHit = object.getBestHit(type);
                    if (bestHit == null)
                        return NULL;
                    return bestHit.getGene().getGeneName();
                }
            });
        }
        // Best family
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new PL_O("-" + Character.toLowerCase(l) + "Family", "Export best " + l + " hit family name (e.g. TRBV12 for TRBV12-3*00)", "Best " + l + " family", "best" + l + "Family") {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit bestHit = object.getBestHit(type);
                    if (bestHit == null)
                        return NULL;
                    return bestHit.getGene().getFamilyName();
                }
            });
        }
        // Best hit score
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new PL_O("-" + Character.toLowerCase(l) + "HitScore", "Export score for best " + l + " hit", "Best " + l + " hit score", "best" + l + "HitScore") {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit bestHit = object.getBestHit(type);
                    if (bestHit == null)
                        return NULL;
                    return String.valueOf(bestHit.getScore());
                }
            });
        }
        // All hits
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new PL_O("-" + Character.toLowerCase(l) + "HitsWithScore", "Export all " + l + " hits with score", "All " + l + " hits", "all" + l + "HitsWithScore") {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit[] hits = object.getHits(type);
                    if (hits.length == 0)
                        return "";
                    StringBuilder sb = new StringBuilder();
                    for (int i = 0; ; i++) {
                        sb.append(hits[i].getGene().getName()).append("(").append(SCORE_FORMAT.format(hits[i].getScore())).append(")");
                        if (i == hits.length - 1)
                            break;
                        sb.append(",");
                    }
                    return sb.toString();
                }
            });
        }
        // All hits without score
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new PL_O("-" + Character.toLowerCase(l) + "Hits", "Export all " + l + " hits", "All " + l + " Hits", "all" + l + "Hits") {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit[] hits = object.getHits(type);
                    if (hits.length == 0)
                        return "";
                    StringBuilder sb = new StringBuilder();
                    for (int i = 0; ; i++) {
                        sb.append(hits[i].getGene().getName());
                        if (i == hits.length - 1)
                            break;
                        sb.append(",");
                    }
                    return sb.toString();
                }
            });
        }
        // All gene names
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new StringExtractor("-" + Character.toLowerCase(l) + "Genes", "Export all " + l + " gene names (e.g. TRBV12-3 for TRBV12-3*00)", "All " + l + " genes", "all" + l + "Genes", type) {

                @Override
                String extractStringForHit(VDJCHit hit) {
                    return hit.getGene().getGeneName();
                }
            });
        }
        // All families
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new StringExtractor("-" + Character.toLowerCase(l) + "Families", "Export all " + l + " gene family anmes (e.g. TRBV12 for TRBV12-3*00)", "All " + l + " families", "all" + l + "Families", type) {

                @Override
                String extractStringForHit(VDJCHit hit) {
                    return hit.getGene().getFamilyName();
                }
            });
        }
        // Best alignment
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new PL_O("-" + Character.toLowerCase(l) + "Alignment", "Export best " + l + " alignment", "Best " + l + " alignment", "best" + l + "Alignment") {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit bestHit = object.getBestHit(type);
                    if (bestHit == null)
                        return NULL;
                    StringBuilder sb = new StringBuilder();
                    for (int i = 0; ; i++) {
                        Alignment<NucleotideSequence> alignment = bestHit.getAlignment(i);
                        if (alignment == null)
                            sb.append(NULL);
                        else
                            sb.append(alignment.toCompactString());
                        if (i == object.numberOfTargets() - 1)
                            break;
                        sb.append(",");
                    }
                    return sb.toString();
                }
            });
        }
        // All alignments
        for (final GeneType type : GeneType.values()) {
            char l = type.getLetter();
            descriptorsList.add(new PL_O("-" + Character.toLowerCase(l) + "Alignments", "Export all " + l + " alignments", "All " + l + " alignments", "all" + l + "Alignments") {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit[] hits = object.getHits(type);
                    if (hits.length == 0)
                        return "";
                    StringBuilder sb = new StringBuilder();
                    for (int j = 0; ; ++j) {
                        for (int i = 0; ; i++) {
                            Alignment<NucleotideSequence> alignment = hits[j].getAlignment(i);
                            if (alignment == null)
                                sb.append(NULL);
                            else
                                sb.append(alignment.toCompactString());
                            if (i == object.numberOfTargets() - 1)
                                break;
                            sb.append(',');
                        }
                        if (j == hits.length - 1)
                            break;
                        sb.append(';');
                    }
                    return sb.toString();
                }
            });
        }
        descriptorsList.add(new FeatureExtractors.NSeqExtractor("-nFeature", "Export nucleotide sequence of specified gene feature", "N. Seq. ", "nSeq") {

            @Override
            public String convert(NSequenceWithQuality seq) {
                return seq.getSequence().toString();
            }
        });
        descriptorsList.add(new FeatureExtractors.NSeqExtractor("-qFeature", "Export quality string of specified gene feature", "Qual. ", "qual") {

            @Override
            public String convert(NSequenceWithQuality seq) {
                return seq.getQuality().toString();
            }
        });
        descriptorsList.add(new FeatureExtractors.WithHeader("-aaFeature", "Export amino acid sequence of specified gene feature", 1, new String[] { "AA. Seq. " }, new String[] { "aaSeq" }) {

            @Override
            protected String extractValue(VDJCObject object, GeneFeature[] parameters) {
                GeneFeature geneFeature = parameters[parameters.length - 1];
                NSequenceWithQuality feature = object.getFeature(geneFeature);
                if (feature == null)
                    return NULL;
                int targetId = object.getTargetContainingFeature(geneFeature);
                TranslationParameters tr = targetId == -1 ? TranslationParameters.FromLeftWithIncompleteCodon : object.getPartitionedTarget(targetId).getPartitioning().getTranslationParameters(geneFeature);
                if (tr == null)
                    return NULL;
                return AminoAcidSequence.translate(feature.getSequence(), tr).toString();
            }
        });
        // descriptorsList.add(new FeatureExtractorDescriptor("-aaFeatureFromLeft", "Export amino acid sequence of " +
        // "specified gene feature starting from the leftmost nucleotide (differs from -aaFeature only for " +
        // "sequences which length are not multiple of 3)", "AA. Seq.", "aaSeq") {
        // @Override
        // public String convert(NSequenceWithQuality seq) {
        // return AminoAcidSequence.translate(seq.getSequence(), FromLeftWithoutIncompleteCodon).toString();
        // }
        // });
        // 
        // descriptorsList.add(new FeatureExtractorDescriptor("-aaFeatureFromRight", "Export amino acid sequence of " +
        // "specified gene feature starting from the rightmost nucleotide (differs from -aaFeature only for " +
        // "sequences which length are not multiple of 3)", "AA. Seq.", "aaSeq") {
        // @Override
        // public String convert(NSequenceWithQuality seq) {
        // return AminoAcidSequence.translate(seq.getSequence(), FromRightWithoutIncompleteCodon).toString();
        // }
        // });
        descriptorsList.add(new FeatureExtractors.NSeqExtractor("-minFeatureQuality", "Export minimal quality of specified gene feature", "Min. qual. ", "minQual") {

            @Override
            public String convert(NSequenceWithQuality seq) {
                return "" + seq.getQuality().minValue();
            }
        });
        descriptorsList.add(new FeatureExtractors.NSeqExtractor("-avrgFeatureQuality", "Export average quality of specified gene feature", "Mean. qual. ", "meanQual") {

            @Override
            public String convert(NSequenceWithQuality seq) {
                return "" + seq.getQuality().meanValue();
            }
        });
        descriptorsList.add(new FeatureExtractors.NSeqExtractor("-lengthOf", "Exports length of specified gene feature.", "Length of ", "lengthOf") {

            @Override
            public String convert(NSequenceWithQuality seq) {
                return "" + seq.size();
            }
        });
        descriptorsList.add(new FeatureExtractors.MutationsExtractor("-nMutations", "Extract nucleotide mutations for specific gene feature; relative to germline sequence.", 1, new String[] { "N. Mutations in " }, new String[] { "nMutations" }) {

            @Override
            String convert(Mutations<NucleotideSequence> mutations, NucleotideSequence seq1, NucleotideSequence seq2, TranslationParameters tr) {
                return mutations.encode(",");
            }
        });
        descriptorsList.add(new FeatureExtractors.MutationsExtractor("-nMutationsRelative", "Extract nucleotide mutations for specific gene feature relative to another feature.", 2, new String[] { "N. Mutations in ", " relative to " }, new String[] { "nMutationsIn", "Relative" }) {

            @Override
            String convert(Mutations<NucleotideSequence> mutations, NucleotideSequence seq1, NucleotideSequence seq2, TranslationParameters tr) {
                return mutations.encode(",");
            }
        });
        final class AAMutations extends FeatureExtractors.MutationsExtractor {

            AAMutations(String command, String description, int nArgs, String[] hPrefix, String[] sPrefix) {
                super(command, description, nArgs, hPrefix, sPrefix);
            }

            @Override
            String convert(Mutations<NucleotideSequence> mutations, NucleotideSequence seq1, NucleotideSequence seq2, TranslationParameters tr) {
                if (tr == null)
                    return "-";
                Mutations<AminoAcidSequence> aaMuts = MutationsUtil.nt2aa(seq1, mutations, tr);
                if (aaMuts == null)
                    return "-";
                return aaMuts.encode(",");
            }
        }
        descriptorsList.add(new AAMutations("-aaMutations", "Extract amino acid mutations for specific gene feature", 1, new String[] { "AA. Mutations in " }, new String[] { "aaMutations" }));
        descriptorsList.add(new AAMutations("-aaMutationsRelative", "Extract amino acid mutations for specific gene feature relative to another feature.", 2, new String[] { "AA. Mutations in ", " relative to " }, new String[] { "aaMutationsIn", "Relative" }));
        final class MutationsDetailed extends FeatureExtractors.MutationsExtractor {

            MutationsDetailed(String command, String description, int nArgs, String[] hPrefix, String[] sPrefix) {
                super(command, description, nArgs, hPrefix, sPrefix);
            }

            @Override
            String convert(Mutations<NucleotideSequence> mutations, NucleotideSequence seq1, NucleotideSequence seq2, TranslationParameters tr) {
                if (tr == null)
                    return "-";
                MutationsUtil.MutationNt2AADescriptor[] descriptors = MutationsUtil.nt2aaDetailed(seq1, mutations, tr, 10);
                if (descriptors == null)
                    return "-";
                StringBuilder sb = new StringBuilder();
                for (int i = 0; i < descriptors.length; i++) {
                    sb.append(descriptors[i]);
                    if (i == descriptors.length - 1)
                        break;
                    sb.append(",");
                }
                return sb.toString();
            }
        }
        String detailedMutationsFormat = "Format <nt_mutation>:<aa_mutation_individual>:<aa_mutation_cumulative>, where <aa_mutation_individual> is an expected amino acid " + "mutation given no other mutations have occurred, and <aa_mutation_cumulative> amino acid mutation is the observed amino acid " + "mutation combining effect from all other. WARNING: format may change in following versions.";
        descriptorsList.add(new MutationsDetailed("-mutationsDetailed", "Detailed list of nucleotide and corresponding amino acid mutations. " + detailedMutationsFormat, 1, new String[] { "Detailed mutations in " }, new String[] { "mutationsDetailedIn" }));
        descriptorsList.add(new MutationsDetailed("-mutationsDetailedRelative", "Detailed list of nucleotide and corresponding amino acid mutations written, positions relative to specified gene feature. " + detailedMutationsFormat, 2, new String[] { "Detailed mutations in ", " relative to " }, new String[] { "mutationsDetailedIn", "Relative" }));
        descriptorsList.add(new ExtractReferencePointPosition());
        descriptorsList.add(new ExtractDefaultReferencePointsPositions());
        descriptorsList.add(new PL_A("-readId", "Export id of read corresponding to alignment", "Read id", "readId") {

            @Override
            protected String extract(VDJCAlignments object) {
                return "" + object.getMinReadId();
            }

            @Override
            public FieldExtractor<VDJCAlignments> create(OutputMode outputMode, String[] args) {
                System.out.println("WARNING: -readId is deprecated. Use -readIds");
                return super.create(outputMode, args);
            }
        });
        descriptorsList.add(new PL_A("-readIds", "Export id of read corresponding to alignment", "Read id", "readId") {

            @Override
            protected String extract(VDJCAlignments object) {
                long[] readIds = object.getReadIds();
                StringBuilder sb = new StringBuilder();
                for (int i = 0; ; i++) {
                    sb.append(readIds[i]);
                    if (i == readIds.length - 1)
                        return sb.toString();
                    sb.append(",");
                }
            }
        });
        descriptorsList.add(new ExtractSequence(VDJCAlignments.class, "-sequence", "Export aligned sequence (initial read), or 2 sequences in case of paired-end reads", "Read(s) sequence", "readSequence"));
        descriptorsList.add(new ExtractSequenceQuality(VDJCAlignments.class, "-quality", "Export initial read quality, or 2 qualities in case of paired-end reads", "Read(s) sequence qualities", "readQuality"));
        descriptorsList.add(new PL_C("-cloneId", "Unique clone identifier", "Clone ID", "cloneId") {

            @Override
            protected String extract(Clone object) {
                return "" + object.getId();
            }
        });
        descriptorsList.add(new PL_C("-count", "Export clone count", "Clone count", "cloneCount") {

            @Override
            protected String extract(Clone object) {
                return "" + object.getCount();
            }
        });
        descriptorsList.add(new PL_C("-fraction", "Export clone fraction", "Clone fraction", "cloneFraction") {

            @Override
            protected String extract(Clone object) {
                return "" + object.getFraction();
            }
        });
        descriptorsList.add(new ExtractSequence(Clone.class, "-sequence", "Export aligned sequence (initial read), or 2 sequences in case of paired-end reads", "Clonal sequence(s)", "clonalSequence"));
        descriptorsList.add(new ExtractSequenceQuality(Clone.class, "-quality", "Export initial read quality, or 2 qualities in case of paired-end reads", "Clonal sequence quality(s)", "clonalSequenceQuality"));
        descriptorsList.add(new PL_A("-descrR1", "Export description line from initial .fasta or .fastq file " + "of the first read (only available if --save-description was used in align command)", "Description R1", "descrR1") {

            @Override
            protected String extract(VDJCAlignments object) {
                List<SequenceRead> reads = object.getOriginalReads();
                if (reads == null)
                    throw new IllegalArgumentException("Error for option \'-descrR1\':\n" + "No description available for read: either re-run align action with -OsaveOriginalReads=true option " + "or don't use \'-descrR1\' in exportAlignments");
                return reads.get(0).getRead(0).getDescription();
            }

            @Override
            public FieldExtractor<VDJCAlignments> create(OutputMode outputMode, String[] args) {
                System.out.println("WARNING: -descrR1 is deprecated. Use -descrsR1");
                return super.create(outputMode, args);
            }
        });
        descriptorsList.add(new PL_A("-descrR2", "Export description line from initial .fasta or .fastq file " + "of the second read (only available if --save-description was used in align command)", "Description R2", "descrR2") {

            @Override
            protected String extract(VDJCAlignments object) {
                List<SequenceRead> reads = object.getOriginalReads();
                if (reads == null)
                    throw new IllegalArgumentException("Error for option \'-descrR1\':\n" + "No description available for read: either re-run align action with -OsaveOriginalReads=true option " + "or don't use \'-descrR1\' in exportAlignments");
                SequenceRead read = reads.get(0);
                if (read.numberOfReads() < 2)
                    throw new IllegalArgumentException("Error for option \'-descrR2\':\n" + "No description available for second read: your input data was single-end");
                return read.getRead(1).getDescription();
            }

            @Override
            public FieldExtractor<VDJCAlignments> create(OutputMode outputMode, String[] args) {
                System.out.println("WARNING: -descrR2 is deprecated. Use -descrsR2");
                return super.create(outputMode, args);
            }
        });
        descriptorsList.add(new PL_A("-descrsR1", "Export description lines from initial .fasta or .fastq file " + "of the first reads (only available if -OsaveOriginalReads=true was used in align command)", "Descriptions R1", "descrsR1") {

            @Override
            protected String extract(VDJCAlignments object) {
                List<SequenceRead> reads = object.getOriginalReads();
                if (reads == null)
                    throw new IllegalArgumentException("Error for option \'-descrR1\':\n" + "No description available for read: either re-run align action with -OsaveOriginalReads option " + "or don't use \'-descrR1\' in exportAlignments");
                StringBuilder sb = new StringBuilder();
                for (int i = 0; ; i++) {
                    sb.append(reads.get(i).getRead(0).getDescription());
                    if (i == reads.size() - 1)
                        return sb.toString();
                    sb.append(",");
                }
            }
        });
        descriptorsList.add(new PL_A("-descrsR2", "Export description lines from initial .fasta or .fastq file " + "of the second reads (only available if -OsaveOriginalReads=true was used in align command)", "Descriptions R2", "descrsR2") {

            @Override
            protected String extract(VDJCAlignments object) {
                List<SequenceRead> reads = object.getOriginalReads();
                if (reads == null)
                    throw new IllegalArgumentException("Error for option \'-descrR1\':\n" + "No description available for read: either re-run align action with -OsaveOriginalReads option " + "or don't use \'-descrR1\' in exportAlignments");
                StringBuilder sb = new StringBuilder();
                for (int i = 0; ; i++) {
                    SequenceRead read = reads.get(i);
                    if (read.numberOfReads() < 2)
                        throw new IllegalArgumentException("Error for option \'-descrsR2\':\n" + "No description available for second read: your input data was single-end");
                    sb.append(read.getRead(1).getDescription());
                    if (i == reads.size() - 1)
                        return sb.toString();
                    sb.append(",");
                }
            }
        });
        descriptorsList.add(new PL_A("-readHistory", "Export read history", "Read history", "readHistory") {

            @Override
            protected String extract(VDJCAlignments object) {
                try {
                    return GlobalObjectMappers.toOneLine(object.getHistory());
                } catch (JsonProcessingException ex) {
                    throw new RuntimeException(ex);
                }
            }
        });
        for (final GeneType type : GeneType.values()) {
            String c = Character.toLowerCase(type.getLetter()) + "IdentityPercents";
            descriptorsList.add(new PL_O("-" + c, type.getLetter() + " alignment identity percents", type.getLetter() + " alignment identity percents", c) {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit[] hits = object.getHits(type);
                    if (hits == null)
                        return NULL;
                    StringBuilder sb = new StringBuilder();
                    sb.append("");
                    for (int i = 0; ; i++) {
                        sb.append(hits[i].getIdentity());
                        if (i == hits.length - 1)
                            return sb.toString();
                        sb.append(",");
                    }
                }
            });
        }
        for (final GeneType type : GeneType.values()) {
            String c = Character.toLowerCase(type.getLetter()) + "BestIdentityPercent";
            descriptorsList.add(new PL_O("-" + c, type.getLetter() + "best alignment identity percent", type.getLetter() + "best alignment identity percent", c) {

                @Override
                protected String extract(VDJCObject object) {
                    VDJCHit hit = object.getBestHit(type);
                    if (hit == null)
                        return NULL;
                    return Float.toString(hit.getIdentity());
                }
            });
        }
        descriptorsList.add(new PL_O("-chains", "Chains", "Chains", "Chains") {

            @Override
            protected String extract(VDJCObject object) {
                return object.commonChains().toString();
            }
        });
        descriptorsList.add(new PL_O("-topChains", "Top chains", "Top chains", "topChains") {

            @Override
            protected String extract(VDJCObject object) {
                return object.commonTopChains().toString();
            }
        });
        descriptors = descriptorsList.toArray(new Field[descriptorsList.size()]);
    }
    return descriptors;
}

Example 4 with SequenceRead

use of com.milaboratory.core.io.sequence.SequenceRead in project mixcr by milaboratory.

the class ActionExportCloneReads method go.

@Override
public void go(ActionHelper helper) throws Exception {
    try (ClnAReader clna = new ClnAReader(parameters.getInputFileName(), VDJCLibraryRegistry.createDefaultRegistry())) {
        VDJCAlignments firstAlignment = clna.readAllAlignments().take();
        if (firstAlignment == null)
            return;
        if (firstAlignment.getOriginalReads() == null)
            throw new ParameterException("Error: original reads were not saved in the .vdjca file: " + "re-run align with '-g' option.");
        int[] cid = parameters.getCloneIds();
        Supplier<IntStream> cloneIds;
        if (cid == null)
            cloneIds = () -> IntStream.range(0, clna.numberOfClones());
        else
            cloneIds = () -> IntStream.of(cid);
        long totalAlignments = cloneIds.get().mapToLong(clna::numberOfAlignmentsInClone).sum();
        AtomicLong alignmentsWritten = new AtomicLong();
        AtomicBoolean finished = new AtomicBoolean(false);
        SmartProgressReporter.startProgressReport("Writing reads", new CanReportProgress() {

            @Override
            public double getProgress() {
                return 1.0 * alignmentsWritten.get() / totalAlignments;
            }

            @Override
            public boolean isFinished() {
                return finished.get();
            }
        });
        boolean paired = firstAlignment.getOriginalReads().get(0).numberOfReads() == 2;
        boolean separate = parameters.doSeparate();
        SequenceWriter globalWriter = separate ? null : createWriter(paired, parameters.getOutputFileName());
        cloneIds.get().forEach(cloneId -> {
            try (SequenceWriter individualWriter = globalWriter == null ? createWriter(paired, cloneFile(parameters.getOutputFileName(), cloneId)) : null) {
                SequenceWriter actualWriter = globalWriter == null ? individualWriter : globalWriter;
                for (VDJCAlignments alignments : CUtils.it(clna.readAlignmentsOfClone(cloneId))) {
                    for (SequenceRead read : alignments.getOriginalReads()) actualWriter.write(read);
                    alignmentsWritten.incrementAndGet();
                }
            } catch (IOException e) {
                throw new RuntimeException(e);
            }
        });
    }
}

Example 5 with SequenceRead

use of com.milaboratory.core.io.sequence.SequenceRead 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;
}