Mercurial > repos > iuc > ngsutils_bam_filter
comparison ngsutils/support/llh.py @ 0:4e4e4093d65d draft
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/ngsutils commit 09194687c74a424732f8b0c017cbb942aad89068
author | iuc |
---|---|
date | Wed, 11 Nov 2015 13:04:07 -0500 |
parents | |
children | 7a68005de299 |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:4e4e4093d65d |
---|---|
1 ''' | |
2 Methods for calculating log-likelihoods for nucleotide frequencies | |
3 ''' | |
4 import math | |
5 import collections | |
6 from ngsutils.support import memoize | |
7 | |
8 _default_background = {'A': 0.3, 'T': 0.3, 'C': 0.2, 'G': 0.2} | |
9 | |
10 NucleotideLogLikelihood = collections.namedtuple('NucleotideLogLikelihood', 'A C G T pseudo') | |
11 | |
12 @memoize | |
13 def pseudo_count(N, bg): | |
14 ''' | |
15 >>> pseudo_count(100, _default_background['A']) | |
16 3 | |
17 >>> pseudo_count(100, _default_background['C']) | |
18 2 | |
19 ''' | |
20 | |
21 return bg * math.sqrt(N) | |
22 | |
23 | |
24 def calc_llh(A, C, G, T, bg=_default_background, pseudo='auto'): | |
25 if pseudo == 'auto': | |
26 N = A + C + G + T | |
27 Ap = float(A) + pseudo_count(N, bg['A']) | |
28 Cp = float(C) + pseudo_count(N, bg['C']) | |
29 Gp = float(G) + pseudo_count(N, bg['G']) | |
30 Tp = float(T) + pseudo_count(N, bg['T']) | |
31 elif pseudo: | |
32 Ap = float(A) + pseudo | |
33 Cp = float(C) + pseudo | |
34 Gp = float(G) + pseudo | |
35 Tp = float(T) + pseudo | |
36 else: | |
37 Ap = float(A) | |
38 Cp = float(C) | |
39 Gp = float(G) | |
40 Tp = float(T) | |
41 | |
42 Np = Ap + Cp + Gp + Tp | |
43 | |
44 freqA = Ap / Np | |
45 freqC = Cp / Np | |
46 freqG = Gp / Np | |
47 freqT = Tp / Np | |
48 | |
49 return NucleotideLogLikelihood(math.log(freqA / bg['A']), math.log(freqC / bg['C']), math.log(freqG / bg['G']), math.log(freqT / bg['T']), pseudo) | |
50 | |
51 | |
52 | |
53 if __name__ == '__main__': | |
54 import doctest | |
55 doctest.testmod() |