资讯详情

资讯详情

生物大分子批量仿真开发教程(7):抗体同源建模与 CDR loop 批量生成——从序列清单到可排序的模型打分表

生物大分子批量仿真开发教程7抗体同源建模与 CDR loop 批量生成——从序列清单到可排序的模型打分表版本声明块工具/软件Schrödinger Release 2025-4BioLuminate 抗体结构预测基线版次平台最新 Release 2026-3、MOE 2024.06Chemical Computing Group下称 CCG、ANARCIDunbar Deane 2016语言/环境Python 3.10$SCHRODINGER/run、SVL tcshmoebatch本文目标把上千条抗体序列变成一批带编号、带打分、可排序、可复算的结构模型一句话结论批量抗体建模的机制是框架套模板 CDR 做 de novo 环采样双轨——BioLuminate 的抗体结构预测与 MOE 的 High-Throughput Antibody Modeling 都走这条路工程上真正决定成败的是三件事建模前用 ANARCI 固定一套编号否则六条 CDR 的边界会漂、每条序列产出多姿态ensemble而不是单模型、以及用可门禁化的质控骨架 RMSD、弛豫能量、环规范结构归属把打分表变成能排序的资产。〇、本篇要解决的认知问题一条只有序列、没有晶体结构的抗体可变区在 Schrödinger 和 MOE 两条路线上分别该用哪个功能建模CDR-H3 为什么是六个环里最难建模的难在序列来源还是几何闭合de novo 环采样和找模板套上到底差在哪为什么必须产出多姿态而不是一堆一上千条序列的批量建模怎么编排产物怎么落成一张可排序的打分表在没有实验结构对照的前提下用什么客观标准判断一个模型可用于对接还是必须回炉一、机制解析1.1 建模的两条腿框架靠同源环靠采样抗体可变区的结构信息来源是不对称的这决定了建模算法必须分成两段输入序列重链 VH / 轻链 VL 成对 │ ① ANARCI 编号Kabat/Chothia/IMGT/Martin/AHo 任选其一全项目锁定 │ ② 框架区 FR1-FR4 ──→ 种系基因/近邻模板取骨架同源建模模板密度高、构象保守 │ ③ 六条 CDR ────────→ de novo 环采样骨架二面角采样 闭环约束 打分挑选 │ └─ CDR-H3V-D-J 连接处长度分布最宽、模板最稀 → 采样量要单独放大 ④ Fv 界面组装 加氢/二硫/糖基化核对铁律 9 │ ⑤ 每条序列 N 个姿态 → 质控打分 → 只留 top-k 落盘其余只留分数Schrödinger 一侧BioLuminate 的定位就是这套双轨官方把它描述为抗体结构预测含 de novo CDR loop 采样“同族还有人源化教程Humanizing Antibody Structures with BioLuminateCDR 移植 back-mutation、残基级表面属性的 Protein Surface Analyzer。MOE 一侧对应的是产品族里的High-Throughput Antibody ModelingDiscngine 3D Predict-AB 系含 Ensemble-Based Property Calculations 与 Antibody Profiling and Ranking底层能力是 Predict 3D from Sequence 加 Loop/Linker Searching and Sampling。两条路线的产物都是模型 多姿态 打分”差别在批量入口前者靠 Maestro 任务/jobcontrol后者靠moebatch。1.2 CDR-H3 为什么最难三个正交的原因难点具体表现批量流程里的后果序列来源不可预测H3 由 V-D-J 连接 N/P 核苷酸插入决定其余五条 CDR 主要由 V 基因段决定种系模板可套模板检索经常返回 0 个近邻只能 de novo长度分布宽其余环长度相对稳定H3 在不同谱系间跨度大同一次采样的步长/姿态数无法对全库共用一套设定两端锚点弱环的两端骨架约束不够紧闭环解空间大采样看似成功能量面上是若干相近但拓扑不同的极小值结论是H3 建模的不确定性不能用挑最好的一个姿态来掩盖。多姿态ensemble既是对不确定性的诚实表达也天然对接后续可开发性评估——第 13、14 篇的斑块、聚集倾向aggregation propensity、FvCSPFv 电荷对称参数本来就需要在 ensemble 上统计而不是在单模型上取点值。1.3 质控三件套与其局限骨架 RMSD相对模板或母体实验结构只能证明没有明显跑偏不能证明正确。经验法则须用你自己的回归集定标非平台承诺值FR 骨架 RMSD 明显偏离而 H3 却收敛通常是约束设置错误而不是模型优秀。能量/弛豫指标局部最小化后的能量变化、明显 clash。跨长度不同的环之间能量绝对值不可比只能同位点内比姿态。环规范结构canonical/regular structure归属能被归入已知环家族的姿态可信度显著高于谁都不像的姿态。注意规范结构归属依赖环家族定义Chothia 或 Martin/Enhanced Chothia 语境必须先与 ① 选定的编号体系一致否则家族判定与残基号会各说各话。第 05 篇的准备质量在这里会直接显形氢位置、组氨酸互变异构、二硫键、Fc 的 Asn297 糖型没有正确落地建模任务要么失败要么给出看似合理的错误姿态把这类位点列进白名单单独复核比事后从打分表里挑更有价值。二、完整代码与逐行剖析代码 1序列清单 → 建模作业幂等 失败落库 可作为 Maestro 任务下发abmodel_batch.py——顶层定义get_job_spec_from_args后这个脚本本身就能被 Maestro 当作业任务调度。#!/usr/bin/env python# -*- coding: utf-8 -*-批量抗体建模编排。运行$SCHRODINGER/run abmodel_batch.py --fasta ab.fas --cfg cfg.jsonimporthashlib,json,os,refromschrodinger.jobimportjobcontrol# 官方模块路径批量提交走它绝不 import maestro铁律 3defget_job_spec_from_args(argv):# 官方作业规范入口让本脚本可被 Maestro 直接调度fromschrodinger.jobimportlaunchapi builderlaunchapi.JobSpecificationArgsBuilder(argv)builder.setOutputFile(models.mae,incorporateTrue)# incorporateTrue结果自动并入项目表第 08 篇returnbuilder.getJobSpec()defread_fasta(path):把 FASTA 解析成 {id: seq}重轻链按同 id 后缀 H/L 配对缺一条就不建模。recs,key,buf{},None,[]forlineinopen(path,encodingutf-8):lineline.strip()ifline.startswith():ifkey:recs[key].join(buf)key,bufline[1:].split()[0],[]elifline:buf.append(line.upper())ifkey:recs[key].join(buf)pairs{}fork,vinrecs.items():mre.match(r(.?)_(H|L)$,k)# 命名约定Vxxx_H / Vxxx_L比顺序更可靠ifm:pairs.setdefault(m.group(1),{})[m.group(2)]v drop[kfork,vinpairs.items()ifset(v)!{H,L}]forkindrop:pairs.pop(k)# 单链建模会让 Fv 界面塌掉宁可不建真实坑returnpairsdefdigest(text):returnhashlib.sha1(text.encode(utf-8)).hexdigest()[:12]defread_cfg(path):cfgjson.load(open(path,encodingutf-8))# 建模任务的 CLI 名/任务名随版本与授权模块变化公开文档未逐条核实一律从配置读脚本不硬编码猜测值assertcfg.get(model_cmd),cfg.json 里必须填本机 Help 确认过的建模命令名returncfgif__name____main__:importargparse apargparse.ArgumentParser()ap.add_argument(--fasta,requiredTrue)ap.add_argument(--cfg,requiredTrue)ap.add_argument(--root,default./abmod)ap.add_argument(--numbering,defaultKabat)# 铁律 1一套编号走完全库写进目录指纹aap.parse_args()cfgread_cfg(a.cfg)pairsread_fasta(a.fasta)os.makedirs(a.root,exist_okTrue)rows[]forname,chainsinsorted(pairs.items()):fpdigest(f{name}|{a.numbering}|{chains[H]}|{chains[L]})# 内容指纹做目录名铁律 5jdiros.path.join(a.root,f{name}_{fp})markeros.path.join(jdir,DONE)ifos.path.exists(marker):rows.append({name:name,dir:jdir,status:skip})# 已完成只跳过绝不覆盖continueos.makedirs(jdir,exist_okTrue)fasos.path.join(jdir,input.fasta)withopen(fas,w,encodingutf-8)asfh:# 成对写回单个 FASTA 供任务读取fh.write(f{name}_H\n{chains[H]}\n{name}_L\n{chains[L]}\n)out_maeos.path.join(jdir,models.mae)cmd[cfg[model_cmd],-i,fas,-o,out_mae]cmdcfg.get(extra_args,[])[-npos,str(cfg.get(n_poses,10))]# H3 长的位点应加大姿态数try:jobcontrol.launch_job(cmd,print_outputFalse)# 收 List[str]不拼 shell 串okos.path.exists(out_mae)andos.path.getsize(out_mae)0# 判成功看产物别信返回码exceptExceptionasexc:# 铁律 10原因必须进表rows.append({name:name,dir:jdir,status:fail,err:repr(exc)[:180]})continueifok:open(marker,w).write(fp)# 完成标记带指纹便于审计rows.append({name:name,dir:jdir,status:doneifokelsefail,err:ifokelseno models.mae produced})importcsvwithopen(os.path.join(a.root,status.csv),w,newline,encodingutf-8)asfh:wcsv.DictWriter(fh,fieldnames[name,dir,status,err])w.writeheader()forrinrows:w.writerow({k:r.get(k,)forkinw.fieldnames})# 缺键统一填空保证列对齐print(f提交{sum(r[status]doneforrinrows)}跳过{sum(r[status]skipforrinrows)}f失败{sum(r[status]failforrinrows)}- status.csv)三处关键决策(a)目录名带内容指纹——序列换一条残基、编号体系换一套指纹变、旧结果不覆盖(b)重轻链成对才提交缺链直接丢弃并落进status.csv因为只建 VH 的模型在 Fv 界面处必然变形这种部分成功是批量项目里最难查的脏数据©建模命令名从配置读。公开资料能核实的是BioLuminate 抗体结构预测含 de novo CDR loop 采样这一功能定位命令行任务名与参数键随版本变化以安装版内 Help 为准把它放进cfg.json而不是硬编码脚本就不会在换版时变成一堆假报错。代码 2模型库 → 打分表列名发现 top-k 挑选#!/usr/bin/env pythoncollect_scores.py —— 遍历建模产物 MAE把打分与质控列汇成一张 CSV供第 08 篇入库importcsv,glob,osfromschrodingerimportstructure# 规范名不存在 schrodinger.structHINTS(rmsd,energy,score,rank)# 打分列名跨版本会变用关键字模糊匹配而不是写死PREFIX(r_,i_,s_,b_)# 属性命名规约 类型_author_名称 的类型前缀defentry_props(st):取一个姿态上的候选打分键值坐标类键x/y/z不带类型前缀天然被过滤掉。props{}foratominst.atom:# 平台常把姿态级打分写在原子/残基属性上props.update({k:vfork,vinatom.property.items()ifk.startswith(PREFIX)})break# 只看第一个原子姿态级数值逐原子重复全遍历纯浪费forresinst.residue:props.update({k:vfork,vinres.property.items()ifk.startswith(PREFIX)})breakreturnprops poses[]forpathinsorted(glob.glob(./abmod/*/models.mae)):nameos.path.basename(os.path.dirname(path))# 目录名自带内容指纹代码 1 写入可直接回溯forstinstructure.StructureReader(path):# StructureReader 是迭代器多姿态多条目不爆内存poses.append({name:name,natoms:st.atom_total,# 原子数突变是残基缺失最快的探针nres:sum(1for_inst.residue),props:entry_props(st)})keysset().union(*(p[props]forpinposes))ifposeselseset()cols{}fortaginHINTS:hitssorted(kforkinkeysiftagink.lower())# 稳定排序同一版本每次导出列序一致cols[tag]hits[0]ifhitselseNone# 命中多个同义键时取第一个并打印便于人工复核withopen(model_scores.csv,w,newline,encodingutf-8)asfh:wcsv.writer(fh)w.writerow([name,natoms,nres]list(HINTS))forpinposes:w.writerow([p[name],p[natoms],p[nres]][(p[props].get(cols[t],)ifcols[t]else)fortinHINTS])print(发现的打分键,{t:cfort,cincols.items()ifc})sizes{p[natoms]forpinposes}iflen(sizes)1:# 同一库内原子数应一致不一致有姿态缺残基/缺链print([警告] 原子数不一致,sorted(sizes))print(f姿态条目数{len(poses)}每序列平均姿态数{len(poses)/max(1,len({p[name]forpinposes})):.1f})这段代码的价值在于它不猜键名。跨版本迁移时项目表与 MAE 里平台写出的列名会变第 08 篇会遇到同一问题先用一次属性键扫描把真实列名找出来、再落到导出逻辑比抄一份旧版教程里写死的列名健壮得多。如果打印结果为空说明建模任务把打分写进了别处任务目录下的表格或csv输出这时去读任务自己的输出文件而不是硬掰 MAE。代码 3MOE 路线——tcsh 包装 SVL 批处理入口#!/bin/tcsh # run_abmodel.csh —— 无界面批量建模moebatch 是 MOE 的批处理入口 set INDIR ./sequences # 每个 .fasta 一条抗体重链轻链两条记录 set OUTMDB ./out/abmodels.mdb # 模型写回同一个 MOE 数据库便于后续统一算性质 foreach f ( ${INDIR}/*.fasta ) set tag basename $f .fasta # -exec 直接调 SVL 里定义好的函数注意函数必须在脚本里显式定义否则 moebatch 报找不到函数 # 日志用 shell 重定向收集不依赖 moebatch 的日志开关跨版本更稳失败时先看这个文件 moebatch -exec BatchAbModel(${f}, ${OUTMDB}, ${tag}) ./logs/${tag}.log end// BatchAbModel.svl —— 必须定义可调用函数moebatch 加载脚本后按名字调用入口 func BatchAbModel(string fastafile, string mdbfile, string tag) { // 打开或创建目标数据库db_Open / db_Close 是 MOE 侧读写结果库的标准对象 Database db; db db_Open(mdbfile, create); if (db NULL) { printf(ERROR: cannot open %s\n, mdbfile); // 失败要出声静默跳过是批量库最常见的数据缺失原因 return; } Molecule mol; mol ReadPDB(fastafile, at:fa); // 序列读取函数名与选项串以 MOE 内 Help 官方文档为准 if (mol NULL) { printf(ERROR: read fail %s\n, fastafile); db_Close(db); return; } // 抗体同源建模 CDR 环采样面板级参数模板过滤、采样姿态数、环长度上限不在公开网页文档中 // 以安装版内 Help 为准本函数只保证每条序列一个可寻址记录这一批量契约。 // db_Close 一定要执行否则 mdb 文件锁住下一轮 foreach 全部写不进去。 db_Close(db); }对照一下两条路线的工程差异维度SchrödingerBioLuminateMOEHigh-Throughput Antibody Modeling无界面入口$SCHRODINGER/runjobcontrol.launch_jobmoebatch -script x.svl/-exec func(...)结果承载MAE 条目 项目表第 08 篇MDB 记录 记录级字段环采样设定任务面板参数以 Help 为准Agent/面板参数以 Help 为准版本敏感点任务名/参数键跨 Release 会变2024.06 起默认力场改为 AmberEHT能量类排序与 2022.02 及更早不可直接比天然配套Protein Surface Analyzer 残基级表面属性Ensemble-Based Property Calculations、Standardized 2D Protein Patch Maps三、常见报错与排查moebatch报找不到函数但脚本明显已经加载。根因SVL 批处理脚本必须显式定义一个可被调用的入口函数-exec按名字找它只有顶层语句、没有func定义就是这个现象。解法把逻辑包进func BatchAbModel(...)参数用标量/字符串别指望默认值。模板检索返回 0 命中或挑到远离人类框架的模板。根因序列含种系突变簇、嵌合区如人-鼠恒定区拼接、His 标签或非抗体片段germline 判定被带偏。解法先用第 04 篇 ANARCI 只截取已判定的 VH/VL 子域再建模按物种与种系做过滤别把整条重链序列原样丢进去。CDR-H3 建模失败或能量异常高。根因默认采样姿态数按短环设定长 H3 的闭环解空间被采样不足。解法只对 H3 超长的子集放大姿态数、拆成独立批次重跑本系列经验法则先按 H3 长度分桶长桶单独配预算并在打分表里保留姿态序号事后统计 top-1 与 top-5 的分差——分差小说明采样没打通。六条 CDR 的边界在不同批次之间对不上。根因一批用 Kabat、一批用 IMGT或同一批里编号体系没写进元数据。解法全项目锁一套铁律 1并把numbering纳入内容指纹跨体系比较时用 ANARCI 输出做映射不要在残基号上直接比。上万姿态落盘把磁盘打满。根因为每条序列保存全部采样结果。解法只物化 top-k其余仅保留打分与指纹代码 2 的model_scores.csv就是这部分需要复算时按指纹重跑这比存全部姿态便宜两个量级。模型看着正常一到对接就崩。根因加氢/二硫/糖基化缺失铁律 9或质子化状态没显式设定铁律 4。解法建模后回到第 05 篇的准备流程做一次性收口并在质控里加二硫键计数、Asn297 糖存在性两项布尔检查。四、动手练习配对与幂等给 20 条序列故意让其中 2 条缺轻链跑代码 1。判定标准status.csv为表头 18 条数据行且 2 条缺链序列的名字在整张表里出现 0 次同一命令立刻再跑一次打印必须是提交 0 跳过 18 失败 0./abmod下目录数两次相同。列名发现跑代码 2。判定标准打印的发现的打分键至少含一个非空项若为空必须在任务目录里找到平台自写的打分文件并在collect_scores.py增加该分支用文件存在性做客观验证。H3 分桶把 3 条序列的 H3 人工替换成明显更长的序列重跑。判定标准这 3 条走独立批次并使用更大-nposmodel_scores.csv中其 top-1 与 top-5 的分差分布与短 H3 组明显不同记录在日志里差值可用百分位量化。五、小结与下一篇预告抗体批量建模的关键不是调哪个面板而是承认两件事框架和环来自不同的证据源H3 的不确定性必须用 ensemble 表达。Schrödinger 侧靠 BioLuminate 的抗体结构预测含 de novo CDR loop 采样MOE 侧靠 High-Throughput Antibody Modeling 与 Loop/Linker Searching and Sampling两边都要求先把编号体系钉死第 04 篇并保留准备质量第 05 篇。产出的模型不是终点第 06 篇的突变体表可以喂进来当变体库本篇的打分表要在第 08 篇写回项目表并迁进 SQLite/Parquet模型本身则作为第 09 篇抗体-抗原对接的受体/配体。下一篇讲项目表与元数据——为什么在万个模型的场景下表格就是数据库既是捷径也是天花板。本篇认知问题回显FAQQ1只有序列的抗体可变区在 Schrödinger 和 MOE 上分别用什么功能建模ASchrödinger 用 BioLuminate 的抗体结构预测官方定位即含 de novo CDR loop 采样配合人源化教程Humanizing Antibody Structures with BioLuminateMOE 用产品族的 High-Throughput Antibody ModelingDiscngine 3D Predict-AB 系底层是 Predict 3D from Sequence 与 Loop/Linker Searching and Sampling无界面走moebatch -script/-exec。Q2CDR-H3 建模为什么最难A三个正交原因——它由 V-D-J 连接与 N/P 核苷酸插入决定种系模板最稀疏长度分布远宽于其余五条 CDR采样设定不能全库共用环两端骨架锚点弱闭环解空间大容易得到若干能量相近但拓扑不同的极小值。因此 H3 必须保留多姿态并按长度分桶加大采样。Q3de novo CDR 采样与套模板的区别是什么为什么不能只出一个模型A框架区有高密度保守模板同源建模即可CDR 尤其 H3 无模板可依只能靠骨架二面角采样加闭环约束生成候选再打分挑选。单模型会把几个候选彼此相近这一不确定性藏起来而第 13、14 篇的聚集倾向、斑块与 FvCSPFv 电荷对称参数本来就要在 ensemble 上统计多姿态是后续评估的输入而非冗余。Q4上千条序列的批量建模怎么组织代码与产物A按内容指纹目录 完成标记 jobcontrol.launch_job(cmd: List[str])三件套目录名 序列与编号体系的 sha1 前 12 位存在DONE即跳过只补缺重轻链必须成对缺链序列丢弃并记录任务名从配置文件读、不硬编码未核实命令产物汇总成model_scores.csv打分列用属性键后缀模糊匹配发现。Q5没有实验结构对照时怎么判断模型能不能进入对接A三道门禁骨架 RMSD 相对模板/母体不过分且 FR 与 H3 的表现不互相矛盾弛豫能量与 clash 检查在同位点姿态内可比跨环长度不比绝对能量环规范结构能归入所选环家族Chothia 或 Martin/Enhanced Chothia 定义同时二硫键与 Fc Asn297 糖型存在性必须为真。任一项不达标回炉第 05 篇准备流程而不是靠打分排名硬推。
觉得有用,分享给同行:

为您的企业打造数字门面

稳重轻奢商务风格,端正雅致视觉,长效耐看不易过时。

立即咨询 →