Search in sources :

Example 1 with FullSiblings

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);
}
Also used : FullSiblings(cz1.breeding.data.FullSiblings)

Example 2 with FullSiblings

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();
}
Also used : FullSiblings(cz1.breeding.data.FullSiblings) HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) DataEntry(cz1.hmm.data.DataEntry) Iterator(java.util.Iterator) ArrayList(java.util.ArrayList) List(java.util.List) HashMap(java.util.HashMap) Map(java.util.Map) HashSet(java.util.HashSet)

Aggregations

FullSiblings (cz1.breeding.data.FullSiblings)2 DataEntry (cz1.hmm.data.DataEntry)1 ArrayList (java.util.ArrayList)1 HashMap (java.util.HashMap)1 HashSet (java.util.HashSet)1 Iterator (java.util.Iterator)1 List (java.util.List)1 Map (java.util.Map)1