hagfish_compplot2

#!/usr/bin/env python

import os
import sys
import math

import numpy as np
import matplotlib as mpl
mpl.use('agg')
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import matplotlib as mpl


import pylab as pl

import subprocess

from hagfishUtils import *

## Arguments: General options 
parser = getHagfishOptparser()
parser.add_option('-S', dest='score', action='store_true',
                  help='plot the hagfish score' )

parser.add_option('-L', dest='library2', action='append',
                  help='Library to compare to')

addPlotParameters(parser)
options, args = parser.parse_args()

l = getLogger('main', options.verbose)

l.info("starting %s" % sys.argv[0])
    
seqId = args[0]

#load the coverage plots
data = hagfishData(options, args)
data2 = hagfishData(options, args, libraries=options.library2)

#prepare the plot
ym1 = quant(data.ok, options.yfrac)
ym2  = quant(data2.ok, options.yfrac)

ymax = max(ym1, ym2)

if options.ymax:
    ymax = float(options.ymax)

l.debug('y1 %s %s' % (ym1, data.median))
l.debug('y2 %s %s' % (ym2, data2.median))
l.debug('ym %s' % (ymax))

title1 = "_".join(options.library)
title2 = "_".join(options.library2)

plot = hagfishPlot(
    options, data, tag='compplot2', 
    data2=data2, ymax = ymax,
    title='comparative plot %s vs %s' % (title1, title2))

class plotter(hagfishPlotBand):
    
    def plotBand(self):
        self.l.debug("Start plotting")
        yc = self.yCorrection

        median1 = self.data.median
        median2 = self.data2.median

        self.ax.fill_between(
            self.locx, 
            self.zero  + yc,
            (self.d1.low) + yc, 
            color = COLDARKBLUE, 
            linewidth=0)

        self.ax.fill_between(
            self.locx, 
            self.d1.low + yc, 
            self.d1.low + self.d1.ok + yc, 
            color = COLDARKGREEN, 
            linewidth=0)

        self.ax.fill_between(
            self.locx, 
            self.d1.low + self.d1.ok + yc, 
            self.d1.low + self.d1.ok + self.d1.high + yc, 
            color = COLDARKRED, 
            linewidth=0)

        self.ax.fill_between(
            self.locx, 
            self.zero  + yc,
            -self.d2.low + yc, 
            color = COLDARKBLUE, 
            linewidth=0)

        self.ax.fill_between(
            self.locx, 
            -self.d2.low + yc, 
            -self.d2.low - self.d2.ok + yc, 
            color = COLDARKGREEN, 
            linewidth=0)

        self.ax.fill_between(
            self.locx, 
            -self.d2.low - self.d2.ok + yc, 
            -self.d2.low - self.d2.ok - self.d2.high + yc, 
            color = COLDARKRED, 
            linewidth=0)

        # self.ax.plot(self.locx, self.d1.low + yc, zorder=5, 
        #              linewidth=0.5, color='black')
        # self.ax.plot(self.locx, -self.d2.low + yc, zorder=5, 
        #              linewidth=0.5, color='black')

        # self.ax.plot(self.locx, self.d1.low + self.d1.ok + yc, zorder=5, 
        #              linewidth=0.5, color='black')
        # self.ax.plot(self.locx, -self.d2.low - self.d2.ok + yc, zorder=5, 
        #              linewidth=0.5, color='black')

        # self.ax.plot(self.locx, self.d1.low + self.d1.ok + self.d1.high + yc, 
        #              zorder=5, linewidth=0.5, color='black')
        # self.ax.plot(self.locx, -self.d2.low - self.d2.ok -self.d2.high + yc, 
        #              zorder=5, linewidth=0.5, color='black')


        if self.__dict__.has_key('nns'):
            self.l.info("plotting N band")
            self.l.info("maxY %s" % self.plot.maxY)
            self.l.info("max min nns %s %s " % (np.max(self.nns), np.min(self.nns)))
            self.l.info("bandbot, bandtop %s %s " % (self.bandBottom, self.bandTop))

            self.ax.fill_between(
                self.locx,
                yc + (self.plot.maxY * self.nns),
                yc - (self.plot.maxY * self.nns),
                linewidth=0, zorder=-10,
                color=COLLIGHTYELLOW)


        if options.score:

            score1 = 0.5 * self.plot.maxY * (
                1 - 2 * np.exp(-1 * (self.d1.ok / median1))
                + np.exp(-1 * ( ( self.d1.ok + self.d1.low + self.d1.high) / median1))
                )

            score2 = 0.5 * self.plot.maxY * (
                1 - 2 * np.exp(-1 * (self.d2.ok / median2))
                + np.exp(-1 * ( ( self.d2.ok + self.d2.low + self.d2.high) / median2))
                )

            self.ax.plot(self.locx, (score1 - score2) + yc, color='black', linewidth=2)

        self.ax.text(0.01 * len(self.locx), yc + (0.9 * self.plot.maxY), title1, 
                     verticalalignment='center',
                     bbox=dict(facecolor='white', alpha=0.6, zorder=20),
                     )
        self.ax.text(0.01 * len(self.locx), yc - (0.9 * self.plot.maxY), title2,
                     verticalalignment='center',
                     bbox=dict(facecolor='white', alpha=0.6, zorder=20),
                     )


        if self.band > 0:
            self.ax.axhline(
                self.yCorrection + self.plot.maxY, linewidth=3,
                alpha=0.7, color=COLDARKYELLOW, zorder=1000)


        
    def setYticks2(self):
        self.plot.yTicks2.append(self.yCorrection - self.data.median)
        self.plot.yTicks2.append(self.yCorrection)
        self.plot.yTicks2.append(self.yCorrection + self.data.median)
        self.plot.yTickLabels2.append("%s" % -self.data.median)
        self.plot.yTickLabels2.append("0")
        self.plot.yTickLabels2.append("%s" % self.data.median)
        self.ax.axhline(
            self.yCorrection - self.data.median,
            alpha=0.3,
            color='black')
        self.ax.axhline(
            self.yCorrection + self.data.median,
            alpha=0.3,
            color='black')
        
plot.plotBands(plotter)
plot.save(tag='%s.%s.compplot2' % (title1,title2))