[med-svn] [Git][med-team/beagle][master] 2 commits: New upstream version 5.1-190921+dfsg
Dylan Aïssi
gitlab at salsa.debian.org
Mon Oct 14 20:28:42 BST 2019
Dylan Aïssi pushed to branch master at Debian Med / beagle
Commits:
17c4ecf0 by Dylan Aïssi at 2019-10-14T19:28:13Z
New upstream version 5.1-190921+dfsg
- - - - -
d8d6175f by Dylan Aïssi at 2019-10-14T19:28:13Z
Update upstream source from tag 'upstream/5.1-190921+dfsg'
Update to upstream version '5.1-190921+dfsg'
with Debian dir 97694f31460058599fb389fb058a6f86e9e611b4
- - - - -
6 changed files:
- imp/ImpStates.java
- main/Main.java
- phase/ImputeBaum.java
- phase/PhaseStates.java
- vcf/AllData.java
- vcf/BitSetGTRec.java
Changes:
=====================================
imp/ImpStates.java
=====================================
@@ -40,13 +40,13 @@ public final class ImpStates {
private final int NIL = -103;
private final ImpIbs ibsHaps;
private final ImpData impData;
- private final int maxStates;
private final int nClusters;
+ private final int maxStates;
private final IntIntMap hapToEnd;
private final PriorityQueue<CompHapSegment> q;
- private final IntList[] compHapToHapList;
- private final IntList[] compHapToEndList;
+ private final IntList[] compositeHapToHap;
+ private final IntList[] compositeHapToEnd;
private final int[] compHapToListIndex;
private final int[] compHapToHap;
@@ -64,10 +64,10 @@ public final class ImpStates {
this.maxStates = impData.par().imp_states();
this.hapToEnd = new IntIntMap(maxStates);
this.q = new PriorityQueue<>(maxStates);
- this.compHapToHapList = IntStream.range(0, maxStates)
+ this.compositeHapToHap = IntStream.range(0, maxStates)
.mapToObj(j -> new IntList())
.toArray(IntList[]::new);
- this.compHapToEndList = IntStream.range(0, maxStates)
+ this.compositeHapToEnd = IntStream.range(0, maxStates)
.mapToObj(j -> new IntList())
.toArray(IntList[]::new);
this.compHapToListIndex = new int[maxStates];
@@ -76,11 +76,11 @@ public final class ImpStates {
}
/**
- * Returns the input data for genotype imputation.
- * @return the input data for genotype imputation
+ * Returns the maximum number of HMM states at a marker.
+ * @return the maximum number of HMM states at a marker
*/
- public ImpData impData() {
- return impData;
+ public int maxStates() {
+ return maxStates;
}
/**
@@ -91,7 +91,7 @@ public final class ImpStates {
* at the {@code m}-th marker in {@code alMatch[m][j]}. The number of
* HMM states states at each marker is returned.
* @param targHap the haplotype index
- * @param hapIndices the two-dimensional array in which
+ * @param haps the two-dimensional array in which
* reference haplotype indices for each HMM state will be stored
* @param alMatch the two-dimensional array in which allele match status
* between the target haplotype and HMM state will be stored
@@ -104,7 +104,7 @@ public final class ImpStates {
* HMM states
* @throws NullPointerException if any array is {@code null}
*/
- public int ibsStates(int targHap, int[][] hapIndices, boolean[][] alMatch) {
+ public int ibsStates(int targHap, int[][] haps, boolean[][] alMatch) {
initializeFields();
for (int j=0, n=ibsHaps.codedSteps().nSteps(); j<n; ++j) {
int[] ibs = ibsHaps.ibsHaps(targHap, j);
@@ -115,93 +115,93 @@ public final class ImpStates {
if (q.isEmpty()) {
fillQWithRandomHaps(targHap);
}
- int numStates = copyData(targHap, hapIndices, alMatch);
+ int numStates = copyData(targHap, haps, alMatch);
return numStates;
}
private void initializeFields() {
hapToEnd.clear();
for (int j=0, n=q.size(); j<n; ++j) {
- compHapToHapList[j].clear();
- compHapToEndList[j].clear();
+ compositeHapToHap[j].clear();
+ compositeHapToEnd[j].clear();
}
q.clear();
}
private void updateFields(int hap, int step) {
if (hapToEnd.get(hap, NIL)==NIL) { // hap not currently in q
+ updateHeadOfQ();
if (q.size()==maxStates) {
- CompHapSegment head = updateHeadAndPoll();
- hapToEnd.remove(head.hap());
+ CompHapSegment head = q.poll();
int modEnd = ibsHaps.codedSteps().stepStart((head.step() + step) >>> 1);
- compHapToHapList[head.compHapIndex()].add(hap); // hap of new segment
- compHapToEndList[head.compHapIndex()].add(modEnd); // end of old segment
-
+ hapToEnd.remove(head.hap());
+ compositeHapToHap[head.compHapIndex()].add(hap); // hap of new segment
+ compositeHapToEnd[head.compHapIndex()].add(modEnd); // end of old segment
head.updateHap(hap);
head.updateStep(step);
q.offer(head);
}
else {
int compHapIndex = q.size();
- compHapToHapList[compHapIndex].add(hap); // hap of new segment
+ compositeHapToHap[compHapIndex].add(hap); // hap of new segment
q.offer(new CompHapSegment(hap, step, compHapIndex));
}
}
hapToEnd.put(hap, step);
}
- private CompHapSegment updateHeadAndPoll() {
- CompHapSegment head = q.poll();
- int latestEnd = hapToEnd.get(head.hap(), NIL);
- while (head.step()!=latestEnd) {
- head.updateStep(latestEnd);
- q.offer(head);
- head = q.poll();
- latestEnd = hapToEnd.get(head.hap(), NIL);
- }
- return head;
- }
-
- private void fillQWithRandomHaps(int hap) {
- assert q.isEmpty();
- int nHaps = impData.nHaps();
- int nStates = Math.min(nHaps-1, maxStates);
- int sample = hap>>1;
- Random rand = new Random(hap);
- for (int i=0; i<nStates; ++i) {
- int h = rand.nextInt(nHaps);
- while ((h>>1)==sample) {
- h = rand.nextInt(nHaps);
+ private void updateHeadOfQ() {
+ CompHapSegment head = q.peek();
+ if (head!=null) {
+ int latestEnd = hapToEnd.get(head.hap(), NIL);
+ while (head.step()!=latestEnd) {
+ head = q.poll();
+ head.updateStep(latestEnd);
+ q.offer(head);
+ head = q.peek();
+ latestEnd = hapToEnd.get(head.hap(), NIL);
}
- q.add(new CompHapSegment(h, nClusters, i));
}
}
private int copyData(int targHap, int[][] hapIndices, boolean[][] alMatch) {
- int nCompHaps = q.size();
+ int nCompositeHaps = q.size();
+ initializeCopy(nCompositeHaps);
int shiftedTargHap = impData.nRefHaps() + targHap;
- initializeCopy(nCompHaps);
for (int m=0; m<nClusters; ++m) {
- int targAllele = impData.allele(m, shiftedTargHap);
- for (int j=0; j<nCompHaps; ++j) {
+ int targAl = impData.allele(m, shiftedTargHap);
+ for (int j=0; j<nCompositeHaps; ++j) {
if (m==compHapToEnd[j]) {
++compHapToListIndex[j];
- compHapToHap[j] = compHapToHapList[j].get(compHapToListIndex[j]);
- compHapToEnd[j] = compHapToEndList[j].get(compHapToListIndex[j]);
+ compHapToHap[j] = compositeHapToHap[j].get(compHapToListIndex[j]);
+ compHapToEnd[j] = compositeHapToEnd[j].get(compHapToListIndex[j]);
+ assert compHapToHap[j] < impData.nRefHaps();
}
hapIndices[m][j] = compHapToHap[j];
- alMatch[m][j] = impData.allele(m, compHapToHap[j])==targAllele;
+ alMatch[m][j] = impData.allele(m, compHapToHap[j])==targAl;
}
}
- return nCompHaps;
+ return nCompositeHaps;
}
private void initializeCopy(int nSlots) {
for (int j=0; j<nSlots; ++j) {
- compHapToEndList[j].add(nClusters); // add missing end of last segment
+ compositeHapToEnd[j].add(nClusters); // add missing end of last segment
compHapToListIndex[j] = 0;
- compHapToHap[j] = compHapToHapList[j].get(0);
- compHapToEnd[j] = compHapToEndList[j].get(0);
+ compHapToHap[j] = compositeHapToHap[j].get(0);
+ compHapToEnd[j] = compositeHapToEnd[j].get(0);
+ }
+ }
+
+ private void fillQWithRandomHaps(int hap) {
+ assert q.isEmpty();
+ int nRefHaps = impData.nRefHaps();
+ int nStates = Math.min(nRefHaps, maxStates);
+ Random rand = new Random(hap);
+ for (int i=0; i<nStates; ++i) {
+ int h = rand.nextInt(nRefHaps);
+ compositeHapToHap[i].add(h); // hap of new segment
+ q.add(new CompHapSegment(h, nClusters, i));
}
}
}
=====================================
main/Main.java
=====================================
@@ -65,8 +65,8 @@ public class Main {
* The program name and version.
*/
public static final String VERSION = "(version 5.1)";
- public static final String PROGRAM = "beagle.24Aug19.3e8.jar";
- public static final String COMMAND = "java -jar beagle.24Aug19.3e8.jar";
+ public static final String PROGRAM = "beagle.21Sep19.ec3.jar";
+ public static final String COMMAND = "java -jar beagle.21Sep19.ec3.jar";
/**
* The copyright string.
@@ -78,7 +78,7 @@ public class Main {
*/
public static final String SHORT_HELP = Main.PROGRAM + " " + VERSION
+ Const.nl + Main.COPYRIGHT
- + Const.nl + "Enter \"java -jar beagle.24Aug19.3e8.jar\" to "
+ + Const.nl + "Enter \"java -jar beagle.21Sep19.ec3.jar\" to "
+ "list command line argument";
private final Par par;
=====================================
phase/ImputeBaum.java
=====================================
@@ -23,6 +23,7 @@ import ints.IntArray;
import ints.IntList;
import java.util.Random;
import vcf.GT;
+import vcf.RefGT;
/**
* <p>Class {@code ImputeBaum} applies the forward and backward algorithms
@@ -48,7 +49,10 @@ public class ImputeBaum {
private final FloatList savedProbs = new FloatList(8);
private final GT hiFreqPhasedGT;
- private final GT allUnphasedGT;
+ private final GT unphTargGT;
+ private final RefGT refGT;
+ private final int nTargHaps;
+ private final int nHaps;
private final int nHiFreqMarkers;
private final HapImputer imputableHaps;
private final IntArray hiFreqIndices;
@@ -74,11 +78,14 @@ public class ImputeBaum {
this.probs = new float[2][nHiFreqMarkers][fwdBwd.maxStates()];
this.hiFreqPhasedGT = phaseData.estPhase().hapsGT();
- this.allUnphasedGT = fpd.targGT();
+ this.unphTargGT = fpd.targGT();
+ this.refGT = fpd.refGT();
+ this.nTargHaps = fpd.targGT().nHaps();
+ this.nHaps = fpd.nHaps();
this.imputableHaps = hapImputer;
this.hiFreqIndices = fpd.hiFreqIndices();
this.rand = new Random(phaseData.seed());
- this.outPhase = new int[2][allUnphasedGT.nMarkers()];
+ this.outPhase = new int[2][unphTargGT.nMarkers()];
}
/**
@@ -109,15 +116,15 @@ public class ImputeBaum {
outPhase[1][end] = hiFreqPhasedGT.allele2(j, sample);
start = end + 1;
}
- imputeInterval(sample, start, allUnphasedGT.nMarkers());
+ imputeInterval(sample, start, unphTargGT.nMarkers());
imputableHaps.setHap(targHap[0], outPhase[0]);
imputableHaps.setHap(targHap[1], outPhase[1]);
}
private void imputeInterval(int sample, int start, int end) {
for (int m=start; m<end; ++m) {
- int a1 = allUnphasedGT.allele1(m, sample);
- int a2 = allUnphasedGT.allele2(m, sample);
+ int a1 = unphTargGT.allele1(m, sample);
+ int a2 = unphTargGT.allele2(m, sample);
if (a1>=0 && a2>=0) {
boolean noFlip = true;
if (a1!=a2) {
@@ -138,7 +145,7 @@ public class ImputeBaum {
}
private float[] unscaledAlProbs(int m, int hapBit, int a1, int a2) {
- float[] alProbs = new float[allUnphasedGT.marker(m).nAlleles()];
+ float[] alProbs = new float[unphTargGT.marker(m).nAlleles()];
boolean rare1 = fpd.isLowFreq(m, a1);
boolean rare2 = fpd.isLowFreq(m, a2);
int mkrA = fpd.prevHiFreqMarker(m);
@@ -148,8 +155,8 @@ public class ImputeBaum {
float[] probsB = probs[hapBit][mkrB];
for (int j=0, n=nStates[hapBit]; j<n; ++j) {
int hap = statesA[j];
- int b1 = allUnphasedGT.allele(m, hap);
- int b2 = allUnphasedGT.allele(m, (hap ^ 0b1));
+ int b1 = allele(m, hap);
+ int b2 = allele(m, (hap ^ 0b1));
if (b1>=0 && b2>=0) {
float wt = fpd.prevWt(m);
float prob = wt*probsA[j] + (1.0f - wt)*probsB[j];
@@ -176,7 +183,7 @@ public class ImputeBaum {
private int imputeAllele(int m, int hapBit) {
savedStates.clear();
savedProbs.clear();
- float[] alProbs = new float[allUnphasedGT.marker(m).nAlleles()];
+ float[] alProbs = new float[unphTargGT.marker(m).nAlleles()];
float unknownAlProb = 0.0f;
int mkrA = fpd.prevHiFreqMarker(m);
int mkrB = Math.min(mkrA + 1, nHiFreqMarkers - 1);
@@ -187,8 +194,8 @@ public class ImputeBaum {
float wt = fpd.prevWt(m);
float prob = wt*probsA[j] + (1.0f - wt)*probsB[j];
int hap = statesA[j];
- int b1 = allUnphasedGT.allele(m, hap);
- int b2 = allUnphasedGT.allele(m, hap ^ 0b1);
+ int b1 = allele(m, hap);
+ int b2 = allele(m, hap ^ 0b1);
if (b1>=0 && b2>=0) {
if (b1==b2) {
alProbs[b1] += prob;
@@ -208,6 +215,18 @@ public class ImputeBaum {
return imputedAllele;
}
+ private int allele(int marker, int hap) {
+ if (hap>=nHaps) {
+ throw new IndexOutOfBoundsException(String.valueOf(hap));
+ }
+ if (hap < nTargHaps) {
+ return unphTargGT.allele(marker, hap);
+ }
+ else {
+ return refGT.allele(marker, hap - nTargHaps);
+ }
+ }
+
private int maxIndex(float[] fa) {
int maxIndex = 0;
for (int j=1; j<fa.length; ++j) {
=====================================
phase/PhaseStates.java
=====================================
@@ -47,8 +47,8 @@ public final class PhaseStates {
private final IntIntMap hapToEnd;
private final PriorityQueue<CompHapSegment> q;
- private final IntList[] compositeHapHap;
- private final IntList[] compositeHapEnd;
+ private final IntList[] compositeHapToHap;
+ private final IntList[] compositeHapToEnd;
private final int[] compHapToListIndex;
private final int[] compHapToHap;
@@ -75,10 +75,10 @@ public final class PhaseStates {
this.minSteps = (int) Math.ceil(defaultMinSteps*scaleFactor);
this.hapToEnd = new IntIntMap(maxStates);
this.q = new PriorityQueue<>(maxStates);
- this.compositeHapHap = IntStream.range(0, maxStates)
+ this.compositeHapToHap = IntStream.range(0, maxStates)
.mapToObj(j -> new IntList())
.toArray(IntList[]::new);
- this.compositeHapEnd = IntStream.range(0, maxStates)
+ this.compositeHapToEnd = IntStream.range(0, maxStates)
.mapToObj(j -> new IntList())
.toArray(IntList[]::new);
this.compHapToListIndex = new int[maxStates];
@@ -200,8 +200,8 @@ public final class PhaseStates {
private void initializeFields() {
hapToEnd.clear();
for (int j=0, n=q.size(); j<n; ++j) {
- compositeHapHap[j].clear();
- compositeHapEnd[j].clear();
+ compositeHapToHap[j].clear();
+ compositeHapToEnd[j].clear();
}
q.clear();
}
@@ -224,15 +224,15 @@ public final class PhaseStates {
CompHapSegment head = q.poll();
int modEnd = phaseData.codedSteps().stepStart((head.step() + step) >>> 1);
hapToEnd.remove(head.hap());
- compositeHapHap[head.compHapIndex()].add(hap); // hap of new segment
- compositeHapEnd[head.compHapIndex()].add(modEnd); // end of old segment
+ compositeHapToHap[head.compHapIndex()].add(hap); // hap of new segment
+ compositeHapToEnd[head.compHapIndex()].add(modEnd); // end of old segment
head.updateHap(hap);
head.updateStep(step);
q.offer(head);
}
else {
int compHapIndex = q.size();
- compositeHapHap[compHapIndex].add(hap); // hap of new segment
+ compositeHapToHap[compHapIndex].add(hap); // hap of new segment
q.offer(new CompHapSegment(hap, step, compHapIndex));
}
}
@@ -264,8 +264,8 @@ public final class PhaseStates {
for (int j=0; j<nCompHaps; ++j) {
if (m==compHapToEnd[j]) {
++compHapToListIndex[j];
- compHapToHap[j] = compositeHapHap[j].get(compHapToListIndex[j]);
- compHapToEnd[j] = compositeHapEnd[j].get(compHapToListIndex[j]);
+ compHapToHap[j] = compositeHapToHap[j].get(compHapToListIndex[j]);
+ compHapToEnd[j] = compositeHapToEnd[j].get(compHapToListIndex[j]);
}
int refAllele = phaseData.allele(m, compHapToHap[j]);
if (isMissing) {
@@ -290,8 +290,8 @@ public final class PhaseStates {
for (int j=0; j<nCompositeHaps; ++j) {
if (m==compHapToEnd[j]) {
++compHapToListIndex[j];
- compHapToHap[j] = compositeHapHap[j].get(compHapToListIndex[j]);
- compHapToEnd[j] = compositeHapEnd[j].get(compHapToListIndex[j]);
+ compHapToHap[j] = compositeHapToHap[j].get(compHapToListIndex[j]);
+ compHapToEnd[j] = compositeHapToEnd[j].get(compHapToListIndex[j]);
}
int refHap = compHapToHap[j];
haps[m][j] = refHap;
@@ -304,10 +304,10 @@ public final class PhaseStates {
private void initializeCopy(int nSlots) {
for (int j=0; j<nSlots; ++j) {
- compositeHapEnd[j].add(nMarkers); // add missing end of last segment
+ compositeHapToEnd[j].add(nMarkers); // add missing end of last segment
compHapToListIndex[j] = 0;
- compHapToHap[j] = compositeHapHap[j].get(0);
- compHapToEnd[j] = compositeHapEnd[j].get(0);
+ compHapToHap[j] = compositeHapToHap[j].get(0);
+ compHapToEnd[j] = compositeHapToEnd[j].get(0);
}
}
@@ -326,7 +326,7 @@ public final class PhaseStates {
while ((h>>1)==sample) {
h = rand.nextInt(nHaps);
}
- compositeHapHap[q.size()].add(h);
+ compositeHapToHap[q.size()].add(h);
q.add(new CompHapSegment(h, nMarkers, i));
}
}
=====================================
vcf/AllData.java
=====================================
@@ -157,7 +157,7 @@ public class AllData implements Data {
.mapToObj(i -> new IntList(16))
.toArray(IntList[]::new);
int nTargSamples = targGT.nSamples();
- int nRefSamples = (refGT!=null) ? refGT.nSamples() : 0;
+ int nRefSamples = (restrictRefGT!=null) ? restrictRefGT.nSamples() : 0;
for (int s=0; s<nTargSamples; ++s) {
int a1 = targGT.allele1(marker, s);
int a2 = targGT.allele2(marker, s);
@@ -169,8 +169,8 @@ public class AllData implements Data {
}
}
for (int s=0; s<nRefSamples; ++s) {
- int a1 = refGT.allele1(marker, s);
- int a2 = refGT.allele2(marker, s);
+ int a1 = restrictRefGT.allele1(marker, s);
+ int a2 = restrictRefGT.allele2(marker, s);
if (a1>=0 && carriers[a1].size()<=maxCarriers) {
carriers[a1].add(nTargSamples + s);
}
=====================================
vcf/BitSetGTRec.java
=====================================
@@ -89,7 +89,7 @@ public final class BitSetGTRec implements GTRec {
private static boolean isRef(int nSamples, long[] isPhased,
long[] isMissing) {
- if (Long.bitCount(isPhased[0])<Long.SIZE) {
+ if (Long.bitCount(isPhased[0])<Math.min(nSamples, Long.SIZE)) {
return false;
}
else {
View it on GitLab: https://salsa.debian.org/med-team/beagle/compare/3509f23f162ed8627c1b8d893411ff88bd5ee82e...d8d6175f10bd2bfb2e409bec23680bc3318cb64e
--
View it on GitLab: https://salsa.debian.org/med-team/beagle/compare/3509f23f162ed8627c1b8d893411ff88bd5ee82e...d8d6175f10bd2bfb2e409bec23680bc3318cb64e
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20191014/59b14569/attachment-0001.html>
More information about the debian-med-commit
mailing list