Inference of CCA1 energy matrix using EMA¶

import pandas
import numpy as np
cca1_pbm = pandas.read_csv('CCA1_RVE1_8mers/CCA1_8mers.txt', sep='\t')

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()]

def X_i_predictions(binned_data, energy_matrix):
    czx = zeros((2,384))
    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,:]))
    
    print likelihood_numerator
    print likelihood_numerator
    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

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 == 43:
        bin_counter = bin_counter + 1
        counter = 0
    counter = counter + 1

#figure(figsize=(10,6))
hist(cca1_pbm_sorted_array[:,1], bins=10)
title("E-score", fontsize=20)
xticks(fontsize=20)
yticks(fontsize=20)
#savefig("CCA1-escore.pdf", format='pdf', dpi=600)

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('cca1_ensemble.npy',ensemble)

ensemble = np.load('ensemble.bp.npy')

normalised_ensemble = []
for i in ensemble:
    normalised_ensemble.append((matrix_normalisation(i)[0],matrix_normalisation(i)[1]))

with open ("TOC1_upstream.txt", "r") as myfile:
    data=myfile.read().replace('\n', '')

with open ("NC_003076.8[25252009..25252439].fa", "r") as myfile:
    data2=myfile.read().replace('\n', '')

from Bio.Seq import *

TOC1_uppstream = Seq(data, alphabet=IUPAC.IUPACUnambiguousDNA())

TOC1 = str(TOC1_uppstream_RC)

def free_energy(seq,energy_matrix):
    """The function gives choose is the sequence is bound by the TF given a energy matrix model.
    """
    energy = 0
    for i in range(8):
        if seq[i] == 'A':
            energy += energy_matrix[0][0,i]
        elif seq[i] == 'T':
            energy += energy_matrix[0][1,i]
        elif seq[i] == 'G':
            energy += energy_matrix[0][2,i]
        elif seq[i] == 'C':
            energy += energy_matrix[0][3,i]
    if energy < energy_matrix[1] :
            return energy
    else:
        return energy

wtCCR2_for = 'GAGGTCAAACCTAGAAAATATCTAAACCTTGAAACCTAG'
wtCCR2_rev = 'CTAGGTTTCAAGGTTTAGATATTTTCTAGGTTTGACCTC'
wtPRR9_for=  'CGATCACAACCACGAAAATATCTTCTCAGAGAAAGAAGA'
wtPRR9_rev=  'TCTTCTTTCTCTGAGAAGATATTTTCGTGGTTGTGATCG'
CBS_PRR9_for='CGATCACAACCACGAAAAAATCTTCTCAGAGAAAGAAGA'
CBS_PRR9_rev='TCTTCTTTCTCTGAGAAGATTTTTTCGTGGTTGTGATCG'

mutCCR2_for = 'GAGGTCAAACCTAGAAAATCGAGAAACCTTGAAACCTAG'
mutCCR2_rev = 'CTAGGTTTCAAGGTTTCTCGATTTTCTAGGTTTGACCTC'
mutPRR9_for = 'CGATCACAACCACGAAAATCGAGTCTCAGAGAAAGAAGA'
mutPRR9_rev = 'TCTTCTTTCTCTGAGACTCGATTTTCGTGGTTGTGATCG'