selective_fixed_labels_py

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyMZkIEcLfjNkvrywRTvzPn8"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "hJgO9_3tBneI"
      },
      "outputs": [],
      "source": [
        " #!/usr/bin/env python\n",
        "\"\"\"\n",
        "Selective FIXED Labels Script - Relative Positioning Version\n",
        "Removes FIXED labels from all positions except those with ≤4.0Å contacts\n",
        "Uses relative positioning within the FIXED region to work across all scaffolds\n",
        "\n",
        "To run:\n",
        "Activate conda base\n",
        "python selective_fixed_labels.py --pdbdir . --verbose --dry_run\n",
        "python selective_fixed_labels.py --pdbdir . --verbose\n",
        "\n",
        "OK this script is specificly for the output of a single RFdiffusion Motif Scaffolding job.\n",
        "PDB files had previpusly been modified to fix the motif that served to scaffold from using /helper_scripts/addFIXEDlabels.py from dl_binder\n",
        "I also has been through analysis using https://github.com/eoinleen/PDB-tools/blob/main/which_residues_to_fix_for_ProteinMPNN_v2.ipynb\n",
        "This is to fix only residues which have heavy atoms within 4.0 of heavy atoms in the target protein.\n",
        "\n",
        "How it works:\n",
        "The script reads existing REMARK PDBinfo-LABEL FIXED lines to identify the motif region,\n",
        "then uses relative positioning within this region to target specific residues. Since\n",
        "RFdiffusion generates variable-length scaffolds, the motif lands at different absolute\n",
        "positions in each structure. By using relative positions (0-indexed within the FIXED\n",
        "region), the same contact residues are identified across all scaffolds regardless of\n",
        "where the motif is positioned. Only residues at pre-defined relative positions\n",
        "corresponding to ≤4.0Å contacts are kept FIXED, while others are removed to allow\n",
        "ProteinMPNN optimization.\n",
        "\"\"\"\n",
        "\n",
        "import os\n",
        "import sys\n",
        "import argparse\n",
        "from collections import defaultdict\n",
        "\n",
        "def parse_args():\n",
        "    parser = argparse.ArgumentParser(description='Apply selective FIXED labels based on contact analysis using relative positioning')\n",
        "    parser.add_argument('--pdbdir', type=str, required=True, help='Directory containing PDB files with FIXED labels')\n",
        "    parser.add_argument('--chain', type=str, default='A', help='Chain to modify')\n",
        "    parser.add_argument('--verbose', action='store_true', help='Verbose output')\n",
        "    parser.add_argument('--dry_run', action='store_true', help='Show what would be changed without modifying files')\n",
        "    return parser.parse_args()\n",
        "\n",
        "def get_critical_relative_positions():\n",
        "    \"\"\"\n",
        "    Returns relative positions within the FIXED region that should remain FIXED\n",
        "    Based on ≤4.0Å contact analysis from the original motif A66-93\n",
        "\n",
        "    FIXED region sequence (28 residues):\n",
        "    Pos:  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27\n",
        "    Orig: 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93\n",
        "    Res:  ?  ?  ?  ?  T  G  L  Q  C  E  M  ?  P  V  ?  Q  T  ?  ?  S  E  ?  ?  ?  ?  ?  ?  ?\n",
        "    Dist: -  -  -  - 2.4 2.8 2.9 3.3 2.8 3.0 2.8  - 3.6 3.7  - 3.4 3.6  -  - 3.9 3.1  -  -  -  -  -  -  -\n",
        "    \"\"\"\n",
        "\n",
        "    # Critical contacts (≤4.0Å) at these relative positions within FIXED region:\n",
        "    critical_relative_positions = [\n",
        "        4,   # Original A70: THR (2.42Å)\n",
        "        5,   # Original A71: GLY (2.80Å)\n",
        "        6,   # Original A72: LEU (2.89Å) ← LEU→PHE target\n",
        "        7,   # Original A73: GLN (3.30Å)\n",
        "        8,   # Original A74: CYS (2.83Å) ← CYS→SER target\n",
        "        9,   # Original A75: GLU (2.95Å)\n",
        "        10,  # Original A76: MET (2.83Å)\n",
        "        12,  # Original A78: PRO (3.64Å)\n",
        "        13,  # Original A79: VAL (3.71Å)\n",
        "        15,  # Original A81: GLN (3.38Å)\n",
        "        16,  # Original A82: THR (3.57Å)\n",
        "        19,  # Original A85: SER (3.89Å)\n",
        "        20,  # Original A86: GLU (3.06Å)\n",
        "    ]\n",
        "\n",
        "    return critical_relative_positions\n",
        "\n",
        "def get_critical_contact_details():\n",
        "    \"\"\"Return detailed information about critical contacts for reporting\"\"\"\n",
        "    return [\n",
        "        (4, 70, 'THR', 2.42),   # Relative pos 4 = Original A70\n",
        "        (5, 71, 'GLY', 2.80),   # Relative pos 5 = Original A71\n",
        "        (6, 72, 'LEU', 2.89),   # Relative pos 6 = Original A72 (LEU→PHE target)\n",
        "        (7, 73, 'GLN', 3.30),   # Relative pos 7 = Original A73\n",
        "        (8, 74, 'CYS', 2.83),   # Relative pos 8 = Original A74 (CYS→SER target)\n",
        "        (9, 75, 'GLU', 2.95),   # Relative pos 9 = Original A75\n",
        "        (10, 76, 'MET', 2.83),  # Relative pos 10 = Original A76\n",
        "        (12, 78, 'PRO', 3.64),  # Relative pos 12 = Original A78\n",
        "        (13, 79, 'VAL', 3.71),  # Relative pos 13 = Original A79\n",
        "        (15, 81, 'GLN', 3.38),  # Relative pos 15 = Original A81\n",
        "        (16, 82, 'THR', 3.57),  # Relative pos 16 = Original A82\n",
        "        (19, 85, 'SER', 3.89),  # Relative pos 19 = Original A85\n",
        "        (20, 86, 'GLU', 3.06),  # Relative pos 20 = Original A86\n",
        "    ]\n",
        "\n",
        "def find_all_fixed_positions(pdb_file):\n",
        "    \"\"\"Extract all FIXED positions from REMARK lines\"\"\"\n",
        "    fixed_positions = []\n",
        "    with open(pdb_file, 'r') as f:\n",
        "        for line in f:\n",
        "            if line.startswith('REMARK PDBinfo-LABEL:') and 'FIXED' in line:\n",
        "                # Extract position number\n",
        "                parts = line.split()\n",
        "                try:\n",
        "                    pos = int(parts[2])\n",
        "                    fixed_positions.append(pos)\n",
        "                except (ValueError, IndexError):\n",
        "                    continue\n",
        "    return sorted(fixed_positions)\n",
        "\n",
        "def convert_relative_to_absolute_positions(fixed_positions, critical_relative_positions):\n",
        "    \"\"\"\n",
        "    Convert relative positions within FIXED region to absolute PDB positions\n",
        "    \"\"\"\n",
        "    if len(fixed_positions) < max(critical_relative_positions) + 1:\n",
        "        raise ValueError(f\"FIXED region too short: {len(fixed_positions)} positions, need at least {max(critical_relative_positions) + 1}\")\n",
        "\n",
        "    critical_absolute_positions = []\n",
        "    for rel_pos in critical_relative_positions:\n",
        "        abs_pos = fixed_positions[rel_pos]\n",
        "        critical_absolute_positions.append(abs_pos)\n",
        "\n",
        "    return sorted(critical_absolute_positions)\n",
        "\n",
        "def remove_non_critical_fixed_labels(pdb_file, critical_positions, dry_run=False):\n",
        "    \"\"\"\n",
        "    Remove FIXED labels from positions not in critical_positions list\n",
        "    \"\"\"\n",
        "    if dry_run:\n",
        "        return True, 0, 0\n",
        "\n",
        "    # Read the file\n",
        "    with open(pdb_file, 'r') as f:\n",
        "        lines = f.readlines()\n",
        "\n",
        "    # Remove non-critical FIXED labels\n",
        "    new_lines = []\n",
        "    removed_count = 0\n",
        "    kept_count = 0\n",
        "\n",
        "    for line in lines:\n",
        "        if line.startswith('REMARK PDBinfo-LABEL:') and 'FIXED' in line:\n",
        "            # Extract position number\n",
        "            parts = line.split()\n",
        "            try:\n",
        "                pos = int(parts[2])\n",
        "                if pos in critical_positions:\n",
        "                    # Keep this FIXED label\n",
        "                    new_lines.append(line)\n",
        "                    kept_count += 1\n",
        "                else:\n",
        "                    # Remove this FIXED label\n",
        "                    removed_count += 1\n",
        "                    # Don't append to new_lines\n",
        "            except (ValueError, IndexError):\n",
        "                # Keep malformed lines as-is\n",
        "                new_lines.append(line)\n",
        "        else:\n",
        "            # Keep all non-FIXED lines\n",
        "            new_lines.append(line)\n",
        "\n",
        "    if removed_count > 0 or kept_count > 0:\n",
        "        # Write back to file\n",
        "        with open(pdb_file, 'w') as f:\n",
        "            f.writelines(new_lines)\n",
        "        return True, removed_count, kept_count\n",
        "\n",
        "    return False, 0, 0\n",
        "\n",
        "def process_scaffold(pdb_file, critical_relative_positions, verbose, dry_run):\n",
        "    \"\"\"Process a single scaffold file\"\"\"\n",
        "    basename = os.path.basename(pdb_file)\n",
        "\n",
        "    if verbose:\n",
        "        print(f\"\\nProcessing: {basename}\")\n",
        "\n",
        "    try:\n",
        "        # Find all current FIXED positions\n",
        "        current_fixed = find_all_fixed_positions(pdb_file)\n",
        "        if verbose:\n",
        "            print(f\"  FIXED region: {len(current_fixed)} positions ({current_fixed[0]}-{current_fixed[-1]})\")\n",
        "\n",
        "        # Convert relative to absolute positions\n",
        "        critical_absolute = convert_relative_to_absolute_positions(current_fixed, critical_relative_positions)\n",
        "\n",
        "        if verbose:\n",
        "            print(f\"  Critical contacts (≤4.0Å): {critical_absolute}\")\n",
        "\n",
        "            # Show the mapping for verification\n",
        "            contact_details = get_critical_contact_details()\n",
        "            print(f\"  Contact mapping:\")\n",
        "            for rel_pos, orig_pos, res_name, distance in contact_details:\n",
        "                abs_pos = current_fixed[rel_pos]\n",
        "                print(f\"    Rel{rel_pos}(A{orig_pos}:{res_name}) → Pos{abs_pos} ({distance}Å)\")\n",
        "\n",
        "        # Determine which will be removed\n",
        "        positions_to_remove = [pos for pos in current_fixed if pos not in critical_absolute]\n",
        "        positions_to_keep = [pos for pos in current_fixed if pos in critical_absolute]\n",
        "\n",
        "        if verbose:\n",
        "            print(f\"  Will remove FIXED from: {positions_to_remove}\")\n",
        "            print(f\"  Will keep FIXED on: {positions_to_keep}\")\n",
        "\n",
        "        if not positions_to_remove:\n",
        "            print(f\"  ✅ No changes needed - all FIXED positions are critical\")\n",
        "            return {'status': 'no_change', 'file': basename, 'fixed_region': f\"{current_fixed[0]}-{current_fixed[-1]}\"}\n",
        "\n",
        "        # Apply changes\n",
        "        if dry_run:\n",
        "            print(f\"  [DRY RUN] Would remove {len(positions_to_remove)} FIXED labels\")\n",
        "            print(f\"  [DRY RUN] Would keep {len(positions_to_keep)} FIXED labels\")\n",
        "            return {'status': 'dry_run', 'file': basename,\n",
        "                    'removed': len(positions_to_remove), 'kept': len(positions_to_keep),\n",
        "                    'fixed_region': f\"{current_fixed[0]}-{current_fixed[-1]}\"}\n",
        "        else:\n",
        "            success, removed_count, kept_count = remove_non_critical_fixed_labels(\n",
        "                pdb_file, critical_absolute, dry_run)\n",
        "\n",
        "            if success:\n",
        "                print(f\"  ✅ Removed {removed_count} FIXED labels, kept {kept_count}\")\n",
        "                return {'status': 'modified', 'file': basename,\n",
        "                        'removed': removed_count, 'kept': kept_count,\n",
        "                        'fixed_region': f\"{current_fixed[0]}-{current_fixed[-1]}\"}\n",
        "            else:\n",
        "                print(f\"  ❌ Failed to modify FIXED labels\")\n",
        "                return {'status': 'failed', 'file': basename}\n",
        "\n",
        "    except Exception as e:\n",
        "        print(f\"  ❌ Error processing file: {e}\")\n",
        "        return {'status': 'error', 'file': basename, 'error': str(e)}\n",
        "\n",
        "def main():\n",
        "    args = parse_args()\n",
        "\n",
        "    # Get critical relative positions based on contact analysis\n",
        "    critical_relative_positions = get_critical_relative_positions()\n",
        "    contact_details = get_critical_contact_details()\n",
        "\n",
        "    print(\"=\"*80)\n",
        "    print(\"SELECTIVE FIXED LABELS SCRIPT - RELATIVE POSITIONING VERSION\")\n",
        "    print(\"=\"*80)\n",
        "    print(f\"Target directory: {args.pdbdir}\")\n",
        "    print(f\"Chain: {args.chain}\")\n",
        "    print(f\"Mode: {'DRY RUN' if args.dry_run else 'LIVE MODIFICATION'}\")\n",
        "    print(f\"Critical relative positions: {critical_relative_positions}\")\n",
        "    print(\"=\"*80)\n",
        "\n",
        "    # Contact analysis summary\n",
        "    print(\"Contact analysis basis (relative positions within FIXED region):\")\n",
        "    for rel_pos, orig_pos, res_name, distance in contact_details:\n",
        "        print(f\"  Rel{rel_pos:2d} (A{orig_pos}:{res_name}) = {distance}Å\")\n",
        "    print(\"=\"*80)\n",
        "\n",
        "    # Find PDB files\n",
        "    pdb_files = [f for f in os.listdir(args.pdbdir) if f.endswith('.pdb')]\n",
        "    if not pdb_files:\n",
        "        print(\"No PDB files found in directory\")\n",
        "        return\n",
        "\n",
        "    print(f\"Found {len(pdb_files)} PDB files\")\n",
        "\n",
        "    # Process each file\n",
        "    results = []\n",
        "    for pdb_file in sorted(pdb_files):\n",
        "        pdb_path = os.path.join(args.pdbdir, pdb_file)\n",
        "        result = process_scaffold(pdb_path, critical_relative_positions, args.verbose, args.dry_run)\n",
        "        results.append(result)\n",
        "\n",
        "    # Summary\n",
        "    print(\"\\n\" + \"=\"*80)\n",
        "    print(\"SUMMARY\")\n",
        "    print(\"=\"*80)\n",
        "\n",
        "    status_counts = defaultdict(int)\n",
        "    total_removed = 0\n",
        "    total_kept = 0\n",
        "\n",
        "    for result in results:\n",
        "        status_counts[result['status']] += 1\n",
        "        if 'removed' in result:\n",
        "            total_removed += result['removed']\n",
        "        if 'kept' in result:\n",
        "            total_kept += result['kept']\n",
        "\n",
        "    for status, count in status_counts.items():\n",
        "        print(f\"{status}: {count} files\")\n",
        "\n",
        "    if total_removed > 0 or total_kept > 0:\n",
        "        print(f\"\\nTotal FIXED labels removed: {total_removed}\")\n",
        "        print(f\"Total FIXED labels kept: {total_kept}\")\n",
        "\n",
        "    # Show FIXED regions found\n",
        "    print(f\"\\nFIXED regions found:\")\n",
        "    for result in results:\n",
        "        if 'fixed_region' in result:\n",
        "            print(f\"  {result['file']}: positions {result['fixed_region']}\")\n",
        "\n",
        "    print(\"\\nNext steps:\")\n",
        "    print(\"1. Run ProteinMPNN with these selective FIXED labels\")\n",
        "    print(\"2. Apply CYS→SER mutations (Step 2)\")\n",
        "    print(\"3. Apply LEU→PHE mutations (Step 3)\")\n",
        "    print(\"=\"*80)\n",
        "\n",
        "if __name__ == \"__main__\":\n",
        "    main()"
      ]
    }
  ]
}