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