Skip to content
Commits on Source (2)
Show whitespace changes
Inline Side-by-side
imp/ImpStates.java
View file @ d8d6175f
......@@ -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();
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.poll();
head = q.peek();
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);
}
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
View file @ d8d6175f
......@@ -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
View file @ d8d6175f
......@@ -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
View file @ d8d6175f
......@@ -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
View file @ d8d6175f
......@@ -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
View file @ d8d6175f
......@@ -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 {
......