# HG changeset patch
# User xuebing
# Date 1333245015 14400
# Node ID 54f508fcdc4015d2ba5cced9c2d3e71eab9a9d85
# Parent  8147a2451531643387a97ca58d521961d7f04afc
Uploaded

diff -r 8147a2451531 -r 54f508fcdc40 seq_to_meme.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/seq_to_meme.py	Sat Mar 31 21:50:15 2012 -0400
@@ -0,0 +1,102 @@
+# -*- coding: iso-8859-1 -*-
+import random,sys,math
+
+#import pylab
+
+def readMotifFile(filename):
+ 
+    try:
+        f=open(filename)
+    except IOError:
+        print "could not open",filename,"Are you sure this file exists?"
+        sys.exit(1)
+    
+    seqs=[]
+    maxL = 0
+    for line in f:
+        if '>' in line or 'N' in line:
+            next
+        else:
+            seq = line.strip().upper()
+            if maxL < len(seq):
+                maxL = len(seq)
+            seqs.append(seq)
+    f.close()
+
+    print len(seqs),'sequences loaded'
+    print 'max seq length:',maxL
+    for i in range(len(seqs)):
+        if len(seqs[i]) < maxL:
+            del seqs[i]
+    print len(seqs),'sequences with length = ',maxL
+    return seqs
+
+
+def createWeightMatrix(seqs,psuedocont):
+
+    motifWidth = len(seqs[0])
+    weightMatrix = []
+    for i in range(motifWidth):
+        weightMatrix.append({'A':psuedocont,'C':psuedocont,'G':psuedocont,'T':psuedocont})
+
+    #Use a for loop to iterate through all the sequences. For each sequence, begin at the start site in starts, and look at motifWidth bases. Count how many times each base appears in each position of the motif
+    #YOUR CODE HERE
+    for seq in seqs:
+        for pos in range(motifWidth):
+            weightMatrix[pos][seq[pos]] = weightMatrix[pos][seq[pos]] + 1.0
+    
+    #Normalize your weight matrix (so that it contains probabilities rather than counts)
+    #Remember the added psuedocounts when you normalize!
+    for pos in range(motifWidth):
+        totalCount = sum(weightMatrix[pos].values())
+        for letter in weightMatrix[pos].keys():
+            weightMatrix[pos][letter] = weightMatrix[pos][letter]/totalCount
+    
+    #Return your weight matrix
+    return weightMatrix
+
+def printMemeFormat(weightMatrix,motifName,filename,nSite,background):
+    f = open(filename,'w')
+    f.write('MEME version 4.4\n\n')
+    
+    f.write('ALPHABET= ACGT\n\n')
+
+    f.write('strands: + -\n\n')
+
+    f.write('Background letter frequencies (from:\n')
+    f.write(background+'\n\n')
+
+    f.write('MOTIF '+motifName+'\n\n') 
+
+    f.write('BL   MOTIF '+motifName+' width='+str(len(weightMatrix))+' seqs='+str(nSite)+'\n')
+    f.write('letter-probability matrix: alength= 4 '+'w= '+str(len(weightMatrix))+' nsites= '+str(nSite)+' E= 0\n')
+    for position in range(len(weightMatrix)):
+        probsThisPosition=weightMatrix[position]
+        f.write('  '+"%.6f" %(probsThisPosition['A'])+'\t  '+"%.6f" %(probsThisPosition['C'])+'\t  '+"%.6f" %(probsThisPosition['G'])+'\t  '+"%.6f" %(probsThisPosition['T'])+'\t\n')
+    f.write('\n\n')
+    f.close()
+    
+#get a two decimal-place string representation of a float f
+def twoDecimal(f):
+    return "%.6f" %(f)
+
+def run():
+
+    #Get file name from command line
+    if len(sys.argv) < 3:
+        print "python seq2meme.py motif_fasta outputfile motifName psuedocont background"
+        sys.exit(1)
+    else:
+        motifFile=sys.argv[1] #
+        outFile=sys.argv[2]
+        motifName=sys.argv[3]
+        psuedocont = float(sys.argv[4])
+        background=' '.join(sys.argv[5].strip().split(','))
+                
+    motifs=readMotifFile(motifFile)
+
+    #Create weight matrix
+    motifWeightMatrix=createWeightMatrix(motifs,psuedocont)
+    printMemeFormat(motifWeightMatrix,motifName,outFile,len(motifs),background)
+run()
+