ANARCI抗体序列编号实战指南：一条命令打通从单条序列到万级批量的全流程-育师

ANARCI抗体序列编号实战指南：一条命令打通从单条序列到万级批量的全流程

【免费下载链接】ANARCIAntibody Numbering and Antigen Receptor ClassIfication项目地址: https://gitcode.com/gh_mirrors/an/ANARCI

如果你的抗体数据来自三个不同的实验室，你会发现同一段 CDR-H3，有人标成 Kabat 的 100A–100H，有人标成 IMGT 的 111–112 附近，还有人干脆不标。论文里"第 104 位半胱氨酸"这种说法，在不同编号体系下根本对不上。这正是 ANARCI（Antibody Numbering and Antigen Receptor ClassIfication）存在的意义：它用隐马尔可夫模型自动识别链类型，再按你指定的国际标准编号方案，给每条抗体序列打上统一的位置标签。本文用一个贯穿始终的案例——把一条真实抗体序列从"裸序列"变成"带编号、带物种、带质量分数的结构化数据"，再扩展到批量场景，带你完整走一遍。

案例起点：把一条 118 个氨基酸的裸序列变成"可对位"的数据

先不做任何安装决策，直接看我们要解决的核心矛盾。下面这条是小鼠重链可变区序列（来自示例数据Example_scripts_and_sequences/12e8.fasta）：

EVQLQQSGAEVVRSGASVKLSCTASGFNIKDYYIHWVKQRPEKGLEWIGWIDPEIGDTEYVPKFQGKATMTADTSSNTAYLQLSSLTSEDTAVYYCNAGHDYDRGRFPYWGQGTLVTVSA

你肉眼能看出它的 CDR1 从哪开始、CDR3 有多长吗？大概率不能。人工数残基既慢又容易错，而且不同人对"框架区/CDR 边界"的理解还不一致。ANARCI 的处理方式是：把序列交给 HMMER 与预建的物种 HMM 库比对，得到最可能的链类型（H/L/A/B）和比对区间，再按编号方案把比对结果翻译成位置号。整个过程的输入输出都很简单，这也是它被大量抗体分析工作流采用的原因。

2 条命令完成安装，第 3 条命令跑出第一个结果

ANARCI 的依赖只有两个硬性要求：Python 环境和 HMMER3。用 conda 可以一次装齐：

conda install -c conda-forge biopython -y conda install -c bioconda hmmer=3.3.2 -y

拿到代码并安装（仓库地址https://gitcode.com/gh_mirrors/an/ANARCI）：

git clone https://gitcode.com/gh_mirrors/an/ANARCI cd ANARCI python setup.py install

这里有个容易被忽略的环节：setup.py install不是简单复制文件，它会自动执行build_pipeline/下的构建流程——从 IMGT 下载生殖系基因序列、用 MUSCLE 比对并训练 HMM。首次安装需要联网，耗时约几分钟，构建产物会被复制到 Python 的anarci包目录里。安装完成后验证：

ANARCI --help

然后直接对刚才那条裸序列出手：

ANARCI -i EVQLQQSGAEVVRSGASVKLSCTASGFNIKDYYIHWVKQRPEKGLEWIGWIDPEIGDTEYVPKFQGKATMTADTSSNTAYLQLSSLTSEDTAVYYCNAGHDYDRGRFPYWGQGTLVTVSA

stdout 会输出编号结果。第一屏是信息头，第二屏是逐位置的编号表。

先读懂输出，再谈优化：编号结果里的五个关键字段

编号文件里每个记录以//结尾。头部的#注释行是你判断结果质量的第一手资料：

# 12e8:H|... # ANARCI numbered # Domain 1 of 1 # Most significant HMM hit #|species|chain_type|e-value|score|seqstart_index|seqend_index| #|mouse|H|8.6e-58|184.9|0|119| # Scheme = imgt H 1 Q H 2 V H 3 Q ...

五个字段各有用途：

字段	含义	判断标准
species / chain_type	命中最优 HMM 的物种与链型	chain_type 决定后续用哪套编号规则
e-value	比对显著程度	远小于 0.001 才值得信任
score	bit 分数	分数越高代表比对越可靠
seqstart_index / seqend_index	编号区段在原始序列中的起止下标	用于定位编号覆盖不到的残基

编号表按链型位置氨基酸三列排列。注意 light chain 统一用L表示，不再区分 kappa/lambda，但链型识别是在 HMM 层面做的，kappa 和 lambda 的模型是分开训练的。如果你把结果喂给下游可视化工具，这个字段可以直接当输入。

顺带提一句反例：对非抗体序列运行 ANARCI，比如示例里的Example_scripts_and_sequences/lysozyme.fasta（溶菌酶），程序只会列出序列名，不产生任何编号。这其实是很好的负向验证手段——确认你的流程不会把非免疫球蛋白乱编成抗体。

六种编号方案怎么选：一张决策表解决 90% 的纠结

-s参数控制方案，命令行使用短代码，Python API 使用全名：

方案	CLI 短代码	适用范围	一句话特性
IMGT	i	所有抗原受体（IG + TCR）	128 个统一位置，CDR3 插入围绕 111/112 对称标注
Kabat	k	仅抗体	经典方案，框架和 CDR 都可能出现插入
Chothia	c	仅抗体	比 Kabat 更贴合结构，CDR-H1 插入位置不同
Martin（增强 Chothia）	m	仅抗体	修正了 Chothia 框架区部分插入位置
AHo	a	所有抗原受体	149 个位置，天然容纳长 CDR，几乎不用插入码
Wolfguy	w	仅抗体	CDR 内不同区段编号，长 CDR 无需插入码

选择逻辑其实很简单：通用交流选 IMGT，对比结构选 Chothia 或 Martin，处理超长 CDR3 优先试 AHo 和 Wolfguy。长 CDR3 是实践中最常见的坑——IMGT 方案里 CDR3 超过一定长度会以插入字母呈现（例如111-ABCD DCBA-112的对称形式），而 AHo 因为预留了 149 个位置，通常能免去插入码的解析麻烦。如果你的研究对象恰好是 CDR-H3 超过 30 个残基的抗体（比如牛、羊、鲨鱼来源），强烈建议同时跑一遍 IMGT 和 AHo 对比。

想要做方案间的一致性检查，可以跑两遍同一输入，把两个方案的编号结果按 IMGT 104 位半胱氨酸对齐，再比对 CDR3 长度标注。示例脚本Example_scripts_and_sequences/run_numbering_benchmark.sh就是同时用 6 种方案处理约一万条 PDB 序列，单轮不到 5 分钟（4 核并行），这个脚本本身也是写批量任务的样板。

批量场景：FASTA 输入、CSV 对齐输出、多核并行一条龙

实际项目里几乎不会只有一条序列。ANARCI 对 FASTA 文件的处理方式和单序列完全一致：

ANARCI -i antibody_sequences.fasta --csv -o results/pdb_imgt --ncpu 4 --assign_germline

这里几个参数值得单独讲：

--csv：把编号结果按链型分别写入独立的 CSV 文件，且各序列按编号方案水平对齐。这意味着你不用自己写代码就能得到"可进 Excel、可做多序列比对"的矩阵，是后续计算 CDR 长度分布、做序列 motif 分析最省力的入口。
--ncpu 4：开启多进程。ANARCI 内部按序列块切分任务，基准测试表明约一万条序列、4 核环境下单方案 < 5 分钟。
--assign_germline：给每条序列追加最相似生殖系基因的注释，形如：

#|mouse|IGHV1-12*01|0.86|IGHJ2*01|0.79|

表示 V 区与 IGHV1-12*01 的序列一致性 0.86、J 区与 IGHJ2*01 一致性 0.79。这对筛选"接近生殖系"的抗体、或者做体细胞高频突变分析非常有用。

另外还有-r ig参数：只保留抗体链、丢弃 TCR。处理混合来源的 PDB 序列时，这是避免 TCR 链污染抗体分析的首选开关。

Python API 进阶：把编号能力塞进你自己的分析管道

CLI 适合人机交互，但自动化流程必须走 API。项目示例Example_scripts_and_sequences/anarci_API_example.py把核心用法全部写清楚了，拆开看只有三件事：

from anarci import anarci # 1. 准备 (名称, 序列) 列表 sequences = [ ("12e8:H", "EVQLQQSGAEVVRSGASVKLSCTASGFNIKDYYIHWVKQRPEKGLEWIGWIDPEIGDTEYVPKFQGKATMTADTSSNTAYLQLSSLTSEDTAVYYCNAGHDYDRGRFPYWGQGTLVTVSAAKTTPPSVYPLAP"), ("lysozyme:A", "KVFGRCELAAAMKRHGLDNYRGYSLGNWVCAAKFESNFNTQATNRNTDGSTDYGILQINSRWWCNDGRTPGSRNLCNIPCSALLSSDITASVNCAKKIVSDGNGMNAWVAWRNRCKGTDVQAWIRGCRL"), ] # 2. 调用，返回三件套 numbering, alignment_details, hit_tables = anarci(sequences, scheme="imgt", output=False) # 3. 逐个解析：未编号的序列为 None for i, name in enumerate([s[0] for s in sequences]): if numbering[i] is None: print(name, "未识别为抗原受体") else: for domain_num, start, end in numbering[i]: print(name, "识别出", len(numbering[i]), "个结构域", "编号区间:", start, "-", end)

返回值的结构值得记住：

numbering[i]：第 i 条序列的编号结果，None表示未识别；否则是结构域列表，每个结构域是(编号表, 起始下标, 结束下标)三元组。scFv 这类含两个可变域的序列会返回两个结构域。
alignment_details[i][j]：第 j 个结构域的比对详情字典，species、chain_type、e-value、score 都在这。
hit_tables[i]：HMMER 对库中所有 HMM 的命中统计，可以看作"这条序列到底像谁的"完整证据链。

如果只想要最轻量的结果，用number()：

from anarci import number numbering, chain_type = number("EVQLQQSGAEVVRSGASVKLSCTASGFNIKDYYIHWVKQRPEKGLEWIGWIDPEIGDTEYVPKFQGKATMTADTSSNTAYLQLSSLTSEDTAVYYCNAGHDYDRGRFPYWGQGTLVTVSAAKTTPPSVYPLAP", scheme="kabat") print(chain_type) # H 或 L

number()只返回第一个结构域和链型，适合快速打标签；处理批量数据时官方明确建议用run_anarci而不是循环调number()，后者会重复加载 HMM 库，慢得多。

另一个值得接入的脚本是Example_scripts_and_sequences/ImmunoPDB.py：它扩展了 Biopython 的 PDBParser，可以直接给 PDB 结构文件里的抗体链重新编号，并按 CDR 区段做注释，输出到xtra字典属性里。做结构分析时，先用它把 PDB 编号统一成 IMGT，再提取 CDR，能省掉大量手工对位。

避坑清单：五个高频翻车点及对症下药

1.number()返回(False, False)不代表序列不是抗体。函数内部有长度保护：序列短于 70 个残基直接返回 False。抗体可变区片段被这样"误杀"时，改用anarci()全参数版本，它不会做这个长度检查。

2. ANARCI 报物种，但它不是物种注释工具。README 里写得很直白：物种判定是编号过程的副产品，官方不建议把它当主要物种注释工具。跨物种嵌合抗体、人源化抗体尤其容易"误判"，此时应看--assign_germline的基因名而不是相信 species 字段。

3. 用 Kabat/Chothia/Martin 编号 TCR 会直接失败。这三个方案只定义了抗体链，对 TCR 序列要么报 AssertionError（API 里返回 False），要么不输出。需要处理 TCR 时换 IMGT 或 AHo，命令行则用-r tcr（配合-s imgt）。

4. 超长 CDR3 编号失败。ANARCI 内部有一个check_for_j机制：当 CDR3 过长导致 V/J 区无法在同一条比对里覆盖时，会尝试把 104 位半胱氨酸之后的序列单独比对 J 区再拼接。如果仍失败，先确认你的方案选的是 IMGT/AHo 而不是位置数较少的方案，再检查序列末端是否缺失 J 区。

5. 安装卡在 HMM 构建阶段。setup.py install需要联网下载 IMGT 生殖系数据。公司内网、代理环境经常在这里超时。解决办法是手动把下载步骤拆出来跑bash build_pipeline/RUN_pipeline.sh，看清楚是哪一步失败，再单独处理网络或镜像问题。

行动清单：按这个顺序练，两天内上手

第 1 小时：装好依赖，用12e8.fasta跑通 CLI，然后换-s k、-s m各跑一遍，肉眼对比 CDR-H1 的编号差异。
第 1 天：把anarci_API_example.py逐行跑通，确保你理解numbering三元组的三个值各自指向什么。
第 1 天下午：用antibody_sequences.fasta+--csv生成对齐矩阵，导入 Excel 或 pandas，统计你手头数据的 CDR3 长度分布。
第 2 天：挑一条你研究里的真实抗体，跑--assign_germline，对比它和最近生殖系基因的一致性，体会高频突变的量级。
进阶目标：把run_numbering_benchmark.sh里的 6 方案并行脚本改造成你自己的批量流程，加入pdb_sequences.fa.txt.gz（约万条 PDB 序列）做压力测试，验证 4 核 5 分钟的性能预期。