dape

#!/usr/bin/env python
"""
DuctApe

Analyze genome(s) and phenome(s)
"""
from ductape import __version__
from ductape.actions import touchProject
from ductape.common.colorlog import ColorFormatter
from ductape.storage.SQLite.database import Organism, Project, Kegg, Biolog
# TODO: handle this imports somewhere else
from matplotlib import cm
from matplotlib import colors
import argparse
import logging.handlers
import os
import sys

__author__ = "Marco Galardini"
__prog__ = "dape"

################################################################################
# Log setup

logger = logging.getLogger('ductape')

################################################################################
# Methods

def dinit(options, wdir, project):
    from ductape.actions import dInit
    if not dInit(project, wdir, options.name, options.descr):
        logger.warning('You can remove or rename the old project file')
        return False
    else:
        return True

def dadd(options, wdir, project):
    from ductape.actions import dAdd
    if not touchProject(project):
        logger.warning('You can setup a new project by running %s init'%
                       __prog__)
        return False
    return dAdd(project, options.orgID, options.n, options.d, options.c)
    
def daddMulti(options, wdir, project):
    from ductape.actions import dAdd
    if not touchProject(project):
        logger.warning('You can setup a new project by running %s init'%
                       __prog__)
        return False
        
    for organism in options.orgID:
        res = dAdd(project, organism)
        if not res:
            logger.warning('Could not add organism (%s)'%organism)
    
    return True

def daddMut(options, wdir, project):
    from ductape.actions import dMutAdd
    if not touchProject(project):
        logger.warning('You can setup a new project by running %s init'%
                       __prog__)
        return False
    return dMutAdd(project, options.mutID, options.m,
                     options.k, options.n, options.d, options.c)

def dstart(options, wdir, project):
    from ductape.actionsterm import fetchKegg
    from ductape.actions import dNet
    from ductape.actions import dCombine
    
    if not touchProject(project):
        logger.warning('You can setup a new project by running %s init'%
                       __prog__)
        return False
    
    # Max value for activity
    biolog = Biolog(project)
    if biolog.isEmpty():
        logger.warning('No phenotypic data is available!')
        return False
    elif biolog.atLeastOneNoParameter():
        logger.warning('Phenome parametrization has not yet been performed!')
        return False
    
    maxAct = biolog.getMaxActivity()
    # Check the activity threshold value
    if options.pthresh >= maxAct:
        logger.warning('The activity threshold is higher than the maximum '+
                       'activity found (%d vs. %d)'%(options.pthresh, maxAct))
        return False
    
    if options.g:
        logger.warning('Skipping mapping to Kegg')
    else:
        # Fetch the Kegg DB?
        if not fetchKegg(project, options.y):
            logger.error('Could not fetch data from KEGG')
            return False
    
    if not dNet(project, options.all, options.paths):
        logger.warning('Combined network analysis failed!')
        return False
    
    return dCombine(project, options.all, options.pthresh)

def dmap(options, wdir, project):
    from ductape.actions import dSetKind, getPathsReacts, getPathsComps
    from ductape.actions import getExclusiveReactions, getExclusiveReactionsMutants
    from ductape.actions import prepareColors, createLegend
    from ductape.kegg.kegg import KeggColor, MapsFetcher
    from ductape.terminal import RunThread
    from ductape.common.utils import rgb_to_hex
    from itertools import combinations
    import numpy as np
    
    if not touchProject(project):
        logger.warning('You can setup a new project by running %s init'%
                       __prog__)
        return False
    
    # Check which maps we have to generate
    # Genome
    proj = Project(project)
    proj.getProject()
    if proj.genome != 'map2kegg':
        logger.warning('Genome mapping to KEGG has not yet been performed!')
    
    # Phenome
    biolog = Biolog(project)
    if biolog.isEmpty():
        logger.warning('No phenotypic data is available!')
        phenome = False
    elif biolog.atLeastOneNoParameter():
        logger.warning('Phenome parametrization has not yet been performed!')
        phenome = False
    else:
        phenome = True
    if phenome and options.s:
        phenome = False
        logger.warning('Skipping the phenomic data')
    
    kegg = Kegg(project)
    
    kind = dSetKind(project)
    
    rpaths = getPathsReacts(project)
    cpaths = getPathsComps(project)
    
    # Check the delta threshold
    maxAct = biolog.getMaxActivity()
    if maxAct is None:
        pass
    elif phenome and maxAct <= options.delta:
        logger.warning('The delta activity threshold is higher than the maximum '+
                       'activity found (%d vs. %d)'%(options.delta,
                                                    maxAct))
        return
    
    if len(options.organisms) == 0:
        # PanGenome or all genomes?
        if kind != 'mutants' and options.all:
            logger.info('Going to generate all the single genomic maps')
            
            organism = Organism(project)
            
            if options.o and phenome:
                if not organism.isOrg(options.o):
                    logger.error('No organism with ID %s is present!'%options.o)
                    return False
                logger.warning('Showing Activity index differences with respect to %s'%options.o)
                
                ref_id = options.o
                
                diffs = [org.org_id
                    for org in organism.getAll()
                    if org.org_id != ref_id]
                
                orgs = set([x.org_id for x in organism.getAll()])
                
                mreacts, ereacts = getExclusiveReactions(project, orgs)
                
                logger.info('Going to generate the metabolic map for reference %s'%
                            ref_id)
                
                # Create the legend
                legend = createLegend('single', project)
                
                if phenome:
                    logger.info('Generating the reactions-only map')
                if not doFetchMaps(project, ref_id, rpaths, cpaths, legend,
                               rorg=mreacts[ref_id], eorg=ereacts[ref_id]):
                    return False
                
                if phenome:
                    for categ in biolog.getCategs():
                        logger.info('Plotting maps for biolog category %s'%categ.category)
                        if not doFetchMaps(project, ref_id, rpaths, cpaths, legend,
                                           categ.category,
                                           rorg=mreacts[ref_id],
                                           eorg=ereacts[ref_id]):
                            return False
                    
                    for org_id in diffs:
                        logger.info('Going to generate the metabolic map for %s (reference %s)'%
                                (org_id,ref_id))
                        
                        # Create the legend
                        legend = createLegend('singlediff', project)
                        
                        if phenome:
                            logger.info('Generating the reactions-only map')
                         
                        # Create the input objects
                        colorPaths = []
                        for path in set(list(rpaths.keys()) + list(cpaths.keys())):
                            KC = KeggColor(path)
                            
                            # Reactions
                            dreact = {}
                            if path in rpaths:
                                for re_id in rpaths[path]:
                                    if re_id in mreacts[org_id]:
                                        dreact[re_id] = colors.cnames['blue']
                                    elif re_id in ereacts[org_id]:
                                        dreact[re_id] = colors.cnames['yellow']
                            KC.setReactions(dreact)
                            KC.setMap(kegg.getPathway(path).html)
                            colorPaths.append(KC)
                            
                        # Go!
                        kmap = MapsFetcher(colorPaths, prefix=org_id, legend=legend)
                
                        if not RunThread(kmap):
                            return False
                        
                        logger.info('%d maps are available in %s'%(len(kmap.pics),
                                                                   kmap._keggroom))
                        
                        if phenome:
                            for categ in biolog.getCategs():
                                logger.info('Plotting maps for biolog category %s'%categ.category)
                                logger.info('Showing differences in the activity >= %d'%options.delta)                            # Get the compounds in the phenomics space
                                corg = {}
                                grey = set()
                                wells = [w for w in biolog.getAllCoByCateg(categ.category)]
                                for well in wells:
                                    mact = biolog.getAvgActivity(well.plate_id, well.well_id, org_id)
                                    ract = biolog.getAvgActivity(well.plate_id, well.well_id, ref_id)
                                    if mact is not None and ract is not None:
                                        diff = ract - mact
                                        if abs(diff) < options.delta:
                                            grey.add(well.co_id)
                                            continue
                                        # Some co_ids are present more than once
                                        if well.co_id not in corg:
                                            corg[well.co_id] = []
                                        corg[well.co_id].append(diff)
                                    else:
                                        grey.add(well.co_id)
                                
                                for k, v in corg.items():
                                    mean = np.array(v).mean()
                                    if mean != np.nan or abs(diff) >= options.delta:
                                        corg[k] = mean
                                    else:
                                        grey.add(k)
                                
                                if len(corg) == 0:
                                    logger.warning('No compounds available for category %s (%s)'%(categ.category, mut_id))
                                    continue
                                
                                # Create the input objects
                                colorPaths = []
                                for path in set(list(rpaths.keys()) + list(cpaths.keys())):
                                    KC = KeggColor(path)

                                    dreact = {}
                                    if path in rpaths:
                                        for re_id in rpaths[path]:
                                            if re_id in mreacts[org_id]:
                                                dreact[re_id] = colors.cnames['blue']
                                            elif re_id in ereacts[org_id]:
                                                dreact[re_id] = colors.cnames['yellow']
                                    
                                    dcomp = {}            
                                    if path in cpaths:
                                        for co_id in cpaths[path]:
                                            if co_id.lstrip('cpd:') in corg:
                                                # We map the values 0-maxAct in a 0-256 window
                                                numcolor = int((corg[co_id.lstrip('cpd:')]*256)/maxAct)
                                                color = cm.PuOr( numcolor )[:3]
                                                color = tuple([int(round(x*255)) for x in color])
                                                dcomp[co_id] = rgb_to_hex(color).upper()
                                            elif co_id.lstrip('cpd:') in grey:
                                                dcomp[co_id] = '#C0C0C0'
                                    
                                    KC.setReactions(dreact)
                                    KC.setCompounds(dcomp)
                                    KC.setMap(kegg.getPathway(path).html)
                                    colorPaths.append(KC)
                                    
                                # Go!
                                prefix = org_id+'_'+categ.category.replace(' ','_').replace('&','and')
                                kmap = MapsFetcher(colorPaths, prefix=prefix, legend=legend)
                        
                                if not RunThread(kmap):
                                    return False
                                
                                logger.info('%d maps are available in %s'%(len(kmap.pics),
                                                                           kmap._keggroom))
                
            elif options.o and not phenome:
                logger.warning('No phenomic data available for the diff mode')
            
            else:
                orgs = [x.org_id for x in organism.getAll()]
                
                mreacts, ereacts = getExclusiveReactions(project, orgs)
        
                for org_id in orgs:
                    logger.info('Going to generate the metabolic map for %s'%org_id)
                
                    # Create the legend
                    legend = createLegend('single', project)
                    
                    if phenome:
                        logger.info('Generating the reactions-only map')
                    if not doFetchMaps(project, org_id, rpaths, cpaths, legend,
                                       rorg=mreacts[org_id], eorg=ereacts[org_id]):
                        return False
                    if phenome:
                        for categ in biolog.getCategs():
                            logger.info('Plotting maps for biolog category %s'%categ.category)
                            if not doFetchMaps(project, org_id, rpaths, cpaths, legend,
                                               categ.category,
                                               rorg=mreacts[org_id],
                                               eorg=ereacts[org_id]):
                                return False
                
        elif proj.isPanGenome() and kind == 'pangenome':
            logger.info('Going to generate a pangenomic metabolic map')
            
            # Some info
            logger.info('Single genome maps can be created using %s map orgID1 orgID2 ...'%
                        __prog__)
            logger.info('All the single genome maps can be created using %s map -a'%
                        __prog__)
            
            conserved = set([r.re_id for r in kegg.getConservedReactions()])
            varlist = [r for r in kegg.getVariableReactions()]
            variable = set([r.re_id for r in varlist])
            nvar = {}
            for r in varlist:
                nvar[r.re_id] = r.orgs
            
            # Go for the colors!
            # Variable
            if len(variable) == 0:
                hexdisp = {}
            else:
                hexdisp = prepareColors(max([r.orgs for r in varlist]),
                                        cm.autumn(np.arange(85,256))[::-1])
            
            # Create the legend
            legend = createLegend('pangenome', project)
            
            if phenome:
                logger.info('Generating the reactions-only map')
            # Create the input objects
            colorPaths = []
            for path in rpaths:
                KC = KeggColor(path)
                dreact = {}
                for re_id in rpaths[path]:
                    if re_id in conserved:
                        dreact[re_id] = colors.cnames['blue']
                    elif re_id in variable:
                        dreact[re_id] = hexdisp[nvar[re_id]]
                KC.setReactions(dreact)
                KC.setMap(kegg.getPathway(path).html)
                colorPaths.append(KC)
                
            # Go!
            kmap = MapsFetcher(colorPaths, prefix='pangenome', legend=legend)
    
            if not RunThread(kmap):
                return False
            
            logger.info('%d maps are available in %s'%(len(kmap.pics),
                                                       kmap._keggroom))

            if phenome:
                for categ in biolog.getCategs():
                    logger.info('Plotting maps for biolog category %s'%categ.category)
                    logger.info(
                            'Showing average differences in the activity >= %d'%options.delta)
                    
                    # Get the compounds in the phenomics space
                    corg = {}
                    # Grey compounds (those below delta)
                    grey = set()
                    # Border compounds (those with at least one difference above delta)
                    border = {}
                    #
                    wells = [w for w in biolog.getAllCoByCateg(categ.category)]
                    for well in wells:
                        acts = [x.avgact 
                                for x in
                                biolog.getAvgActivityEachOrg(well.plate_id,
                                                             well.well_id)]
                        if len(acts) <= 1:
                            continue
                        
                        diffs = []
                        for a, a1 in combinations(acts, 2):
                            diffs.append( abs(a - a1) )
                            
                        avgdiff = np.array(diffs).mean()
                        if avgdiff < options.delta:
                            if max(diffs) >= options.delta:
                                border['cpd:'+well.co_id] = 'red'
                            grey.add(well.co_id)
                            continue
                        
                        # Some co_ids are present more than once
                        if well.co_id not in corg:
                            corg[well.co_id] = []
                        corg[well.co_id].append(avgdiff)
                    
                    toremove = set()
                    for k, v in corg.items():
                        mean = np.array(v).mean()
                        if mean != np.nan:
                            corg[k] = mean
                        else:
                            toremove.add(k)
                    for k in toremove:
                        del corg[k]
                    
                    if len(corg) == 0 and len(grey) == 0:
                        logger.warning('No compounds available for category %s'%(categ.category))
                        continue
                    
                    # Create the input objects
                    colorPaths = []
                    for path in set(list(rpaths.keys()) + list(cpaths.keys())):
                        KC = KeggColor(path)

                        dreact = {}
                        if path in rpaths:
                            for re_id in rpaths[path]:
                                if re_id in conserved:
                                    dreact[re_id] = colors.cnames['blue']
                                elif re_id in variable:
                                    dreact[re_id] = hexdisp[nvar[re_id]]
                        
                        dcomp = {}            
                        if path in cpaths:
                            for co_id in cpaths[path]:
                                if co_id.lstrip('cpd:') in corg:
                                    # We map the values 0-maxAct in a 85-256 window
                                    numcolor = int((corg[co_id.lstrip('cpd:')]*171)/maxAct) + 85
                                    color = cm.Purples( numcolor )[:3]
                                    color = tuple([int(round(x*255)) for x in color])
                                    dcomp[co_id] = rgb_to_hex(color).upper()
                                elif co_id.lstrip('cpd:') in grey:
                                    dcomp[co_id] = '#C0C0C0'
                        
                        KC.setReactions(dreact)
                        KC.setCompounds(dcomp)
                        KC.setBorders(border)
                        KC.setMap(kegg.getPathway(path).html)
                        colorPaths.append(KC)
                        
                    # Go!
                    prefix = 'pangenome_'+categ.category.replace(' ','_').replace('&','and')
                    kmap = MapsFetcher(colorPaths, prefix=prefix, legend=legend)
            
                    if not RunThread(kmap):
                        return False
                    
                    logger.info('%d maps are available in %s'%(len(kmap.pics),
                                                               kmap._keggroom))
                
        elif kind == 'single':
            logger.info('Going to generate a genomic metabolic map')
            
            organism = Organism(project)
            org_id = [x.org_id for x in organism.getAll()][0]
            
            # Create the legend
            legend = createLegend('single', project)
            
            if phenome:
                logger.info('Generating the reactions-only map')
            if not doFetchMaps(project, org_id, rpaths, cpaths, legend):
                return False
            
            if phenome:
                for categ in biolog.getCategs():
                    logger.info('Plotting maps for biolog category %s'%categ.category)
                    if not doFetchMaps(project, org_id, rpaths, cpaths, legend,
                                       categ.category):
                        return False                
        
        elif kind == 'mutants':
            logger.info('Going to generate mutants metabolic map')
        
            organism = Organism(project)
        
            refs = [org.org_id
                    for org in organism.getAll()
                    if not organism.isMutant(org.org_id)]
        
            for ref_id in refs:
                muts = [x for x in organism.getOrgMutants(ref_id)]
                
                mreacts, mix, ereacts = getExclusiveReactionsMutants(project,
                                                                     ref_id,
                                                                     muts)
                
                logger.info('Going to generate the metabolic map for reference %s'%
                            ref_id)
                
                # Create the legend
                legend = createLegend('single', project)
                
                if phenome:
                    logger.info('Generating the reactions-only map')
                if not doFetchMaps(project, ref_id, rpaths, cpaths, legend):
                    return False
                
                if phenome:
                    for categ in biolog.getCategs():
                        logger.info('Plotting maps for biolog category %s'%categ.category)
                        if not doFetchMaps(project, ref_id, rpaths, cpaths, legend,
                                           categ.category):
                            return False
                
                for mut_id in muts:
                    logger.info('Going to generate the metabolic map for mutant %s (parent %s)'%
                            (mut_id,ref_id))
                    
                    # Create the legend
                    legend = createLegend('mutants', project)
                    
                    if phenome:
                        logger.info('Generating the reactions-only map')
                    
                    # Create the input objects
                    colorPaths = []
                    for path in rpaths:
                        KC = KeggColor(path)
                        dreact = {}
                        for re_id in rpaths[path]:
                            if re_id in mix[mut_id]:
                                dreact[re_id] = colors.cnames['green']
                            elif re_id in mreacts[mut_id]:
                                dreact[re_id] = colors.cnames['blue']
                            elif re_id in ereacts[mut_id]:
                                dreact[re_id] = colors.cnames['yellow']
                        KC.setReactions(dreact)
                        KC.setMap(kegg.getPathway(path).html)
                        colorPaths.append(KC)
                        
                    # Go!
                    kmap = MapsFetcher(colorPaths, prefix=mut_id, legend=legend)
            
                    if not RunThread(kmap):
                        return False
                    
                    logger.info('%d maps are available in %s'%(len(kmap.pics),
                                                               kmap._keggroom))
                    
                    if phenome:
                        for categ in biolog.getCategs():
                            logger.info('Plotting maps for biolog category %s'%categ.category)
                            logger.info('Showing differences in the activity >= %d'%options.delta)                            # Get the compounds in the phenomics space
                            corg = {}
                            grey = set()
                            wells = [w for w in biolog.getAllCoByCateg(categ.category)]
                            for well in wells:
                                mact = biolog.getAvgActivity(well.plate_id, well.well_id, mut_id)
                                ract = biolog.getAvgActivity(well.plate_id, well.well_id, ref_id)
                                if mact is not None and ract is not None:
                                    diff = ract - mact
                                    if abs(diff) < options.delta:
                                        grey.add(well.co_id)
                                        continue
                                    # Some co_ids are present more than once
                                    if well.co_id not in corg:
                                        corg[well.co_id] = []
                                    corg[well.co_id].append(diff)
                                else:
                                    grey.add(well.co_id)
                            
                            for k, v in corg.items():
                                mean = np.array(v).mean()
                                if mean != np.nan or abs(diff) >= options.delta:
                                    corg[k] = mean
                                else:
                                    grey.add(k)
                            
                            if len(corg) == 0:
                                logger.warning('No compounds available for category %s (%s)'%(categ.category, mut_id))
                                continue
                            
                            # Create the input objects
                            colorPaths = []
                            for path in set(list(rpaths.keys()) + list(cpaths.keys())):
                                KC = KeggColor(path)
                                
                                dreact = {}
                                if path in rpaths:
                                    for re_id in rpaths[path]:
                                        if re_id in mix[mut_id]:
                                            dreact[re_id] = colors.cnames['green']
                                        elif re_id in mreacts[mut_id]:
                                            dreact[re_id] = colors.cnames['blue']
                                        elif re_id in ereacts[mut_id]:
                                            dreact[re_id] = colors.cnames['yellow']
                                
                                dcomp = {}            
                                if path in cpaths:
                                    for co_id in cpaths[path]:
                                        if co_id.lstrip('cpd:') in corg:
                                            # We map the values 0-maxAct in a 0-256 window
                                            numcolor = int((corg[co_id.lstrip('cpd:')]*256)/maxAct)
                                            color = cm.PuOr( numcolor )[:3]
                                            color = tuple([int(round(x*255)) for x in color])
                                            dcomp[co_id] = rgb_to_hex(color).upper()
                                        elif co_id.lstrip('cpd:') in grey:
                                            dcomp[co_id] = '#C0C0C0'
                                
                                KC.setReactions(dreact)
                                KC.setCompounds(dcomp)
                                KC.setMap(kegg.getPathway(path).html)
                                colorPaths.append(KC)
                                
                            # Go!
                            prefix = mut_id+'_'+categ.category.replace(' ','_').replace('&','and')
                            kmap = MapsFetcher(colorPaths, prefix=prefix, legend=legend)
                    
                            if not RunThread(kmap):
                                return False
                            
                            logger.info('%d maps are available in %s'%(len(kmap.pics),
                                                                       kmap._keggroom))
        else:
            logger.warning('Unattended case %s'%kind)
            return False
            
    elif len(options.organisms) == 1:
        org_id = set(options.organisms).pop()
        
        organism = Organism(project)
        
        if not organism.isOrg(org_id):
            logger.warning('Organism %s is not present!'%org_id)
            return False
        
        logger.info('Going to generate the metabolic map for %s'%org_id)
        
        # Create the legend
        legend = createLegend('single', project)
        
        if phenome:
            logger.info('Generating the reactions-only map')
        if not doFetchMaps(project, org_id, rpaths, cpaths, legend):
            return False
        if phenome:
            for categ in biolog.getCategs():
                logger.info('Plotting maps for biolog category %s'%categ.category)
                if not doFetchMaps(project, org_id, rpaths, cpaths, legend,
                                   categ.category):
                    return False
                    
    else:
        orgs = set(options.organisms)
        
        organism = Organism(project)
        
        # Check
        for org_id in orgs:
            if not organism.isOrg(org_id):
                logger.warning('Organism %s is not present: skipping'%org_id)
                orgs.remove(org_id)
                continue
            
        mreacts, ereacts = getExclusiveReactions(project, orgs)
        
        for org_id in orgs:
            logger.info('Going to generate the metabolic map for %s'%org_id)
        
            # Create the legend
            legend = createLegend('single', project)
            
            if phenome:
                logger.info('Generating the reactions-only map')
            if not doFetchMaps(project, org_id, rpaths, cpaths, legend,
                               rorg=mreacts[org_id], eorg=ereacts[org_id]):
                return False
            if phenome:
                for categ in biolog.getCategs():
                    logger.info('Plotting maps for biolog category %s'%categ.category)
                    if not doFetchMaps(project, org_id, rpaths, cpaths, legend,
                                       categ.category,
                                       rorg=mreacts[org_id],
                                       eorg=ereacts[org_id]):
                        return False
    
    return True

def dimport(options, wdir, project):
    from ductape.actions import dKeggImport
    if not touchProject(project):
        logger.warning('You can setup a new project by running %s init'%
                       __prog__)
        return False
    return dKeggImport(project, options.file)

def dexport(options, wdir, project):
    from ductape.actions import dKeggExport
    if not touchProject(project):
        logger.warning('You can setup a new project by running %s init'%
                       __prog__)
        return False
    return dKeggExport(project)

def dremove(options, wdir, project):
    from ductape.actions import dRemove
    if not touchProject(project):
        logger.warning('Nothing to be removed!\n'+
                       'You can setup a new project by running %s init'%
                       __prog__)
        return False
    return dRemove(project, options.organisms)

def dclear(options, wdir, project):
    from ductape.actions import dClear
    if not touchProject(project):
        logger.warning('Nothing to be cleaned up!\n'+
                       'You can setup a new project by running %s init'%
                       __prog__)
        return False
    return dClear(project, options.keep_org, options.keep_kegg)

def doFetchMaps(project, org_id, rpaths, cpaths, legend=None, category=None,
                rorg=set(), eorg=set()):
    from ductape.kegg.kegg import KeggColor, MapsFetcher
    from ductape.terminal import RunThread
    from ductape.common.utils import rgb_to_hex
    import numpy as np
    
    kegg = Kegg(project)
    biolog = Biolog(project)
    
    # Get the reactions in the organism space
    if len(rorg) == 0:
        rorg = set()
        for oR in kegg.getOrgReact(org_id):
            rorg.add(oR.re_id)
        
    # Get the compounds in the phenomics space
    corg = {}
    if category:
        wells = [w for w in biolog.getAllCoByCateg(category)]
    else:
        wells = []
    for well in wells:
        act = biolog.getAvgActivity(well.plate_id, well.well_id, org_id)
        if act is not None:
            # Some co_ids are present more than once
            if well.co_id not in corg:
                corg[well.co_id] = []
            corg[well.co_id].append(act)
    
    toremove = set()
    for k, v in corg.items():
        mean = np.array(v).mean()
        if mean != np.nan:
            corg[k] = mean
        else:
            toremove.add(k)
    for k in toremove:
        del corg[k]
    
    if category and len(corg) == 0:
        logger.warning('No compounds available for category %s (%s)'%(category, org_id))
        return True
     
    # Create the input objects
    colorPaths = []
    for path in set(list(rpaths.keys()) + list(cpaths.keys())):
        KC = KeggColor(path)
        
        # Reactions
        dreact = {}
        if path in rpaths:
            for re_id in rpaths[path]:
                if re_id in rorg:
                    dreact[re_id] = colors.cnames['blue']
                elif re_id in eorg:
                    dreact[re_id] = colors.cnames['yellow']
                    
        # Compounds
        dcomp = {}
        if path in cpaths:
            for co_id in cpaths[path]:
                if co_id.lstrip('cpd:') in corg:
                    # We map the values 0-maxAct in a 0-256 window
                    numcolor = int((corg[co_id.lstrip('cpd:')]*256)/biolog.getMaxActivity())
                    color = cm.RdYlGn( numcolor )[:3]
                    color = tuple([int(round(x*255)) for x in color])
                    dcomp[co_id] = rgb_to_hex(color).upper()
        
        KC.setReactions(dreact)
        KC.setCompounds(dcomp)
        KC.setMap(kegg.getPathway(path).html)
        colorPaths.append(KC)
        
    # Go!
    if category:
        prefix = org_id + '_' + category.replace(' ','_').replace('&','and')
    else:
        prefix = org_id
    kmap = MapsFetcher(colorPaths, prefix=prefix, legend=legend)

    if not RunThread(kmap):
        return False
    
    logger.info('%d maps are available in %s'%(len(kmap.pics), kmap._keggroom))
    
    return True

################################################################################
# Options

def getOptions():
    description = "Analyze genomes & phenomes"
    parser = argparse.ArgumentParser(description = description,
                                     prog=__prog__)
    parser.add_argument('-p', metavar='project', action='store',
                        dest='project',
                        default='ductape.db',
                        help='Project file')
    parser.add_argument('-w', metavar='workdir', action='store', dest='wdir',
                        default='.',
                        help='Working directory')
    parser.add_argument('-v', action='count',
                        default=0,
                        help='Increase verbosity level')
    parser.add_argument('--version', action='version',
                        version='%(prog)s '+__version__)
    subparsers = parser.add_subparsers()

    parser_init = subparsers.add_parser('init', help='Initialize the project')
    parser_init.add_argument('-n', action="store",
                             dest='name',
                             default = 'Project',
                             help='Project name')
    parser_init.add_argument('-d', metavar='descr', action="store",
                            dest='descr',
                            default = 'DuctApe project',
                            help='Project description')
    parser_init.set_defaults(func=dinit)
    
    parser_add = subparsers.add_parser('add',
                           help='Add an organism')
    parser_add.add_argument('orgID', action='store',
                            help='Organism ID')
    parser_add.add_argument('-n', metavar='name', action="store",
                            default = '',
                            help='Organism name')
    parser_add.add_argument('-d', metavar='descr', action="store",
                            default = '',
                            help='Organism description')
    parser_add.add_argument('-c', metavar='color', action="store",
                            default = '',
                            help='Organism color')
    parser_add.set_defaults(func=dadd)
    
    parser_add_multi = subparsers.add_parser('add-multi',
                       help='Add multiple organisms (providing just the ID)')
    parser_add_multi.add_argument('orgID', action='store', nargs='+',
                            help='Organism ID')
    parser_add_multi.set_defaults(func=daddMulti)
    
    parser_add_mut = subparsers.add_parser('add-mut',
                               help='Add a mutant of an existing organism')
    parser_add_mut.add_argument('mutID', action='store',
                            help='Mutant organism ID')
    parser_add_mut.add_argument('-m', metavar='mutparent', action="store",
                            required = True,
                            help='This mutant parent orgID')
    parser_add_mut.add_argument('-k', metavar='kind', action="store",
                            choices = ['deletion', 'insertion'],
                            default = 'deletion',
                            help='This mutant kind (deletion|insertion)')
    parser_add_mut.add_argument('-n', metavar='name', action="store",
                            default = '',
                            help='Organism name')
    parser_add_mut.add_argument('-d', metavar='descr', action="store",
                            default = '',
                            help='Organism description')
    parser_add_mut.add_argument('-c', metavar='color', action="store",
                            default = '',
                            help='Organism color')
    parser_add_mut.set_defaults(func=daddMut)
    
    parser_start = subparsers.add_parser('start', help='Start the analysis')
    parser_start.add_argument('-g', action="store_true",
                            default=False,
                            help='Skip Kegg mapping')
    parser_start.add_argument('-y', action="store_true",
                            default=False,
                            help='Try to fetch Kegg data even while encountering failures')
    parser_start.add_argument('-a', '--all', action="store_true",
                            default=False,
                            help='Create single organisms net')
    parser_start.add_argument('-s', '--paths', action="store_true",
                            default=False,
                            help='Save each pathway network')
    parser_start.add_argument('-t', '--phenome-threshold', metavar='pthresh',
                              action="store", dest='pthresh',
                            type=float,
                            default=5,
                            help='Activity threshold for the combined matrix '+
                                '(suggested: half of the maximum value)')
    parser_start.set_defaults(func=dstart)
    
    parser_map = subparsers.add_parser('map', help='Fetch metabolic maps')
    parser_map.add_argument('-a', '--all', action="store_true",
                            default=False,
                            help='Plot all the single organisms maps')
    parser_map.add_argument('-o', metavar='difforg', action="store",
                            default = None,
                            help='Diff mode: reference organism ID')
    parser_map.add_argument('-d', metavar='delta', action="store", dest='delta',
                              type=int,
                              default=3,
                              help='Activity delta threshold')
    parser_map.add_argument('-s', action="store_true",
                            default=False,
                            help='Skip the phenomic data')
    parser_map.add_argument('organisms', metavar='orgID', nargs='*',
                            action="store",
                            default=[],
                            help='Organism(s) to be mapped')
    parser_map.set_defaults(func=dmap)
    
    parser_import = subparsers.add_parser('import', help='Import kegg data')
    parser_import.add_argument('file', action="store",
                            help='Kegg dump file')
    parser_import.set_defaults(func=dimport)
    
    parser_export = subparsers.add_parser('export', help='Export kegg data')
    parser_export.set_defaults(func=dexport)
    
    parser_rm = subparsers.add_parser('rm', help='Completely remove an organism')
    parser_rm.add_argument('organisms', metavar='orgID', nargs='+',
                              action="store",
                            help='Organism(s) to be removed')
    parser_rm.set_defaults(func=dremove)
    
    parser_clear = subparsers.add_parser('clear',
                                         help='Clear all the results')
    parser_clear.add_argument('-o', '--keep-org',
                            action="store_true",
                            default=False,
                            help='Keep organisms data')
    parser_clear.add_argument('-k', '--keep-kegg', 
                            action="store_true",
                            default=False,
                            help='Keep KEGG data')
    parser_clear.set_defaults(func=dclear)
    
    return parser.parse_args()

################################################################################

options = getOptions()

logger.setLevel(logging.DEBUG)

ch = logging.StreamHandler()
if options.v == 0:
    ch.setLevel(logging.INFO)
elif options.v >= 1:
    ch.setLevel(logging.DEBUG)
formatter = ColorFormatter('%(asctime)s - $COLOR%(message)s$RESET','%H:%M:%S')
ch.setFormatter(formatter)
logger.addHandler(ch)

fh = logging.handlers.RotatingFileHandler('ductape.log', maxBytes=2000000)
formatter = logging.Formatter('%(asctime)s - %(name)s - [%(levelname)s] - %(message)s',
                            '%Y-%m-%d %H:%M:%S')
fh.setFormatter(formatter)
logger.addHandler(fh)

wdir = os.path.abspath(options.wdir)
if not os.path.exists(wdir):
    try:
        os.mkdir(wdir)
    except:
        logger.error('Could not create working directory %s'%wdir)
    
project = os.path.join(wdir, options.project)

ret = options.func(options, wdir, project)

touchProject(project)

if not ret:
    sys.exit(1)