draw_subnetwork.py for generating a subnetwork consisting of genes that might be responsible for thermomorphogenesis

Use networkx and matplotlib.

Reference:

Quint et al. (2016) Molecular and genetic control of plant thermomorphogenesis.  Nature Plants.

1 files changed, 248 insertions, 0 deletions

diff --git a/Code/draw_subnetwork.py b/Code/draw_subnetwork.py
new file mode 100644
index 0000000..79ab5cd
--- /dev/null
+++ b/Code/draw_subnetwork.py
@@ -0,0 +1,248 @@
+# Usage: python draw_subnetwork.py edges.txt

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

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

+

+import os, sys

+import networkx as nx

+import pylab as plt

+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()    

+

+    for g in gene_lst:

+        G.add_node(gene_dict[g])

+        

+    f = open(edge_fname)

+    for line in f:

+        line = line.strip()

+        lst = line.split('\t')

+        if len(lst) == 10:

+

+            g1 = lst[0].split()[0] # target gene ID

+            g2 = lst[1].split()[0] # source gene ID        

+

+                

+            strength = float(lst[8])

+            method_or_tissue = lst[9]            

+

+

+            g1_label = lst[0].split()[1].split(';')[0] if lst[0].split()[1] != '.' else g1

+            g1_name = lst[0].split()[1] if lst[0].split()[1] != '.' else ''

+                

+

+            g2_label = lst[1].split()[1].split(';')[0] if lst[1].split()[1] != '.' else g2

+            g2_name = lst[1].split()[1] if lst[1].split()[1] != '.' else ''

+

+            if g1 in gene_lst and g2 in gene_lst:

+                G.add_node(gene_dict[g1], full_name=g1_name, label=g1_label) # if g1 is also a TF, then istf='0' will overwrite it in the following for loop

+                G.add_node(gene_dict[g2], full_name=g2_name, label=g2_label) # tf_category contains default TF category code.  It can be modified later given user's input

+                G.add_edge(gene_dict[g2], gene_dict[g1], weight=strength, strength=strength, method=method_or_tissue) # g2 is source, and g1 is target

+

+    f.close()

+

+    return G

+

+

+def compute_total_edge_weight(edges, G):

+    total = 0

+    for e in edges:

+        u = e[0]

+        v = e[1]

+        total += G[u][v]['weight']

+    return total

+

+

+def draw_graph(G, fname):

+    pos=nx.circular_layout(G)

+    tau = 2.5

+    elarge=[(u,v) for (u,v,d) in G.edges(data=True) if d['weight'] >tau]

+    esmall=[(u,v) for (u,v,d) in G.edges(data=True) if d['weight'] <=tau]

+    labels = {}

+    for (n,d) in G.nodes(data=True):

+        if 'label' in d:

+            labels[n] = d['label']

+        else:

+            labels[n] = n

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

+    nx.draw_networkx_edges(G,pos,edgelist=elarge,width=1,alpha=0.2)

+    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.show() # display    

+

+

+def better_date(s):

+    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)

+    all_edges = []

+    all_widths = []

+

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

+        all_edges.append((u, v))

+        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_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.show() # display    

+    

+

+## main

+

+thermomorphogenesis_genes = [

+    'AT4G28720',

+    'AT2G25930',

+    'AT2G40080',

+    'AT3G46640',

+    'AT5G11260',

+    'AT2G43010',

+    'AT3G59060',

+    'AT4G10180',

+    'AT2G32950',

+    'AT3G13550',

+    'AT4G05420',

+    'AT4G21100',

+    'AT2G46340',

+    'AT4G11110',

+    'AT3G15354',

+    'AT1G53090',

+    'AT1G02340',

+    'AT4G08920',

+    'AT4G39950',

+    'AT2G22330',

+    'AT2G42870',

+    'AT5G39860',

+    'AT1G70560',

+    'AT3G62980',

+    'AT4G03190',

+    'AT3G26810',

+    'AT1G12820',

+    'AT4G24390',

+    'AT5G49980',

+    'AT5G01830',

+    'AT5G18010',

+    'AT5G18020',

+    'AT5G18050',

+    'AT5G18060',

+    'AT5G18080',

+    'AT1G29440',

+    'AT1G29510',

+    'AT4G18710',

+    'AT1G75080',

+    'AT1G30330',

+    'AT1G19850',

+    'AT3G33520',

+    'AT4G16280',

+    'AT2G43060',

+    'AT2G18300',

+    'AT4G16780',

+    'AT1G01060',

+    'AT1G22770',

+    'AT4G25420',

+    'AT1G15550',

+    'AT1G78440',

+    'AT5G43700',

+    'AT4G32280',

+    'AT2G38120',

+    'AT1G15580',

+]

+

+

+thermomorphogenesis_genes_small = [

+    'AT2G43010',	#PIF4

+    'AT5G11260',	#HY5

+    'AT2G42870',	#PAR1

+    'AT5G39860',	#PRE1

+    'AT4G16280',	#FCA

+    'AT2G43060',	#IBH1

+    'AT2G18300',	#HBI1

+    'AT4G28720',	#YUC8

+    'AT1G70560',	#TAA1

+    'AT1G30330',	#ARF6

+    'AT1G19850',        #ARF5

+    'AT1G75080',	#BZR1

+    'AT2G25930',	#ELF3

+    'AT2G40080',	#ELF4

+    'AT3G46640',	#LUX

+]

+

+

+gene_dict = {

+    'AT2G43010':'PIF4',

+    'AT5G11260':'HY5',

+    'AT2G42870':'PAR1',

+    'AT5G39860':'PRE1',

+    'AT4G16280':'FCA',

+    'AT2G43060':'IBH1',

+    'AT2G18300':'HBI1',

+    'AT4G28720':'YUC8',

+    'AT1G70560':'TAA1',

+    'AT1G30330':'ARF6',

+    'AT1G19850':'ARF5',

+    'AT1G75080':'BZR1',

+    'AT2G25930':'ELF3',

+    'AT2G40080':'ELF4',

+    'AT3G46640':'LUX'

+}

+

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

+

+graph_lst = []

+graph_names = []

+for fname in sorted(glob.glob('../Analysis/edges.txt.2019*')):

+    if fname == '../Analysis/edges.txt.20190801':

+        continue

+    print(fname)

+    graph_names.append(fname)

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

+

+

+#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)-1):

+    G1 = graph_lst[i]

+    G2 = graph_lst[i+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])

+

+