简介：一份聚焦微生物群落功能预测的专题 PDF，围绕 FAPROTAX 在 16S 扩增子分析中的应用展开，面向微生物生态、环境科学方向的研究生与科研人员，帮助理解环境样本中碳、氢、氮、磷、硫等元素循环的功能注释思路，也适合作为组会汇报与课题设计时的概念参考。内容从数据库的文献构建依据讲起，说明 OTU 分类表与功能注释的关联方式，并延伸到功能—样本矩阵、分组功能丰度热图、环境因子相关性泡泡图等结果呈现，同时对比 ImageGP 一键分析与 Linux 本地运行两种实现路径，也提示了预测准确度与覆盖度的适用边界。压缩包仅含 1 个 PDF 文件，大小约 2.47MB，轻量便于查阅。已有 1148 人学习，读者可据此建立 FAPROTAX 预测流程的整体框架，并通过哺乳动物皮肤微生物、水稻氮利用等案例，把握功能预测在差异挖掘与机制关联中的用法。

1. 拿到 16S 结果之后，功能预测这一步怎么落地

扩增子测序做到最后，ASV 表和物种注释都齐了，合作方或者审稿人往往不会满足于“这批样本里拟杆菌门多、变形菌门少”，他们要的是——这些菌在干哪些活，氮循环、硫循环、碳降解到底有没有差异。FAPROTAX（Functional Annotation of Prokaryotic Taxa）就是为这个问题准备的：它不训练模型，也不依赖参考基因组，而是用一份人工从文献里整理出来的“分类名—功能”对照表，把 16S 的属种注释直接映射成元素循环相关的功能组丰度。它最适合环境、农业、水体沉积物这类关注生物地球化学循环的样本，对肠道、皮肤这类宿主相关样本也能用，但解释要谨慎。想预测整条 KO 通路的人，应该去看 PICRUSt2；想把“这些菌可能参与哪几步氮转化”讲清楚的人，继续往下看。

2. FAPROTAX 在算什么：物种名到元素循环功能的映射表

2.1 数据库不是模型推导，而是文献里“抄”出来的对照表

理解 FAPROTAX 的关键，是先接受它“简单得有点粗暴”的设计。数据库文件 FAPROTAX.txt 本质上是一张两列的映射表：左边是分类单元名（属名或种名），右边是用制表符分隔的若干功能标签，比如nitrogen_fixation、nitrification、sulfate_respiration、methanotrophy、cellulolysis。这些标签不是算法推断出来的，而是整理者逐篇翻阅已发表的纯培养、富集培养和生理实验文献，确认“这个属确实具备这个代谢能力”之后才写进去的。

这样做的好处很直接：每一条注释背后都有一篇可追溯的文献，预测出来的功能不是统计相关性，而是生理学证据。代价同样明显：只有被研究过、被培养过、被验证过的类群才会进表，冷门门类和大量未培养微生物基本是空白。所以它在环境样本上的表现，取决于样本里有多少类群落在“已培养、有文献”的范围里。

这也解释了它为什么和 2016 年那篇海洋微生物组研究绑在一起：那项工作想说明，海洋里功能组成和物种组成并不是一一对应的，同一个功能可以被很多不同类群分担，而 FAPROTAX 提供了把物种翻译成功能的工具，让这种“解耦”能在 16S 数据上被观察。经过几年迭代，数据库在 2020 年 11 月随 1.2.4 版做了一次较大更新，注释比例比早期版本高不少，脚本运行环境也从 Python 2 迁到了 Python 3。

提示：FAPROTAX 的功能标签是“能力”而不是“活性”。一个属被标了denitrification，只说明它有这个基因通路，不代表样本里真的在反硝化。写文章时把措辞放在“潜在功能”“功能潜力”上，比直接写“反硝化强度”稳得多。

2.2 和 PICRUSt2、Tax4Fun 的边界：覆盖度换准确度

同样是 16S 功能预测，几条技术路线的取舍完全不同，选错工具比参数调错更致命。

PICRUSt2 走的是“参考基因组树 + 基因家族推断”的路线，把 ASV 序列放进参考树，推测基因拷贝数，再归到 KO 通路。它的覆盖度很高，几乎每个 ASV 都能落进去，输出的是完整的 KEGG 通路谱，适合做全局的功能差异分析。但它推断的是“基因组潜力”，链条长、假设多，落在具体某一步氮转化上时，中间误差会被放大。

Tax4Fun 系（含 Tax4Fun2）基于 SILVA 注释结果做最近邻匹配，输出的是百分比形式的功能组成，实现简单，但对参考库依赖强，遇到数据库里没有的类群就会失真。

FAPROTAX 的定位正好卡在另一端：覆盖度最低，但每一条映射都有生理学依据，且直接对应到“固氮、硝化、反硝化、硫酸盐还原、甲烷氧化、纤维素降解”这种可以直接和实测理化指标对话的功能组。样本是水稻根际、湖泊水体、海洋沉积物、活性污泥，研究问题是“氮素转化菌群随施肥怎么变”，FAPROTAX 的说服力比通路预测强；样本是人体肠道，研究问题是“整体代谢通路差异”，就该换 PICRUSt2。

我的习惯是两条线并跑：FAPROTAX 出元素循环的功能组，用来和理化数据做关联、画机制图；PICRUSt2 出 KO 通路，用来做无假设的差异筛选。两者结论互相印证时，故事才立得住；只有 FAPROTAX 结果时，讨论部分一定要写清覆盖度这个前提。

2.3 输入的硬性要求：带物种注释的 OTU/ASV 表

FAPROTAX 不吃普通的文本丰度表，它要的是 biom 格式，而且物种注释必须写进 biom 的 observation metadata 里。这一步是绝大多数人第一次翻车的地方。

具体来说，输入需要满足三条：

第一，丰度表是“行 = OTU/ASV，列 = 样本”的矩阵，第一列是 OTU ID，表头是样本名，值是非负整数或抽平后的数。

第二，每个 OTU 必须有分类字符串，且字符串的最后一级落在属或种。FAPROTAX 是靠属名/种名去查表的，如果你的注释只到门（p__Proteobacteria）或者最后一级是g__uncultured，它就只能算作 leftovers。

第三，分类字符串的层级顺序要和数据库预期一致，通常是界门纲目科属种依次排列，用分号分隔。前缀写不写k__、p__无所谓，脚本会做清理，但顺序不能乱，也不能中间插别的字段。

把这三条对齐，后面的命令基本就是流水线。对不齐，就会看到Assigned 12 records to groups这种让人心里一凉的结果——我见过最惨的一次，两千多条 ASV 只匹配上十几条，原因就是上游注释表的属名那一级全带着g__前缀加空格。

3. 本地跑通 1.2.4：从 txt 丰度表到功能矩阵

3.1 建一个干净的 Python3 环境，别和 QIIME2 混用

FAPROTAX 1.2.4 的脚本本身只依赖 numpy，但整条流程要用 biom 命令行工具做格式转换，而 biom 的版本和 QIIME2 环境里的依赖经常打架。血泪经验是：单独建一个环境，别图省事直接在 QIIME2 的 conda 环境里跑。

# 建一个只装必要依赖的干净环境 conda create -n faprotax python=3.8 -y conda activate faprotax # numpy 是脚本硬依赖，biom-format 提供 biom 命令行 conda install -c conda-forge numpy biom-format -y # 验证 biom 命令可用，正常会打印版本号 biom --version

环境建好后，去项目主页下载FAPROTAX_1.2.4.zip并解压，目录里至少要有collapse_table.py和FAPROTAX.txt两个文件。先跑一次帮助，能弹出参数列表就说明解释器路径和依赖都没问题。

# 只看帮助，不做任何计算 python FAPROTAX_1.2.4/collapse_table.py -h

如果这一步报ModuleNotFoundError: No module named 'numpy'，说明 conda 环境没激活对，或者你在别的 shell 里跑；报biom相关错误则是后面章节才会遇到的问题，此处不涉及。先把-h跑通，能省掉后面一半的排查时间。

3.2 txt 转 biom 再挂物种注释，两条命令

假设手上是抽平后的otutab_rare.txt和一份单独的注释文件taxonomy2.txt。注释文件的格式必须是两列、制表符分隔、带表头，表头写作OTUID和taxonomy，顺序不能反：

OTUID taxonomy ASV_1 k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Nitrosomonadales;f__Nitrosomonadaceae;g__Nitrosomonas ASV_2 k__Bacteria;p__Proteobacteria;c__Alphaproteobacteria;o__Rhizobiales;f__Rhizobiaceae;g__Rhizobium

转换分两步：先把文本表转成 biom JSON，再把注释作为 observation metadata 挂上去。

# 第一步：txt 转 biom（JSON 格式，注意表类型写 OTU table） biom convert -i otutab_rare.txt -o otutab_rare.biom \ --table-type="OTU table" --to-json # 第二步：把物种注释写入 observation metadata # --sc-separated 表示 taxonomy 这一列按分号拆成多值 # --observation-header 指定注释文件的表头字段顺序 biom add-metadata -i otutab_rare.biom \ --observation-metadata-fp taxonomy2.txt \ -o otutab_rare_tax.biom \ --sc-separated taxonomy \ --observation-header OTUID,taxonomy # 第三步：抽查前几条 observation，确认 taxonomy 真的进去了 biom summarize-table -i otutab_rare_tax.biom | head -20 biom convert -i otutab_rare_tax.biom -o - --to-tsv | head -5

这三条命令里，--sc-separated是最容易被忽略的参数。它告诉 biom：这一列的值不是普通字符串，而是用分号分隔的多值字段，后续 FAPROTAX 才能按层级取出属名。如果漏写，taxonomy 会被当成一个整串存进去，脚本取最后一段时就会拿到整个字符串，匹配全军覆没。

--observation-header的顺序也要和文件实际列顺序严格对应。写成taxonomy,OTUID而文件是OTUID,taxonomy，biom 不会报错，但会把 ASV ID 塞进 taxonomy 字段，结果是注释全错、匹配为零，属于最隐蔽的一类错误。

3.3 collapse_table.py 的参数逐条拆

输入准备好之后，一条命令出结果：

python FAPROTAX_1.2.4/collapse_table.py \ -i otutab_rare_tax.biom \ -g FAPROTAX_1.2.4/FAPROTAX.txt \ --collapse_by_metadata 'taxonomy' \ -v --force \ -o faprotax.txt \ -r faprotax_report.txt

逐个说清楚：

-i是上一步产出的、带注释的 biom 表。注意不能用原始 txt，也不能用没挂 metadata 的 biom。

-g指向数据库文件FAPROTAX.txt，路径写错会直接报文件不存在。建议用绝对路径或确认当前工作目录。

--collapse_by_metadata 'taxonomy'告诉脚本：按 observation metadata 里名为taxonomy的字段来做归并。这个字段名必须和biom add-metadata时写的一致，大小写敏感。

-v输出详细过程信息，会打印每个功能组匹配到多少条记录、哪些记录被归为 leftovers。排查阶段强烈建议开着。

--force强制覆盖已有输出文件，批处理时必备，否则第二次运行会因为文件存在而中断。

-o是功能丰度矩阵，行是功能组，列是样本，值是该功能组在所有匹配 OTU 上的丰度加和。

-r是明细报告，记录每个功能组由哪些 OTU/ASV 贡献，这是后面做证据链的关键文件。

运行完，屏幕上会出现类似Assigned 979 records to groups, 1924 records were leftovers的统计。这个比例要盯住：匹配上的记录数除以总记录数，就是这次预测的有效覆盖度。经验上环境样本常见三到六成，落到一两成就得回头查注释格式，别急着往下分析。

脚本还支持用--column之类的参数指定从分类字符串里取第几段作为匹配键，不同版本参数名略有差别，跑一次-h确认，比照着旧教程硬套可靠。

3.4 两个输出文件分别看什么

faprotax.txt是主结果，第一列是功能组名，后面每列一个样本。它可以直接读进 R 或 Python 做统计，但要注意两点：一是行数远少于 KO 表，通常几十到一百多行，这是覆盖度决定的；二是同一个 OTU 可能被同时计入多个功能组，所以列求和会大于总 reads 数，这是“一对多”映射的正常现象，不是数据错了。

faprotax_report.txt是黑匣子钥匙。它按功能组组织，每个组下列出归属到该组的 OTU/ASV 及其分类信息。想看“固氮功能到底是谁贡献的”，直接在这里翻；审稿人质疑“你这个固氮功能是不是被某个占主导的属带偏了”，也在这里查。很多人只拿faprotax.txt画图，把 report 丢在一边，等被问细节时才发现答不上来。

4. report 拆成三类表：让结果能进热图和分组比较

4.1 先把 report 清洗成“功能组 + 记录”两行一组

report 文件的结构是：每个功能组先占一行组名，下面跟着若干行以 ASV/OTU ID 开头的记录行，组与组之间用分隔线隔开。这种嵌套结构不适合直接做表，需要先压平成“组名—记录”交替的形式。

# 抓取 ASV 行以及它的上一行（即所属功能组名），去掉纯分隔线 grep 'ASV_' -B 1 faprotax_report.txt | grep -v -P '^--$' > faprotax_report.clean

-B 1的作用是连带输出匹配行的前一行，这样每条记录都跟着它的功能组；-v -P '^--$'把分隔线剔掉，避免后面解析时出现空组。如果你的 OTU ID 不是ASV_开头，把 grep 的模式换成实际前缀，比如OTU_。这一步看着土，但对格式规整的 report 极其有效，比写正则解析整个文件稳。

4.2 用脚本拆出 mat / otu_func / func_otu

清洗完的文件是交替的两行结构，用一段短脚本拆成三张表：OTU×功能的有无矩阵、OTU 到功能的列表、功能到 OTU 的列表。三种表对应三种下游用途，缺一不可。

# report2tables.py # 读入清洗后的 report，输出三类整理表 import sys clean, prefix = sys.argv[1], sys.argv[2] pairs = [] # (功能组, ASV) cur_group = None with open(clean) as fh: for line in fh: line = line.rstrip('\n') if not line.strip(): continue if not line.startswith(prefix): # 不以 ASV_ 开头的行是功能组名 cur_group = line.strip() else: asv = line.split('\t')[0].strip() # 记录行首个字段是 OTU ID if cur_group: pairs.append((cur_group, asv)) # 表一：OTU -> 功能，一行一个 OTU，功能用分号连接 otu2func = {} for g, a in pairs: otu2func.setdefault(a, set()).add(g) with open('faprotax_report.otu_func', 'w') as out: out.write('OTUID\tfunctions\n') for a, gs in otu2func.items(): out.write(f"{a}\t{';'.join(sorted(gs))}\n") # 表二：功能 -> OTU，用于反查某个功能由谁贡献 func2otu = {} for g, a in pairs: func2otu.setdefault(g, set()).add(a) with open('faprotax_report.func_otu', 'w') as out: out.write('function\tOTUIDs\n') for g, as_ in func2otu.items(): out.write(f"{g}\t{';'.join(sorted(as_))}\n") # 表三：OTU x 功能 0/1 矩阵，直接拿去画热图 funcs = sorted({g for g, _ in pairs}) otus = sorted(otu2func) with open('faprotax_report.mat', 'w') as out: out.write('OTUID\t' + '\t'.join(funcs) + '\n') for a in otus: row = ['1' if g in otu2func[a] else '0' for g in funcs] out.write(a + '\t' + '\t'.join(row) + '\n') print(f'{len(otus)} OTUs, {len(funcs)} functions, {len(pairs)} pairs')

跑法是python report2tables.py faprotax_report.clean ASV_，第二个参数就是 grep 用的前缀，两边保持一致。最后打印的三个数字可以当作快速自检：OTU 数、功能组数、配对关系数。配对关系数应该明显大于功能组数，如果两者接近，说明大部分功能组只对应一条记录，样本的功能谱很可能被少数类群撑起来，解读时要收着点。

.mat是 OTU 层级的有无矩阵，适合展示“某个功能由哪些 ASV 分担”；.otu_func适合在补充材料里做逐条溯源；.func_otu适合在正文表格里写“固氮功能由 12 个 ASV 贡献，归属 3 个属”。

4.3 下游三种图的数据准备

主结果faprotax.txt进分析之前要做标准化。因为功能组丰度受测序深度影响，样本间直接比大小没有意义。常见做法是转成相对丰度（每列除以该列总和），或者用 TSS（total sum scaling），也有人用 centered log-ratio 处理组成型数据。抽平后的数据我一般先用相对丰度，热图里再做行方向的 z-score，视觉上更好读。

热图的行是功能组，列是样本，按实验分组排序，配上分组色条，就能看到“某个处理下硝化相关功能整体抬升”这种模式。分组比较用 STAMP 或者 R 里的 Wilcoxon + BH 校正，注意功能组数量少，多重检验校正后很容易全不显著，这时候更该看效应量和具体功能组的生物学含义，而不是死盯 q 值。

和环境因子的关联分析是 FAPROTAX 的高光场景：功能组丰度矩阵与 pH、总氮、有机质这些理化指标做 Spearman 相关，画成气泡图，横纵轴是功能组和环境因子，气泡大小表示相关系数绝对值，颜色表示方向。因为功能组数量有限，一张图能塞下全部，比 KO 通路那种几百行的矩阵清爽得多。做相关时样本量小的项目要谨慎，n 小于 10 时相关系数的可信区间非常宽，图可以画，结论别写死。

5. 常见问题排查：五个把我卡住的现场

5.1 匹配数为零或个位数

现象：命令跑完，屏幕打印Assigned 5 records to groups，两千多条 ASV 几乎全是 leftovers。

原因：八成出在分类字符串上。最常见的是最后一级不是属/种，而是g__uncultured、f__Xanthobacteraceae这类；其次是从上游带下来的g__前缀后面跟了空格或不可见字符；还有一种是分类顺序被打乱，属名跑到了门的前面。

解决：先抽几条看实际字符串，awk -F'\t' 'NR<=3{print $2}' taxonomy2.txt直接肉眼确认；再统计最后一级的分布，如果大量卡在科级，说明上游注释器只给到科，FAPROTAX 无能为力，要么换注释工具重跑，要么在文章里如实写“仅部分类群可归到属级”。空格和特殊字符用sed 's/ //g'清一遍再重新挂注释。

5.2 taxonomy 挂上去了，但脚本报找不到字段

现象：collapse_table.py提示 metadata 里没有taxonomy字段，或者归并结果为空。

原因：biom add-metadata的--observation-header顺序和注释文件列顺序不一致，导致字段名和内容错位；或者--sc-separated拼写成了别的字段名。

解决：用biom convert -i xxx.biom -o - --to-tsv | head -3把 biom 打回文本，直接看 observation metadata 那几列到底叫什么、值长什么样。这一步比猜快得多。

5.3 biom 命令报版本或依赖错误

现象：biom convert抛ImportError，或者与 QIIME2 的命令互相覆盖。

原因：QIIME2 环境里的 biom 是绑定版本，和独立安装的 biom-format 冲突；或者在一个 Python 2 的旧环境里跑 1.2.4 的脚本。

解决：按 3.1 的方式单独建环境，全程只在这个环境里操作。转换和预测分两个环境跑也可以，中间用文件传递，不要试图在一个环境里同时满足所有工具的依赖。

5.4 把功能组丰度当绝对量解读

现象：某样本固氮功能丰度是另一组的 3 倍，直接写“固氮能力增强 3 倍”。

原因：功能组丰度是匹配到的 OTU 丰度加和，既受该类群相对丰度影响，也受测序深度和抽平方式影响，还受覆盖度影响——不同样本匹配上的类群集合可能根本不一样。

解决：统一抽平深度后再比较，报告相对丰度而非绝对 reads；同时打开 report 看两组匹配到的 OTU 是否同源，如果一组匹配到 20 个 ASV、另一组只有 3 个，这个倍数没有意义。表述改成“固氮相关功能组相对丰度显著升高”，并附上贡献 OTU 列表。

5.5 覆盖度低就宣布“样本没有该功能”

现象：结果里没有methanotrophy，直接得出“该环境不存在甲烷氧化菌”的结论。

原因：FAPROTAX 是阳性证据库，没匹配上只代表“没有已知的、被文献收录的类群”，不代表功能不存在。未培养微生物、新属、数据库未收录的类群全都落在 leftovers 里。

解决：先报覆盖度，再谈阴性结果。表述改成“在 FAPROTAX 可注释的范围内未检出”，并在讨论里说明未培养类群的可能贡献。有条件的话用 qPCR 或宏基因组对关键功能基因做一次验证，把话说圆。

6. 让结果经得起追问：抽平、注释库与证据链

功能预测这类工作，最容易被问的不是“你怎么跑的”，而是“你凭什么信”。三个动作能把结果的可信度抬一个台阶。

第一个是抽平的一致性。所有样本必须抽到同一深度再算功能，否则高深度样本匹配到的 OTU 更多，功能组丰度天然偏高，组间差异会掺进测序深度的假信号。抽平深度取样本 reads 数的最小值或去掉极端低值后的分位数，抽三次看功能组排序是否稳定，稳了再定稿。

第二个是注释库的选择。FAPROTAX 不限注释库，但经验上不同库的属名写法、层级完整性差别很大，同一批数据换库重跑，匹配率能差出一大截。我一般用 SILVA 132 打底，它在属级的信息完整度和 FAPROTAX 的命名习惯比较对得上；也见过用 RDP 注释后匹配率明显更高的项目，这跟样本类群有关，没有普适答案。稳妥做法是拿一小批样本分别用两个库跑一遍，比一下 Assigned 比例，再决定用哪个，别直接上全量。

第三个是证据链。每一张功能热图、每一个显著差异的功能组，都要能顺着 report 反查到具体 ASV 和分类。写补充材料时把func_otu表附上，审稿人问“这个功能由谁贡献”，一查就有。下面这张表是我常用的自检清单：

还有一个具体技巧值得单独说：做组间比较前，先把faprotax.txt里在所有样本中相对丰度都低于 0.01% 的功能组过滤掉。这些低丰度功能组本身由一两条 ASV 撑着，稍一抽平就消失，留着只会制造假阳性，让多重检验校正后的结果更难看。过滤阈值别拍脑袋，按你数据的总丰度分布来定，写进方法里就行。

最后提一句和实测数据的交叉验证。如果项目里同时测了总氮、铵态氮、pH，把硝化、反硝化、固氮这几个功能组单独拎出来做散点或相关分析，方向对得上，故事就闭环了；方向对不上，先别急着删数据，回头查抽平和注释库，很多时候问题就出在这两处。从那以后我每次跑 FAPROTAX，都强制先看 Assigned 比例、再抽平三次自检、最后反查一遍 report 里的贡献 OTU，这三步不走完，后面的图一张都不出。希望帮到你。

本文还有配套的精品资源，点击获取