annotate heteroprob.py @ 2:d5ed5c2e25c3 draft

Uploaded
author arkarachai-fungtammasan
date Wed, 22 Apr 2015 12:48:40 -0400
parents 07588b899c13
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
1 ### import libraries ###
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
2 import sys
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
3 import collections, math
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
4 import heapq
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
5 import itertools
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
6
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
7
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
8
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
9 ### basic function ###
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
10 def permuterepeat(n,rlist):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
11 f = math.factorial
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
12 nfac=f(n)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
13 rfaclist=[f(i) for i in rlist]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
14 for rfac in rfaclist:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
15 nfac=nfac/rfac
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
16 return nfac
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
17
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
18 def nCr(n,r):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
19 f = math.factorial
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
20 return f(n) / f(r) / f(n-r)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
21
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
22 def averagelist(a,b,expectedlevelofminor):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
23 product=[]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
24 for i in range(len(a)):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
25 product.append((1-expectedlevelofminor)*a[i]+expectedlevelofminor*b[i])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
26
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
27 return product
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
28
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
29 def complement_base(read):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
30 collect=''
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
31 for i in read:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
32 if i.upper()=='A':
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
33 collect+='T'
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
34 elif i.upper()=='T':
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
35 collect+='A'
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
36 elif i.upper()=='C':
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
37 collect+='G'
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
38 elif i.upper()=='G':
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
39 collect+='C'
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
40 return collect
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
41 def makeallpossible(read):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
42 collect=[]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
43 for i in range(len(read)):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
44 tmp= read[i:]+read[:i]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
45 collect.append(tmp)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
46 collect.append(complement_base(tmp))
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
47 return collect
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
48
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
49 def motifsimplify(base):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
50 '''str--> str
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
51 '''
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
52 motiflength=len(base)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
53 temp=list(set(ALLMOTIF[motiflength]).intersection(set(makeallpossible(base))))
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
54
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
55 return temp[0]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
56
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
57 def majorallele(seq):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
58 binseq=list(set(seq))
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
59 binseq.sort(reverse=True) # highly mutate mode
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
60 #binseq.sort() # majority mode
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
61 storeform=''
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
62 storevalue=0
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
63 for i in binseq:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
64 if seq.count(i)>storevalue:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
65 storeform=i
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
66 storevalue=seq.count(i)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
67
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
68 return int(storeform)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
69
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
70 ### decide global parameter ###
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
71 COORDINATECOLUMN=1
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
72 ALLELECOLUMN=2
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
73 MOTIFCOLUMN=3
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
74 inputname=sys.argv[1]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
75 errorprofile=sys.argv[2]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
76 EXPECTEDLEVELOFMINOR=float(sys.argv[3])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
77 if EXPECTEDLEVELOFMINOR >0.5:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
78 try:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
79 errorexpectcontribution=int('a')
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
80 except Exception, eee:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
81 print eee
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
82 stop_err("Expected contribution of minor allele must be at least 0 and not more than 0.5")
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
83 MINIMUMMUTABLE=0 ###1.2*(1.0/(10**8)) #http://www.ncbi.nlm.nih.gov/pubmed/22914163 Kong et al 2012
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
84
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
85
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
86 ## Fixed global variable
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
87 ALLREPEATTYPE=[1,2,3,4]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
88 ALLREPEATTYPENAME=['mono','di','tri','tetra']
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
89 monomotif=['A','C']
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
90 dimotif=['AC','AG','AT','CG']
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
91 trimotif=['AAC','AAG','AAT','ACC','ACG','ACT','AGC','AGG','ATC','CCG']
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
92 tetramotif=['AAAC','AAAG','AAAT','AACC','AACG','AACT','AAGC','AAGG','AAGT','AATC','AATG','AATT',\
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
93 'ACAG','ACAT','ACCC','ACCG','ACCT','ACGC','ACGG','ACGT','ACTC','ACTG','AGAT','AGCC','AGCG','AGCT',\
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
94 'AGGC','AGGG','ATCC','ATCG','ATGC','CCCG','CCGG','AGTC']
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
95 ALLMOTIF={1:monomotif,2:dimotif,3:trimotif,4:tetramotif}
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
96 monorange=range(5,60)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
97 dirange=range(6,60)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
98 trirange=range(9,60)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
99 tetrarange=range(12,80)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
100 ALLRANGE={1:monorange,2:dirange,3:trirange,4:tetrarange}
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
101
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
102 #########################################
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
103 ######## Prob calculation sector ########
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
104 #########################################
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
105 def multinomial_prob(majorallele,STRlength,motif,probdatabase):
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
106 '''int,int,str,dict-->int
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
107 ### get prob for each STRlength to be generated from major allele
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
108 '''
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
109 #print (majorallele,STRlength,motif)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
110 prob=probdatabase[len(motif)][motif][majorallele][STRlength]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
111 return prob
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
112
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
113 ################################################
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
114 ######## error model database sector ###########
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
115 ################################################
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
116
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
117 ## structure generator
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
118 errormodeldatabase={1:{},2:{},3:{},4:{}}
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
119 sumbymajoralleledatabase={1:{},2:{},3:{},4:{}}
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
120 for repeattype in ALLREPEATTYPE:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
121 for motif in ALLMOTIF[repeattype]:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
122 errormodeldatabase[repeattype][motif]={}
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
123 sumbymajoralleledatabase[repeattype][motif]={}
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
124 for motifsize1 in ALLRANGE[repeattype]:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
125 errormodeldatabase[repeattype][motif][motifsize1]={}
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
126 sumbymajoralleledatabase[repeattype][motif][motifsize1]=0
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
127 for motifsize2 in ALLRANGE[repeattype]:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
128 errormodeldatabase[repeattype][motif][motifsize1][motifsize2]=MINIMUMMUTABLE
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
129 #print errormodeldatabase
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
130 ## read database
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
131
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
132 ## get read count for each major allele
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
133 fd=open(errorprofile)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
134 lines=fd.readlines()
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
135 for line in lines:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
136 temp=line.strip().split('\t')
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
137 t_major=int(temp[0])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
138 t_count=int(temp[2])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
139 motif=temp[3]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
140 sumbymajoralleledatabase[len(motif)][motif][t_major]+=t_count
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
141 fd.close()
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
142 ##print sumbymajoralleledatabase
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
143
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
144 ## get probability
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
145 fd=open(errorprofile)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
146 lines=fd.readlines()
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
147 for line in lines:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
148 temp=line.strip().split('\t')
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
149 t_major=int(temp[0])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
150 t_read=int(temp[1])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
151 t_count=int(temp[2])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
152 motif=temp[3]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
153 if sumbymajoralleledatabase[len(motif)][motif][t_major]>0:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
154 errormodeldatabase[len(motif)][motif][t_major][t_read]=t_count/(sumbymajoralleledatabase[len(motif)][motif][t_major]*1.0)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
155 #errormodeldatabase[repeattype][motif][t_major][t_read]=math.log(t_count/(sumbymajorallele[t_major]*1.0))
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
156
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
157 #else:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
158 # errormodeldatabase[repeattype][motif][t_major][t_read]=0
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
159 fd.close()
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
160 #print errormodeldatabase
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
161 #print math.log(100,10)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
162 #########################################
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
163 ######## input reading sector ###########
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
164 #########################################
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
165
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
166
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
167
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
168 fd = open(inputname)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
169 ##fd=open('sampleinput_C.txt')
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
170 lines=fd.xreadlines()
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
171 for line in lines:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
172 i_read=[]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
173 i2_read=[]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
174 temp=line.strip().split('\t')
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
175 i_coordinate=temp[COORDINATECOLUMN-1]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
176 i_motif=motifsimplify(temp[MOTIFCOLUMN-1])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
177 i_read=temp[ALLELECOLUMN-1].split(',')
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
178 i_read=map(int,i_read)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
179 depth=len(i_read)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
180 heteromajor1=int(temp[6])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
181 heteromajor2=int(temp[7])
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
182
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
183 ### calculate the change to detect combination (using error profile)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
184 heterozygous_collector=0
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
185 alist=[multinomial_prob(heteromajor1,x,i_motif,errormodeldatabase)for x in i_read]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
186 blist=[multinomial_prob(heteromajor2,x,i_motif,errormodeldatabase)for x in i_read]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
187
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
188 ablist=averagelist(alist,blist,EXPECTEDLEVELOFMINOR)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
189
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
190 if 0 in ablist:
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
191 continue
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
192 heterozygous_collector=reduce(lambda y, z: y*z,ablist )
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
193
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
194 ### prob of combination (using multinomial distribution)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
195 frequency_distribution=[len(list(group)) for key, group in itertools.groupby(i_read)]
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
196 ## print frequency_distribution
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
197 expandbypermutation=permuterepeat(depth,frequency_distribution)
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
198
07588b899c13 Uploaded
arkarachai-fungtammasan
parents:
diff changeset
199 print line.strip()+'\t'+str(heterozygous_collector)+'\t'+str(expandbypermutation)+'\t'+str(expandbypermutation*heterozygous_collector)+'\t'+str(depth)