Mercurial > repos > arkarachai-fungtammasan > str_fm
comparison sequencingdepthconversion_G.py @ 0:07588b899c13 draft
Uploaded
| author | arkarachai-fungtammasan |
|---|---|
| date | Wed, 01 Apr 2015 17:05:51 -0400 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:07588b899c13 |
|---|---|
| 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 |