|
4
|
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 |
