Mercurial > repos > arkarachai-fungtammasan > microsatellite_ngs

comparison sequencingdepthconversion_G.py @ 0:20ab85af9505

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author arkarachai-fungtammasan
date Fri, 03 Oct 2014 20:54:30 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:20ab85af9505
1 def stop_err(msg):
2 sys.stderr.write(msg)
3 sys.exit()
4
5 def info2require(X,L,F,r):
6 '''infodepth,readlength,flanksize,repeatlength
7 '''
8 return int(math.ceil((X*L*1.0)/(L-(1*((2*F)+r-1)))))
9
10 def poissondef(meancov,specificcov):
11 nominator=1.0*(meancov**specificcov)*(math.e**(-1*meancov))
12 denominator=math.factorial(specificcov)
13 return nominator/denominator
14
15 def require2recommend(needprob,mindepth):
16 i=mindepth
17 reverseneedprob=1-needprob
18 sumprob=1
19 while sumprob>reverseneedprob: #mean cov
20 sumprob=0
21 for j in range(0,mindepth): #specific cov
22 sumprob+=poissondef(i,j)
23 i+=1
24
25 return i-1
26
27 import sys,math
28
29 repeatlength=int(sys.argv[1])
30 flanksize=int(sys.argv[2])#20
31 readlength=int(sys.argv[3])#100
32 infodepth=int(sys.argv[4])#5
33 probdetection=float(sys.argv[5])#0.90
34
35 if probdetection >1:
36 try:
37 probvalue=int('probvalue')
38 except Exception, eee:
39 print eee
40 stop_err("Proportion of genome to have certain locus specific must be between 0 and 1")
41
42 print 'repeat_length'+'\t'+'read_length'+'\t'+'informative_read_depth''\t'+'=locus_specific_sequencing_depth'+'\t'+'=genome_wide_sequencing_depth'
43 t_requiredepth=info2require(infodepth,readlength,flanksize,repeatlength)
44 t_recomendseq=require2recommend(probdetection,t_requiredepth)
45 preplotlist=[repeatlength,readlength,infodepth,t_requiredepth,t_recomendseq]
46 plotlist=map(str,preplotlist)
47 print '\t'.join(plotlist)
48
49 #print info2require(infodepth,readlength,flanksize,repeatlength)
50 #print poissondef(10,3)
51 #print require2recommend(0.90,80)
52 #informative_read_depth
53 #required_seq_depth
54 #recommend_seq_depth