libprosic 0.4.0

A library for calling of genomic variants using a latent variable model.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119

// use std::f64;
//

// use ordered_float::NotNaN;
// use bio::stats::{LogProb, Prob};
// use bio::stats::combinatorics::combinations;
//
// use model::{Variant, DiscreteAlleleFreqs, AlleleFreq};
// use priors::{Model, InfiniteSitesNeutralVariationModel, PairModel};

use itertools::Itertools;

use model::{AlleleFreq, DiscreteAlleleFreqs};

pub fn allele_freqs(ploidy: u32) -> DiscreteAlleleFreqs {
    (0..ploidy + 1).map(|m| AlleleFreq(m as f64 / ploidy as f64)).collect_vec()
}


/*
pub struct NormalNormalModel {
    inner: InfiniteSitesNeutralVariationModel,
    ploidy: u32,
    allele_freqs: DiscreteAlleleFreqs
}


impl NormalNormalModel {
    pub fn new(ploidy: u32, heterozygosity: Prob) -> Self {
        NormalNormalModel {
            inner: InfiniteSitesNeutralVariationModel::new(2, ploidy, heterozygosity),
            ploidy: ploidy,
            allele_freqs: allele_freqs(ploidy)
        }
    }
}


impl PairModel<DiscreteAlleleFreqs, DiscreteAlleleFreqs> for NormalNormalModel {
    fn prior_prob(&self, af_first: AlleleFreq, af_second: AlleleFreq, _: &Variant) -> LogProb {
        let sample_prior = |af: AlleleFreq| {
            let m = *af * self.ploidy as f64;
            let valid_genotypes = combinations(self.ploidy, m);
            let all_genotypes = combinations()
        };
    }

    fn joint_prob<L, O>(
        &self,
        af_first: &DiscreteAlleleFreqs,
        af_second: &DiscreteAlleleFreqs,
        likelihood_first: &L,
        likelihood_second: &O,
        _: &Variant
    ) -> LogProb where
        L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
        O: Fn(AlleleFreq, AlleleFreq) -> LogProb
    {
        let p_second = LogProb::ln_sum_exp(&af_second.iter().map(|&af_second| {
            likelihood_second(af_second, AlleleFreq(0.0))
        }).collect_vec());

        let prob = LogProb::ln_sum_exp(&af_first.iter().map(|&af_first| {
            let p_first = likelihood_first(af_first, AlleleFreq(0.0));
            let prob = p_first + p_second;

            prob
        }).collect_vec());

        prob
    }

    fn marginal_prob<L, O>(
        &self,
        likelihood_first: &L,
        likelihood_second: &O,
        variant: &Variant
    ) -> LogProb where
        L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
        O: Fn(AlleleFreq, AlleleFreq) -> LogProb
    {
        let p = self.joint_prob(
            self.allele_freqs().0,
            self.allele_freqs().1,
            likelihood_first,
            likelihood_second,
            variant
        );

        p
    }

    fn map<L, O>(
        &self,
        likelihood_first: &L,
        likelihood_second: &O,
        _: &Variant
    ) -> (AlleleFreq, AlleleFreq) where
        L: Fn(AlleleFreq, AlleleFreq) -> LogProb,
        O: Fn(AlleleFreq, AlleleFreq) -> LogProb
    {
        fn calc_map<L: Fn(AlleleFreq, AlleleFreq) -> LogProb>(likelihood: &L, afs: &DiscreteAlleleFreqs) -> AlleleFreq {
            let (_, map) = afs.iter().minmax_by_key(|&af| {
                let p = likelihood(*af, AlleleFreq(0.0));
                NotNaN::new(*p).expect("probability is NaN")
            }).into_option().expect("prior has empty allele frequency spectrum");
            *map
        }

        let map_first = calc_map(likelihood_first, self.allele_freqs().0);
        let map_second = calc_map(likelihood_second, self.allele_freqs().1);

        (map_first, map_second)
    }

    fn allele_freqs(&self) -> (&DiscreteAlleleFreqs, &DiscreteAlleleFreqs) {
        (&self.allele_freqs, &self.allele_freqs)
    }
}*/