use of cz1.breeding.data.FullSiblings in project polyGembler by c-zhou.
the class HiddenMarkovModelBWT method makeWeightingCoeff.
private void makeWeightingCoeff() {
// TODO Auto-generated method stub
this.weights = new double[M][N];
FullSiblings fullSib = new FullSiblings();
double[][] pvals = fullSib.calcDosaConfig(this.de, this.parent_i);
for (int i = 0; i < M; i++) {
if (StatUtils.max(pvals[i]) == 0) {
Arrays.fill(weights[i], 1.0);
continue;
}
final int p_i = Algebra.maxIndex(pvals[i]);
double[] weights_i = weights[i];
// parental samples
if (!miss[i][this.parent_i[0]] && !miss[i][this.parent_i[1]]) {
int[] parentalDosa = fullSib.getParentalDosaByIndex(p_i);
double[] ll0 = this.dp[i][this.parent_i[0]].likelihood, ll1 = this.dp[i][this.parent_i[1]].likelihood;
if (ll0[parentalDosa[0]] + ll1[parentalDosa[1]] < ll0[parentalDosa[1]] + ll1[parentalDosa[0]]) {
int tmp_i = parentalDosa[0];
parentalDosa[0] = parentalDosa[1];
parentalDosa[1] = tmp_i;
}
weights_i[parent_i[0]] = 1.0 + founder_hap_coeff * ll0[parentalDosa[0]] * (N - 2);
weights_i[parent_i[1]] = 1.0 + founder_hap_coeff * ll1[parentalDosa[1]] * (N - 2);
} else {
weights_i[parent_i[0]] = miss[i][this.parent_i[0]] ? 0.1 : 1.0;
weights_i[parent_i[1]] = miss[i][this.parent_i[1]] ? 0.1 : 1.0;
}
boolean[] f1Dosa = fullSib.getF1DosaByIndex(p_i);
// f1 samples
for (int j = 0; j < N; j++) {
if (this.pedigree[j] != -1)
continue;
int dosa = this.getDosage(i, j);
if (miss[i][j])
weights_i[j] = 0.1;
weights_i[j] = f1Dosa[dosa] ? 1.0 : 0.01;
}
Algebra.normalize(weights_i);
Algebra.multiply(weights_i, N);
}
Utils.print(this.weights);
}
use of cz1.breeding.data.FullSiblings in project polyGembler by c-zhou.
the class VCFResampling method run.
@Override
public void run() {
// TODO Auto-generated method stub
FullSiblings fullSib = new FullSiblings(this.ploidy);
final List<String> contigs = DataCollection.getContigList(data_in);
final List<String> samples = DataCollection.getSampleList(data_in);
int[] founder_index = new int[2];
Arrays.fill(founder_index, -1);
for (int i = 0; i != samples.size(); i++) {
if (samples.get(i).equals(founder_haps[0]))
founder_index[0] = i;
else if (samples.get(i).equals(founder_haps[1]))
founder_index[1] = i;
}
for (String contig : contigs) {
DataEntry data = DataCollection.readDataEntry(data_in, new String[] { contig })[0];
double[] position = data.getPosition();
int no = position.length;
if (no <= min_snp) {
int[] out = new int[no];
for (int i = 0; i != no; i++) out[i] = i;
snp_out.put(contig, out);
continue;
}
double[] segs = fullSib.calcSegs(data, founder_index);
int l = 0;
final Map<Integer, List<Integer>> bins = new HashMap<Integer, List<Integer>>();
for (int i = 0; i != no; i++) {
while ((l + 1) * window_size < position[i]) l++;
if (!bins.containsKey(l))
bins.put(l, new ArrayList<Integer>());
bins.get(l).add(i);
}
List<Integer> selected = new ArrayList<Integer>();
final Map<Integer, Iterator<Integer>> bins_it = new HashMap<Integer, Iterator<Integer>>();
for (Map.Entry<Integer, List<Integer>> entry : bins.entrySet()) {
Collections.sort(entry.getValue(), new Comparator<Integer>() {
@Override
public int compare(Integer i0, Integer i1) {
// TODO Auto-generated method stub
return Double.compare(segs[i0], segs[i1]);
}
});
Iterator<Integer> it = entry.getValue().iterator();
selected.add(it.next());
if (it.hasNext())
bins_it.put(entry.getKey(), it);
}
while (selected.size() < min_snp) {
List<Integer> candidates = new ArrayList<Integer>();
Set<Integer> key_rm = new HashSet<Integer>();
for (Map.Entry<Integer, Iterator<Integer>> entry : bins_it.entrySet()) {
Iterator<Integer> it = entry.getValue();
candidates.add(it.next());
if (!it.hasNext())
key_rm.add(entry.getKey());
}
if (candidates.size() == 0)
break;
for (Integer i : key_rm) bins_it.remove(i);
Collections.sort(candidates, new Comparator<Integer>() {
@Override
public int compare(Integer i0, Integer i1) {
// TODO Auto-generated method stub
return Double.compare(segs[i0], segs[i1]);
}
});
for (Integer i : candidates) {
if (selected.size() < min_snp)
selected.add(i);
else
break;
}
}
Collections.sort(selected);
snp_out.put(contig, ArrayUtils.toPrimitive(selected.toArray(new Integer[selected.size()])));
}
int N = 0;
for (Map.Entry<String, int[]> entry : snp_out.entrySet()) N += entry.getValue().length;
System.err.println(N + " SNPs reserved.");
write();
}