draw_subnetwork.py: add a known network (G0) from the thermomorphogenesis paper.

Fixed a bug.

Now I close figure (plt.close()) before creating a new one, to avoid that the current figure is drawn
on top of the old one.

-Hui

1 files changed, 120 insertions, 41 deletions

diff --git a/Code/draw_subnetwork.py b/Code/draw_subnetwork.py
index 79ab5cd..d05f2d4 100644
--- a/Code/draw_subnetwork.py
+++ b/Code/draw_subnetwork.py
@@ -1,7 +1,10 @@
 # Usage: python draw_subnetwork.py edges.txt

-# Purpose: draw a sub-network given a list of genes.  

+# Purpose: draw a sub-network given a list of genes from thermomorphogenesis paper.

+#          The paper Molecular and genetic control of plant thermomorphogenesis. https://doi.org/10.1038/nplants.2015.190

+#

 # Created on 5 December 2019 by Hui Lan (lanhui@zjnu.edu.cn)

 

+

 import os, sys

 import networkx as nx

 import pylab as plt

@@ -9,6 +12,7 @@ import glob
 import math

 from networkx.algorithms.distance_measures import diameter, eccentricity

 

+

 def build_network_from_file(edge_fname, gene_lst, gene_dict):

     G = nx.DiGraph()    

 

@@ -46,6 +50,33 @@ def build_network_from_file(edge_fname, gene_lst, gene_dict):
     return G

 

 

+

+def build_thermomorphogenesis_network(gene_lst, gene_dict):

+    ''' Edges from thermo.png in my asus laptop. '''

+    G = nx.DiGraph()

+

+    for g in gene_lst:

+        G.add_node(gene_dict[g])

+

+    pairs = [('PIF4', 'PAR1'), ('PIF4', 'PRE1'), ('PIF4', 'YUC8'), ('PIF4', 'TAA1'), ('PIF4', 'IAA4'), ('PIF4', 'SAUR21'),

+             ('HY5', 'PAR1'), ('HY5', 'PRE1'), ('HY5', 'YUC8'), ('HY5', 'TAA1'), ('HY5', 'IAA4'), ('HY5', 'SAUR21'),

+             ('FCA', 'PAR1'), ('FCA', 'PRE1'), ('FCA', 'YUC8'), ('FCA', 'TAA1'), ('FCA', 'IAA4'), ('FCA', 'SAUR21'),

+             ('BZR1', 'PAR1'), ('BZR1', 'PRE1'), ('BZR1', 'YUC8'), ('BZR1', 'TAA1'), ('BZR1', 'IAA4'), ('BZR1', 'SAUR21'),

+             ('ARF6', 'PAR1'), ('ARF6', 'PRE1'), ('ARF6', 'YUC8'), ('ARF6', 'TAA1'), ('ARF6', 'IAA4'), ('ARF6', 'SAUR21'),

+             ('PAR1', 'IBH1'), ('PRE1', 'IBH1'),

+             ('PAR1', 'PIF4'), ('PRE1', 'PIF4'),

+             ('IBH1', 'HBI1'),

+             ('IAA4', 'ARF6'),

+             ('ARF6', 'SAUR21')

+    ] # see paper Molecular and genetic control of plant thermomorphogenesis. https://doi.org/10.1038/nplants.2015.190

+

+    for (g2, g1) in pairs:  # g2 is source, g1 is target

+        G.add_edge(g2, g1, weight=2) # g2 is source, and g1 is target

+

+

+    return G

+

+

 def compute_total_edge_weight(edges, G):

     total = 0

     for e in edges:

@@ -71,18 +102,22 @@ def draw_graph(G, fname):
     nx.draw_networkx_edges(G,pos,edgelist=esmall,width=1,alpha=0.1,edge_color='k',style='dashed')

     nx.draw_networkx_labels(G,pos,font_size=8,font_color='k',font_family='sans-serif')

     plt.axis('off')

-    plt.savefig(fname) # save as png

+    plt.savefig(fname) 

+    plt.close()

     #plt.show() # display    

 

 

 def better_date(s):

+    ''' Add a dash between year and month, and a dash between month and day.'''

     if len(s) == 8:

         return '-'.join([s[:4], s[4:6], s[6:]])

     else:

         return s

 

+

 def draw_graph2(G, fname, date):

-    pos=nx.circular_layout(G)

+

+    pos = nx.circular_layout(G)

     all_edges = []

     all_widths = []

 

@@ -91,11 +126,29 @@ def draw_graph2(G, fname, date):
         all_widths.append(math.sqrt(d['weight']))

 

     nx.draw_networkx_nodes(G,pos,alpha=0.05)

-    nx.draw_networkx_edges(G,pos,edgelist=all_edges,width=all_widths,alpha=0.05,edge_color='k',style='dashed')

+    nx.draw_networkx_edges(G,pos,edgelist=all_edges,width=all_widths,alpha=0.2,edge_color='k',style='dashed')

     nx.draw_networkx_labels(G,pos,font_size=11,font_color='b',font_family='sans-serif')

     plt.axis('off')

     plt.title(better_date(date))

-    plt.savefig(fname) # save as png

+    plt.savefig(fname) 

+    plt.close()

+    #plt.show() # display    

+

+

+def draw_graph3(G, fname):

+

+    pos = nx.circular_layout(G)

+    all_edges = []

+

+    for (u, v, d) in G.edges(data=True):

+        all_edges.append((u, v))

+

+    nx.draw_networkx_nodes(G,pos,alpha=0.05)

+    nx.draw_networkx_edges(G,pos,edgelist=all_edges,alpha=0.2,edge_color='k',style='dashed')

+    nx.draw_networkx_labels(G,pos,font_size=11,font_color='b',font_family='sans-serif')

+    plt.axis('off')

+    plt.savefig(fname) 

+    plt.close() # it is important to close the plot before creating another one.  Otherwise, plots will overlap.

     #plt.show() # display    

     

 

@@ -165,6 +218,7 @@ thermomorphogenesis_genes_small = [
     'AT5G11260',	#HY5

     'AT2G42870',	#PAR1

     'AT5G39860',	#PRE1

+    'AT5G43700',        #IAA4

     'AT4G16280',	#FCA

     'AT2G43060',	#IBH1

     'AT2G18300',	#HBI1

@@ -172,18 +226,19 @@ thermomorphogenesis_genes_small = [
     'AT1G70560',	#TAA1

     'AT1G30330',	#ARF6

     'AT1G19850',        #ARF5

-    'AT1G75080',	#BZR1

-    'AT2G25930',	#ELF3

-    'AT2G40080',	#ELF4

-    'AT3G46640',	#LUX

+    'AT5G01830',	#SAUR21    

+    'AT1G75080'   	#BZR1

+#    'AT2G25930',	#ELF3

+#    'AT2G40080',	#ELF4

+#    'AT3G46640',	#LUX

 ]

 

-

 gene_dict = {

     'AT2G43010':'PIF4',

     'AT5G11260':'HY5',

     'AT2G42870':'PAR1',

     'AT5G39860':'PRE1',

+    'AT5G43700':'IAA4',

     'AT4G16280':'FCA',

     'AT2G43060':'IBH1',

     'AT2G18300':'HBI1',

@@ -191,13 +246,18 @@ gene_dict = {
     'AT1G70560':'TAA1',

     'AT1G30330':'ARF6',

     'AT1G19850':'ARF5',

+    'AT5G01830':'SAUR21',

     'AT1G75080':'BZR1',

     'AT2G25930':'ELF3',

     'AT2G40080':'ELF4',

     'AT3G46640':'LUX'

 }

 

-print('Make sub graph...')

+print('Make sub graphs ...')

+

+G0 = build_thermomorphogenesis_network(thermomorphogenesis_genes_small, gene_dict)

+print('Number of edges in the paper thermomorphogenesis is %d' % (len(G0.edges())))

+draw_graph2(G0, '../Data/temp/graph-%s.pdf' % ('20160101'), '')

 

 graph_lst = []

 graph_names = []

@@ -206,43 +266,62 @@ for fname in sorted(glob.glob('../Analysis/edges.txt.2019*')):
         continue

     print(fname)

     graph_names.append(fname)

-    graph_lst.append(build_network_from_file(fname, thermomorphogenesis_genes_small, gene_dict))

+    G = build_network_from_file(fname, thermomorphogenesis_genes_small, gene_dict)

+    graph_lst.append(G)

 

 

-#G1205 = build_network_from_file('../Analysis/edges.txt.20191205', thermomorphogenesis_genes)

-#G1203 = build_network_from_file('../Analysis/edges.txt.20191203', thermomorphogenesis_genes)

-#G1126 = build_network_from_file('../Analysis/edges.txt.20191126', thermomorphogenesis_genes)

-#G1108 = build_network_from_file('../Analysis/edges.txt.20191108', thermomorphogenesis_genes)

+for i in range(len(graph_lst)):

+    G = graph_lst[i]

+    print('Graph from %s' % (graph_names[i]))    

+    e = G.edges(data=True)

+    n = len(e)

+    print('Number of edges is %d' % (n))

+    print('------------------------------------------------------------------------')    

+    draw_graph2(G, '../Data/temp/graph-%s.pdf' % (graph_names[i].split('.')[-1]), graph_names[i].split('.')[-1])

 

-for i in range(len(graph_lst)-1):

-    G1 = graph_lst[i]

-    G2 = graph_lst[i+1]

 

+print('Compute network differences ...')

+Gdiff1 = nx.difference(G0, G)

+draw_graph3(Gdiff1, '../Data/temp/graph-in-G0-not-in-G.pdf')

 

-    print(nx.is_isomorphic(G1, G2))

-    print('Graph 1 from %s' % (graph_names[i]))

-    e1 = G1.edges(data=True)

-    n1 = len(e1)

-    tw1 = compute_total_edge_weight(e1, G1)

-    print('Number of edges is %d' % (n1))

-    if n1 > 0:

-        print('Total edge association strength is %4.2f (avg=%4.2f).' % (tw1, tw1/n1))

-    #print(e1)

+Gdiff2 = nx.difference(G, G0)

+draw_graph3(Gdiff2, '../Data/temp/graph-in-G-not-in-G0.pdf')

 

-    print('Graph 2 from %s' % (graph_names[i+1]))    

-    e2 = G2.edges(data=True)

-    n2 = len(e2)

-    tw2 = compute_total_edge_weight(e2, G2)

-    print('Number of edges is %d' % (n2))

-    if n2 > 0:

-        print('Total edge association strength is %4.2f (avg=%4.2f).' % (tw2, tw2/n2))

-    #print(e2)

+print('Compute network intersection ...')

+Gdiff1 = nx.intersection(G0, G)

+draw_graph3(Gdiff1, '../Data/temp/graph-in-G0-and-in-G.pdf')

 

-    ged = nx.algorithms.similarity.graph_edit_distance(G1, G2)

-    print(ged)

+print('Compute edit distance ...')

+ged = nx.algorithms.similarity.graph_edit_distance(G0, G) # sometimes take very long time to finish

+print('Edit distance is %4.0f' % (ged))

 

-    print('------------------------------------------------------------------------')

+# for i in range(len(graph_lst)-1):

+#     G1 = graph_lst[i]

+#     G2 = graph_lst[i+1]

 

-    draw_graph2(G1, '../Data/temp/graph-%s.png' % (graph_names[i].split('.')[-1]), graph_names[i].split('.')[-1])

 

-    

+#     print(nx.is_isomorphic(G1, G2))

+#     print('Graph 1 from %s' % (graph_names[i]))

+#     e1 = G1.edges(data=True)

+#     n1 = len(e1)

+#     tw1 = compute_total_edge_weight(e1, G1)

+#     print('Number of edges is %d' % (n1))

+#     if n1 > 0:

+#         print('Total edge association strength is %4.2f (avg=%4.2f).' % (tw1, tw1/n1))

+#     #print(e1)

+

+#     print('Graph 2 from %s' % (graph_names[i+1]))    

+#     e2 = G2.edges(data=True)

+#     n2 = len(e2)

+#     tw2 = compute_total_edge_weight(e2, G2)

+#     print('Number of edges is %d' % (n2))

+#     if n2 > 0:

+#         print('Total edge association strength is %4.2f (avg=%4.2f).' % (tw2, tw2/n2))

+#     #print(e2)

+

+#     ged = nx.algorithms.similarity.graph_edit_distance(G1, G2)

+#     print(ged)

+

+#     print('------------------------------------------------------------------------')

+

+#     draw_graph2(G1, '../Data/temp/graph-%s.png' % (graph_names[i].split('.')[-1]), graph_names[i].split('.')[-1])