annotate heteroprob.py @ 1:f265e26ab550

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