# ---
# jupyter:
#   jupytext:
#     formats: ipynb,py
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.9.1+dev
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# ## Preprocess and validate eADAGE network data
#
# eADAGE G-G network/edge list was generated by Alex using the script at [make_GiG_network.R](../generic_expression_patterns_modules/make_GiG_network.R).
#
# Here, we preprocess it to:
# * reweight negative edges
# * add generic gene info from [here](../pseudomonas_analysis/annot_df.tsv)
#
# Output: A new edge list file that can be quickly loaded by `graph-tool` for subsequent analyses.

# +
import os

import numpy as np
import pandas as pd
import graph_tool.all as gt
import matplotlib.pyplot as plt
import seaborn as sns

# +
# relevant file paths
data_dir = "./data"
unprocessed_edge_list = os.path.join(data_dir, "edgeList.csv")

# map of Pa gene names to generic/not generic status, generated by Alex
generic_gene_map = os.path.join("..", "pseudomonas_analysis", "annot_df.tsv")

# save edge list with preprocessed weight information
processed_edge_list = os.path.join(data_dir, "edge_list_processed_unsigned.csv")

# place to save preprocessed graph/attributes, in graph-tool binary format
processed_graph = os.path.join(data_dir, "eadage_generic_graph_unsigned.gt")
# -

# #### Load edge list and handle negative edge weights
#
# Most algorithms for community detection and betweenness don't work with negative weights/correlations. Taking inspiration from WGCNA, we can solve this problem in one of two ways:
# * Unsigned: `weight = abs(corr(g1, g2))`
# * Signed: `weight = abs((1 + corr(g1, g2) / 2)`
#
# We should probably try both eventually, but for now we'll just use the unsigned version (taking absolute value of negative edges). I tend to agree with the rationale in this BioStars post for taking absolute values (unsigned) rather than rescaling (signed approach): https://www.biostars.org/p/144078/#144088

# +
if not os.path.isfile(processed_edge_list):
    edgelist_df = pd.read_csv(unprocessed_edge_list, index_col=0)
    # take absolute value of edge weights
    edgelist_df["weight"] = edgelist_df.weight.abs()
    edgelist_df.to_csv(
        processed_edge_list, columns=["from", "to", "weight"], index=False
    )

edgelist_df = pd.read_csv(processed_edge_list)
edgelist_df.head()
# -

# In this "generic gene map", 1 denotes a generic gene and 0 is all other genes. A gene is considered generic if it had a high percentile from SOPHIE and the manually curated set based on the correlation plot seen [here](../pseudomonas_analysis/2_identify_generic_genes_pathways.ipynb).

annot_df = pd.read_csv(generic_gene_map, sep="\t", index_col=0)
annot_df.head()

G = gt.load_graph_from_csv(
    processed_edge_list,
    skip_first=True,
    directed=False,
    hashed=True,
    eprop_names=["weight"],
    eprop_types=["float"],
)

# +
# add vertex property for generic genes
vprop_generic = G.new_vertex_property("bool")
for ix, v in enumerate(G.vertices()):
    v_name = G.vp["name"][v]
    v_label = annot_df.loc[v_name, "label"]
    vprop_generic[v] = v_label

G.vertex_properties["is_generic"] = vprop_generic
# -

# make sure vertex/edge properties exist
print(G)
print(list(G.vp.keys()))
print(list(G.ep.keys()))

# make sure names/weights from file were loaded properly
for s, t, w in G.iter_edges([G.ep["weight"]]):
    print(G.vp["name"][s], G.vp["name"][t], w)
    if s > 0:
        break

# save graph with attributes to file
G.save(processed_graph, fmt="gt")

# plot generic and non-generic genes, just for fun
# https://stackoverflow.com/a/60462353
red_blue_map = {0: (1, 0, 0, 1), 1: (0, 0, 1, 1)}
plot_color = G.new_vertex_property("vector<double>")
for v in G.vertices():
    plot_color[v] = red_blue_map[G.vp["is_generic"][v]]
gt.graph_draw(G, vertex_fill_color=plot_color)
