{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# 16元环大环内酯分子分析\n",
        "\n",
        "本notebook用于分析 ring16_match_smarts.csv 中的16元环大环内酯分子：\n",
        "1. 计算并统计分子的相对分子质量\n",
        "2. 进行侧链断裂分析\n",
        "3. 统计每个碎片位置的原子数目和分子量分布\n",
        "4. 计算并分析药物性质（LogP、QED等）\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1. 导入必要的库\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "所有模块导入成功!\n"
          ]
        }
      ],
      "source": [
        "import sys\n",
        "from pathlib import Path\n",
        "\n",
        "# 添加项目根目录到 Python 路径\n",
        "notebook_dir = Path().resolve()\n",
        "project_root = notebook_dir.parent\n",
        "sys.path.insert(0, str(project_root))\n",
        "\n",
        "from rdkit import Chem\n",
        "from rdkit.Chem import Descriptors, Crippen, QED\n",
        "from rdkit.Chem import Draw\n",
        "from rdkit.Chem.Draw import rdMolDraw2D\n",
        "from IPython.display import SVG, display\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "\n",
        "# 导入项目模块\n",
        "from src.macrolactone_fragmenter import MacrolactoneFragmenter\n",
        "from src.ring_visualization import get_macrolactone_numbering\n",
        "\n",
        "print(\"所有模块导入成功!\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 2. 读取数据\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "读取数据集: 307 个分子\n",
            "列名: ['Unnamed: 0', 'unique_id', 'IDs', 'molecule_pref_name', 'max_pChEMBL', 'max_pChEMBL_target', '# Target Organisms', 'Target Organisms', '# Known Targets', 'Known Targets', 'target_pref_name', 'smiles', 'has_bridge', 'has_multiple_lactones', 'has_amide', 'parse_failed', 'filter_reason', 'matches_smarts']\n",
            "\n",
            "前5行:\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Unnamed: 0</th>\n",
              "      <th>unique_id</th>\n",
              "      <th>IDs</th>\n",
              "      <th>molecule_pref_name</th>\n",
              "      <th>max_pChEMBL</th>\n",
              "      <th>max_pChEMBL_target</th>\n",
              "      <th># Target Organisms</th>\n",
              "      <th>Target Organisms</th>\n",
              "      <th># Known Targets</th>\n",
              "      <th>Known Targets</th>\n",
              "      <th>target_pref_name</th>\n",
              "      <th>smiles</th>\n",
              "      <th>has_bridge</th>\n",
              "      <th>has_multiple_lactones</th>\n",
              "      <th>has_amide</th>\n",
              "      <th>parse_failed</th>\n",
              "      <th>filter_reason</th>\n",
              "      <th>matches_smarts</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>113</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CHEMBL1354988</td>\n",
              "      <td>NaN</td>\n",
              "      <td>4.85</td>\n",
              "      <td>CHEMBL1293235</td>\n",
              "      <td>3</td>\n",
              "      <td>Homo sapiens, None</td>\n",
              "      <td>3</td>\n",
              "      <td>CHEMBL1293226, CHEMBL612545, CHEMBL1293235</td>\n",
              "      <td>Unchecked, Lysine-specific demethylase 4D-like...</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>NaN</td>\n",
              "      <td>True</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>288</td>\n",
              "      <td>mol_00288</td>\n",
              "      <td>CHEMBL42743</td>\n",
              "      <td>TYLOSIN</td>\n",
              "      <td>8.56</td>\n",
              "      <td>CHEMBL345</td>\n",
              "      <td>1</td>\n",
              "      <td>NaN</td>\n",
              "      <td>1</td>\n",
              "      <td>CHEMBL345</td>\n",
              "      <td>Nucleic Acid</td>\n",
              "      <td>CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>NaN</td>\n",
              "      <td>True</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>313</td>\n",
              "      <td>mol_00313</td>\n",
              "      <td>CHEMBL42743</td>\n",
              "      <td>TYLOSIN</td>\n",
              "      <td>8.56</td>\n",
              "      <td>CHEMBL345</td>\n",
              "      <td>1</td>\n",
              "      <td>NaN</td>\n",
              "      <td>1</td>\n",
              "      <td>CHEMBL345</td>\n",
              "      <td>Nucleic Acid</td>\n",
              "      <td>CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>NaN</td>\n",
              "      <td>True</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>314</td>\n",
              "      <td>mol_00314</td>\n",
              "      <td>CHEMBL42743</td>\n",
              "      <td>TYLOSIN</td>\n",
              "      <td>8.56</td>\n",
              "      <td>CHEMBL345</td>\n",
              "      <td>1</td>\n",
              "      <td>NaN</td>\n",
              "      <td>1</td>\n",
              "      <td>CHEMBL345</td>\n",
              "      <td>Nucleic Acid</td>\n",
              "      <td>CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>NaN</td>\n",
              "      <td>True</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>326</td>\n",
              "      <td>mol_00326</td>\n",
              "      <td>CHEMBL42743</td>\n",
              "      <td>TYLOSIN</td>\n",
              "      <td>8.56</td>\n",
              "      <td>CHEMBL345</td>\n",
              "      <td>1</td>\n",
              "      <td>NaN</td>\n",
              "      <td>1</td>\n",
              "      <td>CHEMBL345</td>\n",
              "      <td>Nucleic Acid</td>\n",
              "      <td>CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>False</td>\n",
              "      <td>NaN</td>\n",
              "      <td>True</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "   Unnamed: 0  unique_id            IDs molecule_pref_name  max_pChEMBL  \\\n",
              "0         113  mol_00113  CHEMBL1354988                NaN         4.85   \n",
              "1         288  mol_00288    CHEMBL42743            TYLOSIN         8.56   \n",
              "2         313  mol_00313    CHEMBL42743            TYLOSIN         8.56   \n",
              "3         314  mol_00314    CHEMBL42743            TYLOSIN         8.56   \n",
              "4         326  mol_00326    CHEMBL42743            TYLOSIN         8.56   \n",
              "\n",
              "  max_pChEMBL_target  # Target Organisms    Target Organisms  # Known Targets  \\\n",
              "0      CHEMBL1293235                   3  Homo sapiens, None                3   \n",
              "1          CHEMBL345                   1                 NaN                1   \n",
              "2          CHEMBL345                   1                 NaN                1   \n",
              "3          CHEMBL345                   1                 NaN                1   \n",
              "4          CHEMBL345                   1                 NaN                1   \n",
              "\n",
              "                                Known Targets  \\\n",
              "0  CHEMBL1293226, CHEMBL612545, CHEMBL1293235   \n",
              "1                                   CHEMBL345   \n",
              "2                                   CHEMBL345   \n",
              "3                                   CHEMBL345   \n",
              "4                                   CHEMBL345   \n",
              "\n",
              "                                    target_pref_name  \\\n",
              "0  Unchecked, Lysine-specific demethylase 4D-like...   \n",
              "1                                       Nucleic Acid   \n",
              "2                                       Nucleic Acid   \n",
              "3                                       Nucleic Acid   \n",
              "4                                       Nucleic Acid   \n",
              "\n",
              "                                              smiles  has_bridge  \\\n",
              "0  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...       False   \n",
              "1  CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...       False   \n",
              "2  CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...       False   \n",
              "3  CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...       False   \n",
              "4  CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...       False   \n",
              "\n",
              "   has_multiple_lactones  has_amide  parse_failed  filter_reason  \\\n",
              "0                  False      False         False            NaN   \n",
              "1                  False      False         False            NaN   \n",
              "2                  False      False         False            NaN   \n",
              "3                  False      False         False            NaN   \n",
              "4                  False      False         False            NaN   \n",
              "\n",
              "   matches_smarts  \n",
              "0            True  \n",
              "1            True  \n",
              "2            True  \n",
              "3            True  \n",
              "4            True  "
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 读取CSV文件\n",
        "csv_file = project_root / 'output' / 'ring16_match_smarts.csv'\n",
        "df = pd.read_csv(csv_file)\n",
        "\n",
        "print(f\"读取数据集: {len(df)} 个分子\")\n",
        "print(f\"列名: {df.columns.tolist()}\")\n",
        "print(f\"\\n前5行:\")\n",
        "df.head()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 3. 计算分子的基本性质\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "正在计算分子性质...\n",
            "\n",
            "成功计算 307 个分子的性质\n",
            "\n",
            "性质统计:\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>mol_weight</th>\n",
              "      <th>logP</th>\n",
              "      <th>num_h_donors</th>\n",
              "      <th>num_h_acceptors</th>\n",
              "      <th>num_rotatable_bonds</th>\n",
              "      <th>tpsa</th>\n",
              "      <th>qed</th>\n",
              "      <th>num_atoms</th>\n",
              "      <th>num_heavy_atoms</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>count</th>\n",
              "      <td>307.000000</td>\n",
              "      <td>307.000000</td>\n",
              "      <td>307.000000</td>\n",
              "      <td>307.000000</td>\n",
              "      <td>307.000000</td>\n",
              "      <td>307.000000</td>\n",
              "      <td>307.000000</td>\n",
              "      <td>307.000000</td>\n",
              "      <td>307.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>mean</th>\n",
              "      <td>786.351365</td>\n",
              "      <td>3.028585</td>\n",
              "      <td>3.781759</td>\n",
              "      <td>13.990228</td>\n",
              "      <td>10.895765</td>\n",
              "      <td>186.006678</td>\n",
              "      <td>0.226296</td>\n",
              "      <td>55.074919</td>\n",
              "      <td>55.074919</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>std</th>\n",
              "      <td>161.791942</td>\n",
              "      <td>1.035552</td>\n",
              "      <td>1.221242</td>\n",
              "      <td>3.663533</td>\n",
              "      <td>3.604946</td>\n",
              "      <td>45.728840</td>\n",
              "      <td>0.129689</td>\n",
              "      <td>11.203885</td>\n",
              "      <td>11.203885</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>min</th>\n",
              "      <td>364.482000</td>\n",
              "      <td>0.525300</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>5.000000</td>\n",
              "      <td>1.000000</td>\n",
              "      <td>80.670000</td>\n",
              "      <td>0.065000</td>\n",
              "      <td>26.000000</td>\n",
              "      <td>26.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>25%</th>\n",
              "      <td>679.938000</td>\n",
              "      <td>2.215600</td>\n",
              "      <td>3.000000</td>\n",
              "      <td>11.000000</td>\n",
              "      <td>8.000000</td>\n",
              "      <td>152.060000</td>\n",
              "      <td>0.146527</td>\n",
              "      <td>48.000000</td>\n",
              "      <td>48.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>50%</th>\n",
              "      <td>811.963000</td>\n",
              "      <td>3.145820</td>\n",
              "      <td>4.000000</td>\n",
              "      <td>15.000000</td>\n",
              "      <td>11.000000</td>\n",
              "      <td>195.150000</td>\n",
              "      <td>0.182802</td>\n",
              "      <td>57.000000</td>\n",
              "      <td>57.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>75%</th>\n",
              "      <td>902.584000</td>\n",
              "      <td>3.717900</td>\n",
              "      <td>5.000000</td>\n",
              "      <td>16.500000</td>\n",
              "      <td>13.000000</td>\n",
              "      <td>218.440000</td>\n",
              "      <td>0.249791</td>\n",
              "      <td>63.000000</td>\n",
              "      <td>63.000000</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>max</th>\n",
              "      <td>1098.334000</td>\n",
              "      <td>5.864900</td>\n",
              "      <td>7.000000</td>\n",
              "      <td>22.000000</td>\n",
              "      <td>18.000000</td>\n",
              "      <td>289.470000</td>\n",
              "      <td>0.725918</td>\n",
              "      <td>78.000000</td>\n",
              "      <td>78.000000</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "        mol_weight        logP  num_h_donors  num_h_acceptors  \\\n",
              "count   307.000000  307.000000    307.000000       307.000000   \n",
              "mean    786.351365    3.028585      3.781759        13.990228   \n",
              "std     161.791942    1.035552      1.221242         3.663533   \n",
              "min     364.482000    0.525300      1.000000         5.000000   \n",
              "25%     679.938000    2.215600      3.000000        11.000000   \n",
              "50%     811.963000    3.145820      4.000000        15.000000   \n",
              "75%     902.584000    3.717900      5.000000        16.500000   \n",
              "max    1098.334000    5.864900      7.000000        22.000000   \n",
              "\n",
              "       num_rotatable_bonds        tpsa         qed   num_atoms  \\\n",
              "count           307.000000  307.000000  307.000000  307.000000   \n",
              "mean             10.895765  186.006678    0.226296   55.074919   \n",
              "std               3.604946   45.728840    0.129689   11.203885   \n",
              "min               1.000000   80.670000    0.065000   26.000000   \n",
              "25%               8.000000  152.060000    0.146527   48.000000   \n",
              "50%              11.000000  195.150000    0.182802   57.000000   \n",
              "75%              13.000000  218.440000    0.249791   63.000000   \n",
              "max              18.000000  289.470000    0.725918   78.000000   \n",
              "\n",
              "       num_heavy_atoms  \n",
              "count       307.000000  \n",
              "mean         55.074919  \n",
              "std          11.203885  \n",
              "min          26.000000  \n",
              "25%          48.000000  \n",
              "50%          57.000000  \n",
              "75%          63.000000  \n",
              "max          78.000000  "
            ]
          },
          "execution_count": 3,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 计算基本性质\n",
        "def calculate_properties(smiles):\n",
        "    \"\"\"计算分子的基本性质\"\"\"\n",
        "    mol = Chem.MolFromSmiles(smiles)\n",
        "    if mol is None:\n",
        "        return None\n",
        "    \n",
        "    props = {\n",
        "        'mol_weight': Descriptors.MolWt(mol),\n",
        "        'logP': Crippen.MolLogP(mol),\n",
        "        'num_h_donors': Descriptors.NumHDonors(mol),\n",
        "        'num_h_acceptors': Descriptors.NumHAcceptors(mol),\n",
        "        'num_rotatable_bonds': Descriptors.NumRotatableBonds(mol),\n",
        "        'tpsa': Descriptors.TPSA(mol),\n",
        "        'qed': QED.qed(mol),\n",
        "        'num_atoms': mol.GetNumAtoms(),\n",
        "        'num_heavy_atoms': mol.GetNumHeavyAtoms()\n",
        "    }\n",
        "    return props\n",
        "\n",
        "# 计算所有分子的性质\n",
        "print(\"正在计算分子性质...\")\n",
        "properties_list = []\n",
        "\n",
        "for idx, row in df.iterrows():\n",
        "    props = calculate_properties(row['smiles'])\n",
        "    if props:\n",
        "        props['unique_id'] = row['unique_id']\n",
        "        properties_list.append(props)\n",
        "\n",
        "# 创建性质DataFrame\n",
        "df_props = pd.DataFrame(properties_list)\n",
        "\n",
        "print(f\"\\n成功计算 {len(df_props)} 个分子的性质\")\n",
        "print(f\"\\n性质统计:\")\n",
        "df_props.describe()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 4. 分子性质分布统计图\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/tmp/ipykernel_3752769/2194066832.py:57: UserWarning: Glyph 20540 (\\N{CJK UNIFIED IDEOGRAPH-503C}) missing from font(s) DejaVu Sans.\n",
            "  plt.tight_layout()\n",
            "/tmp/ipykernel_3752769/2194066832.py:61: UserWarning: Glyph 20540 (\\N{CJK UNIFIED IDEOGRAPH-503C}) missing from font(s) DejaVu Sans.\n",
            "  plt.savefig(output_file, dpi=300, bbox_inches='tight')\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "已保存: /home/zly/project/macro_split/output/ring16_molecular_properties_distribution.png\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "/home/zly/project/macro_split/.pixi/envs/default/lib/python3.14/site-packages/IPython/core/pylabtools.py:170: UserWarning: Glyph 20540 (\\N{CJK UNIFIED IDEOGRAPH-503C}) missing from font(s) DejaVu Sans.\n",
            "  fig.canvas.print_figure(bytes_io, **kw)\n"
          ]
        },
        {
          "data": {
            "image/png": "[encoded data, 174568 characters]",
            "text/plain": [
              "<Figure size 1600x1200 with 4 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# 导入绘图库\n",
        "try:\n",
        "    import seaborn as sns\n",
        "    import matplotlib.pyplot as plt\n",
        "except ImportError:\n",
        "    print(\"需要安装seaborn和matplotlib:\")\n",
        "    print(\"pixi add seaborn matplotlib\")\n",
        "    raise\n",
        "\n",
        "# 设置绘图样式\n",
        "sns.set_style(\"whitegrid\")\n",
        "\n",
        "# 绘制分子质量分布图\n",
        "fig, axes = plt.subplots(2, 2, figsize=(16, 12))\n",
        "fig.suptitle('Distribution of 16-membered Macrocycle Esters Molecular Properties', fontsize=16, fontweight='bold')\n",
        "\n",
        "# 1. 分子量分布\n",
        "ax = axes[0, 0]\n",
        "sns.histplot(df_props['mol_weight'], kde=True, ax=ax, bins=30, color='steelblue', alpha=0.7)\n",
        "ax.set_title(f'Molecular Weight Distribution  (n={len(df_props)})', fontsize=12)\n",
        "ax.set_xlabel('Molecular Weight (Da)', fontsize=10)\n",
        "ax.set_ylabel('Frequency', fontsize=10)\n",
        "ax.axvline(df_props['mol_weight'].mean(), color='red', linestyle='--', label=f'Mean: {df_props[\"mol_weight\"].mean():.1f}')\n",
        "ax.legend()\n",
        "ax.grid(True, alpha=0.3)\n",
        "\n",
        "# 2. LogP分布\n",
        "ax = axes[0, 1]\n",
        "sns.histplot(df_props['logP'], kde=True, ax=ax, bins=30, color='green', alpha=0.7)\n",
        "ax.set_title(f'LogP Distribution (n={len(df_props)})', fontsize=12)\n",
        "ax.set_xlabel('LogP', fontsize=10)\n",
        "ax.set_ylabel('Frequency', fontsize=10)\n",
        "ax.axvline(df_props['logP'].mean(), color='red', linestyle='--', label=f'Mean: {df_props[\"logP\"].mean():.2f}')\n",
        "ax.legend()\n",
        "ax.grid(True, alpha=0.3)\n",
        "\n",
        "# 3. QED分布\n",
        "ax = axes[1, 0]\n",
        "sns.histplot(df_props['qed'], kde=True, ax=ax, bins=30, color='orange', alpha=0.7)\n",
        "ax.set_title(f'QED Distribution (n={len(df_props)})', fontsize=12)\n",
        "ax.set_xlabel('QED值', fontsize=10)\n",
        "ax.set_ylabel('Frequency', fontsize=10)\n",
        "ax.axvline(df_props['qed'].mean(), color='red', linestyle='--', label=f'Mean: {df_props[\"qed\"].mean():.3f}')\n",
        "ax.legend()\n",
        "ax.grid(True, alpha=0.3)\n",
        "\n",
        "# 4. TPSA分布\n",
        "ax = axes[1, 1]\n",
        "sns.histplot(df_props['tpsa'], kde=True, ax=ax, bins=30, color='purple', alpha=0.7)\n",
        "ax.set_title(f'TPSA Distribution (n={len(df_props)})', fontsize=12)\n",
        "ax.set_xlabel('TPSA (Ų)', fontsize=10)\n",
        "ax.set_ylabel('Frequency', fontsize=10)\n",
        "ax.axvline(df_props['tpsa'].mean(), color='red', linestyle='--', label=f'Mean: {df_props[\"tpsa\"].mean():.1f}')\n",
        "ax.legend()\n",
        "ax.grid(True, alpha=0.3)\n",
        "\n",
        "plt.tight_layout()\n",
        "\n",
        "# 保存图片\n",
        "output_file = project_root / 'output' / 'ring16_molecular_properties_distribution.png'\n",
        "plt.savefig(output_file, dpi=300, bbox_inches='tight')\n",
        "print(f\"已保存: {output_file}\")\n",
        "plt.show()\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "MacrolactoneFragmenter 已初始化\n",
            "  环大小: 16\n",
            "\n",
            "准备开始批量处理...\n"
          ]
        }
      ],
      "source": [
        "# 创建16元环大环内酯断裂器\n",
        "fragmenter = MacrolactoneFragmenter(ring_size=16)\n",
        "\n",
        "print(f\"MacrolactoneFragmenter 已初始化\")\n",
        "print(f\"  环大小: {fragmenter.ring_size}\")\n",
        "print(f\"\\n准备开始批量处理...\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "开始批量处理 307 个分子...\n",
            "\n",
            "\n",
            "Processed 307 molecules successfully\n",
            "\n",
            "批处理统计:\n",
            "  成功处理: 307 个分子\n",
            "  总碎片数: 2684\n",
            "  平均碎片数: 8.74\n",
            "  碎片数范围: 7 - 9\n"
          ]
        }
      ],
      "source": [
        "# 批量处理所有分子\n",
        "print(f\"开始批量处理 {len(df)} 个分子...\\n\")\n",
        "\n",
        "batch_results = fragmenter.process_csv(\n",
        "    csv_file,\n",
        "    smiles_column=\"smiles\",\n",
        "    id_column=\"unique_id\",\n",
        "    max_rows=None  # 处理所有分子\n",
        ")\n",
        "\n",
        "print(f\"\\n批处理统计:\")\n",
        "print(f\"  成功处理: {len(batch_results)} 个分子\")\n",
        "\n",
        "total_fragments = sum(len(r.fragments) for r in batch_results)\n",
        "print(f\"  总碎片数: {total_fragments}\")\n",
        "\n",
        "# 统计每个分子的碎片数\n",
        "fragment_counts = [len(r.fragments) for r in batch_results]\n",
        "print(f\"  平均碎片数: {np.mean(fragment_counts):.2f}\")\n",
        "print(f\"  碎片数范围: {np.min(fragment_counts)} - {np.max(fragment_counts)}\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "碎片DataFrame:\n",
            "  形状: (2684, 7)\n",
            "  列: ['fragment_id', 'parent_id', 'parent_smiles', 'cleavage_position', 'fragment_smiles', 'atom_count', 'molecular_weight']\n",
            "\n",
            "前10行:\n"
          ]
        },
        {
          "data": {
            "text/html": [
              "<div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>fragment_id</th>\n",
              "      <th>parent_id</th>\n",
              "      <th>parent_smiles</th>\n",
              "      <th>cleavage_position</th>\n",
              "      <th>fragment_smiles</th>\n",
              "      <th>atom_count</th>\n",
              "      <th>molecular_weight</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>0</th>\n",
              "      <td>mol_00113_frag_0</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>4</td>\n",
              "      <td>*O</td>\n",
              "      <td>2</td>\n",
              "      <td>17.01</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>mol_00113_frag_1</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>5</td>\n",
              "      <td>*C</td>\n",
              "      <td>2</td>\n",
              "      <td>15.04</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>mol_00113_frag_2</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>6</td>\n",
              "      <td>*OC1OC(C)C(O)C(N(C)C)C1O</td>\n",
              "      <td>14</td>\n",
              "      <td>190.22</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>mol_00113_frag_3</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>7</td>\n",
              "      <td>*CCN1CC(C)CC(C)C1</td>\n",
              "      <td>11</td>\n",
              "      <td>140.25</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>mol_00113_frag_4</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>9</td>\n",
              "      <td>*C</td>\n",
              "      <td>2</td>\n",
              "      <td>15.04</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>mol_00113_frag_5</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>10</td>\n",
              "      <td>*=O</td>\n",
              "      <td>2</td>\n",
              "      <td>16.00</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>mol_00113_frag_6</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>13</td>\n",
              "      <td>*C</td>\n",
              "      <td>2</td>\n",
              "      <td>15.04</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>mol_00113_frag_7</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>15</td>\n",
              "      <td>*COC1OC(C)C(O)C(OC)C1OC</td>\n",
              "      <td>15</td>\n",
              "      <td>205.23</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>8</th>\n",
              "      <td>mol_00113_frag_8</td>\n",
              "      <td>mol_00113</td>\n",
              "      <td>CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...</td>\n",
              "      <td>16</td>\n",
              "      <td>*CC</td>\n",
              "      <td>3</td>\n",
              "      <td>29.06</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>mol_00288_frag_0</td>\n",
              "      <td>mol_00288</td>\n",
              "      <td>CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...</td>\n",
              "      <td>4</td>\n",
              "      <td>*O</td>\n",
              "      <td>2</td>\n",
              "      <td>17.01</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "        fragment_id  parent_id  \\\n",
              "0  mol_00113_frag_0  mol_00113   \n",
              "1  mol_00113_frag_1  mol_00113   \n",
              "2  mol_00113_frag_2  mol_00113   \n",
              "3  mol_00113_frag_3  mol_00113   \n",
              "4  mol_00113_frag_4  mol_00113   \n",
              "5  mol_00113_frag_5  mol_00113   \n",
              "6  mol_00113_frag_6  mol_00113   \n",
              "7  mol_00113_frag_7  mol_00113   \n",
              "8  mol_00113_frag_8  mol_00113   \n",
              "9  mol_00288_frag_0  mol_00288   \n",
              "\n",
              "                                       parent_smiles  cleavage_position  \\\n",
              "0  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                  4   \n",
              "1  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                  5   \n",
              "2  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                  6   \n",
              "3  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                  7   \n",
              "4  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                  9   \n",
              "5  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                 10   \n",
              "6  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                 13   \n",
              "7  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                 15   \n",
              "8  CCC1OC(=O)C[C@H](O)C(C)[C@@H](O[C@@H]2OC(C)C(O...                 16   \n",
              "9  CC[C@H]1OC(=O)C[C@@H](O)[C@H](C)[C@@H](O[C@@H]...                  4   \n",
              "\n",
              "            fragment_smiles  atom_count  molecular_weight  \n",
              "0                        *O           2             17.01  \n",
              "1                        *C           2             15.04  \n",
              "2  *OC1OC(C)C(O)C(N(C)C)C1O          14            190.22  \n",
              "3         *CCN1CC(C)CC(C)C1          11            140.25  \n",
              "4                        *C           2             15.04  \n",
              "5                       *=O           2             16.00  \n",
              "6                        *C           2             15.04  \n",
              "7   *COC1OC(C)C(O)C(OC)C1OC          15            205.23  \n",
              "8                       *CC           3             29.06  \n",
              "9                        *O           2             17.01  "
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# 转换为DataFrame便于分析\n",
        "df_fragments = fragmenter.batch_to_dataframe(batch_results)\n",
        "\n",
        "print(f\"碎片DataFrame:\")\n",
        "print(f\"  形状: {df_fragments.shape}\")\n",
        "print(f\"  列: {df_fragments.columns.tolist()}\")\n",
        "print(f\"\\n前10行:\")\n",
        "df_fragments.head(10)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 6. 统计每个位置的碎片分布\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "已准备 14 个位置的碎片数据\n",
            "  位置 3: 0 个碎片\n",
            "  位置 4: 305 个碎片\n",
            "  位置 5: 307 个碎片\n",
            "  位置 6: 307 个碎片\n",
            "  位置 7: 304 个碎片\n",
            "  位置 8: 0 个碎片\n",
            "  位置 9: 307 个碎片\n",
            "  位置 10: 307 个碎片\n",
            "  位置 11: 0 个碎片\n",
            "  位置 12: 0 个碎片\n",
            "  位置 13: 263 个碎片\n",
            "  位置 14: 0 个碎片\n",
            "  位置 15: 277 个碎片\n",
            "  位置 16: 307 个碎片\n"
          ]
        }
      ],
      "source": [
        "# 准备数据：按位置收集碎片\n",
        "position_fragments = {}\n",
        "for position in range(3, 17):  # 位置3-16\n",
        "    position_data = df_fragments[df_fragments['cleavage_position'] == position]\n",
        "    position_fragments[position] = position_data.to_dict('records')\n",
        "\n",
        "print(f\"已准备 {len(position_fragments)} 个位置的碎片数据\")\n",
        "for pos, frags in position_fragments.items():\n",
        "    print(f\"  位置 {pos}: {len(frags)} 个碎片\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 7. 绘制每个位置的碎片分布图（参考test_align_two_molecules.ipynb的绘图逻辑）\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "已保存: /home/zly/project/macro_split/output/atom_count_distribution_ring16.png\n"
          ]
        },
        {
          "data": {
            "image/png": "[encoded data, 195712 characters]",
            "text/plain": [
              "<Figure size 2000x1600 with 16 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "已保存: /home/zly/project/macro_split/output/molecular_weight_distribution_ring16.png\n"
          ]
        },
        {
          "data": {
            "image/png": "[encoded data, 199308 characters]",
            "text/plain": [
              "<Figure size 2000x1600 with 16 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "所有分布图已生成完成!\n"
          ]
        }
      ],
      "source": [
        "# 准备数据\n",
        "properties_to_plot = [\n",
        "    ('atom_count', 'Atom Count', 'Number of Atoms'),\n",
        "    ('molecular_weight', 'Molecular Weight', 'Molecular Weight (Da)')\n",
        "]\n",
        "\n",
        "ring_size = 16\n",
        "positions_to_plot = list(range(3, ring_size + 1))  # 位置3-16\n",
        "\n",
        "# 为每个属性创建图表\n",
        "for prop_key, prop_title, prop_label in properties_to_plot:\n",
        "    fig, axes = plt.subplots(4, 4, figsize=(20, 16))\n",
        "    fig.suptitle(f'{prop_title} Distribution for {ring_size}-Membered Macrolactones', \n",
        "                 fontsize=16, fontweight='bold')\n",
        "    \n",
        "    axes = axes.flatten()\n",
        "    \n",
        "    for idx, position in enumerate(positions_to_plot):\n",
        "        ax = axes[idx]\n",
        "        \n",
        "        if position in position_fragments and len(position_fragments[position]) > 0:\n",
        "            # 提取该位置的所有属性值\n",
        "            values = [frag[prop_key] for frag in position_fragments[position]]\n",
        "            \n",
        "            # 绘制分布图（使用KDE平滑）\n",
        "            sns.histplot(values, kde=True, ax=ax, bins=30, color='steelblue', alpha=0.7)\n",
        "            ax.set_title(f'Position {position}\\n(n={len(values)})', fontsize=10)\n",
        "            ax.set_xlabel(prop_label, fontsize=9)\n",
        "            ax.set_ylabel('Frequency', fontsize=9)\n",
        "            ax.grid(True, alpha=0.3)\n",
        "        else:\n",
        "            ax.text(0.5, 0.5, f'Position {position}\\nNo fragments', \n",
        "                   ha='center', va='center', transform=ax.transAxes, fontsize=10)\n",
        "            ax.set_title(f'Position {position}', fontsize=10)\n",
        "    \n",
        "    # 隐藏多余的子图\n",
        "    for idx in range(len(positions_to_plot), len(axes)):\n",
        "        axes[idx].axis('off')\n",
        "    \n",
        "    plt.tight_layout()\n",
        "    \n",
        "    # 保存图片到output目录\n",
        "    output_file = project_root / 'output' / f'{prop_key}_distribution_ring{ring_size}.png'\n",
        "    plt.savefig(output_file, dpi=300, bbox_inches='tight')\n",
        "    print(f\"已保存: {output_file}\")\n",
        "    plt.show()\n",
        "\n",
        "print(\"\\n所有分布图已生成完成!\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "碎片数据已保存: /home/zly/project/macro_split/output/ring16_fragments_analysis.csv\n",
            "分子性质已保存: /home/zly/project/macro_split/output/ring16_molecular_properties.csv\n",
            "\n",
            "所有分析结果已保存!\n"
          ]
        }
      ],
      "source": [
        "# 保存碎片数据\n",
        "output_fragments_csv = project_root / 'output' / 'ring16_fragments_analysis.csv'\n",
        "df_fragments.to_csv(output_fragments_csv, index=False)\n",
        "print(f\"碎片数据已保存: {output_fragments_csv}\")\n",
        "\n",
        "# 保存分子性质数据\n",
        "output_props_csv = project_root / 'output' / 'ring16_molecular_properties.csv'\n",
        "df_props.to_csv(output_props_csv, index=False)\n",
        "print(f\"分子性质已保存: {output_props_csv}\")\n",
        "\n",
        "print(\"\\n所有分析结果已保存!\")\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 9. 延伸分析建议与总结\n",
        "\n",
        "### 已完成的分析\n",
        "\n",
        "1. **基本性质统计**\n",
        "   - ✅ 分子量、LogP、QED、TPSA等药物性质\n",
        "   - ✅ 氢键供受体数、可旋转键数等\n",
        "   - ✅ 生成了分布图并保存\n",
        "\n",
        "2. **侧链断裂分析**\n",
        "   - ✅ 每个位置（3-16）的碎片数量和类型\n",
        "   - ✅ 碎片的原子数和分子量分布\n",
        "   - ✅ 所有数据保存为CSV格式\n",
        "\n",
        "### 延伸分析建议\n",
        "\n",
        "基于当前分析，以下是值得进一步探索的方向：\n",
        "\n",
        "#### 1. **药物性质深度分析** ⭐⭐⭐（强烈推荐）\n",
        "\n",
        "**LogP分析的价值：**\n",
        "- LogP反映分子的亲脂性，对膜通透性和药物分布至关重要\n",
        "- 大环内酯通常LogP较高（从当前分析可以看到平均值和分布）\n",
        "- 建议：分析不同侧链对LogP的贡献，找出对亲脂性影响最大的位置\n",
        "\n",
        "**QED（定量估计药物相似性）分析：**\n",
        "- QED综合考虑分子量、LogP、氢键等因素，评估\"类药性\"\n",
        "- 大环内酯往往违反\"Lipinski五规则\"（分子量>500），但仍可能是好药\n",
        "- 建议：比较高QED和低QED分子的侧链差异\n",
        "\n",
        "**TPSA（极性表面积）分析：**\n",
        "- TPSA与口服生物利用度密切相关（一般<140Ų为佳）\n",
        "- 建议：分析侧链极性原子（O、N）的分布模式\n",
        "\n",
        "**代码示例：**\n",
        "```python\n",
        "# 分析LogP与侧链的关系\n",
        "# 1. 将分子按LogP分为高、中、低三组\n",
        "# 2. 比较三组在各个位置的侧链类型差异\n",
        "# 3. 找出对LogP贡献最大的位置\n",
        "\n",
        "df_props['logP_group'] = pd.cut(df_props['logP'], bins=3, labels=['Low', 'Medium', 'High'])\n",
        "```\n",
        "\n",
        "#### 2. **侧链特征分析** ⭐⭐⭐\n",
        "\n",
        "**极性侧链vs疏水侧链：**\n",
        "```python\n",
        "def classify_fragment_polarity(smiles):\n",
        "    \\\"\\\"\\\"分类碎片的极性\\\"\\\"\\\"\n",
        "    mol = Chem.MolFromSmiles(smiles)\n",
        "    if mol:\n",
        "        # 统计极性原子(O, N)\n",
        "        polar_atoms = sum(1 for atom in mol.GetAtoms() \n",
        "                         if atom.GetSymbol() in ['O', 'N'])\n",
        "        total_atoms = mol.GetNumAtoms()\n",
        "        return polar_atoms / total_atoms if total_atoms > 0 else 0\n",
        "    return 0\n",
        "\n",
        "# 为每个碎片计算极性比例\n",
        "df_fragments['polarity_ratio'] = df_fragments['fragment_smiles'].apply(\n",
        "    classify_fragment_polarity\n",
        ")\n",
        "\n",
        "# 分析每个位置的极性偏好\n",
        "position_polarity = df_fragments.groupby('cleavage_position')['polarity_ratio'].mean()\n",
        "```\n",
        "\n",
        "#### 3. **结构-活性关系（SAR）分析** ⭐⭐⭐\n",
        "\n",
        "如果您的数据中有生物活性数据（如`max_pChEMBL`）：\n",
        "\n",
        "```python\n",
        "# 读取活性数据\n",
        "df_with_activity = pd.merge(\n",
        "    df_props, \n",
        "    df[['unique_id', 'max_pChEMBL']], \n",
        "    on='unique_id'\n",
        ")\n",
        "\n",
        "# 分析高活性分子的侧链特征\n",
        "high_activity = df_with_activity[df_with_activity['max_pChEMBL'] > 7]\n",
        "low_activity = df_with_activity[df_with_activity['max_pChEMBL'] < 5]\n",
        "\n",
        "# 比较两组在各个位置的碎片分布差异\n",
        "# 找出\"特权侧链\"（privilege side chains）\n",
        "```\n",
        "\n",
        "#### 4. **碎片多样性分析** ⭐⭐\n",
        "\n",
        "```python\n",
        "# 统计每个位置的独特碎片类型\n",
        "unique_fragments_per_position = {}\n",
        "for position in range(3, 17):\n",
        "    position_data = df_fragments[df_fragments['cleavage_position'] == position]\n",
        "    unique_smiles = position_data['fragment_smiles'].unique()\n",
        "    unique_fragments_per_position[position] = {\n",
        "        'total': len(position_data),\n",
        "        'unique': len(unique_smiles),\n",
        "        'diversity_ratio': len(unique_smiles) / len(position_data)\n",
        "    }\n",
        "```\n",
        "\n",
        "#### 5. **聚类分析** ⭐⭐\n",
        "\n",
        "```python\n",
        "from rdkit.Chem import AllChem\n",
        "from sklearn.cluster import KMeans\n",
        "\n",
        "# 基于碎片指纹进行分子聚类\n",
        "def get_fragment_fingerprint(mol_result):\n",
        "    \\\"\\\"\\\"为每个分子生成碎片指纹\\\"\\\"\\\"\n",
        "    # 基于侧链组成创建特征向量\n",
        "    fp = [0] * 14  # 位置3-16\n",
        "    for frag in mol_result.fragments:\n",
        "        pos = frag.cleavage_position\n",
        "        if 3 <= pos <= 16:\n",
        "            fp[pos-3] = frag.molecular_weight\n",
        "    return fp\n",
        "\n",
        "# 聚类找出具有相似侧链模式的分子家族\n",
        "```\n",
        "\n",
        "#### 6. **3D性质分析** ⭐（可选）\n",
        "\n",
        "如果需要更深入的3D分析：\n",
        "```python\n",
        "from rdkit.Chem import AllChem, Descriptors3D\n",
        "\n",
        "# 生成3D构象\n",
        "mol = Chem.MolFromSmiles(smiles)\n",
        "AllChem.EmbedMolecule(mol)\n",
        "AllChem.MMFFOptimizeMolecule(mol)\n",
        "\n",
        "# 计算3D描述符\n",
        "pmi1 = Descriptors3D.PMI1(mol)  # 主惯性矩\n",
        "npr1 = Descriptors3D.NPR1(mol)  # 归一化主惯性比\n",
        "```\n",
        "\n",
        "### 推荐的下一步优先级\n",
        "\n",
        "1. **优先级1（必做）：LogP、QED和TPSA与侧链的关系分析** ⭐⭐⭐\n",
        "   - 对药物设计最直接相关\n",
        "   - 实现简单，计算快速\n",
        "   - 能提供实用的设计指导\n",
        "\n",
        "2. **优先级2（重要）：如果有活性数据，进行SAR分析** ⭐⭐⭐\n",
        "   - 找出活性关键侧链\n",
        "   - 指导后续合成策略\n",
        "\n",
        "3. **优先级3（有价值）：碎片多样性和聚类分析** ⭐⭐\n",
        "   - 发现侧链组合规律\n",
        "   - 可能发现新的设计策略\n",
        "\n",
        "### 其他潜在分析指标\n",
        "\n",
        "**Lipinski规则检查：**\n",
        "```python\n",
        "def check_lipinski(props):\n",
        "    violations = 0\n",
        "    if props['mol_weight'] > 500: violations += 1\n",
        "    if props['logP'] > 5: violations += 1\n",
        "    if props['num_h_donors'] > 5: violations += 1\n",
        "    if props['num_h_acceptors'] > 10: violations += 1\n",
        "    return violations\n",
        "\n",
        "df_props['lipinski_violations'] = df_props.apply(check_lipinski, axis=1)\n",
        "```\n",
        "\n",
        "**Veber规则（口服药物性）：**\n",
        "- TPSA ≤ 140 Ų\n",
        "- 可旋转键数 ≤ 10\n",
        "\n",
        "**立体化学分析：**\n",
        "- 统计手性中心数量\n",
        "- 分析每个位置的立体化学偏好\n",
        "\n",
        "### 总结\n",
        "\n",
        "**已验证的重要指标：**\n",
        "- ✅ LogP：反映亲脂性，影响膜通透性\n",
        "- ✅ QED：综合类药性评分\n",
        "- ✅ TPSA：极性表面积，影响生物利用度\n",
        "\n",
        "**下一步建议：**\n",
        "1. 分析LogP/QED与侧链的关系（最重要）\n",
        "2. 如果有活性数据，进行SAR分析（非常有价值）\n",
        "3. 碎片多样性分析（了解侧链变化空间）\n",
        "4. 聚类分析（发现分子家族）\n",
        "\n",
        "**关键洞察：**\n",
        "大环内酯作为天然产物来源的药物，往往不完全遵循Lipinski规则，但仍能表现出良好的药物性质。\n",
        "理解侧链对这些性质的贡献，可以帮助：\n",
        "- 优化现有分子的成药性\n",
        "- 设计新的类似物\n",
        "- 预测哪些位置的修饰最有价值\n",
        "\n",
        "希望这些分析建议对您的研究有帮助！\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "\n",
            "================================================================================\n",
            "分析完成！\n",
            "================================================================================\n",
            "\n",
            "处理了 307 个分子\n",
            "生成了 2684 个碎片\n",
            "\n",
            "所有结果已保存到: /home/zly/project/macro_split/output\n",
            "\n",
            "生成的文件:\n",
            "  - ring16_molecular_properties_distribution.png (分子性质分布图)\n",
            "  - atom_count_distribution_ring16.png (原子数分布图)\n",
            "  - molecular_weight_distribution_ring16.png (分子量分布图)\n",
            "  - ring16_fragments_analysis.csv (碎片数据)\n",
            "  - ring16_molecular_properties.csv (分子性质数据)\n"
          ]
        }
      ],
      "source": [
        "# 完成！\n",
        "print(\"\\n\" + \"=\"*80)\n",
        "print(\"分析完成！\")\n",
        "print(\"=\"*80)\n",
        "print(f\"\\n处理了 {len(batch_results)} 个分子\")\n",
        "print(f\"生成了 {total_fragments} 个碎片\")\n",
        "print(f\"\\n所有结果已保存到: {project_root / 'output'}\")\n",
        "print(\"\\n生成的文件:\")\n",
        "print(\"  - ring16_molecular_properties_distribution.png (分子性质分布图)\")\n",
        "print(\"  - atom_count_distribution_ring16.png (原子数分布图)\")\n",
        "print(\"  - molecular_weight_distribution_ring16.png (分子量分布图)\")\n",
        "print(\"  - ring16_fragments_analysis.csv (碎片数据)\")\n",
        "print(\"  - ring16_molecular_properties.csv (分子性质数据)\")\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.14.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}
