def main(argv):

    parser = E.OptionParser(version="%prog version: $Id$",
                            usage=globals()["__doc__"])

    parser.add_option("-o", "--filename-ontology", dest="filename_obo",
                      type="string",
                      help="filename with ontology information "
                      "in obo-xml format. "
                      " The latest version can always be"
                      " retrieved at "
                      "`wget http://archive.geneontology.org/"
                      "latest-termdb/go_daily-termdb.obo-xml.gz`."
                      " [default=%default].")

    parser.add_option("-g", "--go-tsv-file", dest="filename_go", type="string",
                      help="filename with gene to go assignments "
                      " This is a tab-separated file with the columns "
                      " <namespace> <gene_id> <goid> <description> <evidence>."
                      " The evidence column is currently ignored. "
                      " [default=%default].")

    parser.add_option("-p", "--pvalues-tsv-file", dest="filename_pvalues",
                      type="string",
                      help="filename with pvalues for each GO term "
                      "This is a tab-separated file with the columns "
                      "<goid> <pvalue>. "
                      "If a third column is present it will be "
                      "interpreted as fold-change. "
                      " [default=%default].")

    parser.add_option("--ontology", dest="ontology", type="choice",
                      action="append",
                      choices=(
                          "biol_process", "cell_location",
                          "mol_function", "all"),
                      help="ontologies to analyze. Ontologies are "
                      "tested separately."
                      " 'all' will test all ontologies. "
                      " [default=%default].")

    parser.add_option("--max-similarity", dest="max_simrel_threshold",
                      type="float",
                      help="cluster until simrel threshold is achieved "
                      "[default=%default].")

    parser.add_option("--min-clusters", dest="min_clusters_threshold",
                      type="float",
                      help="cluster only until at most # clusters remain "
                      "[default=%default].")

    parser.add_option("--palette", dest="palette", type="choice",
                      choices=("rainbow", "gray", "blue-white-red",
                               "autumn", "bone", "cool",
                               "copper", "flag", "gray",
                               "hot", "hsv", "jet", "pink", "prism",
                               "spring", "summer", "winter", "spectral",
                               "RdBu", "RdGy", "BrBG", "BuGn", "Blues",
                               "Greens", "Reds", "Oranges", "Greys"),
                      help="colour palette [default=%default]")

    parser.add_option("--reverse-palette", dest="reverse_palette",
                      action="store_true",
                      help="reverse colour palette [default=%default].")

    parser.set_defaults(filename_obo='go_daily-termdb.obo-xml.gz',
                        filename_go="go.tsv.gz",
                        filename_pvalues=None,
                        filename_bg=None,
                        filename_matrix=None,
                        filename_frequencies=None,
                        ontology=[],
                        # 0.9, 0.7, 0.5 and 0.4
                        # 0.4 very small clusters
                        max_simrel_threshold=0.7,
                        min_frequency=0.05,
                        min_pvalue=0.05,
                        min_child_threshold=0.75,
                        min_clusters_threshold=0,
                        palette="RdBu",
                        reverse_palette=False,
                        )

    options, args = E.Start(parser, argv, add_output_options=True)

    map_ontology2namespace = {
        'biol_process': "biological_process",
        'cell_location': 'cellular_component',
        'mol_function': 'molecular_function'}

    if "all" in options.ontology:
        options.ontology = map_ontology2namespace.keys()

    #########################################
    if options.filename_pvalues is None:
        E.info("reading pvalues from stdin")
        all_term2pvalue, all_term2log2fold = readPValues(options.stdin)
    else:
        E.info("reading pvalues from %s " % options.filename_pvalues)
        with IOTools.openFile(options.filename_pvalues) as infile:
            all_term2pvalue, all_term2log2fold = readPValues(infile)

    E.info("read %i pvalue and %i fold assignments" %
           (len(all_term2pvalue), len(all_term2log2fold)))

    for test_ontology in options.ontology:

        E.info("working on namespace %s" % test_ontology)

        namespace = map_ontology2namespace[test_ontology]

        ################################################################
        # compute or load ancestry information for all nodes
        fn = buildFilename(test_ontology, "ancestors", "tsv.gz")
        if os.path.exists(fn):
            E.info("reading ancestor information from %s" % fn)
            with IOTools.openFile(fn) as infile:
                ancestors = readAncestors(infile)
        else:
            graph = readGOGraph(options.filename_obo, namespace)
            ancestors = computeAncestors(graph, graph.nodes())
            with IOTools.openFile(fn, "w") as outfile:
                outfile.write("goid\tancestors\n")
                outfile.write(
                    "".join(["%s\t%s\n" %
                             (x, ";".join(v))
                             for x, v in ancestors.iteritems()]))

        #########################################
        graph = readGOGraph(options.filename_obo, namespace)

        #########################################
        E.info("reading gene to go assignments from %s " % options.filename_go)
        all_gene2gos, all_go2infos = \
            GO.ReadGene2GOFromFile(
                IOTools.openFile(options.filename_go),
                synonyms=graph.synonyms,
                obsolete=graph.obsoletes)

        #########################################
        # filter pvalues
        synonyms = graph.synonyms
        term2pvalue = dict([(synonyms.get(x, x), y)
                            for x, y in all_term2pvalue.iteritems()])
        term2log2fold = dict([(synonyms.get(x, x), y)
                              for x, y in all_term2log2fold.iteritems()])

        if len(term2pvalue) == 0:
            E.warn("no data - no output produced")
            E.Stop()
            return

        go2info = all_go2infos[test_ontology]
        go2info = collections.defaultdict(
            str, dict([(x.mGOId, x.mDescription) for x in go2info.values()]))
        gene2gos = all_gene2gos[test_ontology]

        if options.filename_bg:
            background = IOTools.readList(
                IOTools.openFile(options.filename_bg))
        else:
            background = list(gene2gos.keys())

        # invert mapping, at the same time insert missing ancestral counts
        go2genes = GO.buildGO2Genes(gene2gos, ancestors=ancestors)

        ################################################################
        # compute or load term frequences
        fn = buildFilename(test_ontology, "termfrequencies", "tsv.gz")
        if os.path.exists(fn):
            E.info("reading term frequencies from %s" % fn)
            with IOTools.openFile(fn) as infile:
                counts, p = readTermFrequencies(infile)
        else:
            counts, p = computeTermFrequencies(graph, go2genes)
            with IOTools.openFile(fn, "w") as outfile:
                outfile.write("goid\tfrequency\tcounts\n")
                for k in counts.keys():
                    outfile.write("%s\t%f\t%s\n" %
                                  (k, p[k], ";".join(counts[k])))

        # compute or load simrel matrix
        fn = buildFilename(test_ontology, "simrel", "matrix.gz")
        if os.path.exists(fn):
            E.info("reading simrel matrix from %s" % fn)
            with IOTools.openFile(fn, "r") as infile:
                matrix, terms = readSimrelMatrix(infile)
        else:
            matrix, terms = buildSimrelMatrix(go2genes, go2info, ancestors, p)
            E.info("writing simrel matrix to file")
            with IOTools.openFile(fn, "w") as outfile:
                IOTools.writeMatrix(outfile, matrix, terms, terms)

        fn = buildFilename(test_ontology, "cluster", "matrix.gz")

        if os.path.exists(fn):
            E.info("reading clustered matrix matrix from %s" % fn)
            with IOTools.openFile(fn, "r") as infile:
                matrix, terms = readSimrelMatrix(infile)
        else:
            matrix, terms, clusters = clusterSimrelMatrix(matrix, terms,
                                                          ancestors, p, counts,
                                                          term2pvalue, go2info,
                                                          options)

            if matrix is None:
                E.warn("%s: no output" % test_ontology)
                continue

            E.info("after clustering: %i x %i matrix" %
                   (len(terms), len(terms)))
            with IOTools.openFile(fn, "w") as outfile:
                IOTools.writeMatrix(outfile, matrix, terms, terms)

            fn = buildFilename(test_ontology, "cluster", "tsv")

            with IOTools.openFile(fn, "w") as outfile:
                outfile.write(
                    "\t".join(("rep",
                               "mem",
                               "pvalue",
                               "frequency",
                               "description")) + "\n")
                for rep, members in clusters.iteritems():
                    for mem in members:
                        outfile.write("%s\t%s\t%f\t%f\t%s\n" %
                                      (rep,
                                       mem,
                                       term2pvalue[mem],
                                       p[mem],
                                       go2info.get(mem, "na")))

        # visualize output
        rows = len(terms)

        # revigo
        #     color = pvalue, size = frequency, egdewith = simrel score
        # me:
        #     size  = pvalue
        #     color = fold (blueish for depletion, reddish for enrichment),
        #     edgewith = can't be done with networkx? (only edgecolor)

        # build graph from matrix
        term_graph = networkx.Graph()
        figure = plt.figure()

        term_graph.add_nodes_from(terms)
        for i, j in itertools.combinations(xrange(rows), 2):
            w = matrix[i][j]
            if w < 0.1:
                continue
            term_graph.add_edge(terms[i], terms[j], weight=w)

        layout = networkx.spring_layout(term_graph)

        if options.reverse_palette:
            color_scheme = eval("matplotlib.cm.%s_r" % options.palette)
        else:
            color_scheme = eval("matplotlib.cm.%s" % options.palette)

        labels = dict([(x, go2info.get(x, x)) for x in terms])
        node_size = [-20.0 * math.log(term2pvalue[x], 10)
                     for x in term_graph.nodes()]
        for x in term2pvalue.keys():
            term2log2fold[x] = numpy.random.normal(0, 5)
        node_color = [term2log2fold[x] for x in term_graph.nodes()]

        networkx.draw_networkx_edges(term_graph, layout, edge_color='0.50')

        vmin = min(node_color)
        vmax = max(node_color)
        vmin = min(vmin, -vmax)
        vmax = max(vmax, -vmin)

        networkx.draw_networkx_nodes(term_graph, layout,
                                     node_size=node_size,
                                     node_color=node_color,
                                     cmap=color_scheme,
                                     alpha=0.8,
                                     vmin=vmin,
                                     vmax=vmax,
                                     linewidths=0.5)

        positions = dict([(x, numpy.array((y[0], y[1] - 0.02)))
                          for x, y in layout.iteritems()])

        networkx.draw_networkx_labels(term_graph,
                                      positions,
                                      labels=labels,
                                      font_size=5,
                                      alpha=0.5)

        # turn of tick marks
        gca = figure.gca()
        gca.set_xticklabels([], visibile=False)
        gca.set_yticklabels([], visibile=False)

        fn = buildFilename(test_ontology, "graph", "svg")
        plt.colorbar()
        plt.savefig(fn)

        fn = buildFilename(test_ontology, "graph", "tsv.gz")
        with IOTools.openFile(fn, "w") as outfile:
            outfile.write("goid\tx\ty\tpvalue\tl2fold\tdescription\n")
            for term in terms:
                outfile.write("%s\t%f\t%f\t%f\t%f\t%s\n" %
                              (term,
                               layout[term][0],
                               layout[term][1],
                               term2pvalue[term],
                               term2log2fold[term],
                               go2info.get(term, "na")))

    E.Stop()
dataset/ETHPy150Open CGATOxford/cgat/scripts/revigo.py/main/original