<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Homework 1
Step1:  How many distinct metabolites are there in the graph?
Step2:  How many reactions?
Step3:  How many edges are there?
Step4:  Calculate the degree k=kin+kout of each of the metabolite vertices (you don’t need to calculate the degree for the vertices that correspond to reactions).
Step5:  What are the top six metabolites in terms of vertex degree in the graph?
Step6:  Plot the distribution of the degrees of these vertices, on log-log scale.
Step7:  Is the degree distribution well-described by a power law?
Yes, the log-log histogram plot decreases linearly.
Calculate the exponent $\alpha$ of the best-fit power-law p_{k} \propto k^{-\alpha}pk∝k−α.
Step8:  How does the \alphaα that you get compare to the estimate of the power-law exponent reported by Jeong et al. in their 2000 article in Nature, “The large-scale organization of metabolic networks” (vol. 407, pp. 651–654)?
The alpha value 2.02 is similar to the gamma value 2.2 in that paper.
Based on structure of the network that you analyzed (bipartite, containing reactions) vs. the structure of the network that they analyzed (network projected to a network containing only metabolites), is it appropriate to compare the exponents? Why or why not?
Comparing the exponents may give a general intuition about the structural similarity between each network, but a direct comparison may be less meaningful.
Calculate the shortest-path-length between all pairs of metabolites (vertices) in the directed graph.
Step9:  Why are some of these path lengths infinite?
In this code, length -1 means there is no path between the vertices.
Calculate the average of the shortest-path-length between all pairs of metabolites in the giant (weakly connected) component of the network, by throwing away any shortest-paths distance value if it is infinite or zero.
Step10:  Calculate the maximum of the shortest-path-length between all pairs of metabolites (throwing away infinite values, as before) in the giant (weakly connected) component of the network (i.e., you are calculating the diameter of the giant component, according to Newman's definition of diameter).
Step11:  Why are the average geodesic distances that we get, roughly twice those reported in Fig. 3b of Jeong et al., 2000?
Our graph is bipartite so there will always be twice as many steps between metabolites than in Jeong.
Calculate the shortest-paths betweenness centrality for all metabolites in the directed network.
Step12:  Plot the scatter plot of betweenness centrality vs. vertex degree for all metabolites, on log-log scale.
Step13:  Among metabolites with degree k=2, what metabolite has highest betweenness centrality in the network?
<ASSISTANT_TASK:>
Python Code:
import pandas
import igraph
edge_list = pandas.read_csv('hsmetnet.txt', sep='\t', header=None).values.tolist()
graph = igraph.Graph.TupleList(edge_list, directed=True)
igraph.summary(graph)
Explanation: Homework 1
End of explanation
flatten = lambda l: sum(l, [])
vertex_set = set(flatten(edge_list))
isreaction = lambda vertex: vertex[0:8]=='REACTION'
metabolite_list = [vertex for vertex in vertex_set if not isreaction(vertex)]
len(metabolite_list)
Explanation: How many distinct metabolites are there in the graph?
End of explanation
reaction_list =  [vertex for vertex in vertex_set if isreaction(vertex)]
len(reaction_list)
Explanation: How many reactions?
End of explanation
len(edge_list)
Explanation: How many edges are there?
End of explanation
nameof = lambda vertex: vertex.attributes()['name']
metabolite_degree = {}
for vertex in graph.vs:
    if not isreaction(nameof(vertex)):
       metabolite_degree[nameof(vertex)] = vertex.degree()
Explanation: Calculate the degree k=kin+kout of each of the metabolite vertices (you don’t need to calculate the degree for the vertices that correspond to reactions).
End of explanation
metabolite_degree_sorted = sorted(metabolite_degree.items(), key=lambda item: item[1], reverse=True)
metabolite_degree_sorted[0:6]
Explanation: What are the top six metabolites in terms of vertex degree in the graph?
End of explanation
import matplotlib.pyplot as plt
import numpy as np
metabolites = [item[0] for item in metabolite_degree.items()]
degrees = [item[1] for item in metabolite_degree.items()]
plt.hist(degrees, bins=np.logspace(np.log10(1),np.log10(1202),50), log=True)
plt.title('Degree Distribution')
plt.xlabel('k')
plt.ylabel('N(k)')
plt.gca().set_xscale('log')
plt.show()
Explanation: Plot the distribution of the degrees of these vertices, on log-log scale.
End of explanation
igraph.statistics.power_law_fit(graph.degree())
Explanation: Is the degree distribution well-described by a power law?
Yes, the log-log histogram plot decreases linearly.
Calculate the exponent $\alpha$ of the best-fit power-law p_{k} \propto k^{-\alpha}pk∝k−α.
End of explanation
vertices = []
for vertex in graph.vs:
    if not isreaction(nameof(vertex)):
        vertices += [vertex]
shortest_paths = [[len(x)-1 for x in graph.get_shortest_paths(v=vertex,to=vertices)] for vertex in vertices]
Explanation: How does the \alphaα that you get compare to the estimate of the power-law exponent reported by Jeong et al. in their 2000 article in Nature, “The large-scale organization of metabolic networks” (vol. 407, pp. 651–654)?
The alpha value 2.02 is similar to the gamma value 2.2 in that paper.
Based on structure of the network that you analyzed (bipartite, containing reactions) vs. the structure of the network that they analyzed (network projected to a network containing only metabolites), is it appropriate to compare the exponents? Why or why not?
Comparing the exponents may give a general intuition about the structural similarity between each network, but a direct comparison may be less meaningful.
Calculate the shortest-path-length between all pairs of metabolites (vertices) in the directed graph.
End of explanation
giant = graph.components(mode=igraph.WEAK).giant()
vertices = []
for vertex in giant.vs:
    if not isreaction(nameof(vertex)):
        vertices += [vertex]
giant_shortest_paths = [[len(x)-1 for x in graph.get_shortest_paths(v=vertex,to=vertices)] for vertex in vertices]
n = 0
m = 0
for vertex_shortest_paths in giant_shortest_paths:
    for path_length in vertex_shortest_paths:
        if path_length > 0:
            n += path_length
            m += 1
n/m
Explanation: Why are some of these path lengths infinite?
In this code, length -1 means there is no path between the vertices.
Calculate the average of the shortest-path-length between all pairs of metabolites in the giant (weakly connected) component of the network, by throwing away any shortest-paths distance value if it is infinite or zero.
End of explanation
diameter = 0
for vertex_shortest_paths in giant_shortest_paths:
    for path_length in vertex_shortest_paths:
        if path_length > 0 and path_length > diameter:
            diameter = path_length
diameter
Explanation: Calculate the maximum of the shortest-path-length between all pairs of metabolites (throwing away infinite values, as before) in the giant (weakly connected) component of the network (i.e., you are calculating the diameter of the giant component, according to Newman's definition of diameter).
End of explanation
vertices = []
for vertex in graph.vs:
    if not isreaction(nameof(vertex)):
        vertices += [vertex]
betweennesses = graph.betweenness(vertices=vertices)
Explanation: Why are the average geodesic distances that we get, roughly twice those reported in Fig. 3b of Jeong et al., 2000?
Our graph is bipartite so there will always be twice as many steps between metabolites than in Jeong.
Calculate the shortest-paths betweenness centrality for all metabolites in the directed network.
End of explanation
degrees = [vertex.degree() for vertex in vertices]
plt.clf()
plt.scatter(degrees, betweennesses)
plt.title('Metabolite Betweenness vs. Degree')
plt.xlabel('Degree')
plt.ylabel('Betweenness')
plt.gca().set_xscale('log')
#plt.gca().set_yscale('log')
plt.show()
Explanation: Plot the scatter plot of betweenness centrality vs. vertex degree for all metabolites, on log-log scale.
End of explanation
highest_betweenness = 0
highest_vertex = None
for degree, betweenness, vertex in zip(degrees, betweennesses, vertices):
    if degree == 2 and betweenness > highest_betweenness:
        highest_betweenness = betweenness
        highest_vertex = vertex
highest_vertex
Explanation: Among metabolites with degree k=2, what metabolite has highest betweenness centrality in the network?
End of explanation
<END_TASK>