[med-svn] [beagle] 01/02: Imported Upstream version 4.1~160727.86a+dfsg
Dylan Aïssi
bob.dybian-guest at moszumanska.debian.org
Fri Sep 30 21:50:40 UTC 2016
This is an automated email from the git hooks/post-receive script.
bob.dybian-guest pushed a commit to branch master
in repository beagle.
commit 6f5d56fdbc37e82721f3abb0916111cc672cf874
Author: Dylan Aïssi <bob.dybian at gmail.com>
Date: Thu Sep 29 23:33:01 2016 +0200
Imported Upstream version 4.1~160727.86a+dfsg
---
blbutil/Utilities.java | 4 +-
dag/MergeableDagLevel.java | 8 +--
main/Main.java | 6 +--
sample/LSHapBaum.java | 126 ++++++++++++++++++---------------------------
vcf/VcfRecBuilder.java | 118 ++++++++++++++++++++----------------------
vcf/VcfWriter.java | 2 +-
6 files changed, 117 insertions(+), 147 deletions(-)
diff --git a/blbutil/Utilities.java b/blbutil/Utilities.java
index b62a575..44df554 100644
--- a/blbutil/Utilities.java
+++ b/blbutil/Utilities.java
@@ -126,7 +126,7 @@ public class Utilities {
*
* @throws NullPointerException if {@code out == null}
*/
- public static void duoPrint(PrintWriter out, String s) {
+ public static void duoPrint(PrintWriter out, CharSequence s) {
System.out.print(s);
out.print(s);
}
@@ -141,7 +141,7 @@ public class Utilities {
*
* @throws NullPointerException if {@code out == null}
*/
- public static void duoPrintln(PrintWriter out, String s) {
+ public static void duoPrintln(PrintWriter out, CharSequence s) {
System.out.println(s);
out.println(s);
}
diff --git a/dag/MergeableDagLevel.java b/dag/MergeableDagLevel.java
index 74c8f4e..0a1135a 100644
--- a/dag/MergeableDagLevel.java
+++ b/dag/MergeableDagLevel.java
@@ -166,16 +166,15 @@ public class MergeableDagLevel {
private void fillArrays(HapsMarker data, float[] weights) {
int parentNode = 0;
+ int nEdges = 0;
for (int hap=0, n=data.nHaps(); hap<n; ++hap) {
int symbol = data.allele(hap);
float count = weights[hap];
int edge = this.outEdges[symbol][parentNode];
if (edge == -1) {
- edge = symbol;
- addEdge(parentNode, symbol, count, edge, hap);
+ addEdge(parentNode, symbol, count, nEdges++, hap);
}
else {
- assert edge == symbol;
assert edge==childNodes[edge];
int child = childNodes[edge];
this.counts[edge] += count;
@@ -183,6 +182,9 @@ public class MergeableDagLevel {
this.child2FirstHap[child] = hap;
}
}
+ if (nEdges < 0.75*nHaps) {
+ reduceEdgeArrayLengths(nEdges);
+ }
}
private void fillArrays(MergeableDagLevel prevLevel, HapsMarker data,
diff --git a/main/Main.java b/main/Main.java
index 79b3dfc..30c0eb3 100644
--- a/main/Main.java
+++ b/main/Main.java
@@ -65,8 +65,8 @@ public class Main {
/**
* The program name and version.
*/
- public static final String program = "beagle.16Jun16.7e4.jar (version 4.1)";
- public static final String command = "java -jar beagle.16Jun16.7e4.jar";
+ public static final String program = "beagle.27Jul16.86a.jar (version 4.1)";
+ public static final String command = "java -jar beagle.27Jul16.86a.jar";
/**
* The copyright string.
@@ -78,7 +78,7 @@ public class Main {
*/
public static final String shortHelp = Main.program
+ Const.nl + Main.copyright
- + Const.nl + "Enter \"java -jar beagle.16Jun16.7e4.jar\" for a "
+ + Const.nl + "Enter \"java -jar beagle.27Jul16.86a.jar\" for a "
+ "summary of command line " + "arguments.";
private final Par par;
diff --git a/sample/LSHapBaum.java b/sample/LSHapBaum.java
index 09dea47..32890b7 100644
--- a/sample/LSHapBaum.java
+++ b/sample/LSHapBaum.java
@@ -42,15 +42,12 @@ public class LSHapBaum {
private final float[] alleleProbs;
private final float[][] fwdVal;
private final float[] bwdVal;
- private final float[] emBwdVal;
private final int[] fwdValueIndex2Marker;
private final RefHapSegs refHapSegs;
private final float[][] fwdHapProbs;
private final float[][] bwdHapProbs;
- private float emBwdValuesSum = 0f;
-
private int windowIndex = -9999;
private int arrayIndex = -9999;
@@ -76,7 +73,6 @@ public class LSHapBaum {
this.fwdValueIndex2Marker = new int[size];
this.fwdVal = new float[size][n];
this.bwdVal = new float[n];
- this.emBwdVal = new float[n];
this.refHapSegs = impData.refHapSegs();
this.fwdHapProbs = new float[impData.nClusters()][];
@@ -101,10 +97,12 @@ public class LSHapBaum {
* {@code hap < 0 || hap >= this.imputationData().targetHapPairs().nHaps()}
*/
public HapAlleleProbs randomHapSample(int hap) {
- int nMarkers = impData.nClusters();
Arrays.fill(alleleProbs, 0f);
- setFwdValues(hap);
- setInitBwdValue(hap);
+ int nMarkers = impData.nClusters();
+ windowIndex = 0;
+ arrayIndex = -1;
+ setForwardValues(0, nMarkers, hap);
+ Arrays.fill(bwdVal, 1.0f/n);
setStateProbs(nMarkers-1, currentIndex());
for (int m=nMarkers-2; m>=0; --m) {
setBwdValue(m, hap);
@@ -123,72 +121,46 @@ public class LSHapBaum {
return impData;
}
- private void setFwdValues(int hap) {
- int nMarkers = impData.nClusters();
- windowIndex = 0;
- arrayIndex = -1;
- for (int m=0; m<nMarkers; ++m) {
- float sum = 0f;
+ private void setForwardValues(int start, int end, int hap) {
+ float lastSum = 1.0f;
+ for (int m=start; m<end; ++m) {
float probRec = impData.pRecomb(m);
+ float probNoRec = 1.0f - probRec;
+ float noErrProb = impData.noErrProb(m);
+ float errProb = impData.errProb(m);
+ float shift = probRec/n;
+ float scale = probNoRec/lastSum;
int prev = currentIndex();
int next = nextIndex();
+ float sum = 0.0f;
fwdValueIndex2Marker[next] = m;
int a = impData.targetAllele(m, hap);
for (int h=0; h<n; ++h) {
- int refAllele = impData.refAllele(m, h);
- float em = (a == refAllele) ? impData.noErrProb(m) : impData.errProb(m);
- float x = m==0 ? 1 : (probRec/n + (1-probRec)*fwdVal[prev][h]);
- fwdVal[next][h] = em*x;
+ float em = (a == impData.refAllele(m, h)) ? noErrProb : errProb;
+ fwdVal[next][h] = m==0 ? em : em*(scale*fwdVal[prev][h] + shift);
sum += fwdVal[next][h];
}
- scale(fwdVal[next], sum);
- }
- }
-
- private static float sum(float[] fa) {
- float sum = 0f;
- for (float f : fa) {
- sum += f;
- }
- return sum;
- }
-
- private static void scale(float[] fa, float divisor) {
- for (int j=0; j<fa.length; ++j) {
- fa[j] /= divisor;
- }
- }
-
- private void setInitBwdValue(int hap) {
- int m = impData.nClusters() - 1;
- float f = 1f/n;
- emBwdValuesSum = 0f;
- int a = impData.targetAllele(m, hap);
- for (int h=0; h<n; ++h) {
- int refAllele = impData.refAllele(m, h);
- float em = (a == refAllele) ? impData.noErrProb(m) : impData.errProb(m);
- bwdVal[h] = f;
- emBwdVal[h] = f*em;
- emBwdValuesSum += emBwdVal[h];
+ lastSum = sum;
}
}
private void setBwdValue(int m, int hap) {
- float bwdValuesSum = 0f;
- float probRec = impData.pRecomb(m + 1);
- float commonTerm = emBwdValuesSum*probRec/n;
+ int mP1 = m + 1;
+ float probRec = impData.pRecomb(mP1);
+ float probNoRec = 1.0f - probRec;
+ float noErrProb = impData.noErrProb(mP1);
+ float errProb = impData.errProb(mP1);
+ float sum = 0f;
+ int al = impData.targetAllele(mP1, hap);
for (int h=0; h<n; ++h) {
- bwdVal[h] = commonTerm + (1-probRec)*emBwdVal[h];
- bwdValuesSum += bwdVal[h];
+ float em = (al==impData.refAllele(mP1, h)? noErrProb : errProb);
+ bwdVal[h] *= em;
+ sum += bwdVal[h];
}
- int a = impData.targetAllele(m, hap);
- emBwdValuesSum = 0f;
+ float scale = probNoRec/sum;
+ float shift = probRec/n;
for (int h=0; h<n; ++h) {
- bwdVal[h] /= bwdValuesSum; // normalize first
- int refAllele = impData.refAllele(m, h);
- float em = (a == refAllele) ? impData.noErrProb(m) : impData.errProb(m);
- emBwdVal[h] = em*bwdVal[h];
- emBwdValuesSum += emBwdVal[h];
+ bwdVal[h] = scale*bwdVal[h] + shift;
}
}
@@ -200,12 +172,27 @@ public class LSHapBaum {
fwdHapProbs[m][refHapSegs.seq(m+1, h)] += stateProbs;
bwdHapProbs[m][refHapSegs.seq(m, h)] += stateProbs;
}
- scale(fwdHapProbs[m], sum(fwdHapProbs[m]));
- scale(bwdHapProbs[m], sum(bwdHapProbs[m]));
+ float sum = sum(fwdHapProbs[m]);
+ scale(fwdHapProbs[m], sum);
+ scale(bwdHapProbs[m], sum);
+ }
+
+ private static float sum(float[] fa) {
+ float sum = 0f;
+ for (float f : fa) {
+ sum += f;
+ }
+ return sum;
+ }
+
+ private static void scale(float[] fa, float divisor) {
+ for (int j=0; j<fa.length; ++j) {
+ fa[j] /= divisor;
+ }
}
private static float threshold(int nSeq) {
- return 0.5f/nSeq;
+ return Math.min(0.005f, 1.0f/nSeq);
}
private void setAlleleProbs(float[] alleleProbs) {
@@ -299,22 +286,7 @@ public class LSHapBaum {
arrayIndex = windowIndex;
int start = fwdValueIndex2Marker[arrayIndex] + 1;
int end = start + ( fwdVal.length - (arrayIndex + 1) );
- for (int m=start; m<end; ++m) {
- float sum = 0f;
- float probRec = impData.pRecomb(m);
- int prev = currentIndex();
- int next = nextIndex();
- fwdValueIndex2Marker[next] = m;
- int a = impData.targetAllele(m, hap);
- for (int h=0; h<n; ++h) {
- int refAllele = impData.refAllele(m, h);
- float em = (a == refAllele) ? impData.noErrProb(m) : impData.errProb(m);
- float x = (probRec/n + (1-probRec)*fwdVal[prev][h]); // since m>0
- fwdVal[next][h] = em*x;
- sum += fwdVal[next][h];
- }
- scale(fwdVal[next], sum);
- }
+ setForwardValues(start, end, hap);
return arrayIndex;
}
else {
diff --git a/vcf/VcfRecBuilder.java b/vcf/VcfRecBuilder.java
index e2356cc..889604b 100644
--- a/vcf/VcfRecBuilder.java
+++ b/vcf/VcfRecBuilder.java
@@ -45,10 +45,10 @@ public final class VcfRecBuilder {
*/
public static final int DEFAULT_INIT_SIZE = 50;
- private static final MathContext mc2 = new MathContext(2);
+ private static final MathContext MC2 = new MathContext(2);
private static final BigDecimal ONE = new BigDecimal(1.0);
- private static final String[] dsProbs = dsProbs();
- private static final String[] r2Probs = r2Probs(dsProbs);
+ private static final String[] DS_VALS = dsProbs();
+ private static final String[] R2_VALS = r2Probs(DS_VALS);
private final StringBuilder sb;
private final R2Estimator r2Est;
@@ -164,63 +164,78 @@ public final class VcfRecBuilder {
/**
* Adds the FORMAT field for a sample to the VCF record for the current
- * marker.
- * @param gtProbs the posterior genotype probabilities
+ * marker. If the specified posterior genotype probabilities do not
+ * sum to 1.0, the specified array will normalized to sum to 1.0.
+ * @param gtypeProbs the posterior genotype probabilities
* @throws IllegalArgumentException if
* {@code gtProbs.length != this.marker().nGenotypes()}
+ * @throws IllegalArgumentException if any element of the specified
+ * array is not a finite non-negative number
* @throws IllegalStateException if {@code this.marker() == null}
* @throws NullPointerException if {@code gtProbs == null}
*/
- public void addSampleData(double[] gtProbs) {
- initializeFieldsAndCheckParameters(gtProbs);
- this.gtProbs[0] = gtProbs[0];
- gt3Probs[0] = gtProbs[0];
- double maxProb = gtProbs[0];
- double sum = gtProbs[0];
- int gt = 1;
- for (int a2=1; a2<nAlleles; ++a2) {
+ public void addSampleData(double[] gtypeProbs) {
+ if (marker==null) {
+ throw new IllegalStateException();
+ }
+ Arrays.fill(gt3Probs, 0.0);
+ Arrays.fill(dose, 0.0);
+ int maxGt = maxIndex(gtypeProbs, nGenotypes);
+ int gt = 0;
+ for (int a2=0; a2<nAlleles; ++a2) {
for (int a1=0; a1<=a2; ++a1) {
- double gtProb = gtProbs[gt];
+ double gtProb = gtypeProbs[gt];
if (gtProb < 0) {
throw new IllegalArgumentException(String.valueOf(gtProb));
}
- this.gtProbs[gt] = gtProb;
- sum += gtProb;
+ gtProbs[gt] = gtProb;
dose[a1] += gtProb;
dose[a2] += gtProb;
- gt3Probs[(a1==0) ? 1 : 2] += gtProb;
- if (gtProb > maxProb) {
+ if (a2==0) {
+ gt3Probs[0] += gtProb;
+ }
+ else {
+ gt3Probs[(a1==0) ? 1 : 2] += gtProb;
+ }
+ if (gt==maxGt) {
allele1 = a1;
allele2 = a2;
- maxProb = gtProb;
}
++gt;
}
}
- if (Math.abs(sum - 1.0) > 1e-6) {
- throw new IllegalArgumentException(String.valueOf(sum));
- }
addToCumAlleleProbs(dose);
r2Est.addSampleData(gt3Probs);
appendFormatData(false);
}
- private void initializeFieldsAndCheckParameters(double[] gtProbs) {
- if (marker==null) {
- throw new IllegalStateException();
+ private int maxIndex(double[] da, int expLength) {
+ if (da.length != expLength) {
+ throw new IllegalArgumentException(String.valueOf(da.length));
+ }
+ int maxIndex = 0;
+ double sum = 0;
+ for (int j=0; j<da.length; ++j) {
+ if (da[j] < 0 || Double.isFinite(da[j])==false) {
+ throw new IllegalArgumentException(String.valueOf(da[j]));
+ }
+ sum += da[j];
+ if (da[j] > da[maxIndex]) {
+ maxIndex = j;
+ }
}
- if (gtProbs.length != nGenotypes) {
- throw new IllegalArgumentException(String.valueOf(gtProbs.length));
+ if (sum != 1.0) {
+ for (int j=0; j<da.length; ++j) {
+ da[j] /= sum;
+ }
}
- Arrays.fill(gt3Probs, 0.0);
- Arrays.fill(dose, 0.0);
- allele1 = 0;
- allele2 = 0;
+ return maxIndex;
}
/**
* Adds the FORMAT field for a sample to the VCF record for the current
- * marker.
+ * marker. If either of the specified posterior allele probabilities
+ * does not sum to 1.0, it will be normalized to sum to 1.0.
* @param alProbs1 the posterior allele probabilities for the individual's
* first allele
* @param alProbs2 the posterior allele probabilities for the individual's
@@ -229,6 +244,8 @@ public final class VcfRecBuilder {
* {@code alProbs1.length != this.marker().nAlleles()}
* @throws IllegalArgumentException if
* {@code alProbs2.length != this.marker().nAlleles()}
+ * @throws IllegalArgumentException if any element of the specified
+ * array is not a finite non-negative number
* @throws IllegalStateException if {@code this.marker() == null}
* @throws NullPointerException if
* {@code alProbs1 == null || alProbs2 == null}
@@ -237,9 +254,9 @@ public final class VcfRecBuilder {
if (marker==null) {
throw new IllegalStateException();
}
- allele1 = maxIndexAndParameterCheck(alProbs1);
- allele2 = maxIndexAndParameterCheck(alProbs2);
Arrays.fill(gt3Probs, 0.0);
+ allele1 = maxIndex(alProbs1, nAlleles);
+ allele2 = maxIndex(alProbs2, nAlleles);
dose[0] = alProbs1[0] + alProbs2[0];
gtProbs[0] = alProbs1[0] * alProbs2[0];
gt3Probs[0] = gtProbs[0];
@@ -261,27 +278,6 @@ public final class VcfRecBuilder {
appendFormatData(isPhased);
}
- private int maxIndexAndParameterCheck(double[] da) {
- if (da.length != nAlleles) {
- throw new IllegalArgumentException(String.valueOf(da.length));
- }
- int maxIndex = 0;
- double sum = 0;
- for (int j=0; j<da.length; ++j) {
- if (da[j] < 0) {
- throw new IllegalArgumentException(String.valueOf(da[j]));
- }
- sum += da[j];
- if (da[j] > da[maxIndex]) {
- maxIndex = j;
- }
- }
- if (Math.abs(sum - 1.0) > 1e-6) {
- throw new IllegalArgumentException(String.valueOf(sum));
- }
- return maxIndex;
- }
-
private void addToCumAlleleProbs(double[] dose) {
for (int j=0; j<dose.length; ++j) {
cumAlleleProbs[j] += dose[j];
@@ -296,13 +292,13 @@ public final class VcfRecBuilder {
if (printDS) {
for (int j=1; j<nAlleles; ++j) {
sb.append( (j==1) ? Const.colon : Const.comma );
- sb.append(dsProbs[(int) Math.rint(100*dose[j])]);
+ sb.append(DS_VALS[(int) Math.rint(100*dose[j])]);
}
}
if (printGP) {
for (int j=0; j<gtProbs.length; ++j) {
sb.append(j==0? Const.colon : Const.comma);
- sb.append(dsProbs[(int) Math.rint(100*gtProbs[j])]);
+ sb.append(DS_VALS[(int) Math.rint(100*gtProbs[j])]);
}
}
}
@@ -337,9 +333,9 @@ public final class VcfRecBuilder {
private void printInfo(PrintWriter out, boolean isImputed) {
if (printDS || printGP) {
out.print("AR2=");
- out.print(r2Probs[(int) Math.rint(100*r2Est.allelicR2())]);
+ out.print(R2_VALS[(int) Math.rint(100*r2Est.allelicR2())]);
out.print(";DR2=");
- out.print(r2Probs[(int) Math.rint(100*r2Est.doseR2())]);
+ out.print(R2_VALS[(int) Math.rint(100*r2Est.doseR2())]);
for (int j=1; j<nAlleles; ++j) {
out.print( (j==1) ? ";AF=" : Const.comma);
out.print(formatProb(cumAlleleProbs[j]/(2*r2Est.nGenotypes())));
@@ -394,10 +390,10 @@ public final class VcfRecBuilder {
private static String formatProb(double p) {
if (p>=0 && p <= 0.5) {
- return new BigDecimal(p).round(mc2).toString();
+ return new BigDecimal(p).round(MC2).toString();
}
else if (p <= 1.0) {
- return new BigDecimal(p-1.0).round(mc2).add(ONE).toString();
+ return new BigDecimal(p-1.0).round(MC2).add(ONE).toString();
}
else {
throw new IllegalArgumentException(String.valueOf(p));
diff --git a/vcf/VcfWriter.java b/vcf/VcfWriter.java
index 31a1639..b410514 100644
--- a/vcf/VcfWriter.java
+++ b/vcf/VcfWriter.java
@@ -29,7 +29,7 @@ public final class VcfWriter {
private static final String dr2Info = "##INFO=<ID=DR2,Number=1,Type=Float,"
+ "Description=\"Dosage R-Squared: estimated squared correlation between "
+ "estimated REF dose [P(RA) + 2*P(RR)] and true REF dose\">";
- private static final String impInfo = "##INFO=<ID=IMP,Number=1,Type=Flag,"
+ private static final String impInfo = "##INFO=<ID=IMP,Number=0,Type=Flag,"
+ "Description=\"Imputed marker\">";
private static final String gtFormat = "##FORMAT=<ID=GT,Number=1,Type=String,"
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/beagle.git
More information about the debian-med-commit
mailing list