Search in sources :

Example 11 with GeneType

use of io.repseq.core.GeneType in project mixcr by milaboratory.

the class TargetMerger method merge.

/**
 * @param targetLeft       left sequence
 * @param targetRight      right sequence
 * @param trySequenceMerge whether to try merging using sequence overlap (if alignment overlap failed)
 */
public TargetMergingResult merge(AlignedTarget targetLeft, AlignedTarget targetRight, boolean trySequenceMerge) {
    for (GeneType geneType : GeneType.VJC_REFERENCE) {
        if (!hasAlignments(targetLeft, geneType) || !hasAlignments(targetRight, geneType))
            continue;
        List<HitMappingRecord> als = extractSortedHits(targetLeft, targetRight, geneType);
        Alignment<NucleotideSequence>[] topHits = als.get(0).alignments;
        if (topHits[0] != null && topHits[1] != null) {
            final Alignment<NucleotideSequence> left = topHits[0];
            final Alignment<NucleotideSequence> right = topHits[1];
            final int from = Math.max(left.getSequence1Range().getFrom(), right.getSequence1Range().getFrom());
            final int to = Math.min(left.getSequence1Range().getTo(), right.getSequence1Range().getTo());
            if (to <= from)
                continue;
            int delta = left.convertToSeq2Position(from) - right.convertToSeq2Position(from);
            if (delta != left.convertToSeq2Position(to) - right.convertToSeq2Position(to))
                continue;
            int seq1Offset = delta > 0 ? delta : 0;
            int seq2Offset = delta > 0 ? 0 : -delta;
            int overlap = Math.min(targetLeft.getTarget().size() - seq1Offset, targetRight.getTarget().size() - seq2Offset);
            int mismatches = SequencesUtils.mismatchCount(targetLeft.getTarget().getSequence(), seq1Offset, targetRight.getTarget().getSequence(), seq2Offset, overlap);
            double identity = MismatchOnlyPairedReadMerger.identity(identityType, targetLeft.getTarget(), seq1Offset, targetRight.getTarget(), seq2Offset, overlap);
            if (identity < minimalAlignmentMergeIdentity)
                return new TargetMergingResult(geneType);
            final AlignedTarget merge = merge(targetLeft, targetRight, delta, OverlapType.AlignmentOverlap, mismatches);
            return new TargetMergingResult(true, null, merge, PairedReadMergingResult.MATCH_SCORE * (overlap - mismatches) + PairedReadMergingResult.MISMATCH_SCORE * mismatches, overlap, mismatches, delta);
        }
    }
    if (!trySequenceMerge)
        return new TargetMergingResult();
    final PairedReadMergingResult merge = merger.merge(targetLeft.getTarget(), targetRight.getTarget());
    if (!merge.isSuccessful())
        return new TargetMergingResult();
    return new TargetMergingResult(false, null, merge(targetLeft, targetRight, merge.getOffset(), OverlapType.SequenceOverlap, merge.getErrors()), merge.score(), merge.getOverlap(), merge.getErrors(), merge.getOffset());
}
Also used : Alignment(com.milaboratory.core.alignment.Alignment) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) PairedReadMergingResult(com.milaboratory.core.merger.PairedReadMergingResult) GeneType(io.repseq.core.GeneType)

Example 12 with GeneType

use of io.repseq.core.GeneType in project mixcr by milaboratory.

the class TargetMerger method merge.

@SuppressWarnings("unchecked")
public AlignedTarget merge(AlignedTarget targetLeft, AlignedTarget targetRight, int offset, OverlapType overlapType, int nMismatches) {
    if (offset < 0)
        return merge(targetRight, targetLeft, -offset, overlapType, nMismatches);
    final NSequenceWithQuality mergedTarget = merger.overlap(targetLeft.getTarget(), targetRight.getTarget(), offset);
    EnumMap<GeneType, VDJCHit[]> result = new EnumMap<>(GeneType.class);
    for (GeneType geneType : GeneType.VJC_REFERENCE) {
        final BatchAlignerWithBaseParameters bp = ((KGeneAlignmentParameters) alignerParameters.getGeneAlignerParameters(geneType)).getParameters();
        final VDJCHit[] leftHits = targetLeft.getAlignments().getHits(geneType);
        final VDJCHit[] rightHits = targetRight.getAlignments().getHits(geneType);
        GeneFeature alignedFeature = leftHits.length == 0 ? rightHits.length == 0 ? null : rightHits[0].getAlignedFeature() : leftHits[0].getAlignedFeature();
        Map<VDJCGeneId, HitMappingRecord> map = extractHitsMapping(targetLeft, targetRight, geneType);
        ArrayList<VDJCHit> resultingHits = new ArrayList<>();
        for (Map.Entry<VDJCGeneId, HitMappingRecord> mE : map.entrySet()) {
            final VDJCGene gene = mE.getValue().gene;
            Alignment<NucleotideSequence> mergedAl = merge(bp.getScoring(), extractBandedWidth(bp), mergedTarget.getSequence(), offset, mE.getValue().alignments[0], mE.getValue().alignments[1]);
            resultingHits.add(new VDJCHit(gene, mergedAl, alignedFeature));
        }
        Collections.sort(resultingHits);
        // final float relativeMinScore = extractRelativeMinScore(bp);
        // int threshold = (int) (resultingHits.size() > 0 ? resultingHits.get(0).getScore() * relativeMinScore : 0);
        // for (int i = resultingHits.size() - 1; i > 0; --i)
        // if (resultingHits.get(i).getScore() < threshold)
        // resultingHits.remove(i);
        result.put(geneType, resultingHits.toArray(new VDJCHit[resultingHits.size()]));
    }
    VDJCAlignments alignments = new VDJCAlignments(result, new NSequenceWithQuality[] { mergedTarget }, new SequenceHistory[] { new SequenceHistory.Merge(overlapType, targetLeft.getHistory(), targetRight.getHistory(), offset, nMismatches) }, VDJCAlignments.mergeOriginalReads(targetLeft.getAlignments(), targetRight.getAlignments()));
    AlignedTarget resultTarget = new AlignedTarget(alignments, 0);
    for (BPoint bPoint : BPoint.values()) {
        int leftPoint = targetLeft.getBPoint(bPoint);
        int rightPoint = targetRight.getBPoint(bPoint);
        if (leftPoint != -1 && rightPoint != -1)
            throw new IllegalArgumentException("Same bPoint defined in both input targets.");
        else if (leftPoint != -1)
            resultTarget = resultTarget.setBPoint(bPoint, leftPoint);
        else if (rightPoint != -1)
            resultTarget = resultTarget.setBPoint(bPoint, offset + rightPoint);
    }
    return resultTarget;
}
Also used : GeneFeature(io.repseq.core.GeneFeature) SequenceHistory(com.milaboratory.mixcr.basictypes.SequenceHistory) VDJCGeneId(io.repseq.core.VDJCGeneId) NSequenceWithQuality(com.milaboratory.core.sequence.NSequenceWithQuality) VDJCAlignments(com.milaboratory.mixcr.basictypes.VDJCAlignments) BatchAlignerWithBaseParameters(com.milaboratory.core.alignment.batch.BatchAlignerWithBaseParameters) KGeneAlignmentParameters(com.milaboratory.mixcr.vdjaligners.KGeneAlignmentParameters) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) VDJCGene(io.repseq.core.VDJCGene) GeneType(io.repseq.core.GeneType) VDJCHit(com.milaboratory.mixcr.basictypes.VDJCHit)

Example 13 with GeneType

use of io.repseq.core.GeneType in project mixcr by milaboratory.

the class VDJCAlignmentsDifferenceReader method compare.

private Diff compare(VDJCAlignments first, VDJCAlignments second) {
    // fixme correct in 2.2
    if (first == null && second == null)
        return null;
    else if (first == null)
        return new Diff(first, second, DiffStatus.AlignmentPresentOnlyInSecond);
    else if (second == null)
        return new Diff(first, second, DiffStatus.AlignmentPresentOnlyInFirst);
    else if (first.getMinReadId() < second.getMinReadId())
        return new Diff(first, second, DiffStatus.AlignmentPresentOnlyInFirst);
    else if (first.getMinReadId() > second.getMinReadId())
        return new Diff(first, second, DiffStatus.AlignmentPresentOnlyInSecond);
    else {
        boolean diffGeneFeature = true;
        if (featureToCompare != null)
            diffGeneFeature = !Objects.equals(first.getFeature(featureToCompare), second.getFeature(featureToCompare));
        EnumMap<GeneType, Boolean> diffHits = new EnumMap<>(GeneType.class);
        boolean same = !diffGeneFeature;
        for (GeneType geneType : GeneType.VDJC_REFERENCE) {
            boolean b = sameHits(first.getHits(geneType), second.getHits(geneType), hitsCompareLevel);
            diffHits.put(geneType, !b);
            if (!b)
                same = false;
        }
        if (same)
            return new Diff(first, second, DiffStatus.AlignmentsAreSame);
        else
            return new Diff(first, second, DiffStatus.AlignmentsAreDifferent, new DiffReason(diffGeneFeature, diffHits));
    }
}
Also used : GeneType(io.repseq.core.GeneType) EnumMap(java.util.EnumMap)

Example 14 with GeneType

use of io.repseq.core.GeneType 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 15 with GeneType

use of io.repseq.core.GeneType in project mixcr by milaboratory.

the class VDJCAlignerAbstract method init.

@Override
protected void init() {
    DAlignerParameters dAlignerParameters = parameters.getDAlignerParameters();
    List<VDJCGene> dGenes = genesToAlign.get(GeneType.Diversity);
    if (dAlignerParameters != null && dGenes.size() != 0)
        singleDAligner = new SingleDAligner(dAlignerParameters, genesToAlign.get(GeneType.Diversity));
    vAligner = createKAligner(GeneType.Variable);
    jAligner = createKAligner(GeneType.Joining);
    cAligner = createKAligner(GeneType.Constant);
    Chains chains = new Chains();
    for (VDJCGene gene : getUsedGenes()) chains = chains.merge(gene.getChains());
    filters = new EnumMap<>(GeneType.class);
    for (GeneType geneType : GeneType.VJC_REFERENCE) {
        HashMap<String, BitArray> f = new HashMap<>();
        for (final String chain : chains) {
            BatchAlignerWithBaseWithFilter<NucleotideSequence, VDJCGene, AlignmentHit<NucleotideSequence, VDJCGene>> aligner = getAligner(geneType);
            if (aligner != null)
                f.put(chain, aligner.createFilter(new Filter<VDJCGene>() {

                    @Override
                    public boolean accept(VDJCGene object) {
                        return object.getChains().contains(chain);
                    }
                }));
        }
        filters.put(geneType, f);
    }
}
Also used : HashMap(java.util.HashMap) Chains(io.repseq.core.Chains) AlignmentHit(com.milaboratory.core.alignment.batch.AlignmentHit) NucleotideSequence(com.milaboratory.core.sequence.NucleotideSequence) VDJCGene(io.repseq.core.VDJCGene) GeneType(io.repseq.core.GeneType) BitArray(com.milaboratory.util.BitArray)

Aggregations

GeneType (io.repseq.core.GeneType)25 NucleotideSequence (com.milaboratory.core.sequence.NucleotideSequence)14 GeneFeature (io.repseq.core.GeneFeature)8 VDJCHit (com.milaboratory.mixcr.basictypes.VDJCHit)7 VDJCGene (io.repseq.core.VDJCGene)5 PairedRead (com.milaboratory.core.io.sequence.PairedRead)4 NSequenceWithQuality (com.milaboratory.core.sequence.NSequenceWithQuality)4 EnumMap (java.util.EnumMap)4 Test (org.junit.Test)4 Alignment (com.milaboratory.core.alignment.Alignment)3 VDJCAlignments (com.milaboratory.mixcr.basictypes.VDJCAlignments)3 Chains (io.repseq.core.Chains)3 ReferencePoint (io.repseq.core.ReferencePoint)3 VDJCGeneId (io.repseq.core.VDJCGeneId)3 ArrayList (java.util.ArrayList)3 MultiAlignmentHelper (com.milaboratory.core.alignment.MultiAlignmentHelper)2 AlignmentHit (com.milaboratory.core.alignment.batch.AlignmentHit)2 VDJCAlignerParameters (com.milaboratory.mixcr.vdjaligners.VDJCAlignerParameters)2 HashMap (java.util.HashMap)2 Map (java.util.Map)2