--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/heteroprob.py	Fri Oct 03 20:54:30 2014 -0400
@@ -0,0 +1,199 @@
+### import libraries ###
+import sys
+import collections, math
+import heapq
+import itertools
+
+
+
+### basic function ###
+def permuterepeat(n,rlist):
+    f = math.factorial
+    nfac=f(n)
+    rfaclist=[f(i) for i in rlist]
+    for rfac in rfaclist:
+        nfac=nfac/rfac
+    return nfac
+
+def nCr(n,r):
+    f = math.factorial
+    return f(n) / f(r) / f(n-r)
+    
+def averagelist(a,b,expectedlevelofminor):
+    product=[]
+    for i in range(len(a)):
+        product.append((1-expectedlevelofminor)*a[i]+expectedlevelofminor*b[i])
+  
+    return product
+        
+def complement_base(read):
+    collect=''
+    for i in read:
+        if i.upper()=='A':
+            collect+='T'
+        elif i.upper()=='T':
+            collect+='A'
+        elif i.upper()=='C':
+            collect+='G'
+        elif i.upper()=='G':
+            collect+='C'
+    return collect
+def makeallpossible(read):
+    collect=[]
+    for i in range(len(read)):
+        tmp= read[i:]+read[:i]
+        collect.append(tmp)
+        collect.append(complement_base(tmp))
+    return collect
+
+def motifsimplify(base):
+    '''str--> str
+    '''
+    motiflength=len(base)
+    temp=list(set(ALLMOTIF[motiflength]).intersection(set(makeallpossible(base))))
+    
+    return temp[0]
+
+def majorallele(seq):
+    binseq=list(set(seq))  
+    binseq.sort(reverse=True)   # highly mutate mode
+    #binseq.sort()              # majority mode
+    storeform=''
+    storevalue=0
+    for i in binseq:
+        if seq.count(i)>storevalue:
+            storeform=i
+            storevalue=seq.count(i)
+            
+    return int(storeform)
+
+### decide global parameter ###
+COORDINATECOLUMN=1
+ALLELECOLUMN=2
+MOTIFCOLUMN=3
+inputname=sys.argv[1]
+errorprofile=sys.argv[2]
+EXPECTEDLEVELOFMINOR=float(sys.argv[3])
+if EXPECTEDLEVELOFMINOR >0.5:
+	try:
+		errorexpectcontribution=int('a')
+	except Exception, eee:
+		print eee
+		stop_err("Expected contribution of minor allele must be at least 0 and not more than 0.5")
+MINIMUMMUTABLE=0 ###1.2*(1.0/(10**8))  #http://www.ncbi.nlm.nih.gov/pubmed/22914163 Kong et al 2012
+
+
+## Fixed global variable
+ALLREPEATTYPE=[1,2,3,4]
+ALLREPEATTYPENAME=['mono','di','tri','tetra']
+monomotif=['A','C']
+dimotif=['AC','AG','AT','CG']
+trimotif=['AAC','AAG','AAT','ACC','ACG','ACT','AGC','AGG','ATC','CCG']
+tetramotif=['AAAC','AAAG','AAAT','AACC','AACG','AACT','AAGC','AAGG','AAGT','AATC','AATG','AATT',\
+'ACAG','ACAT','ACCC','ACCG','ACCT','ACGC','ACGG','ACGT','ACTC','ACTG','AGAT','AGCC','AGCG','AGCT',\
+'AGGC','AGGG','ATCC','ATCG','ATGC','CCCG','CCGG','AGTC']
+ALLMOTIF={1:monomotif,2:dimotif,3:trimotif,4:tetramotif}
+monorange=range(5,60)
+dirange=range(6,60)
+trirange=range(9,60)
+tetrarange=range(12,80)
+ALLRANGE={1:monorange,2:dirange,3:trirange,4:tetrarange}
+
+#########################################
+######## Prob calculation sector ########
+#########################################
+def multinomial_prob(majorallele,STRlength,motif,probdatabase):
+    '''int,int,str,dict-->int
+    ### get prob for each STRlength to be generated from major allele
+    '''
+    #print (majorallele,STRlength,motif)
+    prob=probdatabase[len(motif)][motif][majorallele][STRlength]
+    return prob
+
+################################################
+######## error model database sector ###########
+################################################
+
+## structure generator
+errormodeldatabase={1:{},2:{},3:{},4:{}}
+sumbymajoralleledatabase={1:{},2:{},3:{},4:{}}
+for repeattype in ALLREPEATTYPE:
+    for motif in ALLMOTIF[repeattype]:
+        errormodeldatabase[repeattype][motif]={}
+        sumbymajoralleledatabase[repeattype][motif]={}
+        for motifsize1 in ALLRANGE[repeattype]:
+            errormodeldatabase[repeattype][motif][motifsize1]={}
+            sumbymajoralleledatabase[repeattype][motif][motifsize1]=0
+            for motifsize2 in ALLRANGE[repeattype]:
+                errormodeldatabase[repeattype][motif][motifsize1][motifsize2]=MINIMUMMUTABLE
+#print errormodeldatabase
+## read database
+
+## get read count for each major allele
+fd=open(errorprofile)
+lines=fd.readlines()
+for line in lines:
+    temp=line.strip().split('\t')
+    t_major=int(temp[0])
+    t_count=int(temp[2])
+    motif=temp[3]
+    sumbymajoralleledatabase[len(motif)][motif][t_major]+=t_count
+fd.close()
+##print sumbymajoralleledatabase
+
+## get probability
+fd=open(errorprofile)
+lines=fd.readlines()
+for line in lines:
+    temp=line.strip().split('\t')
+    t_major=int(temp[0])
+    t_read=int(temp[1])
+    t_count=int(temp[2])
+    motif=temp[3]
+    if sumbymajoralleledatabase[len(motif)][motif][t_major]>0:
+        errormodeldatabase[len(motif)][motif][t_major][t_read]=t_count/(sumbymajoralleledatabase[len(motif)][motif][t_major]*1.0)
+        #errormodeldatabase[repeattype][motif][t_major][t_read]=math.log(t_count/(sumbymajorallele[t_major]*1.0))
+        
+    #else:
+    #    errormodeldatabase[repeattype][motif][t_major][t_read]=0
+fd.close()
+#print errormodeldatabase    
+#print math.log(100,10)
+#########################################
+######## input reading sector ###########
+#########################################
+
+
+
+fd = open(inputname)
+##fd=open('sampleinput_C.txt')
+lines=fd.xreadlines()
+for line in lines:
+    i_read=[]
+    i2_read=[]
+    temp=line.strip().split('\t')
+    i_coordinate=temp[COORDINATECOLUMN-1]
+    i_motif=motifsimplify(temp[MOTIFCOLUMN-1])
+    i_read=temp[ALLELECOLUMN-1].split(',')
+    i_read=map(int,i_read)
+    depth=len(i_read)
+    heteromajor1=int(temp[6])
+    heteromajor2=int(temp[7])
+
+### calculate the change to detect combination (using error profile)
+    heterozygous_collector=0  
+    alist=[multinomial_prob(heteromajor1,x,i_motif,errormodeldatabase)for x in i_read]
+    blist=[multinomial_prob(heteromajor2,x,i_motif,errormodeldatabase)for x in i_read]
+      
+    ablist=averagelist(alist,blist,EXPECTEDLEVELOFMINOR)
+       
+    if 0 in ablist:
+        continue
+    heterozygous_collector=reduce(lambda y, z: y*z,ablist )
+
+### prob of combination (using multinomial distribution)
+    frequency_distribution=[len(list(group)) for key, group in itertools.groupby(i_read)]
+    ## print frequency_distribution
+    expandbypermutation=permuterepeat(depth,frequency_distribution)
+
+    print line.strip()+'\t'+str(heterozygous_collector)+'\t'+str(expandbypermutation)+'\t'+str(expandbypermutation*heterozygous_collector)+'\t'+str(depth)