[med-svn] [Git][med-team/beagle][upstream] New upstream version 5.1-190921+dfsg

Dylan Aïssi gitlab at salsa.debian.org
Mon Oct 14 20:28:54 BST 2019



Dylan Aïssi pushed to branch upstream at Debian Med / beagle


Commits:
17c4ecf0 by Dylan Aïssi at 2019-10-14T19:28:13Z
New upstream version 5.1-190921+dfsg
- - - - -


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/commit/17c4ecf0446518fbdf7063b525bdd875d10535b1

-- 
View it on GitLab: https://salsa.debian.org/med-team/beagle/commit/17c4ecf0446518fbdf7063b525bdd875d10535b1
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/6b1604a6/attachment-0001.html>


More information about the debian-med-commit mailing list