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

Example 16 with NSequenceWithQuality

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

the class MergerTest method mAssert.

public static void mAssert(String originalSeq, String seq1, String seq2, String expectedSequence, String expectedQuality) {
    NucleotideSequence originalSequence = new NucleotideSequence(originalSeq);
    NSequenceWithQuality target1 = new NSequenceWithQuality(seq1, lets('A', seq1.length()));
    NSequenceWithQuality target2 = new NSequenceWithQuality(seq2, lets('B', seq2.length()));
    LinearGapAlignmentScoring<NucleotideSequence> scoring = new LinearGapAlignmentScoring<>(NucleotideSequence.ALPHABET, 4, -5, -9);
    Alignment<NucleotideSequence> alignment1 = Aligner.alignLocalLinear(scoring, originalSequence, target1.getSequence());
    Alignment<NucleotideSequence> alignment2 = Aligner.alignLocalLinear(scoring, originalSequence, target2.getSequence());
    NSequenceWithQuality result = Merger.merge(new Range(0, originalSequence.size()), new Alignment[] { alignment1, alignment2 }, new NSequenceWithQuality[] { target1, target2 });
    Assert.assertEquals(expectedSequence, result.getSequence().toString());
    Assert.assertEquals(expectedQuality, result.getQuality().toString());
}

Example 17 with NSequenceWithQuality

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

the class PartialAlignmentsAssembler method searchOverlaps.

@SuppressWarnings("unchecked")
private OverlapSearchResult searchOverlaps(final VDJCAlignments rightAl, final boolean allowChimeras) {
    final Chains jChains = rightAl.getAllChains(GeneType.Joining);
    int rightTargetId = getRightPartitionedSequence(rightAl);
    if (rightTargetId == -1)
        return null;
    rightParts.incrementAndGet();
    final VDJCPartitionedSequence rightTarget = rightAl.getPartitionedTarget(rightTargetId);
    NSequenceWithQuality rightSeqQ = rightTarget.getSequence();
    NucleotideSequence rightSeq = rightSeqQ.getSequence();
    int stop = rightTarget.getPartitioning().getPosition(ReferencePoint.JBeginTrimmed);
    if (stop == -1)
        stop = rightTarget.getSequence().size();
    else
        stop -= kOffset;
    // black list of left parts failed due to inconsistent overlapped alignments (failed AMerge)
    TLongHashSet blackList = new TLongHashSet();
    SEARCH_LEFT_PARTS: while (true) {
        int maxOverlap = -1, maxDelta = -1, maxOverlapIndexInList = -1, maxBegin = -1, maxEnd = -1;
        List<KMerInfo> maxOverlapList = null;
        boolean isMaxOverOverlapped = false;
        for (int rFrom = 0; rFrom < stop && rFrom + kValue < rightSeqQ.size(); rFrom++) {
            long kMer = kMer(rightSeqQ.getSequence(), rFrom, kValue);
            List<KMerInfo> match = kToIndexLeft.get(kMer);
            if (match == null)
                continue;
            out: for (int i = 0; i < match.size(); i++) {
                boolean isOverOverlapped = false;
                final VDJCAlignments leftAl = match.get(i).getAlignments();
                if (blackList.contains(leftAl.getAlignmentsIndex()))
                    continue;
                if (// You shall not merge with yourself
                leftAl.getAlignmentsIndex() == rightAl.getAlignmentsIndex() || alreadyMergedIds.contains(leftAl.getAlignmentsIndex()))
                    continue;
                // Checking chains compatibility
                if (!allowChimeras && jChains != null && !leftAl.getAllChains(GeneType.Variable).intersects(jChains))
                    continue;
                // Check for the same V
                if (leftAl.getBestHit(GeneType.Variable) != null && rightAl.getBestHit(GeneType.Variable) != null && !leftAl.hasCommonGenes(GeneType.Variable, rightAl))
                    continue;
                final NucleotideSequence leftSeq = leftAl.getPartitionedTarget(getLeftPartitionedSequence(leftAl)).getSequence().getSequence();
                int lFrom = match.get(i).kMerPositionFrom;
                // begin and end in the coordinates of left
                int delta, begin = delta = lFrom - rFrom;
                if (begin < 0) {
                    begin = 0;
                    isOverOverlapped = true;
                }
                int end = leftSeq.size();
                if (end >= rightSeq.size() + delta) {
                    end = rightSeq.size() + delta;
                    isOverOverlapped = true;
                }
                for (int j = begin; j < end; j++) if (leftSeq.codeAt(j) != rightSeq.codeAt(j - delta))
                    continue out;
                int overlap = end - begin;
                if (maxOverlap < overlap) {
                    maxOverlap = overlap;
                    maxOverlapList = match;
                    maxOverlapIndexInList = i;
                    maxDelta = delta;
                    maxBegin = begin;
                    maxEnd = end;
                    isMaxOverOverlapped = isOverOverlapped;
                }
            }
        }
        if (maxOverlapList == null)
            return null;
        if (maxOverlap < minimalAssembleOverlap)
            return null;
        final KMerInfo left = maxOverlapList.remove(maxOverlapIndexInList);
        VDJCAlignments leftAl = left.alignments;
        // final long readId = rightAl.getReadId();
        ArrayList<AlignedTarget> leftTargets = extractAlignedTargets(leftAl);
        ArrayList<AlignedTarget> rightTargets = extractAlignedTargets(rightAl);
        AlignedTarget leftCentral = leftTargets.get(left.targetId);
        AlignedTarget rightCentral = rightTargets.get(rightTargetId);
        AlignedTarget central = targetMerger.merge(leftCentral, rightCentral, maxDelta, SequenceHistory.OverlapType.CDR3Overlap, 0);
        // Setting overlap position
        central = AlignedTarget.setOverlapRange(central, maxBegin, maxEnd);
        final List<AlignedTarget> leftDescriptors = new ArrayList<>(2), rightDescriptors = new ArrayList<>(2);
        for (int i = 0; i < left.targetId; ++i) leftDescriptors.add(leftTargets.get(i));
        for (int i = left.targetId + 1; i < leftAl.numberOfTargets(); ++i) rightDescriptors.add(leftTargets.get(i));
        for (int i = 0; i < rightTargetId; ++i) leftDescriptors.add(rightTargets.get(i));
        for (int i = rightTargetId + 1; i < rightAl.numberOfTargets(); ++i) rightDescriptors.add(rightTargets.get(i));
        // Merging to VJ junction
        List<AlignedTarget>[] allDescriptors = new List[] { leftDescriptors, rightDescriptors };
        TargetMerger.TargetMergingResult bestResult = TargetMerger.FAILED_RESULT;
        int bestI;
        // Trying to merge left and right reads to central one
        for (List<AlignedTarget> descriptors : allDescriptors) do {
            bestI = -1;
            for (int i = 0; i < descriptors.size(); i++) {
                TargetMerger.TargetMergingResult result = targetMerger.merge(descriptors.get(i), central);
                if (result.failedDueInconsistentAlignments()) {
                    // Inconsistent alignments -> retry
                    blackList.add(leftAl.getAlignmentsIndex());
                    continue SEARCH_LEFT_PARTS;
                }
                if (bestResult.getScore() < result.getScore()) {
                    bestResult = result;
                    bestI = i;
                }
            }
            if (bestI != -1) {
                assert bestResult != TargetMerger.FAILED_RESULT;
                central = bestResult.getResult();
                descriptors.remove(bestI);
            }
        } while (bestI != -1);
        // Merging left+left / right+right
        outer: for (int d = 0; d < allDescriptors.length; d++) {
            List<AlignedTarget> descriptors = allDescriptors[d];
            for (int i = 0; i < descriptors.size(); i++) for (int j = i + 1; j < descriptors.size(); j++) {
                TargetMerger.TargetMergingResult result = targetMerger.merge(descriptors.get(i), descriptors.get(j));
                if (result.failedDueInconsistentAlignments()) {
                    // Inconsistent alignments -> retry
                    blackList.add(leftAl.getAlignmentsIndex());
                    continue SEARCH_LEFT_PARTS;
                }
                if (result.isSuccessful()) {
                    descriptors.set(i, result.getResult());
                    descriptors.remove(j);
                    --d;
                    continue outer;
                }
            }
        }
        if (isMaxOverOverlapped)
            overoverlapped.incrementAndGet();
        // Creating pre-list of resulting targets
        List<AlignedTarget> result = new ArrayList<>();
        result.addAll(leftDescriptors);
        result.add(central);
        result.addAll(rightDescriptors);
        // Ordering and filtering final targets
        return new OverlapSearchResult(maxOverlapList, left, AlignedTarget.orderTargets(result));
    }
}

Example 18 with NSequenceWithQuality

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

the class VDJCAlignmentsFormatter method getTargetAsMultiAlignment.

public static MultiAlignmentHelper getTargetAsMultiAlignment(VDJCObject vdjcObject, int targetId, boolean addHitScore, boolean addReads) {
    if (addReads && !(vdjcObject instanceof VDJCAlignments))
        throw new IllegalArgumentException("Read alignments supported only for VDJCAlignments.");
    NSequenceWithQuality target = vdjcObject.getTarget(targetId);
    NucleotideSequence targetSeq = target.getSequence();
    SequencePartitioning partitioning = vdjcObject.getPartitionedTarget(targetId).getPartitioning();
    List<Alignment<NucleotideSequence>> alignments = new ArrayList<>();
    List<String> alignmentLeftComments = new ArrayList<>();
    List<String> alignmentRightComments = new ArrayList<>();
    for (GeneType gt : GeneType.values()) for (VDJCHit hit : vdjcObject.getHits(gt)) {
        Alignment<NucleotideSequence> alignment = hit.getAlignment(targetId);
        if (alignment == null)
            continue;
        alignment = alignment.invert(targetSeq);
        alignments.add(alignment);
        alignmentLeftComments.add(hit.getGene().getName());
        alignmentRightComments.add(" " + (int) (hit.getAlignment(targetId).getScore()) + (addHitScore ? " (" + (int) (hit.getScore()) + ")" : ""));
    }
    // Adding read information
    if (addReads) {
        VDJCAlignments vdjcAlignments = (VDJCAlignments) vdjcObject;
        SequenceHistory history = vdjcAlignments.getHistory(targetId);
        List<SequenceHistory.RawSequence> reads = history.rawReads();
        for (SequenceHistory.RawSequence read : reads) {
            NucleotideSequence seq = vdjcAlignments.getOriginalSequence(read.index).getSequence();
            int offset = history.offset(read.index);
            Alignment<NucleotideSequence> alignment = Aligner.alignOnlySubstitutions(targetSeq, seq, offset, seq.size(), 0, seq.size(), AffineGapAlignmentScoring.IGBLAST_NUCLEOTIDE_SCORING);
            alignments.add(alignment);
            alignmentLeftComments.add(read.index.toString());
            alignmentRightComments.add("");
        }
    }
    MultiAlignmentHelper helper = MultiAlignmentHelper.build(MultiAlignmentHelper.DEFAULT_SETTINGS, new Range(0, target.size()), targetSeq, alignments.toArray(new Alignment[alignments.size()]));
    if (!alignments.isEmpty())
        drawPoints(helper, partitioning, POINTS_FOR_REARRANGED);
    drawAASequence(helper, partitioning, targetSeq);
    helper.addSubjectQuality("Quality", target.getQuality());
    helper.setSubjectLeftTitle("Target" + targetId);
    helper.setSubjectRightTitle(" Score" + (addHitScore ? " (hit score)" : ""));
    for (int i = 0; i < alignmentLeftComments.size(); i++) {
        helper.setQueryLeftTitle(i, alignmentLeftComments.get(i));
        helper.setQueryRightTitle(i, alignmentRightComments.get(i));
    }
    return helper;
}

Example 19 with NSequenceWithQuality

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

the class VDJCAlignments method getOriginalSequence.

public NSequenceWithQuality getOriginalSequence(SequenceHistory.FullReadIndex index) {
    if (originalReads == null)
        throw new IllegalStateException("Original reads are not saved for the alignment object.");
    SequenceRead read = null;
    for (SequenceRead originalRead : originalReads) if (originalRead.getId() == index.readId) {
        read = originalRead;
        break;
    }
    if (read == null)
        throw new IllegalArgumentException("No such read index: " + index);
    NSequenceWithQuality seq = read.getRead(index.mateIndex).getData();
    return index.isReverseComplement ? seq.getReverseComplement() : seq;
}

Example 20 with NSequenceWithQuality

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

the class VDJCObject method getTargetContainingFeature.

public final int getTargetContainingFeature(GeneFeature feature) {
    NSequenceWithQuality tmp;
    int targetIndex = -1, quality = -1;
    for (int i = 0; i < targets.length; ++i) {
        tmp = getPartitionedTarget(i).getFeature(feature);
        if (tmp != null && quality < tmp.getQuality().minValue())
            targetIndex = i;
    }
    return targetIndex;
}