# Usage: python draw_subnetwork.py edges.txt # 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 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 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: 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) 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) 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.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) 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 ## 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 'AT5G43700', #IAA4 'AT4G16280', #FCA 'AT2G43060', #IBH1 'AT2G18300', #HBI1 'AT4G28720', #YUC8 'AT1G70560', #TAA1 'AT1G30330', #ARF6 'AT1G19850', #ARF5 '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', 'AT4G28720':'YUC8', 'AT1G70560':'TAA1', 'AT1G30330':'ARF6', 'AT1G19850':'ARF5', 'AT5G01830':'SAUR21', 'AT1G75080':'BZR1', 'AT2G25930':'ELF3', 'AT2G40080':'ELF4', 'AT3G46640':'LUX' } 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 = [] for fname in sorted(glob.glob('../Analysis/edges.txt.2019*')): if fname == '../Analysis/edges.txt.20190801': continue print(fname) graph_names.append(fname) G = build_network_from_file(fname, thermomorphogenesis_genes_small, gene_dict) graph_lst.append(G) 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]) print('Compute network differences ...') Gdiff1 = nx.difference(G0, G) draw_graph3(Gdiff1, '../Data/temp/graph-in-G0-not-in-G.pdf') Gdiff2 = nx.difference(G, G0) draw_graph3(Gdiff2, '../Data/temp/graph-in-G-not-in-G0.pdf') print('Compute network intersection ...') Gdiff1 = nx.intersection(G0, G) draw_graph3(Gdiff1, '../Data/temp/graph-in-G0-and-in-G.pdf') 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)) # 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])