Search in sources :

Example 1 with Range

use of com.milaboratory.core.Range in project mixcr by milaboratory.

the class FullSeqAssembler method buildClone.

/* ================================= Re-align and build final clone ====================================== */
private Clone buildClone(double count, BranchSequences targets) {
    Alignment<NucleotideSequence>[] vTopHitAlignments = new Alignment[targets.ranges.length], jTopHitAlignments = new Alignment[targets.ranges.length];
    VDJCHit hit = clone.getBestHit(Variable);
    if (hit == null)
        throw new UnsupportedOperationException("No V hit.");
    NucleotideSequence vTopReferenceSequence = hit.getGene().getFeature(hit.getAlignedFeature());
    hit = clone.getBestHit(Joining);
    if (hit == null)
        throw new UnsupportedOperationException("No J hit.");
    NucleotideSequence jTopReferenceSequence = hit.getGene().getFeature(hit.getAlignedFeature());
    // Excessive optimization
    CachedIntArray cachedIntArray = new CachedIntArray();
    int assemblingFeatureLength = targets.assemblingFeatureLength;
    for (int i = 0; i < targets.ranges.length; i++) {
        Range range = targets.ranges[i];
        NucleotideSequence sequence = targets.sequences[i].getSequence();
        // Asserts
        if (range.getFrom() < N_LEFT_DUMMIES + lengthV && range.getTo() >= N_LEFT_DUMMIES + lengthV && i != targets.assemblingFeatureTargetId)
            throw new RuntimeException();
        if (range.getTo() >= rightAssemblingFeatureBound && range.getFrom() < rightAssemblingFeatureBound && i != targets.assemblingFeatureTargetId)
            throw new RuntimeException();
        if (range.getTo() < N_LEFT_DUMMIES + lengthV) {
            boolean floatingLeftBound = !parameters.alignedRegionsOnly && i == 0 && alignerParameters.getVAlignerParameters().getParameters().isFloatingLeftBound();
            // Can be reduced to a single statement
            if (range.getFrom() < N_LEFT_DUMMIES)
                // This target contain extra non-V nucleotides on the left
                vTopHitAlignments[i] = alignLinearSeq1FromRight(((LinearGapAlignmentScoring<NucleotideSequence>) alignerParameters.getVAlignerParameters().getScoring()), vTopReferenceSequence, sequence.getSequence(), 0, range.getTo() - N_LEFT_DUMMIES, 0, sequence.size(), !floatingLeftBound, cachedIntArray);
            else if (floatingLeftBound)
                vTopHitAlignments[i] = alignLinearSeq1FromRight(((LinearGapAlignmentScoring<NucleotideSequence>) alignerParameters.getVAlignerParameters().getScoring()), vTopReferenceSequence, sequence.getSequence(), range.getFrom() - N_LEFT_DUMMIES, range.length(), 0, sequence.size(), false, cachedIntArray);
            else
                vTopHitAlignments[i] = Aligner.alignGlobal(alignerParameters.getVAlignerParameters().getScoring(), vTopReferenceSequence, sequence, range.getFrom() - N_LEFT_DUMMIES, range.length(), 0, sequence.size());
        } else if (i == targets.assemblingFeatureTargetId) {
            /*
                 *  V gene
                 */
            boolean vFloatingLeftBound = !parameters.alignedRegionsOnly && i == 0 && alignerParameters.getVAlignerParameters().getParameters().isFloatingLeftBound();
            // Can be reduced to a single statement
            if (range.getFrom() < N_LEFT_DUMMIES)
                // This target contain extra non-V nucleotides on the left
                vTopHitAlignments[i] = alignLinearSeq1FromRight(((LinearGapAlignmentScoring<NucleotideSequence>) alignerParameters.getVAlignerParameters().getScoring()), vTopReferenceSequence, sequence.getSequence(), 0, lengthV, 0, targets.assemblingFeatureOffset, !vFloatingLeftBound, cachedIntArray);
            else if (vFloatingLeftBound)
                vTopHitAlignments[i] = alignLinearSeq1FromRight(((LinearGapAlignmentScoring<NucleotideSequence>) alignerParameters.getVAlignerParameters().getScoring()), vTopReferenceSequence, sequence.getSequence(), range.getFrom() - N_LEFT_DUMMIES, lengthV - (range.getFrom() - N_LEFT_DUMMIES), 0, targets.assemblingFeatureOffset, false, cachedIntArray);
            else
                vTopHitAlignments[i] = Aligner.alignGlobal(alignerParameters.getVAlignerParameters().getScoring(), vTopReferenceSequence, sequence, range.getFrom() - N_LEFT_DUMMIES, lengthV - (range.getFrom() - N_LEFT_DUMMIES), 0, targets.assemblingFeatureOffset);
            /*
                 *  J gene
                 */
            boolean jFloatingRightBound = !parameters.alignedRegionsOnly && i == targets.ranges.length - 1 && alignerParameters.getJAlignerParameters().getParameters().isFloatingRightBound();
            if (range.getTo() >= rightAssemblingFeatureBound + jLength)
                // This target contain extra non-J nucleotides on the right
                jTopHitAlignments[i] = alignLinearSeq1FromLeft(((LinearGapAlignmentScoring<NucleotideSequence>) alignerParameters.getJAlignerParameters().getScoring()), jTopReferenceSequence, sequence.getSequence(), jOffset, jLength, targets.assemblingFeatureOffset + assemblingFeatureLength, sequence.size() - (targets.assemblingFeatureOffset + assemblingFeatureLength), !jFloatingRightBound, cachedIntArray);
            else if (jFloatingRightBound)
                jTopHitAlignments[i] = alignLinearSeq1FromLeft(((LinearGapAlignmentScoring<NucleotideSequence>) alignerParameters.getJAlignerParameters().getScoring()), jTopReferenceSequence, sequence.getSequence(), jOffset, range.getTo() - rightAssemblingFeatureBound, targets.assemblingFeatureOffset + assemblingFeatureLength, sequence.size() - (targets.assemblingFeatureOffset + assemblingFeatureLength), false, cachedIntArray);
            else
                jTopHitAlignments[i] = Aligner.alignGlobal(alignerParameters.getJAlignerParameters().getScoring(), jTopReferenceSequence, sequence, jOffset, range.getTo() - rightAssemblingFeatureBound, targets.assemblingFeatureOffset + assemblingFeatureLength, sequence.size() - (targets.assemblingFeatureOffset + assemblingFeatureLength));
        } else if (range.getFrom() > rightAssemblingFeatureBound) {
            boolean floatingRightBound = !parameters.alignedRegionsOnly && i == targets.ranges.length - 1 && alignerParameters.getJAlignerParameters().getParameters().isFloatingRightBound();
            if (range.getTo() >= rightAssemblingFeatureBound + jLength)
                // This target contain extra non-J nucleotides on the right
                jTopHitAlignments[i] = alignLinearSeq1FromLeft(((LinearGapAlignmentScoring<NucleotideSequence>) alignerParameters.getJAlignerParameters().getScoring()), jTopReferenceSequence, sequence.getSequence(), jOffset + (range.getFrom() - rightAssemblingFeatureBound), jLength - (range.getFrom() - rightAssemblingFeatureBound), 0, sequence.size(), !floatingRightBound, cachedIntArray);
            else if (floatingRightBound)
                jTopHitAlignments[i] = alignLinearSeq1FromLeft(((LinearGapAlignmentScoring<NucleotideSequence>) alignerParameters.getJAlignerParameters().getScoring()), jTopReferenceSequence, sequence.getSequence(), jOffset + (range.getFrom() - rightAssemblingFeatureBound), range.length(), 0, sequence.size(), false, cachedIntArray);
            else
                jTopHitAlignments[i] = Aligner.alignGlobal(alignerParameters.getJAlignerParameters().getScoring(), jTopReferenceSequence, sequence, jOffset + (range.getFrom() - rightAssemblingFeatureBound), range.length(), 0, sequence.size());
        } else
            throw new RuntimeException();
    }
    NSequenceWithQuality assemblingFeatureSeq = targets.sequences[targets.assemblingFeatureTargetId].getRange(targets.assemblingFeatureOffset, targets.assemblingFeatureOffset + targets.assemblingFeatureLength);
    Clone clone = cloneFactory.create(0, count, geneScores, new NSequenceWithQuality[] { assemblingFeatureSeq });
    vTopHitAlignments[targets.assemblingFeatureTargetId] = mergeTwoAlignments(vTopHitAlignments[targets.assemblingFeatureTargetId], clone.getBestHit(Variable).getAlignment(0).move(targets.assemblingFeatureOffset));
    jTopHitAlignments[targets.assemblingFeatureTargetId] = mergeTwoAlignments(clone.getBestHit(Joining).getAlignment(0).move(targets.assemblingFeatureOffset), jTopHitAlignments[targets.assemblingFeatureTargetId]);
    EnumMap<GeneType, VDJCHit[]> hits = new EnumMap<>(GeneType.class);
    for (GeneType gt : GeneType.VDJC_REFERENCE) hits.put(gt, Arrays.stream(clone.getHits(gt)).map(h -> moveHitTarget(h, targets.assemblingFeatureTargetId, targets.assemblingFeatureOffset, targets.ranges.length)).toArray(VDJCHit[]::new));
    VDJCHit[] tmp = hits.get(Variable);
    tmp[0] = substituteAlignments(tmp[0], vTopHitAlignments);
    tmp = hits.get(Joining);
    tmp[0] = substituteAlignments(tmp[0], jTopHitAlignments);
    return new Clone(targets.sequences, hits, count, 0);
}
Also used : IntStream(java.util.stream.IntStream) VDJCPartitionedSequence(com.milaboratory.mixcr.basictypes.VDJCPartitionedSequence) TIntObjectHashMap(gnu.trove.map.hash.TIntObjectHashMap) java.util(java.util) TObjectFloatHashMap(gnu.trove.map.hash.TObjectFloatHashMap) Clone(com.milaboratory.mixcr.basictypes.Clone) Supplier(java.util.function.Supplier) TIntIntIterator(gnu.trove.iterator.TIntIntIterator) com.milaboratory.core.alignment(com.milaboratory.core.alignment) VDJCAlignerParameters(com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters) CloneFactory(com.milaboratory.mixcr.assembler.CloneFactory) MutationsBuilder(com.milaboratory.core.mutations.MutationsBuilder) Constants(gnu.trove.impl.Constants) io.repseq.core(io.repseq.core) Range(com.milaboratory.core.Range) TObjectIntHashMap(gnu.trove.map.hash.TObjectIntHashMap) CUtils(cc.redberry.pipe.CUtils) TIntIntHashMap(gnu.trove.map.hash.TIntIntHashMap) OutputPort(cc.redberry.pipe.OutputPort) Collectors(java.util.stream.Collectors) com.milaboratory.core.sequence(com.milaboratory.core.sequence) TIntHashSet(gnu.trove.set.hash.TIntHashSet) Joining(io.repseq.core.GeneType.Joining) Stream(java.util.stream.Stream) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit) VDJCGenes(io.repseq.gen.VDJCGenes) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) Variable(io.repseq.core.GeneType.Variable) Range(com.milaboratory.core.Range) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit) Clone(com.milaboratory.mixcr.basictypes.Clone)

Example 2 with Range

use of com.milaboratory.core.Range in project mixcr by milaboratory.

the class FullSeqAssembler method alignLinearSeq1FromRight.

// fixme write docs
static Alignment<NucleotideSequence> alignLinearSeq1FromRight(LinearGapAlignmentScoring<NucleotideSequence> scoring, NucleotideSequence seq1, NucleotideSequence seq2, int offset1, int length1, int offset2, int length2, boolean global, CachedIntArray cache) {
    MutationsBuilder<NucleotideSequence> mutations = new MutationsBuilder<>(NucleotideSequence.ALPHABET);
    BandedSemiLocalResult result;
    int width = 2 * length1;
    if (global) {
        int seq2added = width;
        if (length2 > length1) {
            length2 = Math.min(length2, length1 + width);
            seq2added = Math.min(length2, width + (length2 - length1));
        }
        result = BandedLinearAligner.alignLeftAdded0(scoring, seq1, seq2, offset1, length1, 0, offset2, length2, seq2added, width, mutations, cache);
    } else
        result = BandedLinearAligner.alignSemiLocalRight0(scoring, seq1, seq2, offset1, length1, offset2, length2, width, Integer.MIN_VALUE, mutations, cache);
    return new Alignment<>(seq1, mutations.createAndDestroy(), new Range(result.sequence1Stop, offset1 + length1), new Range(result.sequence2Stop, offset2 + length2), result.score);
}
Also used : MutationsBuilder(com.milaboratory.core.mutations.MutationsBuilder) Range(com.milaboratory.core.Range)

Example 3 with Range

use of com.milaboratory.core.Range in project mixcr by milaboratory.

the class FullSeqAssembler method alignLinearSeq1FromLeft.

// fixme write docs
static Alignment<NucleotideSequence> alignLinearSeq1FromLeft(LinearGapAlignmentScoring<NucleotideSequence> scoring, NucleotideSequence seq1, NucleotideSequence seq2, int offset1, int length1, int offset2, int length2, boolean global, CachedIntArray cache) {
    MutationsBuilder<NucleotideSequence> mutations = new MutationsBuilder<>(NucleotideSequence.ALPHABET);
    BandedSemiLocalResult result;
    int width = 2 * length1;
    if (global) {
        int seq2added = width;
        if (length2 > length1) {
            length2 = Math.min(length2, length1 + width);
            seq2added = Math.min(length2, width + (length2 - length1));
        }
        result = BandedLinearAligner.alignRightAdded0(scoring, seq1, seq2, offset1, length1, 0, offset2, length2, seq2added, width, mutations, cache);
    } else
        result = BandedLinearAligner.alignSemiLocalLeft0(scoring, seq1, seq2, offset1, length1, offset2, length2, width, Integer.MIN_VALUE, mutations, cache);
    return new Alignment<>(seq1, mutations.createAndDestroy(), new Range(offset1, result.sequence1Stop + 1), new Range(offset2, result.sequence2Stop + 1), result.score);
}
Also used : MutationsBuilder(com.milaboratory.core.mutations.MutationsBuilder) Range(com.milaboratory.core.Range)

Example 4 with Range

use of com.milaboratory.core.Range in project mixcr by milaboratory.

the class FullSeqAssembler method toPointSequences.

List<PointSequence> toPointSequences(VDJCAlignments alignments, int iTarget) {
    // 
    // N_LEFT_DUMMIES     assemblingFeatureLength
    // ------|--------------|--------------|------------------------>
    // ↓              ↓              ↓
    // 0000000vvvvvvvvvvvvvvCDR3CDR3CDR3CDR3jjjjjjjjjjjjjjjjCCCCCCCCC
    VDJCPartitionedSequence target = alignments.getPartitionedTarget(iTarget);
    NSequenceWithQuality targetSeq = alignments.getTarget(iTarget);
    VDJCHit vHit = alignments.getBestHit(Variable);
    Alignment<NucleotideSequence> vAlignment = (vHit == null || vHit.getAlignment(iTarget) == null || !Objects.equals(genes.v, vHit.getGene()) || vHit.getAlignment(iTarget).getSequence1Range().getFrom() > lengthV) ? null : vHit.getAlignment(iTarget);
    VDJCHit jHit = alignments.getBestHit(Joining);
    Alignment<NucleotideSequence> jAlignment = (jHit == null || jHit.getAlignment(iTarget) == null || !Objects.equals(genes.j, jHit.getGene()) || jHit.getAlignment(iTarget).getSequence1Range().getTo() < jOffset) ? null : jHit.getAlignment(iTarget);
    List<PointSequence> points = new ArrayList<>();
    if (target.getPartitioning().isAvailable(assemblingFeature.getFirstPoint())) {
        // This target contains left edge of the assembling feature
        int leftStop = target.getPartitioning().getPosition(assemblingFeature.getFirstPoint());
        if (hasV) {
            if (vAlignment != null)
                toPointSequencesByAlignments(points, vAlignment, targetSeq, new Range(parameters.alignedRegionsOnly ? vAlignment.getSequence2Range().getFrom() : 0, leftStop), N_LEFT_DUMMIES);
        } else if (!parameters.alignedRegionsOnly)
            toPointSequencesNoAlignments(points, targetSeq, new Range(0, leftStop), N_LEFT_DUMMIES - leftStop);
    } else if (hasV && vAlignment != null)
        // This target ends before beginning (left edge) of the assembling feature
        toPointSequencesByAlignments(points, vAlignment, targetSeq, new Range(parameters.alignedRegionsOnly ? vAlignment.getSequence2Range().getFrom() : 0, vAlignment.getSequence2Range().getTo()), N_LEFT_DUMMIES);
    if (target.getPartitioning().isAvailable(assemblingFeature.getLastPoint())) {
        // This target contains right edge of the assembling feature
        int rightStart = target.getPartitioning().getPosition(assemblingFeature.getLastPoint());
        if (hasJ) {
            if (jAlignment != null)
                toPointSequencesByAlignments(points, jAlignment, targetSeq, new Range(rightStart, parameters.alignedRegionsOnly ? jAlignment.getSequence2Range().getTo() : targetSeq.size()), N_LEFT_DUMMIES + lengthV + assemblingFeatureLength - jOffset);
        } else
            toPointSequencesNoAlignments(points, targetSeq, new Range(rightStart, targetSeq.size()), N_LEFT_DUMMIES + lengthV + assemblingFeatureLength - rightStart);
    } else if (hasJ && jAlignment != null)
        // This target starts after the end (right edge) of the assembling feature
        toPointSequencesByAlignments(points, jAlignment, targetSeq, new Range(jAlignment.getSequence2Range().getFrom(), parameters.alignedRegionsOnly ? jAlignment.getSequence2Range().getTo() : targetSeq.size()), N_LEFT_DUMMIES + lengthV + assemblingFeatureLength - jOffset);
    return points;
}
Also used : VDJCPartitionedSequence(com.milaboratory.mixcr.basictypes.VDJCPartitionedSequence) Range(com.milaboratory.core.Range) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit)

Example 5 with Range

use of com.milaboratory.core.Range in project mixcr by milaboratory.

the class FullSeqAssembler method toPointSequencesByAlignments.

void toPointSequencesByAlignments(List<PointSequence> points, Alignment<NucleotideSequence> alignment, NSequenceWithQuality seq2, Range seq2Range, int offset) {
    // if (seq2Range.length() == 0)
    // return;
    alignment = AlignmentUtils.shiftIndelsAtHomopolymers(alignment);
    Range alSeq2Range = alignment.getSequence2Range(), alSeq2RangeIntersection = alSeq2Range.intersectionWithTouch(seq2Range), alSeq1RangeIntersection = convertToSeq1Range(alignment, alSeq2RangeIntersection);
    assert alSeq1RangeIntersection != null;
    int shift;
    // left
    shift = offset + alignment.getSequence1Range().getFrom() - alignment.getSequence2Range().getFrom();
    for (int i = seq2Range.getFrom(); i < alSeq2RangeIntersection.getFrom(); ++i) points.add(createPointSequence(i + shift, seq2, i, i + 1, alignment.getSequence2Range()));
    // central
    for (int i = alSeq1RangeIntersection.getFrom(); i < alSeq1RangeIntersection.getTo(); ++i) points.add(createPointSequence(i + offset, seq2, alignment.convertToSeq2Range(new Range(i, i + 1)), alignment.getSequence2Range()));
    // right
    shift = offset + alignment.getSequence1Range().getTo() - alignment.getSequence2Range().getTo();
    for (int i = alSeq2RangeIntersection.getTo(); i < seq2Range.getTo(); ++i) points.add(createPointSequence(i + shift, seq2, i, i + 1, alignment.getSequence2Range()));
}
Also used : Range(com.milaboratory.core.Range)

Aggregations

Range (com.milaboratory.core.Range)45 Test (org.junit.Test)23 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)19 VDJCHit (com.milaboratory.mixcr.basictypes.VDJCHit)4 ArrayList (java.util.ArrayList)4 SingleFastqReaderTest (com.milaboratory.core.io.sequence.fastq.SingleFastqReaderTest)3 MutationsBuilder (com.milaboratory.core.mutations.MutationsBuilder)3 ReferencePoint (io.repseq.core.ReferencePoint)3 Path (java.nio.file.Path)3 NSequenceWithQuality (com.milaboratory.core.sequence.NSequenceWithQuality)2 Clone (com.milaboratory.mixcr.basictypes.Clone)2 VDJCPartitionedSequence (com.milaboratory.mixcr.basictypes.VDJCPartitionedSequence)2 VDJCAlignerParameters (com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters)2 CUtils (cc.redberry.pipe.CUtils)1 OutputPort (cc.redberry.pipe.OutputPort)1 com.milaboratory.core.alignment (com.milaboratory.core.alignment)1 Alignment (com.milaboratory.core.alignment.Alignment)1 LinearGapAlignmentScoring (com.milaboratory.core.alignment.LinearGapAlignmentScoring)1 MultiAlignmentHelper (com.milaboratory.core.alignment.MultiAlignmentHelper)1 com.milaboratory.core.sequence (com.milaboratory.core.sequence)1