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

Example 1 with TranslationParameters

use of com.milaboratory.core.sequence.TranslationParameters 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 2 with TranslationParameters

use of com.milaboratory.core.sequence.TranslationParameters in project repseqio by repseqio.

the class InferAnchorPointsAction method go.

public void go(GeneFeature geneFeature, String inputFile, String outputFile) throws Exception {
    VDJCLibraryRegistry.resetDefaultRegistry();
    VDJCLibraryRegistry reg = VDJCLibraryRegistry.getDefault();
    // Map of library names
    Map<String, String> libraryNameToAddress = new HashMap<>();
    // Registering reference library
    int i = 0;
    for (String refAddress : params.getReference()) {
        String name = REFERENCE_LIBRARY_PREFIX + (i++);
        reg.registerLibraries(refAddress, name);
        libraryNameToAddress.put(name, refAddress);
    }
    if (params.getReference().isEmpty()) {
        reg.loadAllLibraries("default");
        for (VDJCLibrary library : reg.getLoadedLibraries()) libraryNameToAddress.put(library.getName(), "built-in");
    }
    // Registering target library
    reg.registerLibraries(inputFile, TARGET_LIBRARY_NAME);
    // Compile gene filter
    Pattern namePattern = params.name == null ? null : Pattern.compile(params.name);
    List<VDJCLibrary> refLibraries = reg.getLoadedLibrariesByNamePattern(REFERENCE_LIBRARY_PATTERN);
    SimpleBatchAlignerParameters<AminoAcidSequence> aParams = new SimpleBatchAlignerParameters<>(5, 0.4f, params.getAbsoluteMinScore(geneFeature), true, AffineGapAlignmentScoring.getAminoAcidBLASTScoring(BLASTMatrix.BLOSUM62, -10, -1));
    SimpleBatchAligner<AminoAcidSequence, Ref> aligner = new SimpleBatchAligner<>(aParams);
    int dbSize = 0;
    for (VDJCLibrary lib : refLibraries) {
        for (VDJCGene gene : lib.getGenes()) {
            NucleotideSequence nSeq = gene.getFeature(geneFeature);
            if (nSeq == null)
                continue;
            ReferencePoint frameReference = GeneFeature.getFrameReference(geneFeature);
            ReferencePoints partitioning = gene.getPartitioning();
            if (frameReference == null)
                continue;
            int relativePosition = partitioning.getRelativePosition(geneFeature, frameReference);
            if (relativePosition < 0)
                continue;
            TranslationParameters frame = withIncompleteCodon(relativePosition);
            AminoAcidSequence aaSequence = AminoAcidSequence.translate(nSeq, frame);
            aligner.addReference(aaSequence, new Ref(gene, frame, nSeq.size()));
            ++dbSize;
        }
    }
    System.out.println("DB size: " + dbSize);
    System.out.println();
    // Checking that db is not empty
    if (dbSize == 0)
        throw new RuntimeException("No reference genes.");
    ArrayList<VDJCLibraryData> result = new ArrayList<>();
    ByteArrayOutputStream bos = new ByteArrayOutputStream();
    PrintStream bufferPS = new PrintStream(bos);
    // Iteration over target genes
    for (VDJCLibrary lib : reg.getLoadedLibrariesByName(TARGET_LIBRARY_NAME)) {
        ArrayList<VDJCGeneData> genes = new ArrayList<>();
        for (VDJCGene targetGene : lib.getGenes()) {
            bos.reset();
            PrintStream ps = params.outputOnlyModified() ? bufferPS : System.out;
            if (namePattern != null && !namePattern.matcher(targetGene.getName()).matches()) {
                if (!params.outputOnlyModified())
                    genes.add(targetGene.getData());
                continue;
            }
            ps.println("Processing: " + targetGene.getName() + " (" + (targetGene.isFunctional() ? "F" : "P") + ") " + targetGene.getChains());
            // Getting gene feature sequence from target gene
            NucleotideSequence nSeq = targetGene.getFeature(geneFeature);
            if (nSeq == null) {
                ps.println("Failed to extract " + GeneFeature.encode(geneFeature));
                ps.println("================");
                ps.println();
                if (!params.outputOnlyModified())
                    genes.add(targetGene.getData());
                continue;
            }
            // Alignment result
            AlignmentResult<AlignmentHit<AminoAcidSequence, Ref>> bestAResult = null;
            TranslationParameters bestFrame = null;
            // Searching for best alignment
            for (TranslationParameters frame : TRANSLATION_PARAMETERS) {
                AminoAcidSequence aaSeq = AminoAcidSequence.translate(nSeq, frame);
                AlignmentResult<AlignmentHit<AminoAcidSequence, Ref>> r = aligner.align(aaSeq);
                if (r != null && r.hasHits() && (bestAResult == null || bestAResult.getBestHit().getAlignment().getScore() < r.getBestHit().getAlignment().getScore())) {
                    bestAResult = r;
                    bestFrame = frame;
                }
            }
            if (bestFrame == null) {
                ps.println("No alignments found.");
                if (!params.outputOnlyModified())
                    genes.add(targetGene.getData());
                continue;
            }
            List<AlignmentHit<AminoAcidSequence, Ref>> hits = bestAResult.getHits();
            VDJCGeneData targetGeneData = targetGene.getData().clone();
            boolean anyPointChanged = false;
            for (int ai = 0; ai < hits.size(); ai++) {
                // Accumulate output
                ByteArrayOutputStream localBos = new ByteArrayOutputStream();
                PrintStream localPS = new PrintStream(localBos);
                Alignment<AminoAcidSequence> bestAlignment = hits.get(ai).getAlignment();
                Ref bestRef = hits.get(ai).getRecordPayload();
                VDJCGene bestReferenceGene = bestRef.gene;
                localPS.println("Aligned with " + bestReferenceGene.getName() + " from " + libraryNameToAddress.get(bestReferenceGene.getParentLibrary().getName()) + " ; Score = " + bestAlignment.getScore());
                AlignmentHelper alignmentHelper = bestAlignment.getAlignmentHelper();
                for (AlignmentHelper h : alignmentHelper.split(150)) localPS.println(h + "\n");
                ReferencePoints targetPartitioning = targetGene.getPartitioning();
                ReferencePoints referencePartitioning = bestReferenceGene.getPartitioning();
                for (GeneFeature.ReferenceRange range : geneFeature) for (ReferencePoint point : range.getIntermediatePoints()) {
                    localPS.print(point + ": ");
                    boolean isAvailable = targetPartitioning.isAvailable(point);
                    if (isAvailable) {
                        localPS.println("already set");
                    }
                    if (!referencePartitioning.isAvailable(point)) {
                        if (!isAvailable)
                            localPS.println("not set in reference gene");
                        continue;
                    }
                    int ntPositionInReference = referencePartitioning.getRelativePosition(geneFeature, point);
                    // Projecting position
                    AminoAcidSequence.AminoAcidSequencePosition aaPositionInReferece = AminoAcidSequence.convertNtPositionToAA(ntPositionInReference, bestRef.ntSeqLength, bestRef.frame);
                    if (aaPositionInReferece == null) {
                        if (!isAvailable)
                            localPS.println("failed to convert to aa position in ref");
                        continue;
                    }
                    int aaPositionInTarget = Alignment.aabs(bestAlignment.convertToSeq2Position(aaPositionInReferece.aminoAcidPosition));
                    if (aaPositionInTarget == -1) {
                        if (!isAvailable)
                            localPS.println("failed to project using alignment");
                        continue;
                    }
                    int ntPositionInTarget = AminoAcidSequence.convertAAPositionToNt(aaPositionInTarget, nSeq.size(), bestFrame);
                    if (ntPositionInTarget == -1) {
                        if (!isAvailable)
                            localPS.println("failed");
                        continue;
                    }
                    ntPositionInTarget += aaPositionInReferece.positionInTriplet;
                    ntPositionInTarget = targetPartitioning.getAbsolutePosition(geneFeature, ntPositionInTarget);
                    if (ntPositionInTarget == -1) {
                        if (!isAvailable)
                            localPS.println("failed");
                        continue;
                    }
                    if (isAvailable) {
                        int existingPosition = targetPartitioning.getPosition(point);
                        if (existingPosition != ntPositionInTarget) {
                            localPS.println("inferred position differs from existing value.  existing: " + existingPosition + ", inferred: " + ntPositionInTarget);
                        }
                        continue;
                    }
                    localPS.println(ntPositionInTarget);
                    targetGeneData.getAnchorPoints().put(point, (long) ntPositionInTarget);
                    anyPointChanged = true;
                }
                if (!anyPointChanged) {
                    if (!params.outputOnlyModified() && ai == 0)
                        ps.write(localBos.toByteArray());
                } else {
                    ps.write(localBos.toByteArray());
                    break;
                }
            }
            ps.println("================");
            ps.println();
            if (anyPointChanged && params.outputOnlyModified())
                System.out.write(bos.toByteArray());
            genes.add(targetGeneData);
        }
        result.add(new VDJCLibraryData(lib.getData(), genes));
    }
    VDJCDataUtils.writeToFile(result, outputFile, false);
}

Example 3 with TranslationParameters

use of com.milaboratory.core.sequence.TranslationParameters in project repseqio by repseqio.

the class FromFastaActionAbstract method go.

@Override
public void go(ActionHelper helper) throws Exception {
    Pattern functionalityRegexp = params.getFunctionalityRegexp();
    GeneType geneType = params.getGeneType();
    Map<String, VDJCGeneData> genes = new HashMap<>();
    Path libraryPath = Paths.get(params.getOutputJSON()).toAbsolutePath();
    // Parsing -P or --gene-feature parameters
    Map<ReferencePoint, Integer> points = new HashMap<>();
    if (params.geneFeature != null) {
        GeneFeature gf = params.getGeneFeature();
        points.put(gf.getFirstPoint(), 0);
        points.put(gf.getLastPoint(), -1);
    } else
        for (Map.Entry<String, String> p : params.points.entrySet()) {
            ReferencePoint anchorPoint = ReferencePoint.getPointByName(p.getKey());
            if (anchorPoint == null)
                throw new IllegalArgumentException("Unknown anchor point: " + p.getKey());
            int position = Integer.decode(p.getValue());
            points.put(anchorPoint, position);
        }
    // Check
    for (Map.Entry<ReferencePoint, Integer> entry : points.entrySet()) if (entry.getKey().getGeneType() != null && entry.getKey().getGeneType() != geneType)
        throw new IllegalArgumentException("Incompatible anchor point and gene type: " + entry.getKey() + " / " + geneType);
    try (FastaReader reader = new FastaReader<>(params.getInput(), null);
        SequenceStorage storage = params.doEmbedSequences() ? new EmbeddedWriter() : params.getOutputFasta() == null ? new ExistingFileWriter(libraryPath, Paths.get(params.getInput()).toAbsolutePath()) : new FastaSequenceStorage(libraryPath, Paths.get(params.getOutputFasta()).toAbsolutePath())) {
        for (FastaReader.RawFastaRecord record : CUtils.it((OutputPortCloseable<FastaReader.RawFastaRecord>) reader.asRawRecordsPort())) {
            StringWithMapping swm = StringWithMapping.removeSymbol(record.sequence, params.getPaddingCharacter());
            NucleotideSequence seq = new NucleotideSequence(swm.getModifiedString());
            if (seq.containsWildcards()) {
                System.out.println("Sequence dropped because contain wildcards: " + record.description);
                continue;
            }
            String[] fields = record.description.split("\\|");
            String geneName = fields[params.nameIndex];
            boolean functionality = true;
            if (params.functionalityIndex != null)
                functionality = functionalityRegexp.matcher(fields[params.functionalityIndex]).matches();
            SortedMap<ReferencePoint, Long> anchorPoints = new TreeMap<>();
            for (Map.Entry<String, String> p : params.patterns.entrySet()) {
                ReferencePoint anchorPoint = ReferencePoint.getPointByName(p.getKey());
                if (anchorPoint == null)
                    throw new IllegalArgumentException("Unknown anchor point: " + p.getKey());
                if (anchorPoint.getGeneType() != null && anchorPoint.getGeneType() != geneType)
                    throw new IllegalArgumentException("Incompatible anchor point and gene type: " + anchorPoint + " / " + geneType);
                Pattern pattern = Pattern.compile(p.getValue());
                int position = -1;
                for (boolean stops : new boolean[] { false, true }) for (int f = 0; f < 3; f++) {
                    if (position != -1)
                        continue;
                    TranslationParameters tp = TranslationParameters.withoutIncompleteCodon(f);
                    AminoAcidSequence aa = AminoAcidSequence.translate(seq, tp);
                    if (!stops && aa.containStops())
                        continue;
                    String str = aa.toString();
                    Matcher matcher = pattern.matcher(str);
                    if (matcher.find()) {
                        int aaPosition = matcher.start(1);
                        position = AminoAcidSequence.convertAAPositionToNt(aaPosition, seq.size(), tp);
                    }
                }
                if (position == -1)
                    continue;
                anchorPoints.put(anchorPoint, (long) position);
            }
            for (Map.Entry<ReferencePoint, Integer> p : points.entrySet()) {
                // AA patterns have priority over positional anchor points
                if (anchorPoints.containsKey(p.getKey()))
                    continue;
                // Converting position using
                int position = swm.convertPosition(p.getValue());
                // Can't be converted (e.g. position of padding symbol) skipping
                if (position == -1)
                    continue;
                anchorPoints.put(p.getKey(), (long) position);
            }
            if (genes.containsKey(geneName)) {
                if (params.getIgnoreDuplicates()) {
                    System.out.println("Ignored: Duplicate records for " + geneName);
                    continue;
                } else
                    throw new IllegalArgumentException("Duplicate records for " + geneName);
            }
            BaseSequence baseSequence = storage.storeSequence(seq, geneName, record.description);
            VDJCGeneData gene = new VDJCGeneData(baseSequence, geneName, geneType, functionality, new Chains(params.chain), new TreeMap<String, SortedSet<String>>(), anchorPoints).addMetaValue(KnownVDJCGeneMetaFields.COMMENTS, record.description);
            genes.put(geneName, gene);
        }
        VDJCLibraryData library = new VDJCLibraryData(params.taxonId, params.speciesNames, new ArrayList<>(genes.values()), new TreeMap<String, SortedSet<String>>(), storage.getBase()).addMetaValue(KnownVDJCLibraryMetaFields.COMMENTS, "Imported from: " + params.getInput());
        VDJCDataUtils.writeToFile(new VDJCLibraryData[] { library }, params.getOutputJSON(), false);
    }
}