diff --git a/scripts/analyze_validation_fragment_library.py b/scripts/analyze_validation_fragment_library.py index 98de199..596ad48 100644 --- a/scripts/analyze_validation_fragment_library.py +++ b/scripts/analyze_validation_fragment_library.py @@ -3,6 +3,7 @@ from __future__ import annotations import argparse from math import ceil from pathlib import Path +import sqlite3 import matplotlib @@ -448,6 +449,48 @@ def format_position_mapping(positions: list[int], ring_size: int) -> str: return ", ".join(f"{position} → {mirror_ring_position(position, ring_size)}" for position in positions) +def build_zero_fragment_parent_table(db_path: str | Path, ring_size: int) -> pd.DataFrame: + with sqlite3.connect(db_path) as connection: + tables = { + row[0] + for row in connection.execute( + "SELECT name FROM sqlite_master WHERE type='table'" + ) + } + if "fragment_library_entries" in tables: + query = """ + SELECT ml_id, num_sidechains, cleavage_positions + FROM parent_molecules + WHERE classification = 'standard_macrolactone' + AND ring_size = ? + AND processing_status = 'success' + AND NOT EXISTS ( + SELECT 1 + FROM fragment_library_entries fle + WHERE fle.source_parent_ml_id = parent_molecules.ml_id + AND fle.source_type = 'validation_extract' + AND fle.splice_ready = 1 + ) + ORDER BY ml_id + """ + else: + # The lightweight test database only includes parent_molecules. + query = """ + SELECT ml_id, num_sidechains, cleavage_positions + FROM parent_molecules + WHERE classification = 'standard_macrolactone' + AND ring_size = ? + AND processing_status = 'success' + AND num_sidechains = 0 + ORDER BY ml_id + """ + return pd.read_sql_query( + query, + connection, + params=[ring_size], + ) + + def build_markdown_report( output_dir: Path, analysis_df: pd.DataFrame, @@ -655,6 +698,8 @@ def build_markdown_report_zh( ring_df: pd.DataFrame, filtered_ring_df: pd.DataFrame, filter_candidates: pd.DataFrame, + standard_success_parent_count: int, + zero_fragment_ring_parents: pd.DataFrame, diversity_gt3: pd.DataFrame, position_counts: pd.DataFrame, ring_sensitivity_table: pd.DataFrame, @@ -688,8 +733,20 @@ def build_markdown_report_zh( "", f"- 当前验证后的可拼接碎片库包含 **{len(analysis_df):,}** 条片段记录，来源于 **{analysis_df['source_parent_ml_id'].nunique():,}** 个母体分子。", f"- 其中 {ring_size} 元环子集包含 **{len(ring_df):,}** 条片段记录，来源于 **{ring_df['source_parent_ml_id'].nunique():,}** 个母体分子。", + f"- 这里的“可拼接碎片”指验证库中 `source_type='validation_extract'` 且 `splice_ready=1` 的单锚点侧链片段。桥环、稠环或任何具有多个环连接点的侧链都已经在生成阶段被排除，不会进入这份库。", + f"- 16 元环子集是按母体元数据里的 `ring_size=16` 直接过滤出来的，不是从片段反推 ring size。", f"- 用于设计相关位点分析的严格子集定义为：片段重原子数 **>= {design_min_atoms}**。", "", + "## 16 元环母体计数口径说明", + "", + f"- 在数据库里，`standard_macrolactone + ring_size={ring_size} + processing_status=success` 一共有 **{standard_success_parent_count:,}** 个母体。", + f"- 其中 **{ring_df['source_parent_ml_id'].nunique():,}** 个母体至少产出过 1 条可拼接片段，所以进入了当前报告的 16 元环片段统计。", + f"- 剩余 **{len(zero_fragment_ring_parents):,}** 个母体没有任何可拼接片段；它们的共同特征是 `num_sidechains=0`、`cleavage_positions=[]`。", + "", + "```text", + zero_fragment_ring_parents.to_string(index=False) if not zero_fragment_ring_parents.empty else "No zero-fragment parents.", + "```", + "", "## 全库碎片大小结论", "", f"- 默认清洗阈值建议使用 `<= {design_min_atoms - 2}` 重原子删除。该阈值会删除 **{int(conservative_filter.removed_rows):,}** 条记录（{conservative_filter.removed_row_fraction:.1%}），但仅删除 **{int(conservative_filter.removed_unique_fragments):,}** 个唯一片段（{conservative_filter.removed_unique_fraction:.1%}）。", @@ -730,7 +787,7 @@ def build_markdown_report_zh( "", "## 桥环 / 稠环干扰的敏感性分析", "", - "桥连或双锚点侧链不会进入当前片段库，因为断裂逻辑只保留与主环存在 **1 个连接点** 的侧链组件。也就是说，真正的 bridge / fused multi-anchor components 已被代码层面排除。", + "桥连或双锚点侧链不会进入当前片段库，因为断裂逻辑只保留与主环存在 **1 个连接点** 的侧链组件。也就是说，真正的 bridge / fused multi-anchor components 已被代码层面排除。上面那 6 个 16 元环母体并不是这类“被误收进来的桥环碎片”，而是根本没有任何可拼接外侧链，所以不会产生 fragment 行。", "", "但是，需要额外区分另一类情况：**cyclic single-anchor side chains**。这类片段虽然只在一个位置连到主环，因此会被保留下来，但片段自身可能包含糖环、杂环或其他环状骨架，仍然会显著影响位点多样性排名。", "", @@ -905,6 +962,21 @@ def main(argv: list[str] | None = None) -> None: ].copy() ring_df = analysis_df[analysis_df["ring_size"] == args.ring_size].copy() filtered_ring_df = ring_df[ring_df["fragment_atom_count"] >= args.design_min_atoms].copy() + zero_fragment_ring_parents = build_zero_fragment_parent_table(args.db, args.ring_size) + with sqlite3.connect(args.db) as connection: + standard_success_parent_count = int( + pd.read_sql_query( + """ + SELECT COUNT(*) AS count + FROM parent_molecules + WHERE classification = 'standard_macrolactone' + AND ring_size = ? + AND processing_status = 'success' + """, + connection, + params=[args.ring_size], + ).iloc[0]["count"] + ) if analysis_df.empty: raise ValueError("No splice-ready standard macrolactone fragments available for analysis.") @@ -1003,6 +1075,8 @@ def main(argv: list[str] | None = None) -> None: ring_df, filtered_ring_df, filter_candidates, + standard_success_parent_count, + zero_fragment_ring_parents, diversity_gt3, position_counts, ring_sensitivity, diff --git a/validation_output/fragment_library_analysis/fragment_library_analysis_report_zh.md b/validation_output/fragment_library_analysis/fragment_library_analysis_report_zh.md index 9c4e383..f4b2071 100644 --- a/validation_output/fragment_library_analysis/fragment_library_analysis_report_zh.md +++ b/validation_output/fragment_library_analysis/fragment_library_analysis_report_zh.md @@ -1,11 +1,205 @@ # 大环内酯碎片库分析报告（中文） +## 数据筛选流程图 + +下面这张图把本报告里的数据是如何一步步筛出来的串起来。它分成两段： + +- 第一段是 `validation` 阶段，负责把原始 MacrolactoneDB 母体变成 `fragment_library.csv` +- 第二段是 `fragment_library_analysis` 阶段，负责从 `fragment_library.csv` 里统计出本报告里的 34,829、4,451、8,108、1,105 等数字 + +```mermaid +flowchart TD + A["原始输入 CSV
/Users/lingyuzeng/project/macro_split/data/MacrolactoneDB/ring12_20/temp.csv"] --> B["stratified_sample_by_ring_size()
按 12-20 元环分层抽样 10%"] + B --> C["classify_macrocycle()
只保留 standard_macrolactone"] + C -->|non_standard_macrocycle / not_macrolactone| C1["processing_status = skipped
不进入裂解"] + C -->|standard_macrolactone| D["find_macrolactone_candidates()
确认主环候选"] + D --> E["number_macrolactone()
建立 canonical numbering"] + E --> F["遍历 position > 2 的环位点"] + F --> G["跳过 ring atom
跳过 intrinsic lactone neighbor"] + G --> H["collect_fragmentable_side_chain_atoms()
只保留可拼接单锚点侧链"] + H -->|None| H1["丢弃该邻居/该裂解尝试"] + H -->|atoms returned| I["build_fragment_with_isotope()
生成同位素标记碎片"] + I --> J["Chem.MolFromSmiles(plain_smiles)
SMILES 必须有效"] + J -->|invalid| J1["丢弃该碎片"] + J -->|valid| K["写入 SideChainFragment"] + K --> L["写入 FragmentLibraryEntry
source_type='validation_extract'
splice_ready=True"] + L --> M["parent.processing_status = success
num_sidechains += 1"] + M --> N["导出 fragment_library.csv / summary.csv / fragments.db"] + N --> O["load_fragment_library_dataset()
把 fragment_library.csv 与 parent_molecules 合并"] + O --> P["annotate_fragment_atom_counts()
按 plain SMILES 计算重原子数"] + P --> Q["build_fragment_atom_count_summary()
得到 34,829 / 4,451 / 1,852"] + P --> R["build_position_diversity_table()
得到位点统计表"] + P --> S["ring_size=16 过滤
得到 8,108 条片段、1,105 个母体"] + S --> T["> 3 重原子过滤
得到设计相关子集与位点多样性结论"] +``` + +这张流程图对应的是一条从“原始母体”到“最终碎片统计”的两段式筛选链路。这里的“原始输入 CSV”指的是传入 `validator.run(input_csv)` 的 MacrolactoneDB ring12_20 数据集文件，在当前仓库里的默认入口脚本 [scripts/validate_macrolactone_db.py](/Users/lingyuzeng/project/macro_split/scripts/validate_macrolactone_db.py) 中，默认值就是 `data/MacrolactoneDB/ring12_20/temp.csv`，对应的绝对路径是 `/Users/lingyuzeng/project/macro_split/data/MacrolactoneDB/ring12_20/temp.csv`。第一段发生在 `validation` 阶段：程序先读取这个输入 CSV，对所有分子做 `classify_macrocycle()` 预分类，并按 12 到 20 元环分层抽样。这里的“抽 10%”不是对整个数据集一次性随机抽样，而是先按 `_ring_size` 把分子分到 12、13、14、15、16、17、18、19、20 这九个层级，再对每个层级分别抽取 `max(1, int(len(group) * sample_ratio))` 条记录；`sample_ratio` 的默认值是 `0.1`，`random_state` 的默认值是 `42`。抽样之后，只有被识别为 `standard_macrolactone` 的分子才会继续处理；`non_standard_macrocycle` 和 `not_macrolactone` 会直接被标记为 `skipped`，不会进入裂解步骤。对标准大环内酯，程序先用 `find_macrolactone_candidates()` 确认主环候选，再用 `number_macrolactone()` 建立 canonical numbering，之后只遍历 `position > 2` 的环位点。对于每个位点，代码会跳过环内原子和 intrinsic lactone neighbor，只对真正位于外侧链上的邻居进行处理；随后通过 `collect_fragmentable_side_chain_atoms()` 判断这段侧链是否满足“单锚点、可拼接”的要求。如果该函数返回 `None`，说明这条侧链不符合拼接条件，整次裂解尝试就会被丢弃。若侧链可拼接，则继续用 `build_fragment_with_isotope()` 生成同位素标记碎片，并用 `Chem.MolFromSmiles(plain_smiles)` 做一次有效性检查；只有 SMILES 能被正确解析的碎片才会被写入 `SideChainFragment` 和 `FragmentLibraryEntry`。在写入统一碎片库时，这些记录会被固定标记为 `source_type='validation_extract'` 和 `splice_ready=True`，表示它们是后续拼接设计可以直接使用的碎片。母体记录则会被回写为 `processing_status = success`，并记录该分子实际产生了多少条侧链碎片以及对应的裂解位点。 + +第二段发生在 `fragment_library_analysis` 阶段：分析脚本只读取已经生成好的 `fragment_library.csv`，再从 SQLite 数据库中读取 `parent_molecules` 表，把 `source_parent_ml_id` 和 `ml_id` 对齐后合并母体元数据。这个合并步骤的作用是把每条碎片重新挂回到它来自哪个母体、属于哪种分类、对应哪个环大小以及母体是否成功处理等信息上；如果有任何碎片在数据库里找不到对应母体，分析会直接报错，而不会静默继续。合并完成后，脚本会根据 `fragment_smiles_plain` 计算重原子数，形成后续所有统计的基础字段。全库的 `34,829` 条记录和 `4,451` 个去重母体，就是在这一步对 `fragment_library_entries` 做总量统计得到的；随后再按 `ring_size=16` 过滤，就得到 16 元环子集的 `8,108` 条碎片记录和 `1,105` 个母体。接下来如果进一步应用 `>= 4` 重原子的设计相关过滤，就会得到用于位点多样性分析的严格子集；`build_position_diversity_table()` 就是在这个子集上按 `cleavage_position` 分组，统计每个位点的总碎片数、唯一碎片数、香农熵、Tanimoto 距离和原子数分布，从而得出本报告里 16 元环位点排序和热点结论。 + +### 标准大环内酯的识别规则 + +这里再单独把“标准环怎么判”的逻辑说清楚，因为这是整个筛选链路的根。 + +- 入口是 [analyzer.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/analyzer.py#L19-L26) 调用 [\_core.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/_core.py#L89-L137) 里的 `classify_macrolactone()` +- `find_macrolactone_candidates()` 会遍历 `mol.GetRingInfo().AtomRings()`，只保留环长在 `12..20` 的候选环 +- 每个候选环都必须能找到一个乳内酯特征： + - 1 个羰基碳 + - 1 个酯氧 + - 能组成 lactone ring 的闭环 +- 如果找不到任何 12-20 元 lactone ring，就会被判为 `not_macrolactone` +- 如果找到多个彼此重叠的候选环，就会被判为 `non_standard_macrocycle`，这是桥环、稠环或多环重叠结构的主要排除门槛 +- 如果候选环唯一，再调用 `build_numbering_result()` 建立编号，随后检查 `3..N` 位是否存在非碳原子： + - 只要环上除 1、2 位以外出现了 O、N 等非 C 原子，就会被判为 `non_standard_macrocycle` + - 对应的原因码是 `contains_non_carbon_ring_atoms_outside_positions_1_2` +- 只有“候选环唯一 + 3..N 全为碳 + 环长在 12-20”同时成立时，才会被判为 `standard_macrolactone` + +这也解释了你提到的两类结构： + +- **环肽/非标准大环**：如果环上 3..N 出现 N、O 等非碳原子，虽然可能也是大环，但不会进入标准大环内酯集合 +- **桥环/稠环**：如果一个分子里存在多个重叠的 lactone 候选环，程序会直接把它归为 `non_standard_macrocycle`，不会进入后续裂解 + +### 口径对照 + +- `34,829` 是最终 `fragment_library_entries` 中的总行数，不是原始输入总分子数 +- `4,451` 是这 34,829 条片段对应的去重母体数 +- `4,482` 是 `summary.csv` / 数据库里所有 `standard_macrolactone + success` 母体数 +- 两者之差 `31` 个母体，是已经成功识别为标准大环内酯、但 `num_sidechains=0`，因此没有产生任何可拼接片段的分子 +- `1,111` 是 `standard_macrolactone + ring_size=16 + success` 母体数 +- `1,105` 是其中至少产出 1 条可拼接片段的 16 元环母体数 +- 两者之差 `6` 个母体，也是没有任何可拼接外侧链的 16 元环分子 + +## 分步筛选说明（对应代码） + +### 1. 输入与分层抽样 + +- 入口在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L38-L67) 的 `run()` +- 原始输入先读成 `DataFrame` +- 然后调用 [sampling.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/sampling.py#L8-L60) 里的 `stratified_sample_by_ring_size(df, self.sample_ratio, self.smiles_col)` +- 这一步先对所有分子做 `classify_macrocycle()` 预分类，再按 `_ring_size` 在 `12..20` 各层分别抽样 +- 默认 `sample_ratio=0.1`，`random_state=42` +- 这一步的意义是：按 12-20 元环大小分层抽样，而不是对整个库做一次性随机抽样 + +### 2. 母体分类 + +- 入口在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L72-L146) 的 `_process_molecule()` +- 用 `self.analyzer.classify_macrocycle(smiles)` 把分子分成三类： + - `standard_macrolactone` + - `non_standard_macrocycle` + - `not_macrolactone` +- 只有 `standard_macrolactone` 会进入后续裂解 +- 其他两类直接设为 `processing_status = skipped` + +### 3. 标准大环内酯的候选确认与编号 + +- 入口在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L148-L186) +- 先调用 `find_macrolactone_candidates(mol, ring_size=parent.ring_size)` 找主环候选 +- 如果候选为空，或同一分子里出现多个无法唯一确定的候选，会报错或终止 +- 然后调用 `number_macrolactone(mol, ring_size=parent.ring_size)` 建立 canonical numbering +- 这里的编号规则就是本仓库统一规则：`1 = 内酯羰基碳`，`2 = 相邻酯氧`，`3..N` 顺着环唯一遍历 + +### 4. 逐位点裂解 + +- 入口仍在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L187-L264) +- 循环只看 `position > 2` 的环位点，`1` 和 `2` 不参与裂解 +- 对每个位点： + - 先取该位点的环原子 + - 再遍历其邻居 + - 若邻居还是环内原子，则跳过 + - 若是 intrinsic lactone neighbor，也跳过 +- 剩下的邻居才进入 `collect_fragmentable_side_chain_atoms(...)` + +### 5. 侧链可拼接性过滤 + +- 这一步是最关键的筛选门槛，代码仍在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L206-L214) +- `collect_fragmentable_side_chain_atoms(...)` 必须返回一个有效的原子集合 +- 如果返回 `None`，说明该侧链不满足“可拼接单锚点侧链”的条件，直接丢弃 +- 这也是桥环、稠环、双锚点或其他不可拼接外侧链被排除的核心位置 + +### 6. 碎片构建与有效性检查 + +- 入口在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L216-L228) +- 用 `build_fragment_with_isotope(...)` 生成： + - `fragment_smiles_labeled` + - `fragment_smiles_plain` + - `original_bond_type` +- 接着用 `Chem.MolFromSmiles(plain_smiles)` 再做一次有效性检查 +- 如果 plain SMILES 无法解析，这条碎片也会被丢弃 + +### 7. 写入片段库 + +- 入口在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L230-L264) +- 每条通过筛选的碎片都会写入两个表： + - `SideChainFragment` + - `FragmentLibraryEntry` +- 其中 `FragmentLibraryEntry` 固定写入： + - `source_type='validation_extract'` + - `splice_ready=True` + - `dummy_atom_count=1` +- 也就是说，进入本报告的都是已经被判定为“可直接用于单锚点拼接”的碎片 + +### 8. 母体状态回写 + +- 入口在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L273-L279) +- 处理完成后，母体会被标记为 `processing_status = success` +- 同时写回： + - `num_sidechains = len(fragments)` + - `cleavage_positions = [...]` +- 如果标准大环内酯没有任何可拼接侧链，母体仍然会是 `success`，但 `num_sidechains=0` + +### 9. 统一导出 + +- 入口在 [validator.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/validator.py#L280-L383) +- 这一步导出三类关键文件： + - `summary.csv` + - `summary_statistics.json` + - `fragment_library.csv` +- 其中本报告使用的核心输入是 `fragment_library.csv` 和数据库里的 `parent_molecules` + +### 10. 报告统计 + +- 入口在 [fragment_library_analysis.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/fragment_library_analysis.py#L44-L69) +- `load_fragment_library_dataset()` 会： + - 读取 `fragment_library.csv` + - 从 SQLite 里读取 `parent_molecules` + - 按 `source_parent_ml_id = ml_id` 合并 + - 校验母体元数据不能缺失 +- 然后 `annotate_fragment_atom_counts()` 会按 `fragment_smiles_plain` 计算重原子数 + +### 11. 汇总输出 + +- 入口在 [fragment_library_analysis.py](/Users/lingyuzeng/project/macro_split/src/macro_lactone_toolkit/validation/fragment_library_analysis.py#L72-L183) +- `build_fragment_atom_count_summary()` 生成全库统计，得到本报告中的： + - `rows = 34,829` + - `unique_parent_molecules = 4,451` + - `unique_fragment_smiles = 1,852` +- `build_filter_candidate_table()` 生成阈值删除表，因此报告里能看到 `<= 2`、`<= 3`、`<= 4` 等过滤候选 +- `build_position_diversity_table()` 生成位点多样性表，因此报告里能看到 16 元环的位点排序与 `> 3` 重原子子集结论 + ## 数据范围 - 当前验证后的可拼接碎片库包含 **34,829** 条片段记录，来源于 **4,451** 个母体分子。 - 其中 16 元环子集包含 **8,108** 条片段记录，来源于 **1,105** 个母体分子。 +- 这里的“可拼接碎片”指验证库中 `source_type='validation_extract'` 且 `splice_ready=1` 的单锚点侧链片段。桥环、稠环或任何具有多个环连接点的侧链都已经在生成阶段被排除，不会进入这份库。 +- 16 元环子集是按母体元数据里的 `ring_size=16` 直接过滤出来的，不是从片段反推 ring size。 - 用于设计相关位点分析的严格子集定义为：片段重原子数 **>= 4**。 +## 16 元环母体计数口径说明 + +- 在数据库里，`standard_macrolactone + ring_size=16 + processing_status=success` 一共有 **1,111** 个母体。 +- 其中 **1,105** 个母体至少产出过 1 条可拼接片段，所以进入了当前报告的 16 元环片段统计。 +- 剩余 **6** 个母体没有任何可拼接片段；它们的共同特征是 `num_sidechains=0`、`cleavage_positions=[]`。 + +```text + ml_id num_sidechains cleavage_positions +ML00006860 0 [] +ML00007029 0 [] +ML00007030 0 [] +ML00007031 0 [] +ML00007032 0 [] +ML00008015 0 [] +``` + ## 全库碎片大小结论 - 默认清洗阈值建议使用 `<= 2` 重原子删除。该阈值会删除 **28,069** 条记录（80.6%），但仅删除 **26** 个唯一片段（1.4%）。 @@ -41,7 +235,7 @@ ## 桥环 / 稠环干扰的敏感性分析 -桥连或双锚点侧链不会进入当前片段库，因为断裂逻辑只保留与主环存在 **1 个连接点** 的侧链组件。也就是说，真正的 bridge / fused multi-anchor components 已被代码层面排除。 +桥连或双锚点侧链不会进入当前片段库，因为断裂逻辑只保留与主环存在 **1 个连接点** 的侧链组件。也就是说，真正的 bridge / fused multi-anchor components 已被代码层面排除。上面那 6 个 16 元环母体并不是这类“被误收进来的桥环碎片”，而是根本没有任何可拼接外侧链，所以不会产生 fragment 行。 但是，需要额外区分另一类情况：**cyclic single-anchor side chains**。这类片段虽然只在一个位置连到主环，因此会被保留下来，但片段自身可能包含糖环、杂环或其他环状骨架，仍然会显著影响位点多样性排名。 @@ -110,4 +304,3 @@ - 16 元环位点多样性图：`ring16_position_diversity_gt3.png` - 16 元环桥环/带环侧链敏感性图：`ring16_position_ring_sensitivity.png` - 16 元环药化热点对比图：`ring16_medchem_hotspot_comparison.png` - diff --git a/validation_output/macrolactone_validation_fragmentation_enumeration_report_zh.md b/validation_output/macrolactone_validation_fragmentation_enumeration_report_zh.md new file mode 100644 index 0000000..682dd15 --- /dev/null +++ b/validation_output/macrolactone_validation_fragmentation_enumeration_report_zh.md @@ -0,0 +1,380 @@ +# MacrolactoneDB 验证、碎片库分析与 Tylosin 枚举报告 + +本文档把 `temp.csv` 的验证筛选、侧链碎片库生成、碎片大小分析、16 元环位点提取，以及 tylosin 枚举结果串成一条完整数据链路。目标是让没有上下文的读者也能直接理解： + +1. 原始数据如何从 11,036 个分子收缩到标准大环内酯集合。 +2. 为什么最终碎片库只有 4,451 个母体进入可拼接片段分析。 +3. 为什么后续侧链筛选采用 `>3` 重原子作为下限。 +4. 16 元环的 `3/4/12/13` 位点碎片是如何导出并用于 tylosin 枚举的。 + +## 结果文件一览 + +| 文件 | 作用 | +|---|---| +| [validation_output/summary.csv](/Users/lingyuzeng/project/macro_split/validation_output/summary.csv) | 逐分子验证结果与 `processing_status` | +| [validation_output/summary_statistics.json](/Users/lingyuzeng/project/macro_split/validation_output/summary_statistics.json) | 验证统计汇总 | +| [validation_output/fragments.db](/Users/lingyuzeng/project/macro_split/validation_output/fragments.db) | 验证阶段 SQLite 数据库 | +| [validation_output/fragment_library.csv](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library.csv) | 统一碎片库导出 | +| [validation_output/fragment_library_analysis/fragment_atom_count_summary.csv](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/fragment_atom_count_summary.csv) | 碎片原子数概览，含最大值 | +| [validation_output/fragment_library_analysis/fragment_atom_count_frequency.csv](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/fragment_atom_count_frequency.csv) | 原子数频率分布 | +| [validation_output/fragment_library_analysis/fragment_atom_count_filter_candidates.csv](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/fragment_atom_count_filter_candidates.csv) | 不同下限阈值的删减效果 | +| [validation_output/fragment_library_analysis/analysis_summary.txt](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/analysis_summary.txt) | 文字版分析摘要 | +| [validation_output/ring16_position_fragment_exports/pos3_fragments_dedup.csv](/Users/lingyuzeng/project/macro_split/validation_output/ring16_position_fragment_exports/pos3_fragments_dedup.csv) | 16 元环 3 位碎片库 | +| [validation_output/ring16_position_fragment_exports/pos4_fragments_dedup.csv](/Users/lingyuzeng/project/macro_split/validation_output/ring16_position_fragment_exports/pos4_fragments_dedup.csv) | 16 元环 4 位碎片库 | +| [validation_output/ring16_position_fragment_exports/pos12_fragments_dedup.csv](/Users/lingyuzeng/project/macro_split/validation_output/ring16_position_fragment_exports/pos12_fragments_dedup.csv) | 16 元环 12 位碎片库 | +| [validation_output/ring16_position_fragment_exports/pos13_fragments_dedup.csv](/Users/lingyuzeng/project/macro_split/validation_output/ring16_position_fragment_exports/pos13_fragments_dedup.csv) | 16 元环 13 位碎片库 | +| [validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/tylosin_scheme_b_unique_products.csv](/Users/lingyuzeng/project/macro_split/validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/tylosin_scheme_b_unique_products.csv) | tylosin 枚举唯一产物 | +| [validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/tylosin_scheme_b_fix_pos13.sqlite](/Users/lingyuzeng/project/macro_split/validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/tylosin_scheme_b_fix_pos13.sqlite) | tylosin 枚举 provenance 数据库 | + +## 数据流总览 + +```mermaid +flowchart TD + A["data/MacrolactoneDB/ring12_20/temp.csv
11,036 个分子"] --> B["classify_macrocycle()
standard / non-standard / not"] + B -->|standard_macrolactone| C["逐个标准大环内酯做 canonical numbering"] + B -->|non_standard_macrocycle / not_macrolactone| D["processing_status = skipped"] + C --> E["遍历 position > 2 的环位点"] + E --> F["collect_fragmentable_side_chain_atoms()
只保留单锚点可拼接侧链"] + F -->|None| G["丢弃该裂解尝试"] + F -->|atoms returned| H["build_fragment_with_isotope()
生成 labeled / plain SMILES"] + H --> I["写入 side_chain_fragments"] + H --> J["写入 fragment_library_entries
source_type = validation_extract"] + J --> K["fragment_library.csv / fragments.db"] + K --> L["按 plain SMILES 计算重原子数"] + L --> M["全库阈值分析 + 位点多样性分析"] + M --> N["导出 ring16_position_fragment_exports/"] + N --> O["tylosin scheme_b_fix_pos13 枚举"] +``` + +## 1. 术语与编号规则 + +| 术语 | 含义 | +|---|---| +| `canonical numbering` | `1 = 内酯羰基碳`，`2 = 相邻酯氧`，`3..N` 按从 2 位出发沿环唯一遍历顺序编号 | +| `mirror mapping` | 16 元环固定镜像关系，`3→16`、`4→15`、`5→14`、`6→13`、`7→12`、`8→11`、`9→10` | +| `standard_macrolactone` | 只有一个可唯一确认的 12-20 元 lactone 环，且 `3..N` 全部为碳 | +| `non_standard_macrocycle` | 有重叠候选环，或者 `3..N` 位置出现非碳原子，属于桥环 / 稠环 / 非标准大环 | +| `not_macrolactone` | 找不到有效的 12-20 元 lactone 环 | +| `cleavage_position` | 侧链连接到环上的位置编号 | +| `splice_ready` | 该碎片是单锚点、可直接用于重新拼接的碎片 | +| `fragment_smiles_plain` | 去掉 isotope 但保留 dummy 原子 `*` 的碎片 SMILES | +| `fragment_smiles_labeled` | 带 isotope 的碎片 SMILES，例如 `[13*]...` | + +这里特别强调一点：碎片去重时考虑了 dummy 原子 `*`，但**不**把 dummy 上的 isotope 作为唯一性依据。也就是说，`*C`、`*O` 这类结构会按 `plain SMILES` 去重，但 `[`13*`]...` 这类位置标签不会把同一 plain 结构拆成多个独立 chemotype。 + +## 2. 从 `temp.csv` 到标准大环内酯 + +原始输入文件是 [data/MacrolactoneDB/ring12_20/temp.csv](/Users/lingyuzeng/project/macro_split/data/MacrolactoneDB/ring12_20/temp.csv)，共 **11,036** 个分子。 + +### 2.1 分类结果 + +| 分类 | 数量 | 占总输入比例 | +|---|---:|---:| +| `non_standard_macrocycle` | 6,336 | 57.41% | +| `standard_macrolactone` | 4,482 | 40.61% | +| `not_macrolactone` | 218 | 1.98% | + +这一步的关键不是“有没有大环”本身，而是“是否满足标准大环内酯定义”： + +1. 必须存在 12-20 元 lactone 候选环。 +2. 候选环必须唯一，不能是多个重叠候选。 +3. 环上 `3..N` 必须全部是碳。 + +因此，很多结构虽然是宏环，但会因为桥环、稠环、杂原子插入或候选不唯一而被排除在 `standard_macrolactone` 之外。 + +### 2.2 标准大环内酯的 ring size 分布 + +在 **4,482** 个标准大环内酯中，ring size 分布如下： + +| ring size | 数量 | 占标准集比例 | +|---|---:|---:| +| 12 | 426 | 9.50% | +| 13 | 198 | 4.42% | +| 14 | 2,478 | 55.29% | +| 15 | 43 | 0.96% | +| 16 | 1,111 | 24.79% | +| 17 | 42 | 0.94% | +| 18 | 140 | 3.12% | +| 19 | 6 | 0.13% | +| 20 | 38 | 0.85% | + +这个分布说明，当前验证集中最主要的标准大环内酯集中在 **14 元环** 和 **16 元环**，其中 14 元环超过一半，16 元环接近四分之一。 + +### 2.3 为什么 4,482 会进一步变成 4,451 + +这 4,482 个标准大环内酯里，有 **31 个分子没有任何可拼接侧链碎片**，因此不会进入 `fragment_library.csv`。 + +换句话说： + +| 阶段 | 数量 | 占标准集比例 | +|---|---:|---:| +| `standard_macrolactone` 总数 | 4,482 | 100.00% | +| 至少产生 1 条侧链碎片 | 4,451 | 99.31% | +| `num_sidechains = 0` | 31 | 0.69% | + +这 31 个分子不是“分类失败”，而是“分类成功但没有通过侧链可拼接过滤”。它们的 `cleavage_positions` 为空，说明在 `position > 2` 的环位点上，没有找到满足规则的外侧链。 + +### 2.4 实际进入碎片库的 ring size 分布 + +在最终 **4,451** 个 fragment-bearing 母体中，ring size 分布如下： + +| ring size | 数量 | 占 fragment-bearing 母体比例 | +|---|---:|---:| +| 12 | 422 | 9.48% | +| 13 | 191 | 4.29% | +| 14 | 2,475 | 55.61% | +| 15 | 41 | 0.92% | +| 16 | 1,105 | 24.83% | +| 17 | 36 | 0.81% | +| 18 | 139 | 3.12% | +| 19 | 4 | 0.09% | +| 20 | 38 | 0.85% | + +这张表才是后续碎片库真正的母体分布，因为它排除了那 31 个“标准但无侧链”的分子。 + +## 3. 碎片库如何生成 + +碎片库的生成逻辑在验证流程中完成，结果写入： + +| 输出 | 含义 | +|---|---| +| [validation_output/fragment_library.csv](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library.csv) | 统一碎片库导出 | +| [validation_output/fragments.db](/Users/lingyuzeng/project/macro_split/validation_output/fragments.db) | `parent_molecules`、`side_chain_fragments`、`fragment_library_entries` 等表 | + +对每个 `standard_macrolactone`，程序会： + +1. 先做 canonical numbering。 +2. 只遍历 `position > 2` 的位点。 +3. 跳过环内原子和 intrinsic lactone neighbor。 +4. 用 `collect_fragmentable_side_chain_atoms()` 识别单锚点、可拼接侧链。 +5. 用 `build_fragment_with_isotope()` 生成 `fragment_smiles_labeled` 和 `fragment_smiles_plain`。 +6. 只有 `Chem.MolFromSmiles(plain_smiles)` 可解析时，碎片才会写入数据库和 CSV。 + +因此，这份库并不是“把所有侧链都随意裂开”，而是“只保留单锚点、可重拼接、SMILES 合法的侧链碎片”。 + +## 4. 统一碎片库的规模与去重 + +`fragment_library.csv` 的核心统计如下： + +| 指标 | 数值 | +|---|---:| +| 行数 | 34,829 | +| 去重母体数 | 4,451 | +| 唯一 `fragment_smiles_plain` 数 | 1,852 | +| 唯一 `fragment_smiles_labeled` 数 | 2,051 | + +这里最容易误解的一点是：`1,852` 不是“把 dummy 原子去掉以后才算出来的”，而是按 `fragment_smiles_plain` 去重得到的。`plain` 仍然保留 dummy 原子 `*`，只是去掉了 isotope 标签，因此： + +1. `*C` 和 `*O` 这样的结构会被当成不同碎片。 +2. `[13*]C` 和 `[16*]C` 会因为 isotope 标签被抹掉而收敛成同一个 plain 结构。 + +换言之，`1,852` 说明的是“**去掉位置标签后的 chemotype 数量**”，而不是“把 dummy 原子也删掉后的数量”。 + +### 4.1 高频碎片 + +当前库里最常见的 plain 碎片是： + +| 碎片 | 计数 | 占总行数比例 | +|---|---:|---:| +| `*C` | 15,927 | 45.73% | +| `*O` | 4,722 | 13.56% | +| `*=O` | 3,167 | 9.09% | +| `*CC` | 2,261 | 6.49% | + +这四类最简单的碎片合计已经占到总行数的绝大部分。它们非常适合解释为什么后续要使用 `>3` 重原子作为设计相关下限，因为库里真正占量的是大量一到三原子的“噪声型”片段，而不是高信息密度的大碎片。 + +### 4.2 为何 `1852` 仍然不算“小” + +`34,829` 条行数去重后剩 `1,852` 个 unique plain SMILES，并不反常，原因是： + +1. 很多碎片在不同母体、不同位置上会重复出现。 +2. 小碎片的复用率极高，尤其是 `*C`、`*O`、`*=O`、`*CC` 这类。 +3. 只要保留 dummy 原子 `*`，就会把“裸碎片”与“带连接点的片段”区分开来，因此 unique 数不会无限压缩。 + +## 5. 侧链筛选阈值：为什么下限是 `>3` 重原子 + +关于你关心的“上限和下限”： + +1. **下限是明确的**：后续设计相关侧链分析采用 `>3` 重原子。 +2. **上限不是硬编码规则**：当前代码和报告没有设置统一的碎片上限。 +3. **如果问当前库的经验上沿**：`fragment_atom_count_summary.csv` 显示最大重原子数是 **48**，`p95 = 14`，`p99 = 27`。 + +### 5.1 结果文件位置 + +如果你要找“上限是多少”的依据，应该看这几个文件： + +| 文件 | 作用 | +|---|---| +| [fragment_atom_count_summary.csv](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/fragment_atom_count_summary.csv) | 给出 `max_atom_count = 48`、`p95 = 14`、`p99 = 27` 等摘要 | +| [fragment_atom_count_frequency.csv](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/fragment_atom_count_frequency.csv) | 给出每个原子数的频次分布 | +| [fragment_atom_count_filter_candidates.csv](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/fragment_atom_count_filter_candidates.csv) | 给出不同下限阈值的删减效果 | +| [analysis_summary.txt](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/analysis_summary.txt) | 文字版摘要，最适合快速引用 | + +### 5.2 为什么下限选 `>3` + +阈值候选表显示： + +| 下限规则 | 删除行数 | 删除行比例 | 删除 unique fragments | +|---|---:|---:|---:| +| `<=1` | 23,994 | 68.89% | 10 | +| `<=2` | 28,069 | 80.59% | 26 | +| `<=3` | 28,550 | 81.97% | 52 | +| `<=4` | 29,045 | 83.39% | 88 | +| `<=5` | 29,272 | 84.04% | 141 | + +这组数据给出的结论非常清楚： + +1. `<=3` 已经去掉了 **82%** 左右的记录。 +2. 但它只去掉了 **2.8%** 的 unique fragments。 +3. 换句话说，低于等于 3 重原子的碎片主要贡献的是“行数噪声”，不是“结构多样性”。 + +所以，`>3` 是一个很自然的设计相关下限，因为它把大量单原子、双原子、三原子碎片排除掉，同时尽可能保留更像“真实取代基”的片段。 + +### 5.3 上限为什么不设成硬规则 + +当前库的重原子数分布是长尾的，最大值到 48，但： + +1. 大于 14 的片段已经属于少数。 +2. 大于 27 的片段更是极少。 +3. 这些长尾片段中仍可能包含糖环、杂环、稠环等有意义的单锚点片段。 + +因此，当前项目更适合把上限当作“任务依赖的经验裁剪”，而不是“全局硬上限”。 + +如果你必须写一个论文里的操作性描述，可以这样写： + +> 本研究采用 `>3` 重原子作为设计相关侧链下限；上限不做全局硬截断，仅在具体下游任务中按库规模与化学可解释性进行任务特异性裁剪。当前库的经验分布显示重原子数最大为 48，95% 分位数为 14，99% 分位数为 27。 + +这个表述比“固定上限 = 某个整数”更稳妥，也更符合当前代码和数据实际。 + +## 6. 16 元环的位点碎片提取 + +16 元环是本项目里最重要的 ring size 之一，后续位点导出固定在 canonical 位置 `3/4/12/13`。 + +### 6.1 位点导出数量 + +`validation_output/ring16_position_fragment_exports/` 下四个导出文件的 unique fragment 数量如下： + +| 位点 | unique fragment 数 | +|---|---:| +| 3 | 121 | +| 4 | 70 | +| 12 | 99 | +| 13 | 198 | + +所以，你前面问的“3、4、12 位置是不是 121、70、99？”答案是：**对，正确**。如果把 13 位也算上，那么 13 位是 **198** 条 unique fragments。 + +### 6.2 对应的全量与设计相关子集 + +| 位点 | 全量碎片数 | unique plain SMILES | `>3` 重原子后碎片数 | `>3` unique plain SMILES | +|---|---:|---:|---:|---:| +| 3 | 1,048 | 121 | 269 | 117 | +| 4 | 595 | 70 | 269 | 63 | +| 12 | 930 | 99 | 177 | 83 | +| 13 | 876 | 198 | 709 | 193 | + +这说明： + +1. 13 位在天然库里最丰富。 +2. 3 位和 4 位的 unique chemotype 也不低。 +3. 12 位虽然总量不如 13 位，但设计相关子集仍然相当可观。 + +### 6.3 位置编号与文献标签的关系 + +由于项目统一使用 canonical numbering，16 元环的这四个位点在文献视角下可以理解为镜像位置： + +| canonical position | literature-style mirrored label | +|---|---| +| 3 | 16 | +| 4 | 15 | +| 12 | 7 | +| 13 | 6 | + +因此，当前数据里最富集的 canonical 位点 `13, 3, 4, 12`，对应文献常说的 `6, 16, 15, 7` 一组方向镜像标签。 + +## 7. tylosin 的枚举结果 + +本项目里 tylosin 的枚举结果保存在： + +| 文件 | 含义 | +|---|---| +| [validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/tylosin_scheme_b_unique_products.csv](/Users/lingyuzeng/project/macro_split/validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/tylosin_scheme_b_unique_products.csv) | 唯一产物表 | +| [validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/tylosin_scheme_b_fix_pos13.sqlite](/Users/lingyuzeng/project/macro_split/validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/tylosin_scheme_b_fix_pos13.sqlite) | 组合 provenance 与运行元数据 | +| [validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/README.md](/Users/lingyuzeng/project/macro_split/validation_output/enumeration/per_scaffold/tylosin/scheme_b_fix_pos13/README.md) | 运行说明 | + +### 7.1 枚举设定 + +该结果对应 `scheme_b_fix_pos13`，核心设定是： + +| 参数 | 值 | +|---|---:| +| `reference_slug` | `tylosin` | +| `replace_positions` | `[3, 4, 12, 13]` | +| `fixed_positions` | `[13]` | +| `candidate_counts` | `{3: 121, 4: 70, 12: 99}` | + +这意味着： + +1. 13 位的糖基保持不变。 +2. 只对 3、4、12 位进行枚举。 +3. 拼接空间是这三个位点片段库的笛卡尔积。 + +### 7.2 枚举规模 + +| 指标 | 数值 | +|---|---:| +| 尝试组合数 | 838,530 | +| 成功拼接数 | 838,530 | +| 去重后唯一产物数 | 810,810 | +| 合并掉的重复组合数 | 27,720 | +| 重复率 | 3.31% | + +这里最关键的一点是： + +`838,530 = 121 × 70 × 99` + +所以原始组合数就是 3 个位点片段库的直积。去重后变成 810,810，说明有一小部分不同组合最终坍缩成了同一个 canonical product SMILES。 + +### 7.3 唯一产物分布 + +`occurrence_count` 的分布是： + +| occurrence_count | 唯一产物数 | +|---|---:| +| 1 | 790,020 | +| 2 | 13,860 | +| 3 | 6,930 | + +这说明大多数 unique product 都只对应单一路径，但也有一部分产物由多个组合收敛而来。 + +## 8. 可以直接写进论文的结论 + +### 8.1 验证与碎片库生成 + +在 MacrolactoneDB `ring12_20` 验证集中，11,036 个输入分子经分类后得到 4,482 个 `standard_macrolactone`，其中 31 个标准大环没有任何可拼接侧链碎片，最终形成 4,451 个 fragment-bearing 母体和 34,829 条可拼接碎片记录。 + +### 8.2 侧链筛选阈值 + +碎片大小分布极度右偏，`*C`、`*O`、`*=O`、`*CC` 四类最简单碎片已占据绝大多数记录。基于当前库的频率分布，建议将 `>3` 重原子作为后续设计相关侧链的下限；同时，当前库没有全局硬上限，最高观察值为 48 重原子，95% 分位数为 14，99% 分位数为 27。 + +### 8.3 16 元环位点 + +在 canonical 16 元环编号下，`3/4/12/13` 位导出的 unique fragments 数分别为 121、70、99、198。若按文献镜像标签表达，则对应 `16/15/7/6` 位。 + +### 8.4 tylosin 枚举 + +以 tylosin 为参考骨架、固定 13 位糖基不变时，`3/4/12` 位枚举共产生 838,530 个组合，去重后得到 810,810 个 unique products，重复折叠率为 3.31%。 + +## 9. 推荐的图表引用 + +| 图 | 文件 | +|---|---| +| 全库碎片原子数分布 | [fragment_atom_count_distribution.png](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/fragment_atom_count_distribution.png) | +| 16 元环位点数量对比 | [ring16_position_count_comparison.png](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/ring16_position_count_comparison.png) | +| 16 元环设计相关碎片 boxplot | [ring16_position_atom_count_boxplot_gt3.png](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/ring16_position_atom_count_boxplot_gt3.png) | +| 16 元环位点多样性图 | [ring16_position_diversity_gt3.png](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/ring16_position_diversity_gt3.png) | +| 16 元环药化热点对比 | [ring16_medchem_hotspot_comparison.png](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/ring16_medchem_hotspot_comparison.png) | +| 16 元环环状侧链敏感性 | [ring16_position_ring_sensitivity.png](/Users/lingyuzeng/project/macro_split/validation_output/fragment_library_analysis/ring16_position_ring_sensitivity.png) | + +## 10. 一句话总结 + +这条链路的核心结论是：`temp.csv` 中 11,036 个分子经过标准大环内酯筛选后，真正进入可拼接碎片库的是 4,451 个母体、34,829 条碎片；碎片大小分布强烈偏向 1-3 重原子小片段，因此后续设计分析采用 `>3` 重原子作为下限，而 16 元环的 3、4、12、13 位碎片与 tylosin 的 fixed-pos13 枚举结果一起构成了后续论文中最重要的参数和设计依据。