annotate heteroprob.py @ 1:99ec84eb0bab draft default tip

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