Graph.py

#!/usr/bin/env python
import graph_tool.all as gt
import csv
import re
import sys, getopt
import os 
import argparse
from shutil import copyfile

#arguments
parser = argparse.ArgumentParser()
parser.add_argument('--output','-o',default='./', help="ouput directory [default: ./]")
parser.add_argument('--input','-i',default='./singlenode.csv', help="input directory [default: ./singlenode.csv]")
args = parser.parse_args()

#load matrix
g=gt.load_graph_from_csv(args.input,skip_first=True,string_vals=True,csv_options={"delimiter":","},directed=False)

#g, pos = gt.triangulation(random((500, 2)) * 4, type="delaunay")

tree = gt.min_spanning_tree(g)
gt.graph_draw(g, edge_color=tree, output=args.output+"min_tree.svg")

vp, ep = gt.betweenness(g)
gt.graph_draw(g, vertex_fill_color=vp,output=args.output+"betweenness.pdf",edge_pen_width=gt.prop_to_size(ep, mi=0.5, ma=5),vorder=vp)

#minimize blockmodel with and without overlap
state_sbm_overlap_ndc = gt.minimize_blockmodel_dl(g, overlap=True, deg_corr=False)
state_sbm_overlap_dc  = gt.minimize_blockmodel_dl(g, overlap=True, deg_corr=True)

#equilibration
gt.mcmc_equilibrate(state_sbm_overlap_ndc)
gt.mcmc_equilibrate(state_sbm_overlap_dc)

if state_sbm_overlap_ndc.entropy() < state_sbm_overlap_dc.entropy():
	state_sbm_overlap_ndc.draw(output=args.output+"sbm_overlap_ndc.svg")
	state_sbm_overlap=state_sbm_overlap_ndc
else:
	state_sbm_overlap_dc.draw(output=args.output+"sbm_overlap_dc.svg")
	state_sbm_overlap=state_sbm_overlap_dc
	
#minimize blockmodel
state_sbm_ndc = gt.minimize_blockmodel_dl(g, overlap=True, deg_corr=False)
state_sbm_dc  = gt.minimize_blockmodel_dl(g, overlap=True, deg_corr=True)

#equilibration
gt.mcmc_equilibrate(state_sbm_ndc)
gt.mcmc_equilibrate(state_sbm_dc)

if state_sbm_ndc.entropy() < state_sbm_dc.entropy():
	state_sbm_ndc.draw(output=args.output+"sbm_ndc.svg")
	state_sbm = state_sbm_ndc
else:
	state_sbm_dc.draw(output=args.output+"sbm_dc.svg")
	state_sbm = state_sbm_dc

#nested blockmodel with and without equilibration
#without equlibration
state_nested_ndc = gt.minimize_nested_blockmodel_dl(g, deg_corr=False)
state_nested_dc = gt.minimize_nested_blockmodel_dl(g, deg_corr=True)

if state_nested_ndc.entropy() < state_nested_dc.entropy():
	state_nested_ndc.draw(output=args.output+"nested_ndc.svg")
	state_nested = state_nested_ndc
else:
	state_nested_dc.draw(output=args.output+"nested_dc.svg")
	state_nested = state_nested_dc

#with equilibration
state_nested_equil_ndc = state_nested_ndc.copy(sampling=True)
gt.mcmc_equilibrate(state_nested_equil_ndc)

state_nested_equil_dc = state_nested_dc.copy(sampling=True)
gt.mcmc_equilibrate(state_nested_equil_dc)

if state_nested_ndc.entropy() < state_nested_dc.entropy():
	state_nested_ndc.draw(output=args.output+"nested_equil_ndc.svg")
	state_nested_equil = state_nested_equil_ndc
else:
	state_nested_dc.draw(output=args.output+"nested_equil_dc.svg")
	state_nested_equil = state_nested_equil_dc

#central edges view
bv, be = gt.betweenness(g)
u=gt.GraphView(g,directed=False)
gt.graph_draw(u,vertex_fill_color=bv, output=args.output+"central-edges-view.svg")


#write files

#sbm with overlap
b = state_sbm_overlap.get_majority_blocks()
d = {}
for i in range(0,len(list(b))):
    d[g.vp.name[i]]=b[i]

with open(args.output+'block_member_sbm_overlap.csv', 'w') as csv_file:
    writer = csv.writer(csv_file)
    for gene, group in d.items():
       writer.writerow([gene, group])

#sbm
s = state_sbm.get_blocks()
e = {}
for i in range(0,len(list(g.vp.name))):
    e[g.vp.name[i]]=s[i]

with open(args.output+'block_member_sbm.csv', 'w') as csv_file:
    writer = csv.writer(csv_file)
    for gene, group in e.items():
       writer.writerow([gene, group])

#nested sbm
lstate = state_nested.levels[0]
n=lstate.get_blocks()
f = {}
for i in range(0,len(list(n))):
    f[g.vp.name[i]]=n[i]

with open(args.output+'block_member_nested.csv', 'w') as csv_file:
    writer = csv.writer(csv_file)
    for gene, group in f.items():
       writer.writerow([gene, group])

#nested sbm with equilibration
lstate2 = state_nested_equil.levels[0]
eq=lstate2.get_blocks()
h = {}
for i in range(0,len(list(eq))):
    h[g.vp.name[i]]=eq[i]

with open(args.output+'block_member_nested_equil.csv', 'w') as csv_file:
    writer = csv.writer(csv_file)
    for gene, group in h.items():
       writer.writerow([gene, group])

#writing best model to file 
#nested models tend to be under estimated
states_value = [(state_nested.entropy()-(state_sbm.get_N() * 5)),(state_nested_equil.entropy()-(state_sbm.get_N() * 3)),state_sbm.entropy(),state_sbm_overlap.entropy()]
states = ['nested','nested_equil','sbm','sbm_overlap']
state_best = states[states_value.index(min(states_value))]

copyfile('block_member_'+state_best+".csv","./block_member.csv")
print(state_best+" seems to be the best fitting algorithm and the memberships of it are saved to block_member.csv")

#save graph to file

g.save(args.output+"graph.xml.gz")