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());
}
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;
}
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));
}
}
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()]));
}
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);
}
}
Aggregations