Skip to content

Commit

Permalink
Score of C and D alignments optionally included in total score. For #18.
Browse files Browse the repository at this point in the history
  • Loading branch information
@dbolotin
dbolotin committed Sep 7, 2015
1 parent 8837e8d commit 2973ffa
Show file tree
Hide file tree
Showing 3 changed files with 117 additions and 63 deletions.
129 changes: 81 additions & 48 deletions src/main/java/com/milaboratory/mixcr/vdjaligners/VDJCAlignerPVFirst.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,16 +38,10 @@
import com.milaboratory.core.sequence.NucleotideSequence;
import com.milaboratory.mixcr.basictypes.VDJCAlignments;
import com.milaboratory.mixcr.basictypes.VDJCHit;
import com.milaboratory.mixcr.reference.Allele;
import com.milaboratory.mixcr.reference.GeneFeature;
import com.milaboratory.mixcr.reference.GeneType;
import com.milaboratory.mixcr.reference.ReferencePoint;
import com.milaboratory.mixcr.reference.*;
import gnu.trove.map.hash.TIntObjectHashMap;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Comparator;
import java.util.List;
import java.util.*;

public final class VDJCAlignerPVFirst extends VDJCAlignerAbstract<PairedRead> {
public VDJCAlignerPVFirst(VDJCAlignerParameters parameters) {
Expand Down Expand Up @@ -92,23 +86,19 @@ public VDJCAlignmentResult<PairedRead> process(PairedRead input) {
return new VDJCAlignmentResult<>(input);
}

// Performing alignment of C and D genes if corresponding parameters are set to include their scores to
// the total score value
bestHelper.performCDAlignment();

// Calculates if this score is bigger then the threshold
if (bestHelper.score() < parameters.getMinSumScore()) {
onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
return new VDJCAlignmentResult<>(input);
}

//if (bestHelper.hasHits())

// Finally filtering hits inside this helper to meet minSumScore and maxHits limits
bestHelper.filterHits(parameters.getMinSumScore(), parameters.getMaxHits());

//else {
// onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
// // If hits for V or J are missing
// return new VDJCAlignmentResult<>(input);
//}

// If hits for V or J are missing after filtration
if (!bestHelper.isGood()) {
onFailedAlignment(input, VDJCAlignmentFailCause.LowTotalScore);
Expand Down Expand Up @@ -148,6 +138,7 @@ final class PAlignmentHelper {
final KAlignmentResult<?>[] vResults;
KAlignmentResult<?>[] jResults;
PairedHit[] vHits, jHits;
VDJCHit[] dHits = null, cHits = null;
PairedHit bestVHits;

PAlignmentHelper(PairedTarget target, KAlignmentResult... vResults) {
Expand Down Expand Up @@ -242,8 +233,22 @@ PairedHit[] extractDoubleHits(KAlignmentResult<?>... results) {
* Returns sum score for this targets.
*/
float score() {
return (vHits.length > 0 ? vHits[0].sumScore : 0.0f) +
(jHits != null && jHits.length > 0 ? jHits[0].sumScore : 0.0f);
// Adding V score
float score = vHits.length > 0 ? vHits[0].sumScore : 0.0f;

// Adding J score
if (jHits != null && jHits.length > 0)
score += jHits[0].sumScore;

// Adding C score
if (parameters.doIncludeCScore() && cHits != null && cHits.length > 0)
score += cHits[0].getScore();

// Adding D score
if (parameters.doIncludeDScore() && dHits != null && dHits.length > 0)
score += dHits[0].getScore();

return score;
}

void addHits(TIntObjectHashMap<PairedHit> hits, KAlignmentResult<?> result, int index) {
Expand All @@ -270,23 +275,42 @@ void addHits(TIntObjectHashMap<PairedHit> hits, KAlignmentResult<?> result, int
VDJCAlignments createResult(long readId, VDJCAlignerPVFirst aligner) {
VDJCHit[] vHits = convert(this.vHits, GeneType.Variable, aligner);
VDJCHit[] jHits = convert(this.jHits, GeneType.Joining, aligner);
VDJCHit[] dHits = null, cHits = null;

return new VDJCAlignments(readId, vHits, dHits, jHits, cHits, target.targets);
}

/**
* Preforms J alignment after V alignments are built.
*/
void performJAlignment() {
jHits = extractDoubleHits(jResults = new KAlignmentResult[]{
performJAlignment(0),
performJAlignment(1)
});

calculateScoreAndSort(jHits);
}

/**
* Perform final alignment of D and C genes on fully marked-up reads (with by V and J alignments).
*/
void performCDAlignment() {
PairedHit bestVHit = vHits[0];
PairedHit bestJHit = jHits[0];

//Alignment of D gene
if (singleDAligner != null) {
PreVDJCHit[][] preDHits = new PreVDJCHit[2][];
Arrays.fill(preDHits, zeroArray);

for (int i = 0; i < 2; ++i) {
//for target0
Alignment<NucleotideSequence> vAlignment = vHits[0].getAlignment(i);
Alignment<NucleotideSequence> jAlignment = jHits[0].getAlignment(i);
Alignment<NucleotideSequence> vAlignment = bestVHit.get(i).getAlignment();
Alignment<NucleotideSequence> jAlignment = bestJHit.get(i).getAlignment();
if (vAlignment == null || jAlignment == null)
continue;
int from = vAlignment.getSequence2Range().getTo(),
to = jAlignment.getSequence2Range().getFrom();
if (from > to)
if (from >= to)
continue;
List<PreVDJCHit> temp = singleDAligner.align0(target.targets[i].getSequence(),
getPossibleDLoci(vHits, jHits), from, to);
Expand All @@ -302,31 +326,17 @@ VDJCAlignments createResult(long readId, VDJCAlignerPVFirst aligner) {
KAlignmentHit[][] results = new KAlignmentHit[2][];
Arrays.fill(results, zeroKArray);
for (int i = 0; i < 2; ++i) {
//for target0
Alignment<NucleotideSequence> jAlignment = jHits[0].getAlignment(i);
Alignment<NucleotideSequence> jAlignment = bestJHit.get(i).getAlignment();
if (jAlignment == null)
continue;
int from = jAlignment.getSequence2Range().getTo();
List<KAlignmentHit> temp = cAligner.align(target.targets[i].getSequence(), from, target.targets[i].size()).getHits();
List<KAlignmentHit> temp = cAligner.align(target.targets[i].getSequence(), from,
target.targets[i].size()).getHits();
results[i] = temp.toArray(new KAlignmentHit[temp.size()]);
}
cHits = combine(getCAllelesToAlign(),
parameters.getFeatureToAlign(GeneType.Constant), results);
}

return new VDJCAlignments(readId, vHits, dHits, jHits, cHits, target.targets);
}

/**
* Preforms J alignment after V alignments are built.
*/
void performJAlignment() {
jHits = extractDoubleHits(jResults = new KAlignmentResult[]{
performJAlignment(0),
performJAlignment(1)
});

calculateScoreAndSort(jHits);
}

/**
Expand Down Expand Up @@ -358,12 +368,17 @@ KAlignmentResult performJAlignment(int index) {
* Filters hit to finally meet maxHit and minScore limits.
*/
public void filterHits(float minTotalScore, int maxHits) {
final float minVScore = Math.max(parameters.getRelativeMinVScore() * vHits[0].sumScore,
minTotalScore - jHits[0].sumScore);
this.vHits = extractHits(minVScore, vHits, maxHits);
// Calculate this value once to use twice in the code below
float totalMScore = minTotalScore - score();

float minScore = Math.max(
parameters.getRelativeMinVScore() * vHits[0].sumScore,
totalMScore + vHits[0].sumScore
);
this.vHits = extractHits(minScore, vHits, maxHits);

if (vHits.length > 0)
this.jHits = extractHits(minTotalScore - vHits[0].sumScore, jHits, maxHits);
if (vHits.length > 0 && jHits.length > 0)
this.jHits = extractHits(totalMScore + jHits[0].sumScore, jHits, maxHits);
}

/**
Expand All @@ -382,6 +397,15 @@ private PairedHit[] extractHits(float minScore, PairedHit[] result, int maxHits)
}
}

public Set<Locus> getPossibleDLoci(PairedHit[] vHits, PairedHit[] jHits) {
EnumSet<Locus> loci = EnumSet.noneOf(Locus.class);
for (PairedHit vHit : vHits)
loci.add(getAllele(GeneType.Variable, vHit.getId()).getLocus());
for (PairedHit jHit : jHits)
loci.add(getAllele(GeneType.Variable, jHit.getId()).getLocus());
return loci;
}

/**
* Converts array of "internal" PairedHits to a double array of KAlignmentHits to pass this value to a VDJAlignment
* constructor (VDJAlignmentImmediate).
Expand All @@ -407,13 +431,14 @@ static void calculateScoreAndSort(PairedHit[] hits) {
* read.
*/
static final class PairedHit {
KAlignmentHit<?> hit0, hit1;
KAlignmentHit<NucleotideSequence> hit0, hit1;
float sumScore = -1, vEndScore = -1;

PairedHit() {
}

PairedHit(KAlignmentHit hit0, KAlignmentHit hit1) {
PairedHit(KAlignmentHit<NucleotideSequence> hit0, KAlignmentHit<NucleotideSequence> hit1) {
assert hit0 == null || hit1 == null || hit0.getId() == hit1.getId();
this.hit0 = hit0;
this.hit1 = hit1;
}
Expand Down Expand Up @@ -450,12 +475,13 @@ void set(int i, KAlignmentHit hit) {
this.hit0 = hit;
else
this.hit1 = hit;
assert hit0 == null || hit1 == null || hit0.getId() == hit1.getId();
}

/**
* To use this hit as an array of two single hits.
*/
KAlignmentHit get(int i) {
KAlignmentHit<NucleotideSequence> get(int i) {
assert i == 0 || i == 1;

if (i == 0)
Expand All @@ -464,6 +490,13 @@ KAlignmentHit get(int i) {
return hit1;
}

/**
* Returns id of reference sequence
*/
int getId() {
return hit0 == null ? hit1.getId() : hit0.getId();
}

/**
* Converts this object to a VDJCHit
*/
Expand Down
49 changes: 34 additions & 15 deletions src/main/java/com/milaboratory/mixcr/vdjaligners/VDJCAlignerParameters.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,32 +50,22 @@ public final class VDJCAlignerParameters implements HasFeatureToAlign, java.io.S
@JsonIgnore
protected final EnumMap<GeneType, GeneAlignmentParameters> alignmentParameters;
protected VJAlignmentOrder vjAlignmentOrder;
protected boolean includeDScore, includeCScore;
protected float minSumScore;
protected int maxHits;
protected float relativeMinVFR3CDR3Score;
protected float relativeMinVScore;
protected PairedEndReadsLayout readsLayout;
protected MergerParameters mergerParameters;

protected VDJCAlignerParameters(EnumMap<GeneType, GeneAlignmentParameters> alignmentParameters, float minSumScore,
int maxHits, float relativeMinVScore, float relativeMinVFR3CDR3Score,
PairedEndReadsLayout readsLayout,
MergerParameters mergerParameters) {
this.alignmentParameters = alignmentParameters;
this.minSumScore = minSumScore;
this.maxHits = maxHits;
this.relativeMinVScore = relativeMinVScore;
this.relativeMinVFR3CDR3Score = relativeMinVFR3CDR3Score;
this.readsLayout = readsLayout;
this.mergerParameters = mergerParameters;
}

@JsonCreator
public VDJCAlignerParameters(@JsonProperty("vParameters") KGeneAlignmentParameters vParameters,
@JsonProperty("dParameters") DAlignerParameters dParameters,
@JsonProperty("jParameters") KGeneAlignmentParameters jParameters,
@JsonProperty("cParameters") KGeneAlignmentParameters cParameters,
@JsonProperty("vjAlignmentOrder") VJAlignmentOrder vjAlignmentOrder,
@JsonProperty("includeDScore") boolean includeDScore,
@JsonProperty("includeCScore") boolean includeCScore,
@JsonProperty("minSumScore") float minSumScore,
@JsonProperty("maxHits") int maxHits,
@JsonProperty("relativeMinVScore") float relativeMinVScore,
Expand All @@ -88,6 +78,8 @@ public VDJCAlignerParameters(@JsonProperty("vParameters") KGeneAlignmentParamete
setGeneAlignerParameters(GeneType.Joining, jParameters);
setGeneAlignerParameters(GeneType.Constant, cParameters);
this.vjAlignmentOrder = vjAlignmentOrder;
this.includeDScore = includeDScore;
this.includeCScore = includeCScore;
this.minSumScore = minSumScore;
this.maxHits = maxHits;
this.relativeMinVScore = relativeMinVScore;
Expand Down Expand Up @@ -177,6 +169,22 @@ public void setVjAlignmentOrder(VJAlignmentOrder vjAlignmentOrder) {
this.vjAlignmentOrder = vjAlignmentOrder;
}

public boolean doIncludeDScore() {
return includeDScore;
}

public void setIncludeDScore(boolean includeDScore) {
this.includeDScore = includeDScore;
}

public boolean doIncludeCScore() {
return includeCScore;
}

public void setIncludeCScore(boolean includeCScore) {
this.includeCScore = includeCScore;
}

public VDJCAlignerParameters setMinSumScore(float minSumScore) {
this.minSumScore = minSumScore;
return this;
Expand Down Expand Up @@ -244,6 +252,9 @@ public MergerParameters getMergerParameters() {
public String toString() {
return "VDJCAlignerParameters{" +
"alignmentParameters=" + alignmentParameters +
", vjAlignmentOrder=" + vjAlignmentOrder +
", includeDScore=" + includeDScore +
", includeCScore=" + includeCScore +
", minSumScore=" + minSumScore +
", maxHits=" + maxHits +
", relativeMinVFR3CDR3Score=" + relativeMinVFR3CDR3Score +
Expand All @@ -261,6 +272,9 @@ public boolean equals(Object o) {
VDJCAlignerParameters that = (VDJCAlignerParameters) o;

if (maxHits != that.maxHits) return false;
if (that.vjAlignmentOrder == vjAlignmentOrder) return false;
if (that.includeDScore == includeDScore) return false;
if (that.includeCScore == includeCScore) return false;
if (Float.compare(that.minSumScore, minSumScore) != 0) return false;
if (Float.compare(that.relativeMinVFR3CDR3Score, relativeMinVFR3CDR3Score) != 0) return false;
if (Float.compare(that.relativeMinVScore, relativeMinVScore) != 0) return false;
Expand All @@ -277,6 +291,10 @@ public int hashCode() {
int result = alignmentParameters.hashCode();
result = 31 * result + (minSumScore != +0.0f ? Float.floatToIntBits(minSumScore) : 0);
result = 31 * result + maxHits;
result = 31 * result + vjAlignmentOrder.hashCode();
result = 31 * result + (includeDScore ? 17 : 0);
result = 31 * result + (includeCScore ? 17 : 0);
result = 31 * result + maxHits;
result = 31 * result + (relativeMinVFR3CDR3Score != +0.0f ? Float.floatToIntBits(relativeMinVFR3CDR3Score) : 0);
result = 31 * result + (relativeMinVScore != +0.0f ? Float.floatToIntBits(relativeMinVScore) : 0);
result = 31 * result + (readsLayout != null ? readsLayout.hashCode() : 0);
Expand All @@ -286,7 +304,8 @@ public int hashCode() {

@Override
public VDJCAlignerParameters clone() {
return new VDJCAlignerParameters(getVAlignerParameters(), getDAlignerParameters(), getJAlignerParameters(), getCAlignerParameters(),
vjAlignmentOrder, minSumScore, maxHits, relativeMinVFR3CDR3Score, relativeMinVScore, readsLayout, mergerParameters);
return new VDJCAlignerParameters(getVAlignerParameters(), getDAlignerParameters(), getJAlignerParameters(),
getCAlignerParameters(), vjAlignmentOrder, includeDScore, includeCScore, minSumScore, maxHits,
relativeMinVFR3CDR3Score, relativeMinVScore, readsLayout, mergerParameters);
}
}
2 changes: 2 additions & 0 deletions src/main/resources/parameters/vdjcaligner_parameters.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,8 @@
}
},
"vjAlignmentOrder": "VThenJ",
"includeDScore": false,
"includeCScore": false,
"mergerParameters": {
"minimalOverlap": 17,
"minimalIdentity": 0.9
Expand Down

0 comments on commit 2973ffa

Please sign in to comment.