Examples with AlignmentHelper com.milaboratory.core.alignment.AlignmentHelper used on opensource projects

Example 1 with AlignmentHelper

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

use of com.milaboratory.core.alignment.AlignmentHelper 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;
}

Example 3 with AlignmentHelper

use of com.milaboratory.core.alignment.AlignmentHelper 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);
}