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

Example 6 with Alignment

use of com.milaboratory.core.alignment.Alignment in project mixcr by milaboratory.

the class Merger method merge.

public static NSequenceWithQuality merge(Range range, Alignment<NucleotideSequence>[] alignments, NSequenceWithQuality[] targets) {
    // Checking arguments
    if (alignments.length != targets.length)
        throw new IllegalArgumentException();
    // Extracting reference sequence
    NucleotideSequence sequence = alignments[0].getSequence1();
    for (int i = 1; i < alignments.length; i++) if (sequence != alignments[i].getSequence1())
        throw new IllegalArgumentException("Different reference sequences.");
    int position = range.getFrom();
    int localPosition = 0;
    int[] mPointers = new int[alignments.length];
    for (int i = 0; i < alignments.length; i++) mPointers[i] = firstIndexAfter(alignments[i].getAbsoluteMutations(), position);
    // Aggregator of quality information
    SequenceQualityBuilder qualityBuilder = new SequenceQualityBuilder().ensureCapacity(range.length());
    // Aggregator of mutations
    MutationsBuilder<NucleotideSequence> mutationsBuilder = new MutationsBuilder<>(NucleotideSequence.ALPHABET);
    do {
        int winnerIndex = -1;
        byte bestQuality = -1;
        int winnerMPointer = -1, winnerMLength = -1;
        for (int i = 0; i < alignments.length; i++) {
            Alignment<NucleotideSequence> al = alignments[i];
            if (!al.getSequence1Range().contains(position))
                continue;
            // Current mutation index
            int mPointer = mPointers[i];
            Mutations<NucleotideSequence> mutations = al.getAbsoluteMutations();
            // Number of mutations with the same position
            int length;
            if (mPointer < mutations.size() && mutations.getPositionByIndex(mPointer) == position) {
                length = 1;
                while (mPointer + length < mutations.size() && mutations.getPositionByIndex(mPointer + length) == position) ++length;
            } else
                length = 0;
            byte pointQuality = getQuality(position, mPointer, length, al, targets[i].getQuality());
            if (bestQuality < pointQuality) {
                winnerIndex = i;
                winnerMPointer = mPointer;
                winnerMLength = length;
                bestQuality = pointQuality;
            }
        }
        Mutations<NucleotideSequence> winnerMutations = alignments[winnerIndex].getAbsoluteMutations();
        byte sumQuality = 0;
        // Second pass to calculate score
        OUT: for (int i = 0; i < alignments.length; i++) {
            Alignment<NucleotideSequence> al = alignments[i];
            if (!al.getSequence1Range().contains(position))
                continue;
            // Current mutation index
            int mPointer = mPointers[i];
            Mutations<NucleotideSequence> mutations = al.getAbsoluteMutations();
            // Number of mutations with the same position
            int length;
            if (mPointer < mutations.size() && mutations.getPositionByIndex(mPointer) == position) {
                length = 1;
                while (mPointer + length < mutations.size() && mutations.getPositionByIndex(mPointer + length) == position) ++length;
            } else
                length = 0;
            // Advancing pointer
            mPointers[i] += length;
            if (length != winnerMLength)
                continue;
            for (int k = 0; k < length; ++k) if (mutations.getMutation(mPointer + k) != winnerMutations.getMutation(winnerMPointer + k))
                continue OUT;
            byte pointQuality = getQuality(position, mPointer, length, al, targets[i].getQuality());
            sumQuality += pointQuality;
        }
        // todo ??
        if (winnerIndex == -1)
            return null;
        sumQuality = (byte) Math.min(sumQuality, MAX_QUALITY);
        qualityBuilder.append(sumQuality);
        if (winnerMLength != 0)
            mutationsBuilder.append(winnerMutations, winnerMPointer, winnerMLength);
        ++localPosition;
    } while (++position < range.getTo());
    Mutations<NucleotideSequence> mutations = mutationsBuilder.createAndDestroy();
    NSequenceWithQuality toMutate = new NSequenceWithQuality(sequence.getRange(range), qualityBuilder.createAndDestroy());
    return MutationsUtil.mutate(toMutate, mutations.move(-range.getFrom()));
}