Inference of LUX energy matrix using EMA¶

def X_i_predictions(binned_data, energy_matrix):
    czx = zeros((2,32))
    for k in range(shape(binned_data)[0]):
        energy = 0
        for i in range(8):
            if binned_data[k,0][i] == 'A':
                energy += energy_matrix[0][0,i]
            elif binned_data[k,0][i] == 'T':
                energy += energy_matrix[0][1,i]
            elif binned_data[k,0][i] == 'G':
                energy += energy_matrix[0][2,i]
            elif binned_data[k,0][i] == 'C':
                energy += energy_matrix[0][3,i]
        if energy < energy_matrix[1] :
            czx[0,binned_data[k,2]] += 1
        else:
            czx[1,binned_data[k,2]] += 1
    return czx


def likelihood(czx):
    likelihood_numerator = 0
    for j in range(shape(czx)[0]):
        for k in  range(shape(czx)[1]):
            likelihood_numerator += math.lgamma(czx[j,k]+1)
            
    likelihood_denominator = 0
    for j in range(shape(czx)[0]):
        likelihood_denominator += math.lgamma(31+sum(czx[j,:])+1)
    
    return likelihood_numerator - likelihood_denominator


def matrix_perturbation(matrix):
    choose = numpy.random.choice(3)
    perturbed = copy(matrix[0])
    
    if choose == 0:
        x = numpy.random.choice(4)
        y = numpy.random.choice(8)
        perturbed[x,y] = perturbed[x,y]+numpy.random.normal(0.0,.1)
        return perturbed, matrix[1]
        
    elif choose == 1:
        mu2=matrix[1]+numpy.random.normal(0.0,.1)
        
        return perturbed, mu2
    
    elif choose == 2:
        temp = numpy.random.normal(0.0,.1)
        mu2=matrix[1]+temp
        temp2 = numpy.random.choice(8)
        for k in range(4):
            perturbed[k,temp2] = perturbed[k,temp2] + temp        
        
        return perturbed, mu2

def matrix_normalisation(energy_matrix_preturbed):
    temp = copy(energy_matrix_preturbed[0])
    x= 0
    for i in range(8):
        x+=min(temp[:,i]) 
    
    for i in range(8):
        for j in range(4):
            energy_matrix_preturbed[0][j,i] = (energy_matrix_preturbed[0][j,i] - min(temp[:,i]))/(energy_matrix_preturbed[1]-x)
    
    energy_matrix_preturbed[1] = 1.0
    return copy(energy_matrix_preturbed)


def seed_matrix(initial_sequences):
    """ Takes a set of squences and returns a matrix with starting energies for matrix inference using
    kinney et al method
    """
    matrix = zeros((4,8))
    for eightmere in initial_sequences:
        for base in range(len(eightmere[0])):
            if 'A' in eightmere[0][base]:
                matrix[0,base]=1+matrix[0,base]
            elif 'T' in eightmere[0][base]:
                matrix[1,base]=1+matrix[1,base]
            elif 'G' in eightmere[0][base]:
                matrix[2,base]=1+matrix[2,base]
            elif 'C' in eightmere[0][base]:
                matrix[3,base]=1+matrix[3,base]

    for i in range(matrix.shape[1]):
        matrix[:,i]= matrix[:,i]/sum(matrix[:,i])
        matrix[:,i]= matrix[:,i]/max(matrix[:,i])
        matrix[:,i] = abs(matrix[:,i]-max(matrix[:,i]))
    
    return matrix

import pandas

forward = dict([])
reverse = dict([])
forward[cca1_pbm['8-mer'][0]] = 0
for i in range(len(cca1_pbm['8-mer'])-1):
    if cca1_pbm['8-mer'][i+1][::-1] in forward:
        reverse[cca1_pbm['8-mer'][i+1]] = i+1
    else:
        forward[cca1_pbm['8-mer'][i+1]] = i+1
    
cca1_pbm=cca1_pbm.loc[forward.values()]

bin_counter = 0
counter = 0
for i in range(shape(cca1_pbm_sorted_array)[0]):
    cca1_pbm_sorted_array[i,2] = bin_counter
    if counter == 516:
        bin_counter = bin_counter + 1
        counter = 0
    counter = counter + 1

initial_matrix= cca1_pbm_sorted_array[-10:]

for i in range(steps):    
    
    energy_matrix_preturbed2 = matrix_perturbation(energy_matrix_preturbed)
    

    if  energy_matrix_preturbed2[1] >= 0.0  and energy_matrix_preturbed2[1] <=  8. and energy_matrix_preturbed2[0].min() >= 0. and energy_matrix_preturbed2[0].max() <=1. :  
        M2_predictions = X_i_predictions(cca1_pbm_sorted_array,energy_matrix_preturbed2)
        L2 = likelihood(M2_predictions)
        #print L2-L1
        
        if numpy.random.binomial(1,min(1.,exp(L2-L1))):
            #print energy_matrix_preturbed2
            if i % 5 == 0:
                print 'Accepted '+str(len(ensemble))
                print 'Ratio '+ str(acceptance/(rejection+acceptance))
            L1 = L2
            
            energy_matrix_preturbed = energy_matrix_preturbed2
            
            chain.append(L2)
            acceptance +=1.
            ensemble.append(copy(energy_matrix_preturbed2))
            #print "MCMC steps " + str(len(chain))

    else:
        #print 'Rejected '+str(i)
        rejection +=1.
    

save('lux_ensemble.npy',ensemble)