#!/usr/bin/env catnip
# Site de liaison protéine-ligand : poche de l'indinavir dans la protéase du VIH-1 (PDB 1HSG)
#
# 1HSG est le complexe de la protéase du VIH-1 avec l'indinavir (Crixivan), un
# inhibiteur de protéase utilisé en clinique. On cherche la POCHE de liaison :
# l'ensemble des résidus protéiques dont au moins un atome approche un atome du
# ligand sous un seuil documenté.
#
# - Détection de poche : pour chaque résidu protéique, on prend la distance
# minimale entre l'un de ses atomes et l'un des 45 atomes du ligand. Un résidu
# est « dans la poche » si cette distance est <= SEUIL. C'est un critère
# géométrique de proximité, pas une preuve d'interaction chimique.
#
# - Classification heuristique du contact (par distance, pas par chimie) :
# < 3.5 Å suggère une liaison H possible (distance donneur-accepteur typique),
# < 4.5 Å reste dans la gamme des contacts de van der Waals. Ce sont des
# étiquettes indicatives ; établir une vraie liaison H exige les angles, les
# hydrogènes et le type d'atome, pas seulement une distance.
#
# gemmi lit le mmCIF ; on extrait tôt name/seqid/chaîne/coordonnées dans des
# structs et arrays Catnip. scipy calcule les distances min par cdist ; gemmi
# fournit un second calcul indépendant (pos.dist) pour le contrôle numérique.
#
# Provenance des données :
# PDB 1HSG — protéase du VIH-1 (homodimère, chaînes A/B, 99 résidus chacune)
# complexée à l'indinavir (ligand MK1, 45 atomes, chaîne B seqid 902)
# Résolution 2.0 Å — groupe d'espace P 21 21 2
# DOI de la structure : 10.2210/pdb1HSG/pdb
# Source : RCSB PDB / wwPDB — format PDBx/mmCIF — licence CC0 (domaine public)
#
# Faits de la littérature (présentés comme tels, non dérivés ici) : la protéase du
# VIH-1 est un homodimère dont la dyade catalytique Asp25/Asp25' porte la coupure
# peptidique ; l'indinavir est un inhibiteur cliniquement utilisé. Si Asp25 (chaîne
# A ou B) ressort dans la poche, c'est cohérent avec ce site actif connu — une
# cohérence, pas une preuve tirée des distances.
#
# Portée / exclusions : seuls les résidus protéiques (ATOM) sont analysés. Sont
# exclus le ligand lui-même (MK1) et les 127 eaux cristallographiques (HOH) ;
# aucun autre HETATM non-eau n'est présent. Aucune symétrie cristallographique
# n'est appliquée : les distances portent sur les coordonnées de la maille
# asymétrique déposée (les images de symétrie proches ne sont pas comptées).
#
# DEPS: gemmi numpy scipy matplotlib
gemmi = import('gemmi')
numpy = import('numpy')
sp_dist = import('scipy.spatial.distance')
mpl = import('matplotlib')
mpl.use('Agg') # backend headless : rendu fichier, aucune fenêtre
plt = import('matplotlib.pyplot')
import('pathlib', 'Path')
script_dir = Path(META.file).parent
output_dir = script_dir / 'output'
output_dir.mkdir(exist_ok=True)
# Paramètres géométriques (documentés, unités explicites)
THRESHOLD = 4.5 # seuil de poche : distance min atome-ligand <= 4.5 Å (van der Waals)
HBOND_CUT = 3.5 # sous ce seuil, une liaison H est envisageable (heuristique)
PDB_ID = '1HSG'
LIGAND_LABEL = "indinavir, inhibiteur de protéase"
PROVENANCE = "RCSB PDB / wwPDB, PDBx/mmCIF, licence CC0"
DOI = "10.2210/pdb1HSG/pdb"
# Le ligand réduit à son identité et son nombre d'atomes (FAITS).
struct Ligand {
name: str; chain: str; seqid: int; n_atoms: int;
}
# Gamme de contact déduite de la seule distance minimale. Étiquette HEURISTIQUE :
# une distance ne prouve pas une interaction, elle la rend plausible ou non.
union ContactClass {
hbond_range; vdw_range
label(self): str => {
match self {
ContactClass.hbond_range => { "liaison H possible (< 3.5 Å, heuristique)" }
ContactClass.vdw_range => { "van der Waals (< 4.5 Å, heuristique)" }
}
}
}
class_code = (c: ContactClass): int => {
match c {
ContactClass.hbond_range => { 0 }
ContactClass.vdw_range => { 1 }
}
}
class_color = (c: ContactClass): str => {
match c {
ContactClass.hbond_range => { '#d62728' }
ContactClass.vdw_range => { '#1f77b4' }
}
}
classify_contact = (d: float): ContactClass => {
if d < HBOND_CUT { ContactClass.hbond_range } else { ContactClass.vdw_range }
}
# Un résidu scanné : distance min au ligand par deux voies (cd = scipy cdist sur
# coordonnées extraites, gd = gemmi pos.dist indépendant) + drapeau protéique.
# gd n'est calculé que pour les résidus protéiques (le contrôle ne porte que sur
# eux) : la branche else des eaux/ligand n'évalue pas la boucle imbriquée.
struct ResidueScan {
chain: str; name: str; seqid: int; cd: float; gd: float; is_protein: int;
}
# Un résidu de poche : identité + distance min + classe de contact heuristique.
struct PocketResidue {
chain: str; name: str; seqid: int; min_dist: float; contact;
}
# Rapport des grandeurs CALCULÉES (distinctes des faits extraits et des
# classifications heuristiques imprimés autour).
struct NeighborhoodReport {
n_pocket: int; n_protein: int; threshold: float;
closest: str; mean_dist: float; max_ctrl_err: float;
display(self) => {
print("⇒ Mesures calculées (géométrie, pas d'interprétation biologique)")
print(f" Résidus de poche : {self.n_pocket} sur {self.n_protein} résidus protéiques")
print(f" critère : distance min atome-résidu ↔ atome-ligand <= {self.threshold} Å")
print(f" Résidu le plus proche : {self.closest}")
print(f" Distance min moyenne sur la poche : {round(self.mean_dist, 3)} Å")
print(f" Contrôle numérique (scipy cdist vs gemmi pos.dist) : erreur max = {self.max_ctrl_err} Å")
}
}
# --- Chargement et repérage du ligand (FAITS lus dans le fichier) ---
cif_path = script_dir / 'data/1hsg.cif'
st = gemmi.read_structure(str(cif_path))
model = st[0]
# Le ligand = l'unique résidu HETATM non-eau. Balayage while (pas de for fiable).
lig_res = None
ci = 0
while ci < len(model) {
scan_ch = model[ci]
ri = 0
while ri < len(scan_ch) {
scan_r = scan_ch[ri]
if scan_r.het_flag == 'H' and scan_r.name != 'HOH' {
lig_res = scan_r
}
ri = ri + 1
}
ci = ci + 1
}
ligand = Ligand(lig_res.name, model[len(model) - 1].name, lig_res.seqid.num, len(lig_res))
# Coordonnées du ligand extraites tôt dans un array numpy (n_atoms, 3).
lig_coords = numpy.arange(len(lig_res)).tolist().[(i) => {
at = lig_res[int(i)]
list(float(at.pos.x), float(at.pos.y), float(at.pos.z))
}]
lig_arr = numpy.array(lig_coords)
# --- Scan de tous les résidus (MESURES) ---
# Pour un résidu : distance min au ligand par deux implémentations indépendantes.
scan_residue = (chain_name: str, r): ResidueScan => {
n_at = len(r)
coords = numpy.arange(n_at).tolist().[(ai) => {
at = r[int(ai)]
list(float(at.pos.x), float(at.pos.y), float(at.pos.z))
}]
res_arr = numpy.array(coords)
cd = float(numpy.min(sp_dist.cdist(res_arr, lig_arr)))
is_prot = if r.het_flag != 'H' { 1 } else { 0 }
# Contrôle gemmi : min sur (atome résidu × atome ligand) de pos.dist (C++).
gd = if is_prot == 1 {
per_atom = numpy.arange(n_at).tolist().[(ai) => {
at = r[int(ai)]
to_lig = numpy.arange(len(lig_res)).tolist().[(li) => {
float(at.pos.dist(lig_res[int(li)].pos))
}]
float(numpy.min(numpy.array(to_lig)))
}]
float(numpy.min(numpy.array(per_atom)))
} else { 0.0 }
ResidueScan(chain_name, r.name, r.seqid.num, cd, gd, is_prot)
}
# Balayage plat de tous les résidus de toutes les chaînes (nesté puis aplati).
per_chain = numpy.arange(len(model)).tolist().[(cix) => {
ch = model[int(cix)]
numpy.arange(len(ch)).tolist().[(rix) => {
scan_residue(ch.name, ch[int(rix)])
}]
}]
all_res = reduce(per_chain, (a, b) => { a + b })
# Sélection de la poche : protéique ET sous le seuil. Filtrage par masque numpy
# (Catnip ne filtre pas sur des appels de fonction), tri par distance croissante.
cds = numpy.array(all_res.[(rd) => { rd.cd }])
prot = numpy.array(all_res.[(rd) => { rd.is_protein }])
mask = numpy.logical_and(numpy.equal(prot, 1), numpy.less_equal(cds, THRESHOLD))
sel_idx = numpy.arange(len(all_res))[mask].tolist()
order = numpy.argsort(cds[mask]).tolist()
pocket = order.[(o) => {
rd = all_res[int(sel_idx[int(o)])]
PocketResidue(rd.chain, rd.name, rd.seqid, rd.cd, classify_contact(rd.cd))
}]
n_protein = int(numpy.sum(numpy.equal(prot, 1)))
mean_dist = float(numpy.mean(numpy.array(pocket.[(p) => { p.min_dist }])))
# Contrôle numérique : sur les résidus protéiques, écart entre cdist (scipy) et
# pos.dist (gemmi). Deux normes euclidiennes indépendantes sur les mêmes atomes :
# l'écart max doit être ~0 (valide l'extraction des coordonnées).
scan_err = sel_idx.[(k) => {
rd = all_res[int(k)]
abs(rd.cd - rd.gd)
}]
max_ctrl_err = float(numpy.max(numpy.array(scan_err)))
closest = pocket[0]
closest_str = f"{closest.chain} {closest.name}{closest.seqid} à {round(closest.min_dist, 2)} Å"
# --- Rapport à trois niveaux ---
print(f"⇒ Faits extraits du fichier ({PDB_ID})")
print(f" Provenance : {PROVENANCE}")
print(f" DOI structure : {DOI}")
print(f" Structure : {st.name}, résolution {st.resolution} Å, groupe d'espace {st.spacegroup_hm}")
print(f" Chaînes protéiques : {model[0].name} et {model[len(model) - 1].name} (homodimère)")
print(f" Ligand : {ligand.name} — {LIGAND_LABEL}, chaîne {ligand.chain}, seqid {ligand.seqid}, " +
f"{ligand.n_atoms} atomes")
print(" Exclusions : ligand MK1, 127 eaux (HOH), aucun autre HETATM ; symétrie non appliquée")
print(" Fait de littérature : protéase du VIH-1 = homodimère, dyade catalytique Asp25/Asp25'")
print()
report = NeighborhoodReport(len(pocket), n_protein, THRESHOLD, closest_str, mean_dist, max_ctrl_err)
report.display()
print()
print("⇒ Classification heuristique des contacts (par distance, indicative)")
codes = numpy.array(pocket.[(p) => { class_code(p.contact) }])
n_hbond = int(numpy.sum(numpy.equal(codes, 0)))
n_vdw = int(numpy.sum(numpy.equal(codes, 1)))
print(f" {ContactClass.hbond_range.label()} : {n_hbond} résidus")
print(f" {ContactClass.vdw_range.label()} : {n_vdw} résidus")
# Cohérence avec le site actif connu : présence d'Asp25 (chaîne A ou B) dans la poche.
asp25 = pocket.[(p) => { if p.name == 'ASP' and p.seqid == 25 { 1 } else { 0 } }]
n_asp25 = int(numpy.sum(numpy.array(asp25)))
if n_asp25 > 0 {
print(f" Note : Asp25 présent dans la poche ({n_asp25} copie(s)) — cohérent avec la dyade catalytique")
} else {
print(" Note : Asp25 absent de la poche au seuil courant")
}
print()
# Liste détaillée de la poche (résidu, distance min, classe heuristique).
print("⇒ Résidus de poche (triés par distance croissante)")
pocket.[(p) => {
print(f" {p.chain} {p.name}{p.seqid} : {round(p.min_dist, 2)} Å — {p.contact.label()}")
}]
print()
# --- Figure : profil de distance min par résidu de poche ---
labels = pocket.[(p) => { f"{p.chain} {p.name}{p.seqid}" }]
dists = numpy.array(pocket.[(p) => { p.min_dist }])
colors = pocket.[(p) => { class_color(p.contact) }]
ypos = numpy.arange(len(pocket))
fig, ax = plt.subplots(figsize=tuple(8, 10))
ax.barh(ypos, dists, color=colors, edgecolor='black', linewidth=0.4)
ax.axvline(HBOND_CUT, color='#d62728', linestyle='--', linewidth=1.0, label="seuil liaison H (3.5 Å)")
ax.axvline(THRESHOLD, color='#1f77b4', linestyle='--', linewidth=1.0, label="seuil poche (4.5 Å)")
ax.set_yticks(ypos.tolist())
ax.set_yticklabels(labels)
ax.invert_yaxis() # le plus proche en haut
ax.set_xlabel("distance min atome-résidu ↔ atome-ligand (Å)")
ax.set_title(f"Poche de liaison de l'indinavir ({PDB_ID}) — {len(pocket)} résidus")
ax.legend(loc='lower right')
fig.tight_layout()
output_path = output_dir / 'protein_ligand_neighborhood.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>Protein-ligand neighborhood</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)
}