New upstream version 5.1-190921+dfsg (17c4ecf0) · Commits · Debian Med / beagle

imp/ImpStates.java

+54 −54

Original line number	Diff line number	Diff line
		@@ -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

+3 −3

Original line number	Diff line number	Diff line
		@@ -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

+31 −12

Original line number	Diff line number	Diff line
		@@ -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 = wtprobsA[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 = wtprobsA[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

+17 −17

Original line number	Diff line number	Diff line
		@@ -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

+3 −3

Original line number	Diff line number	Diff line
		@@ -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);
		}