From 9f70ca557368dcdc2e890b80f34459e9f6b29dec Mon Sep 17 00:00:00 2001 From: Hui Lan Date: Thu, 5 Dec 2019 22:36:15 +0800 Subject: 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. --- Code/draw_subnetwork.py | 248 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 248 insertions(+) create mode 100644 Code/draw_subnetwork.py (limited to 'Code') 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]) + + -- cgit v1.2.1