#!/usr/bin/env catnip
# Réseau d'interactions de résidus (RIN) et centralités de la crambine (PDB 1CRN)
#
# La crambine est une petite protéine hydrophobe de 46 résidus (chaîne A). On la
# décrit ici comme un GRAPHE de contacts, pas comme une séquence :
#
#  - NŒUD : un résidu, étiqueté par son nom et son numéro (ex. THR1, CYS3).
#  - ARÊTE : une paire de résidus « en contact » spatial. Le contact se mesure
#    entre carbones Cα (carbone central du squelette peptidique) : deux résidus
#    sont reliés si leur distance Cα-Cα est < SEUIL et s'ils sont séparés d'au
#    moins SEQ_SEP le long de la chaîne. Le filtre |i - j| >= SEQ_SEP écarte les
#    voisins de séquence triviaux (liaison peptidique), pour ne garder que les
#    rapprochements dus au repliement tridimensionnel.
#
# Sur ce graphe on calcule trois centralités (Freeman 1978, Wikipedia
# « Centrality ») : degré (nombre de contacts, normalisé), intermédiarité
# (betweenness : fraction de plus courts chemins passant par le nœud) et
# proximité (closeness : inverse de la distance moyenne aux autres nœuds). Elles
# repèrent des HUBS DE CONTACT STRUCTURAUX — une propriété topologique du graphe,
# pas une importance biologique démontrée (voir la lecture heuristique en fin).
#
# gemmi lit le mmCIF ; on extrait tôt name/seqid/x/y/z dans des structs Catnip.
# numpy/scipy calculent la matrice de distances ; networkx porte le graphe et les
# centralités ; Catnip organise la sélection des arêtes, le rapport et la figure.
#
# Provenance des données :
#   PDB 1CRN — crambine, Crambe abyssinica, résolution 0.54 Å (chaîne A, 46 résidus)
#   Citation primaire : Teeter M.M. (1984) PNAS 81:6014-6018
#   Source : RCSB PDB / wwPDB — format PDBx/mmCIF — licence CC0 (domaine public)
#
# Portée de l'exemple : seuls les Cα de la chaîne A sont analysés. Sont exclus les
# eaux, hétéro-atomes et ligands (absents ici : la chaîne A est 100 % acides
# aminés). find_atom('CA', '*') retourne le conformère primaire ; ce dépôt ne
# contient aucun altloc secondaire (colonne altloc entièrement vide). Aucune
# symétrie cristallographique n'est appliquée : distances sur l'unité asymétrique.
#
# DEPS: gemmi numpy scipy networkx matplotlib

gemmi = import('gemmi')
numpy = import('numpy')
sp_dist = import('scipy.spatial.distance')
nx = import('networkx')
mpl = import('matplotlib')
mpl.use('Agg')  # backend headless : rendu fichier, aucune fenêtre
plt = import('matplotlib.pyplot')
import('pathlib', 'Path')
import('builtins', 'sorted')

script_dir = Path(META.file).parent
output_dir = script_dir / 'output'
output_dir.mkdir(exist_ok=True)

# Paramètres du graphe (documentés, unités explicites)
SEQ_SEP = 2      # séparation minimale en séquence : |i - j| >= 2 (exclut i, i+1)
THRESHOLD = 8.0  # seuil de contact Cα-Cα en Ångström (convention type CASP)

PDB_ID = '1CRN'
CITATION = "Teeter M.M. (1984) PNAS 81:6014-6018"
PROVENANCE = "RCSB PDB / wwPDB, PDBx/mmCIF, licence CC0"

# Un résidu réduit à son identité et à la position de son Cα. On ne broadcaste
# jamais sur les objets gemmi bruts : tout est copié ici dans des scalaires.
struct Residue {
    name: str; seqid: int; x: float; y: float; z: float;

    label(self): str => { f"{self.name}{self.seqid}" }
}

# Une arête du RIN : deux index de résidus en contact et leur distance Cα-Cα (Å).
struct Edge {
    i: int; j: int; dist: float;
}

# Centralités d'un nœud (toutes normalisées dans [0, 1]).
struct Centrality {
    label: str; degree: float; betweenness: float; closeness: float;
}

# Rapport des grandeurs CALCULÉES (distinctes des faits extraits et de la lecture
# heuristique). Le contrôle numérique est le lemme des poignées de main
# (handshake lemma) : la somme des degrés vaut exactement deux fois le nombre
# d'arêtes, puisque chaque arête est comptée à ses deux extrémités.
struct NetworkReport {
    n_nodes: int; n_edges: int; avg_degree: float;
    deg_sum: int; handshake_ok: bool; seq_sep: int; threshold: float;

    display(self) => {
        print("⇒ Mesures calculées (topologie du graphe, pas d'interprétation biologique)")
        print(f"  Nœuds (résidus)   : {self.n_nodes}")
        print(f"  Arêtes (contacts) : {self.n_edges}")
        print(f"    définition : |i - j| >= {self.seq_sep} et distance Cα-Cα < {self.threshold} Å")
        print(f"  Degré moyen : {round(self.avg_degree, 3)} contacts par résidu")
        print(f"  Contrôle (lemme des poignées de main) : Σdeg = {self.deg_sum}, 2·|E| = {2 * self.n_edges} → égalité {self.handshake_ok}")
    }
}

# --- Chargement et extraction (FAITS lus dans le fichier) ---

cif_path = script_dir / 'data/1crn.cif'
st = gemmi.read_structure(str(cif_path))
chain = st[0][0]  # modèle 0, chaîne A
n = len(chain)
idx = numpy.arange(n).tolist()

# Broadcast sur des index entiers (pas sur les résidus gemmi) : chaque itération
# lit son Cα et copie name/seqid/x/y/z dans un struct Catnip.
residues = idx.[(i) => {
    r = chain[int(i)]
    ca = r.find_atom('CA', '*')  # conformère primaire (aucun altloc ici)
    Residue(r.name, r.seqid.num, float(ca.pos.x), float(ca.pos.y), float(ca.pos.z))
}]

labels = residues.[(r) => { r.label() }]

# --- Matrice de distances et arêtes (MESURES) ---

coords = numpy.array(residues.[(r) => { list(r.x, r.y, r.z) }])
dmat = sp_dist.cdist(coords, coords)  # (n, n) distances Cα-Cα en Å

# triu_indices(n, SEQ_SEP) : couples (i, j) avec j - i >= SEQ_SEP, donc i < j et
# séparation de séquence suffisante. Indexation à plat car Catnip ne parse pas
# l'indexation numpy à deux axes dmat[ii, jj].
tri = numpy.triu_indices(n, SEQ_SEP)
ii = tri[0]
jj = tri[1]
dv = dmat.ravel()[ii * n + jj]
sel = numpy.less(dv, THRESHOLD)
sel_i = ii[sel].tolist()
sel_j = jj[sel].tolist()
sel_d = dv[sel].tolist()

edges = numpy.arange(len(sel_i)).tolist().[(k) => {
    Edge(int(sel_i[int(k)]), int(sel_j[int(k)]), float(sel_d[int(k)]))
}]

# --- Construction du graphe networkx ---

g = nx.Graph()
for lab in labels {
    g.add_node(lab)
}
for e in edges {
    g.add_edge(labels[e.i], labels[e.j], weight=e.dist)
}

# --- Centralités (MESURES) ---

deg = nx.degree_centrality(g)
bet = nx.betweenness_centrality(g)
clo = nx.closeness_centrality(g)

cents = labels.[(l) => { Centrality(l, deg[l], bet[l], clo[l]) }]
top_deg = sorted(cents, key=(c) => { -c.degree })[:5]
top_bet = sorted(cents, key=(c) => { -c.betweenness })[:5]
top_clo = sorted(cents, key=(c) => { -c.closeness })[:5]

# Contrôle numérique : lemme des poignées de main. La somme des degrés (comptée
# nœud par nœud sur le graphe networkx) doit égaler 2·|E|.
degs = numpy.array(labels.[(l) => { g.degree(l) }])
deg_sum = int(numpy.sum(degs))
n_edges = g.number_of_edges()
handshake_ok = deg_sum == 2 * n_edges
avg_degree = float(2 * n_edges) / n

# --- Rapport à trois niveaux ---

print(f"⇒ Faits extraits du fichier ({PDB_ID})")
print(f"  Provenance : {PROVENANCE}")
print(f"  Citation primaire : {CITATION}")
print(f"  Structure : {st.name}, chaîne {chain.name}, {n} résidus (Cα)")
print()

report = NetworkReport(n, n_edges, avg_degree, deg_sum, handshake_ok, SEQ_SEP, THRESHOLD)
report.display()
print()

print("⇒ Top 5 hubs de contact par centralité (topologique, normalisée)")
print("  Degré (nombre de contacts) :")
for c in top_deg {
    print(f"    {c.label} : {round(c.degree, 3)}")
}
print("  Intermédiarité (betweenness) :")
for c in top_bet {
    print(f"    {c.label} : {round(c.betweenness, 3)}")
}
print("  Proximité (closeness) :")
for c in top_clo {
    print(f"    {c.label} : {round(c.closeness, 3)}")
}
print()

print("⇒ Lecture heuristique (À NE PAS confondre avec une preuve biologique)")
print("  La centralité de graphe repère des hubs de contact STRUCTURAUX : des")
print("  résidus au carrefour de nombreux contacts non-locaux. C'est une propriété")
print("  TOPOLOGIQUE du réseau de distances, pas une importance fonctionnelle ou")
print("  catalytique (qui exige des données biochimiques). Repère indicatif : la")
print("  crambine porte 3 ponts disulfure (CYS 3-40, 4-32, 16-26) qui rigidifient")
print("  le cœur ; voir des cystéines parmi les hubs est COHÉRENT avec cette")
print("  rigidité, mais la centralité ne le DÉMONTRE pas.")
print()

# --- Figure : réseau dessiné, nœuds colorés/taillés par centralité de degré ---

node_deg = numpy.array(labels.[(l) => { deg[l] }])
node_size = labels.[(l) => { 320.0 * deg[l] + 40.0 }]

pos = nx.spring_layout(g, seed=42)  # layout déterministe (graine fixée)

fig, ax = plt.subplots(figsize=tuple(9, 8))
nx.draw_networkx_edges(g, pos, ax=ax, edge_color='#cccccc', width=0.6)
nodes = nx.draw_networkx_nodes(
    g,
    pos,
    ax=ax,
    nodelist=labels,
    node_color=node_deg,
    node_size=node_size,
    cmap='plasma',
)
fig.colorbar(nodes, ax=ax, label="centralité de degré", fraction=0.046, pad=0.04)
nx.draw_networkx_labels(g, pos, ax=ax, font_size=5)
ax.set_title(f"Réseau d'interactions de résidus ({PDB_ID}) — seuil {THRESHOLD} Å, |i-j| >= {SEQ_SEP}")
ax.axis('off')

fig.tight_layout()
output_path = output_dir / 'residue_interaction_network.png'
fig.savefig(str(output_path), dpi=120, bbox_inches='tight')

print(f"⇒ Figure → {output_path}")

# Aperçu navigateur : sert la figure une fois puis rend la main.
# --no-browser garde un chemin headless (le PNG reste écrit ci-dessus).
if '--no-browser' not in import('sys').argv {
    http = import('http')
    b64 = import('base64').b64encode(output_path.read_bytes()).decode('ascii')
    http.serve(f'<!doctype html><meta charset="utf-8"><title>Residue interaction network</title><body style="margin:0;background:#0d1117"><img style="max-width:100%;display:block;margin:0 auto" src="data:image/png;base64,{b64}">',
        0, 'text/html; charset=utf-8', True)
}