use of com.milaboratory.core.sequence.NSequenceWithQuality 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()));
}
Aggregations