getDownstream

#!/usr/bin/python
'''
From a regulatory network get all the downstream genes
(including nested regulators and operons)

The network should be in GML format
'''

def getOptions():
    import argparse

    # create the top-level parser
    description = ("From a regulatory network get all the downstream genes")
    parser = argparse.ArgumentParser(description = description)
    parser.add_argument('-N', '--name', action="store_true",
            default=False,
            dest='name',
            help='The regulator name is provided [Default: locus_tag]')
    parser.add_argument('-G', '--genes', action="store_true",
            default=False,
            dest='genes',
            help='Print the number of genes and operons')
    parser.add_argument('gmlfile', action='store', 
                            help='Regulatory network GML file')
    parser.add_argument('regulator', action='store',
                            help='Regulator to be searched')
    return parser.parse_args()

options = getOptions()

import networkx as nx

net = nx.parse_gml(open(options.gmlfile))

import sys
# Is the regulator really inside the network
if options.regulator not in net.nodes() and not options.name:
    if options.genes:
        print '\t'.join([options.regulator, 'NA', 'NA'])
    sys.exit(0)
    #raise KeyError('%s not present in network %s'%(options.regulator, options.gmlfile))
if options.regulator not in set(net.node[x].get('name', None) for x in net) and options.name:
    if options.genes:
        print '\t'.join([options.regulator, 'NA', 'NA'])
    sys.exit(0)
    #raise KeyError('%s not present in network %s'%(options.regulator, options.gmlfile))

if options.name:
    for n in net:
        name = net.node[n].get('name', None)
        if name == options.regulator:
            snode = n
else:
    snode = options.regulator

downstream = set()

genes = set()
operons = set()

nodes = nx.depth_first_search.dfs_successors(net, snode)
for n in nodes:
    if net.node[n]['kind'] == 'operon':
        operons.add(n)
    else:
        genes.add(n)

    downstream.add(n)
    for n1 in nodes[n]:
        if net.node[n1]['kind'] == 'operon':
            operons.add(n1)
        else:
            genes.add(n1)

        downstream.add(n1)

if options.genes:
    print '\t'.join([options.regulator, str(len(genes)), str(len(operons))])
else:
    for n in downstream:
        print n