Search in sources :

Example 16 with VDJCHit

use of com.milaboratory.mixcr.basictypes.VDJCHit in project mixcr by milaboratory.

the class ActionExportParameters method getFilter.

@SuppressWarnings("unchecked")
public Filter<T> getFilter() {
    List<Filter<T>> filters = new ArrayList<>();
    final Chains chains = getChains();
    filters.add(new Filter<T>() {

        @Override
        public boolean accept(T object) {
            for (GeneType gt : GeneType.VJC_REFERENCE) {
                VDJCHit bestHit = object.getBestHit(gt);
                if (bestHit != null && chains.intersects(bestHit.getGene().getChains()))
                    return true;
            }
            return false;
        }
    });
    if (filters.isEmpty())
        return ACCEPT_ALL;
    if (filters.size() == 1)
        return filters.get(0);
    return and(filters.toArray(new Filter[filters.size()]));
}
Also used : Filter(cc.redberry.primitives.Filter) Chains(io.repseq.core.Chains) GeneType(io.repseq.core.GeneType) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit)

Example 17 with VDJCHit

use of com.milaboratory.mixcr.basictypes.VDJCHit in project mixcr by milaboratory.

the class VDJCAlignerPVFirst method combine.

@SuppressWarnings("unchecked")
public static VDJCHit[] combine(final GeneFeature feature, final AlignmentHit<NucleotideSequence, VDJCGene>[][] hits) {
    for (int i = 0; i < hits.length; i++) Arrays.sort(hits[i], GENE_ID_COMPARATOR);
    ArrayList<VDJCHit> result = new ArrayList<>();
    // Sort-join-like algorithm
    int i;
    VDJCGene minGene;
    Alignment<NucleotideSequence>[] alignments;
    final int[] pointers = new int[hits.length];
    while (true) {
        minGene = null;
        for (i = 0; i < pointers.length; ++i) if (pointers[i] < hits[i].length && (minGene == null || minGene.getId().compareTo(hits[i][pointers[i]].getRecordPayload().getId()) > 0))
            minGene = hits[i][pointers[i]].getRecordPayload();
        // All pointers > hits.length
        if (minGene == null)
            break;
        // Collecting alignments for minAllele
        alignments = new Alignment[hits.length];
        for (i = 0; i < pointers.length; ++i) if (pointers[i] < hits[i].length && minGene == hits[i][pointers[i]].getRecordPayload()) {
            alignments[i] = hits[i][pointers[i]].getAlignment();
            ++pointers[i];
        }
        // Collecting results
        result.add(new VDJCHit(minGene, alignments, feature));
    }
    VDJCHit[] vdjcHits = result.toArray(new VDJCHit[result.size()]);
    Arrays.sort(vdjcHits);
    return vdjcHits;
}
Also used : VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit)

Example 18 with VDJCHit

use of com.milaboratory.mixcr.basictypes.VDJCHit in project mixcr by milaboratory.

the class VDJCAlignerPVFirst method processPartial.

private VDJCAlignmentResult<PairedRead> processPartial(PairedRead input, PAlignmentHelper[] helpers) {
    // Calculates which PTarget was aligned with the highest score
    helpers[0].performCDAlignment();
    PAlignmentHelper bestHelper = helpers[0];
    for (int i = 1; i < helpers.length; ++i) {
        helpers[i].performCDAlignment();
        if (bestHelper.score() < helpers[i].score())
            bestHelper = helpers[i];
    }
    if (!bestHelper.hasVOrJHits()) {
        onFailedAlignment(input, VDJCAlignmentFailCause.NoHits);
        return new VDJCAlignmentResult<>(input);
    }
    // Calculates if this score is bigger then the threshold
    if (bestHelper.score() < parameters.getMinSumScore()) {
        onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
        return new VDJCAlignmentResult<>(input);
    }
    // Finally filtering hits inside this helper to meet minSumScore and maxHits limits
    bestHelper.filterHits(parameters.getMinSumScore(), parameters.getMaxHits());
    // TODO do we really need this ?
    if (!bestHelper.hasVOrJHits()) {
        onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
        return new VDJCAlignmentResult<>(input);
    }
    VDJCAlignments alignments = bestHelper.createResult(input.getId(), this, input);
    // Final check
    if (!parameters.getAllowNoCDR3PartAlignments()) {
        // CDR3 Begin / End
        boolean containCDR3Parts = false;
        final VDJCHit bestV = alignments.getBestHit(GeneType.Variable);
        final VDJCHit bestJ = alignments.getBestHit(GeneType.Joining);
        for (int i = 0; i < 2; i++) {
            if ((bestV != null && bestV.getAlignment(i) != null && bestV.getAlignment(i).getSequence1Range().getTo() >= bestV.getGene().getPartitioning().getRelativePosition(parameters.getFeatureToAlign(GeneType.Variable), reqPointL)) || (bestJ != null && bestJ.getAlignment(i) != null && bestJ.getAlignment(i).getSequence1Range().getFrom() <= bestJ.getGene().getPartitioning().getRelativePosition(parameters.getFeatureToAlign(GeneType.Joining), reqPointR))) {
                containCDR3Parts = true;
                break;
            }
        }
        if (!containCDR3Parts) {
            onFailedAlignment(input, VDJCAlignmentFailCause.NoCDR3Parts);
            return new VDJCAlignmentResult<>(input);
        }
    }
    // Read successfully aligned
    onSuccessfulAlignment(input, alignments);
    if (bestHelper.vChimera)
        onSegmentChimeraDetected(GeneType.Variable, input, alignments);
    if (bestHelper.jChimera)
        onSegmentChimeraDetected(GeneType.Joining, input, alignments);
    return new VDJCAlignmentResult<>(input, alignments);
}
Also used : VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit)

Example 19 with VDJCHit

use of com.milaboratory.mixcr.basictypes.VDJCHit in project mixcr by milaboratory.

the class CloneAccumulator method accumulate.

public synchronized void accumulate(ClonalSequence data, VDJCAlignments alignment, boolean mapped) {
    if (!mapped) {
        // Core sequence accumulation
        coreCount += alignment.getNumberOfReads();
        // Accumulate information about V-D-J alignments only for strictly clustered reads
        // (only for core clonotypes members)
        float score;
        // Accumulate information about all genes
        for (GeneType geneType : GeneType.VJC_REFERENCE) {
            TObjectFloatHashMap<VDJCGeneId> geneScores = this.geneScores.get(geneType);
            VDJCHit[] hits = alignment.getHits(geneType);
            if (hits.length == 0)
                continue;
            if (geneScores == null)
                this.geneScores.put(geneType, geneScores = new TObjectFloatHashMap<>());
            for (VDJCHit hit : hits) {
                // Calculating sum of natural logarithms of scores
                score = hit.getScore();
                geneScores.adjustOrPutValue(hit.getGene().getId(), score, score);
            }
        }
        aggregator.aggregate(data.getConcatenated().getQuality());
    } else
        // Mapped sequence accumulation
        mappedCount += alignment.getNumberOfReads();
}
Also used : GeneType(io.repseq.core.GeneType) VDJCGeneId(io.repseq.core.VDJCGeneId) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit)

Example 20 with VDJCHit

use of com.milaboratory.mixcr.basictypes.VDJCHit in project mixcr by milaboratory.

the class CloneFactory method create.

public Clone create(int id, double count, EnumMap<GeneType, TObjectFloatHashMap<VDJCGeneId>> geneScores, NSequenceWithQuality[] targets) {
    EnumMap<GeneType, VDJCHit[]> hits = new EnumMap<>(GeneType.class);
    for (GeneType geneType : GeneType.VJC_REFERENCE) {
        VJCClonalAlignerParameters vjcParameters = parameters.getVJCParameters(geneType);
        if (vjcParameters == null)
            continue;
        GeneFeature featureToAlign = featuresToAlign.get(geneType);
        TObjectFloatHashMap<VDJCGeneId> gtGeneScores = geneScores.get(geneType);
        if (gtGeneScores == null)
            continue;
        GeneFeature[] intersectingFeatures = new GeneFeature[assemblingFeatures.length];
        for (int i = 0; i < assemblingFeatures.length; ++i) {
            intersectingFeatures[i] = GeneFeature.intersection(featureToAlign, assemblingFeatures[i]);
            if (intersectingFeatures[i] != null)
                switch(geneType) {
                    case Variable:
                        if (!intersectingFeatures[i].getFirstPoint().equals(assemblingFeatures[i].getFirstPoint()))
                            throw new IllegalArgumentException();
                        break;
                    case Joining:
                        if (!intersectingFeatures[i].getLastPoint().equals(assemblingFeatures[i].getLastPoint()))
                            throw new IllegalArgumentException();
                        break;
                }
        }
        VDJCHit[] result = new VDJCHit[gtGeneScores.size()];
        int pointer = 0;
        TObjectFloatIterator<VDJCGeneId> iterator = gtGeneScores.iterator();
        while (iterator.hasNext()) {
            iterator.advance();
            VDJCGene gene = usedGenes.get(iterator.key());
            Alignment<NucleotideSequence>[] alignments = new Alignment[assemblingFeatures.length];
            for (int i = 0; i < assemblingFeatures.length; ++i) {
                if (intersectingFeatures[i] == null)
                    continue;
                NucleotideSequence referenceSequence = gene.getFeature(featureToAlign);
                Range rangeInReference = gene.getPartitioning().getRelativeRange(featureToAlign, intersectingFeatures[i]);
                if (rangeInReference == null || referenceSequence == null)
                    continue;
                Boolean leftSide;
                switch(geneType) {
                    case Variable:
                        leftSide = intersectingFeatures[i].getLastPoint().isTrimmable() ? true : null;
                        break;
                    case Joining:
                        leftSide = intersectingFeatures[i].getFirstPoint().isTrimmable() ? false : null;
                        break;
                    case Constant:
                        leftSide = null;
                        break;
                    default:
                        throw new RuntimeException();
                }
                BandedAlignerParameters<NucleotideSequence> alignmentParameters = vjcParameters.getAlignmentParameters();
                int referenceLength = rangeInReference.length();
                NucleotideSequence target = targets[i].getSequence();
                if (alignmentParameters.getScoring() instanceof LinearGapAlignmentScoring) {
                    if (leftSide == null) {
                        alignments[i] = BandedLinearAligner.align((LinearGapAlignmentScoring<NucleotideSequence>) alignmentParameters.getScoring(), referenceSequence, target, rangeInReference.getFrom(), referenceLength, 0, target.size(), alignmentParameters.getWidth());
                    } else if (leftSide) {
                        assert rangeInReference.getFrom() + referenceLength == referenceSequence.size();
                        alignments[i] = BandedLinearAligner.alignSemiLocalLeft((LinearGapAlignmentScoring<NucleotideSequence>) alignmentParameters.getScoring(), referenceSequence, target, rangeInReference.getFrom(), referenceLength, 0, target.size(), alignmentParameters.getWidth(), alignmentParameters.getStopPenalty());
                    } else {
                        assert rangeInReference.getFrom() == 0;
                        // int offset2 = Math.max(0, target.size() - referenceLength);
                        alignments[i] = BandedLinearAligner.alignSemiLocalRight((LinearGapAlignmentScoring<NucleotideSequence>) alignmentParameters.getScoring(), referenceSequence, target, rangeInReference.getFrom(), referenceLength, 0, target.size(), alignmentParameters.getWidth(), alignmentParameters.getStopPenalty());
                    }
                } else {
                    if (leftSide == null) {
                        alignments[i] = BandedAffineAligner.align((AffineGapAlignmentScoring<NucleotideSequence>) alignmentParameters.getScoring(), referenceSequence, target, rangeInReference.getFrom(), referenceLength, 0, target.size(), alignmentParameters.getWidth());
                    } else if (leftSide) {
                        assert rangeInReference.getFrom() + referenceLength == referenceSequence.size();
                        alignments[i] = BandedAffineAligner.semiLocalRight((AffineGapAlignmentScoring<NucleotideSequence>) alignmentParameters.getScoring(), referenceSequence, target, rangeInReference.getFrom(), referenceLength, 0, target.size(), alignmentParameters.getWidth());
                    } else {
                        assert rangeInReference.getFrom() == 0;
                        // int offset2 = Math.max(0, target.size() - referenceLength);
                        alignments[i] = BandedAffineAligner.semiLocalLeft((AffineGapAlignmentScoring<NucleotideSequence>) alignmentParameters.getScoring(), referenceSequence, target, rangeInReference.getFrom(), referenceLength, 0, target.size(), alignmentParameters.getWidth());
                    }
                }
            }
            result[pointer++] = new VDJCHit(gene, alignments, featureToAlign, iterator.value());
        }
        Arrays.sort(result, 0, pointer);
        hits.put(geneType, pointer < result.length ? Arrays.copyOf(result, pointer) : result);
    }
    // D
    NucleotideSequence sequenceToAlign = targets[indexOfAssemblingFeatureWithD].getSequence();
    int from = 0;
    int to = sequenceToAlign.size();
    VDJCHit[] hs = hits.get(GeneType.Variable);
    if (hs != null && hs.length > 0) {
        int p = hs[0].getPartitioningForTarget(indexOfAssemblingFeatureWithD).getPosition(ReferencePoint.VEndTrimmed);
        if (p != -1) {
            if (p < 0)
                p = -2 - p;
            from = p;
        }
    }
    hs = hits.get(GeneType.Joining);
    if (hs != null && hs.length > 0) {
        int p = hs[0].getPartitioningForTarget(indexOfAssemblingFeatureWithD).getPosition(ReferencePoint.JBeginTrimmed);
        if (p != -1) {
            if (p < 0)
                p = -2 - p;
            to = p;
        }
    }
    if (from < to)
        hits.put(GeneType.Diversity, dAligner.align(sequenceToAlign, VDJCAligner.getPossibleDLoci(hits.get(GeneType.Variable), hits.get(GeneType.Joining)), from, to, indexOfAssemblingFeatureWithD, assemblingFeatures.length));
    else
        hits.put(GeneType.Diversity, new VDJCHit[0]);
    return new Clone(targets, hits, count, id);
}
Also used : Range(com.milaboratory.core.Range) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit) Clone(com.milaboratory.mixcr.basictypes.Clone)

Aggregations

VDJCHit (com.milaboratory.mixcr.basictypes.VDJCHit)22 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)8 VDJCAlignments (com.milaboratory.mixcr.basictypes.VDJCAlignments)8 GeneType (io.repseq.core.GeneType)7 VDJCGene (io.repseq.core.VDJCGene)7 Range (com.milaboratory.core.Range)4 ArrayList (java.util.ArrayList)4 Alignment (com.milaboratory.core.alignment.Alignment)3 NSequenceWithQuality (com.milaboratory.core.sequence.NSequenceWithQuality)3 Clone (com.milaboratory.mixcr.basictypes.Clone)3 ReferencePoint (io.repseq.core.ReferencePoint)3 VDJCGeneId (io.repseq.core.VDJCGeneId)3 PairedRead (com.milaboratory.core.io.sequence.PairedRead)2 PairedFastqReader (com.milaboratory.core.io.sequence.fastq.PairedFastqReader)2 VDJCPartitionedSequence (com.milaboratory.mixcr.basictypes.VDJCPartitionedSequence)2 VDJCAlignerParameters (com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters)2 Chains (io.repseq.core.Chains)2 GeneFeature (io.repseq.core.GeneFeature)2 Test (org.junit.Test)2 CUtils (cc.redberry.pipe.CUtils)1