hagfish_wiggle

#!/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()