generate wiggle formatted plots

Author: Mark Fiers

hagfish_wiggle

@@ -0,0 +1,153 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import pickle
+
+import numpy as np
+
+import logging
+import optparse
+
+from hagfish_file_util import *
+
+## Arguments: General options 
+parser = optparse.OptionParser()
+parser.set_defaults(binSize=500)
+parser.add_option('-v', dest='verbose', action="count", 
+                  help='Show debug information')
+parser.add_option('-b', dest='binSize',  type='int',
+                  help='binSize')
+parser.add_option('-s', dest='seqId',  
+                  help='seqId to use in the output..')
+
+options, args = parser.parse_args()
+
+l = logging.getLogger('hagfish')
+handler = logging.StreamHandler()
+logmark = chr(27) + '[0;37;44mHAGFISH' + \
+          chr(27) + '[0m ' 
+
+formatter = logging.Formatter(
+    logmark + '%(levelname)-6s %(message)s')
+
+handler.setFormatter(formatter)
+l.addHandler(handler)
+
+if options.verbose >= 2:
+    l.setLevel(logging.DEBUG)
+elif options.verbose == 1:
+    l.setLevel(logging.INFO)
+else:
+    l.setLevel(logging.WARNING)
+
+
+def generate_histogram(F,
+                       cptype,
+                       category,
+                       seqId,
+                       group,
+                       data, 
+                       bin):
+    
+    if category == 'low':
+        color = '001,086,102'
+    elif category == 'ok':
+        color = '000,125,028'
+    else:
+        color = '166,16,0'
+
+    F.write('track name="Hf_%s_%s_%s" description="Hagfish wiggle plot %s %s for %s" color=%s\n' % 
+            (cptype, category, seqId, cptype, category, seqId, color))
+    F.write("fixedStep chrom=%s start=1 step=%d\n"% (
+            seqId, bin))
+
+	
+    for r in range(0, len(data)+bin-1,bin):
+        score = np.average(data[r:r+bin-1])
+        if not np.isnan(score):
+            F.write("%.2f\n" % score)
+            
+if __name__ == '__main__':
+    
+    #read an arbitrary seqId file
+    for f in os.listdir('seqInfo'):
+        if '.seqinfo' in f:
+            seqInfoFile = os.path.join('seqInfo', f)
+            break
+        else:
+            l.critical("cannot find a seqInfo file")
+            sys.exit(-1)
+
+    l.info("reading %s to get seqinfo" % seqInfoFile)
+    with open(seqInfoFile) as F:
+        seqInfo = pickle.load(F)
+        
+    l.info("discovered %d sequences" % len(seqInfo))
+
+    binSize = options.binSize
+    if len(args) > 0:
+        seqIds = args
+    else:
+        seqIds = seqInfo.keys()
+
+    if not os.path.exists('wiggle'):
+        os.mkdir('wiggle')
+
+    if options.seqId: assert(len(seqIds) == 1)
+        
+    for seqId in seqIds:
+        l.info("Processing %s" % seqId)
+        
+        if options.seqId:
+            outSeqId = options.seqId
+        else:
+            outSeqId = seqId
+
+        l.info("Outputting %s" % outSeqId)
+        
+        FG = open(os.path.join('wiggle', '%s.icp.ok.WIG' % (outSeqId)), 'w')
+        FS = open(os.path.join('wiggle', '%s.icp.low.WIG' % (outSeqId)), 'w')
+        FL = open(os.path.join('wiggle', '%s.icp.high.WIG' % (outSeqId)), 'w')
+        FC = open(os.path.join('wiggle', '%s.icp.score.WIG' % (outSeqId)), 'w')
+
+        FGE = open(os.path.join('wiggle', '%s.ecp.ok.WIG' % (outSeqId)), 'w')
+        FSE = open(os.path.join('wiggle', '%s.ecp.low.WIG' % (outSeqId)), 'w')
+        FLE = open(os.path.join('wiggle', '%s.ecp.high.WIG' % (outSeqId)), 'w')
+        
+        file_base = os.path.join('combined', seqId)
+        r_ok = np_load(file_base, 'r_ok')
+        r_low = np_load(file_base, 'r_low')
+        r_high = np_load(file_base, 'r_high')
+        r_ok_ends = np_load(file_base, 'r_ok_ends')
+        r_low_ends = np_load(file_base, 'r_low_ends')
+        r_high_ends = np_load(file_base, 'r_high_ends')
+
+        l.info("read %d datapoints" % len(r_ok))
+        median = np.median(r_ok)
+        score = 1 - 2 * np.exp(-1 * (r_ok / median)) \
+                + np.exp(-1 * ((r_ok + r_low + r_high)/median))
+
+        seqLen = len(r_ok)
+        l.info("sequence length of %s is %d" % (seqId, seqLen))
+
+       
+        generate_histogram(FG,'ICP', 'ok', outSeqId, 'ok', r_ok, binSize)
+        generate_histogram(FS,'ICP', 'low', outSeqId, 'short', r_low, binSize)
+        generate_histogram(FL,'ICP', 'high', outSeqId, 'long', r_high, binSize)
+
+        generate_histogram(FGE,'ECP', 'ok', outSeqId, 'ok_ends', r_ok_ends, binSize)
+        generate_histogram(FSE,'ECP', 'low', outSeqId, 'short_ends', r_low_ends, binSize)
+        generate_histogram(FLE,'ECP', 'high', outSeqId, 'long_ends', r_high_ends, binSize)
+
+        generate_histogram(FC,'ICP', 'Score', outSeqId, 'score', score, binSize)
+
+        FG.close()
+        FS.close()
+        FL.close()
+        FC.close()
+
+        FGE.close()
+        FSE.close()
+        FLE.close()
+