Pyoverdine BGC Analysis

from rdkit import Chem
from rdkit.Chem import rdCoordGen
from rdkit.Chem.Draw import rdMolDraw2D
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
from PIL import Image
import io

# 1. CORRECTED SMILES
# - Removed extra 'N' at start
# - Changed 'CCCNO' to 'CCCN(O)' for proper Hydroxyornithine structure
smiles_stereo = (
    "N[C@@H](CCC(=O)O)C(=O)"       # L-Glu (atoms 0-9)
    "N[C@H](Cc1ccc(O)cc1)C(=O)"    # D-Tyr (atoms 10-21)
    "N[C@@H](CCN)C(=O)"            # L-Dab (atoms 22-28)
    "N[C@@H](CC(=O)O)C(=O)"        # L-Asp (atoms 29-36)
    "N[C@@H](C)C(=O)"              # L-Ala (atoms 37-41)
    "N[C@@H](CC(=O)O)C(=O)"        # L-Asp (atoms 42-49)
    "N[C@H](CCCN(O))C(=O)"         # D-OHOrn (atoms 50-58)
    "NCC(=O)"                      # Gly (atoms 59-62)
    "N[C@@H](CCCN(O))C(=O)O"       # L-OHOrn (atoms 63-72)
)

mol = Chem.MolFromSmiles(smiles_stereo)
rdCoordGen.AddCoords(mol)

# 2. UPDATED Atom Ranges (Shifted by -1 due to removal of extra N)
residue_indices = {
    'L-Glu': range(0, 9),      # 9 atoms
    'D-Tyr': range(9, 21),     # 12 atoms
    'L-Dab': range(21, 28),    # 7 atoms
    'L-Asp₁': range(28, 36),   # 8 atoms
    'L-Ala': range(36, 41),    # 5 atoms
    'L-Asp₂': range(41, 49),   # 8 atoms
    'D-OHOrn': range(49, 58),  # 9 atoms
    'Gly': range(58, 62),      # 4 atoms
    'L-OHOrn': range(62, 72)   # 10 atoms
}

# 3. Assign colors (Same palette)
colors = [
    (0.9, 0.3, 0.3),   # Red - Glu
    (0.9, 0.6, 0.2),   # Orange - Tyr
    (0.9, 0.9, 0.2),   # Yellow - Dab
    (0.4, 0.8, 0.4),   # Green - Asp1
    (0.3, 0.7, 0.7),   # Teal - Ala
    (0.3, 0.5, 0.9),   # Blue - Asp2
    (0.5, 0.3, 0.9),   # Purple - D-OHOrn
    (0.8, 0.4, 0.8),   # Magenta - Gly
    (0.9, 0.5, 0.6),   # Pink - L-OHOrn
]

highlight_atoms = {}
highlight_colors = {}

for i, (res_name, atom_range) in enumerate(residue_indices.items()):
    for atom_idx in atom_range:
        highlight_atoms[atom_idx] = colors[i]
        highlight_colors[atom_idx] = colors[i]

# 4. Draw the molecule
drawer = rdMolDraw2D.MolDraw2DCairo(1400, 800)
drawer.drawOptions().addAtomIndices = False
drawer.drawOptions().padding = 0.05
# Optional: Make bonds thicker for better visibility
drawer.drawOptions().bondLineWidth = 2
drawer.DrawMolecule(mol, highlightAtoms=list(highlight_atoms.keys()),
                    highlightAtomColors=highlight_colors)
drawer.FinishDrawing()

# 5. Combine with Legend in Matplotlib
mol_png = drawer.GetDrawingText()
mol_img = Image.open(io.BytesIO(mol_png))

fig, (ax_mol, ax_legend) = plt.subplots(1, 2, figsize=(18, 8), 
                                         gridspec_kw={'width_ratios': [2.5, 1]})

# Display Molecule
ax_mol.imshow(mol_img)
ax_mol.axis('off')
ax_mol.set_title('Predicted Pyoverdine Precursor (Ferribactin)\nPseudomonas putida (Houseplant Root Isolate)', 
                  fontsize=16, fontweight='bold', pad=20)

# Create Legend
legend_elements = [Patch(facecolor=c, label=l) for c, l in zip(colors, [
    '1. L-Glu (becomes side chain)',
    '2. D-Tyr (chromophore core)',
    '3. L-Dab (chromophore core)',
    '4. L-Asp',
    '5. L-Ala',
    '6. L-Asp',
    '7. D-OHOrn (Fe³⁺ ligand)',
    '8. Gly (predicted)',
    '9. L-OHOrn (Fe³⁺ ligand)'
])]

ax_legend.axis('off')
ax_legend.legend(handles=legend_elements, loc='upper right', fontsize=12, 
                  title='Peptide Sequence', title_fontsize=14, frameon=True,
                  fancybox=True, shadow=True)

# Add Annotation Box
annotations = """
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
STRUCTURE PREDICTION
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
► Side chain: Succinamide
  (inferred from pvdN presence)

► Chromophore: Dihydroxyquinoline
  (formed from residues 1-3 cyclization)

► Fe³⁺ Binding:
  Via 2× Hydroxamate groups (residues 7,9)
  and Chromophore catechol.

► Status: NOVEL VARIANT
  Sequence distinct from known P. putida types.
  "Gly" at pos 8 inferred from "NH2" prediction.
  LC-MS/MS required for final confirmation.
"""
ax_legend.text(0.5, 0.0, annotations, transform=ax_legend.transAxes, 
               fontsize=11, verticalalignment='bottom', horizontalalignment='center',
               family='monospace', bbox=dict(boxstyle='round', facecolor='#f8f9fa', alpha=1.0))

plt.tight_layout()
plt.savefig('pyoverdine_final_structure.png', dpi=150, bbox_inches='tight')
plt.show()