Can now adjust the genotype likelihoods of a variant returned from the rod. This automatically causes the lodBtr, lodBtnb, and genotype to be recomputed.

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@1041 348d0f76-0448-11de-a6fe-93d51630548a
This commit is contained in:
kiran 2009-06-18 07:26:37 +00:00
parent 9a7cec7d2e
commit 7a921c908c
2 changed files with 39 additions and 1 deletions

View File

@ -5,6 +5,8 @@ import org.broadinstitute.sting.utils.GenomeLoc;
import java.io.IOException;
public class rodVariants extends BasicReferenceOrderedDatum {
private enum Genotype { AA, AC, AG, AT, CC, CG, CT, GG, GT, TT }
private GenomeLoc loc;
private char refBase = 'N';
private int depth;
@ -75,4 +77,40 @@ public class rodVariants extends BasicReferenceOrderedDatum {
public float getLodBtnb() { return lodBtnb; }
public float[] getGenotypeLikelihoods() { return genotypeLikelihoods; }
public void adjustLikelihoods(double[] likelihoods) {
for (int likelihoodIndex = 0; likelihoodIndex < likelihoods.length; likelihoodIndex++) {
genotypeLikelihoods[likelihoodIndex] += likelihoods[likelihoodIndex];
}
String bestGenotype = "NN";
double bestLikelihood = Double.NEGATIVE_INFINITY;
double nextBestLikelihood = Double.NEGATIVE_INFINITY;
double refLikelihood = Double.NEGATIVE_INFINITY;
for (int likelihoodIndex = 0; likelihoodIndex < likelihoods.length; likelihoodIndex++) {
if (genotypeLikelihoods[likelihoodIndex] > bestLikelihood) {
bestLikelihood = genotypeLikelihoods[likelihoodIndex];
bestGenotype = Genotype.values()[likelihoodIndex].toString();
}
}
for (int likelihoodIndex = 0; likelihoodIndex < likelihoods.length; likelihoodIndex++) {
if (genotypeLikelihoods[likelihoodIndex] > nextBestLikelihood && genotypeLikelihoods[likelihoodIndex] < bestLikelihood) {
nextBestLikelihood = genotypeLikelihoods[likelihoodIndex];
}
}
for (int likelihoodIndex = 0; likelihoodIndex < likelihoods.length; likelihoodIndex++) {
if (refBase == Genotype.values()[likelihoodIndex].toString().charAt(0) &&
refBase == Genotype.values()[likelihoodIndex].toString().charAt(1)) {
refLikelihood = genotypeLikelihoods[likelihoodIndex];
}
}
this.bestGenotype = bestGenotype;
this.lodBtr = (float) (bestLikelihood - refLikelihood);
this.lodBtnb = (float) (bestLikelihood - nextBestLikelihood);
}
}

View File

@ -47,7 +47,7 @@ public class VariantFiltrationWalker extends LocusWalker<Integer, Integer> {
if (feature.equalsIgnoreCase("binomial")) { ivf = new IVFBinomialStrand(); }
else { throw new StingException(String.format("Cannot understand feature '%s'\n", feature)); }
variant.multiplyLikelihoods(ivf.compute(ref, context));
variant.adjustLikelihoods(ivf.compute(ref, context));
vwriter.println(variant);
}