Binding energy statistics in the yeast genome (intergenic regions). Here we analyse signatures of selection from genomic data using transcription factor binding sites as an example. Directory data/ contains a file alignment.txt.NaN.removed where intergenic sequences for four yeast species (S.cer, S.par, S.mik, S.bay) are given. Note we have removed insertions and deletions from the alignment so that analyses which in this assignment focus on a single species aspects are a bit simpler.

Python code that is given:

from setup import * %matplotlib inline

#import data to a dataframe datafile="data/alignment.txt.NaN.removed"; seq_igs=pd.read_csv(datafile,sep="\s+"); tfBindingFile="data/Emat.abf1.kinney"; # imports an energy matrix modelling transcription factor to DNA binding for factor Abf1 # the dimensions are 4 x 20, correspoding to the four nucleotides A,C,G,T and the binding sites motif length 20 # smaller energy values indicate better binding Emat=np.array(pd.read_csv(tfBindingFile,header=None,sep="\s+"))

def getEnergy(seq,Emat): s1=list(seq.replace('A','0').replace('C','1').replace('G','2').replace('T','3').replace('N','4')) Lmat=Emat.shape[1] Lseq=len(s1) Ev=[] for i in range(0,Lseq-Lmat+1): E=0.0; k=0; flag=0; eps=0.0; for j in range(i,i+Lmat): nuc=int(s1[j]) if nuc < 4: eps=Emat[nuc,k]; else: #remove sequences with missing data flag=1; E+=eps k+=1; if flag==0: Ev.append(E) return Ev

def randomise(seq): s1=list(seq.replace('A','0').replace('C','1').replace('G','2').replace('T','3').replace('N','4')) random.shuffle(s1) return "".join(s1)

Question:

Visualise the binding energy matrix Emat with elements ()k(a) where k denotes column (binding site position) and ,,,aA,C,G,T nucleotides. What is the best possible binding sequence according to this model?

Explain what the functions getEnergy() and randomise() do, and add some comments to the code to help a reader understand how they work.