FAPROTAX 16S元素循环功能预测：输入制备、参数与验收

简介：这份PDF讲义面向微生物组16S扩增子分析的学习者与研究人员，聚焦FAPROTAX原核物种功能注释及碳、氢、氮、磷、硫等元素循环的功能预测，适合已具备OTU表与基础生信背景、希望把分类结果推进到生态功能层面的读者。文件共1个PDF，大小约2.47MB，篇幅集中，便于按主题通读或检索。内容围绕FAPROTAX数据库的构建依据与适用边界展开，说明其基于已发表可培养菌文献、预测准确度较好但覆盖度有限，并涉及OTU分类表经Python脚本转换为功能矩阵、ImageGP一键分析与Linux本地分析两条路径。结果部分展示功能描述与样本对应的矩阵，以及功能丰度热图、环境因子相关性泡泡图等应用，并引用海洋微生物组、哺乳动物皮肤微生物和水稻氮利用等案例，帮助读者理解如何解读功能差异、评估方法局限并形成后续统计分析思路。已有1148人学习。

1. 元素循环功能预测：一份 16S 表能回答到什么程度

有位做湿地沉积物的同学拿来一份抽平后的 OTU 表和一张采样信息表，想知道这批样本在碳氮磷硫的元素循环功能上到底有没有分层。项目里没有宏基因组预算，测序只做到 16S V3–V4。FAPROTAX 就是为这种处境准备的工具：数据库把物种分类名和已发表文献里验证过的功能一一对上，再用一个 Python 脚本把带物种注释的 OTU 分类表折叠成“功能 × 样本”矩阵，元素循环相关的通路能直接落到表上，后面接统计和热图就行。它的边界同样清楚——条目来自可培养菌的已验证文献，命中的部分准确度较好，覆盖度却偏低，一份表里往往只有一部分 OTU 能进功能表，其余全被记成 leftovers。所以这套方案适合环境、农业、水体这类以生物地球化学循环为关注点的课题，不适合拿它当宏基因组的平替。

2. 技术选型与数据结构：FAPROTAX 在四类预测工具里的位置

选型这一步翻车的成本最高，因为工具选错，后面所有图表都要重做。16S 功能预测常见的有四条路：KO 通路预测、元素循环预测、表型预测、以及 R 包实现的参考库映射。它们的输入长得几乎一样，输出含义却完全不同，混用会直接毁掉讨论部分。

2.1 什么时候选 FAPROTAX、什么时候换别的

判断标准其实就一条：你的科学问题落在哪个层级。

如果问题是“这批样本的氮循环潜力如何”“稻田根际的硫代谢功能是否富集”“不同深度水体的碳固定途径有无差异”，FAPROTAX 是对口的，因为它的功能条目本身就是按元素循环过程组织的，产甲烷、甲烷氧化、固氮、硝化、反硝化、硫酸盐还原这些都能直接读出来，写讨论时不用再做二次翻译。

如果问题是“整体代谢通路谱系差异”“某个 KO 模块是否富集”，那属于 PICRUSt2 的战场，参考库更大、覆盖更广，但代价是推断链条更长，结果里混着大量未验证的推测。

如果只关心革兰氏染色、氧气利用、生物膜形成这类表型，那是表型预测工具的活，FAPROTAX 表里根本没有这些字段。

还有一个容易忽略的差别：不同工具输出量纲不同。有的输出百分比，有的输出预测计数，FAPROTAX 输出的是每个功能在各样本中的丰度矩阵，量纲跟输入 OTU 表一致。这意味着两套结果不能塞进同一张热图里做 z-score 比较，量纲不同，颜色再好看也是错的。

我一般的做法是先写一句能被答辩老师追问三次的问题陈述，再反查工具能不能回答。问题陈述里出现“元素循环”“碳氮磷硫”“生物地球化学”这类词，选 FAPROTAX；出现“通路”“KEGG”“代谢模块”，换工具。

2.2 数据库与映射逻辑：一张文献表和一段查表脚本

FAPROTAX 的构成比很多人想的简单：一张人工整理的对应表，加一个折叠脚本。

对应表是制表符分隔的文本文件，每行给出一个功能条目和与之相关的物种分类名。关系是多对多——一个功能对应多个物种，一个物种也可能同时挂多个功能，比如某类芽孢杆菌既可能出现在固氮条目下，也可能出现在其他代谢条目下。所以折叠不是一对一替换，而是把同一功能下所有命中 OTU 的丰度加起来。

脚本的匹配逻辑是沿分类层级逐级向下查表：拿到 OTU 的 taxonomy 字符串，按分号切层，从门到属逐级比对，最细一级能命中数据库就计入，命中不了就往上一层退，一直退到底还没有，这条 OTU 就进 leftovers。关键点在于——它只查表，不做任何推断。数据库里没有的物种，无论 sequence 相似度多高，都不会被赋功能。这也是覆盖度低的根本原因：大量环境样本里的未培养菌，分类注释只到门或纲，根本没机会命中。

版本上，1.2.4 版在 2020 年 11 月更新过数据库，同时把脚本运行环境切到 Python 3，注释比例比早期版本有明显提升。物种注释数据库它不限定，但经验上 SILVA 132 与新版库的配合更好，而 RDP 注释的结果匹配率往往最高，这一点在做预实验时值得花半天对比一下，比后面反复解释“为什么注释率只有三成”省事得多。

2.3 输入输出数据结构：从 OTU 表到功能矩阵

输入侧只有两个东西：一张带物种注释的 OTU 表（biom 格式，物种注释挂在 observation metadata 上），以及数据库文件。输出侧是两张表，结构差异很大，用途也完全不同。

功能矩阵的行名是功能描述，不是 KO 号，也不是通路 ID，所以后续做富集分析时没法直接对接 KEGG 注释体系，这一点在方法部分要写明白，别让审稿人误以为你做了通路富集。

report 文件的信息密度被严重低估。它记录的是每个功能条目下的来源 OTU 及贡献量，等于给你留了一条从“功能显著富集”倒推回“具体是哪些菌在贡献”的证据链。答辩时被问“你凭什么说固氮功能升高了”，能当场翻出这张表，比任何解释都有说服力。

3. 输入表制备：环境搭建与 biom 格式改造

这一步看着是纯操作，实际是整个流程里报错最密集的地方。绝大多数失败不是脚本的问题，而是输入的 biom 少了元数据、或者注释列名跟参数对不上。

3.1 运行环境与依赖的确定

FAPROTAX 是纯脚本方案，没有网络服务，全部在本地跑。需要的环境很朴素：Python 3、numpy、以及处理 biom 格式的 biom-format 包。

我一般不去动系统自带的 Python，直接复用一个已有的 conda 环境（做过扩增子分析的机器上通常都有 Qiime2 环境，里面 numpy 和 biom 都是现成的），省掉一堆编译问题。

# 解压发布包（压缩包从项目官方发布页获取，注意不要用几年前的旧版本） unzip FAPROTAX_1.2.4.zip -d ~/biosoft/ cd ~/biosoft/FAPROTAX_1.2.4 # 复用已有分析环境，避免污染系统 Python conda activate qiime2-2020.6 # 冒烟测试：能打印帮助信息说明脚本与依赖都正常 python collapse_table.py -h | head -20

冒烟测试这一步别省。它的作用是把“环境问题”和“数据问题”提前分开，否则后面报错时你分不清是 numpy 没装好还是 biom 表结构不对。

如果上一条命令报ModuleNotFoundError: No module named 'numpy'，说明当前环境缺依赖，补装即可；如果是 biom 相关命令找不到，装的是 biom-format 而不是 biom：

conda install -y numpy biom-format

参数说明：-y是免交互确认，批量部署时必须有；环境名按你自己机器上的实际情况改，这里只是示例。

3.2 OTU 表转 biom 并写入物种注释

脚本只认 biom 格式，所以文本 OTU 表要先转，而且必须把物种注释挂上去。整个改造是两步：格式转换 + 元数据添加。

# 第一步：文本表转 biom json 格式 # 输入表必须是抽平后的（文件名里的 rare 就是抽平标记），未抽平会让功能丰度不可比 biom convert -i otutab_rare.txt -o otutab_rare.biom \ --table-type="OTU table" --to-json # 第二步：把物种注释作为 observation metadata 挂到 biom 上 # --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

注释文件taxonomy2.txt是两列制表符分隔的纯文本，第一行是表头，后面每行一个 OTU：

OTUID taxonomy ASV_1 d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Enterobacterales;f__Enterobacteriaceae;g__Escherichia ASV_2 d__Bacteria;p__Bacteroidetes;c__Bacteroidia;o__Bacteroidales;f__Bacteroidaceae;g__Bacteroides

几个参数值得逐个说清：-i和-o是输入输出，别写反；--table-type对 OTU 表固定写OTU table；--to-json输出可读的 JSON 而非二进制 HDF5，调试时方便直接打开看结构；--sc-separated taxonomy告诉 biom 这个字段是用分号分隔的列表，脚本后面就靠这个分隔符切层级；--observation-header里的两个列名必须与文件表头逐字一致，多一个空格都会导致元数据挂空。

这里有个隐含要求：taxonomy 字符串里要保留属名（g__那一段）。只到门或纲的注释，脚本只能往上层退，命中概率大幅下降。所以做注释时尽量别在属级提前截断。

3.3 输入表自检：三个必须先过的检查点

在跑预测之前，我习惯先做三项检查，把能提前发现的问题挡在门外，比事后从 leftovers 里猜原因高效得多。

第一项，OTU ID 是否一一对应。OTU 表和注释文件的 ID 集合必须完全一致，不能有注释表里多出来的、也不能有 OTU 表里没注释的。用一条命令就能验：

# 比较两边 OTU 数量，数字对不上就是 ID 或表头出了问题 grep -c '^ASV_' taxonomy2.txt head -2 otutab_rare.txt | tail -1 | tr '\t' '\n' | grep -c 'ASV_'

第二项，taxonomy 分隔符。有些流程输出的分类注释用竖线或逗号分隔，直接喂进去会整条字符串当一级处理，结果全部 leftovers。跑之前head看一眼实际分隔符，跟--sc-separated对齐。

第三项，是否抽平。功能预测的输出是绝对丰度折叠，样本测序深度差两三倍时，功能表的样本间比较毫无意义。抽平（或至少做相对丰度转换）要在这一步之前完成，不能等出结果再补。

4. 功能预测执行与结果整理

脚本跑起来只要一条命令，真正花时间的是参数确认和结果拆分。这一步的输出直接决定后面能画什么图。

4.1 collapse_table.py 的参数怎么配

# -i 带物种注释的 biom；-g 数据库文件路径 # --collapse_by_metadata 指定按哪一列元数据折叠，必须与 add-metadata 里的列名一致 # -v 输出过程信息；--force 允许覆盖已有结果（重复跑时必需） # -o 功能丰度矩阵；-r 功能到 OTU 的来源报告 python collapse_table.py \ -i otutab_rare_tax.biom \ -g FAPROTAX.txt \ --collapse_by_metadata 'taxonomy' \ -v --force \ -o faprotax.txt \ -r faprotax_report.txt

参数里最容易错的是--collapse_by_metadata。它的值不是随便起的名字，而是 biom 里那一列的列名，也就是上一步--observation-header里声明的taxonomy。写成taxonomy_name、Taxonomy、或者带引号多打了空格，脚本不会报“列名错误”，而是安静地一个都匹配不上，最后你会看到一行刺眼的统计：Assigned 0 records to groups。

跑完终端会打印一行汇总，格式类似Assigned 979 records to groups, 1924 records were leftovers。这两个数字是后面所有解释的起点，务必记下来写进分析日志：Assigned 是命中数据库的 OTU 条数，leftovers 是没命中的。经验上环境样本的命中比例本就不高，但低到个位数百分比就值得停下来查原因，而不是硬着头皮往下做。

--force这个参数也值得提一句，第一次跑无所谓，第二次改参数重跑时如果没有它，脚本会拒绝覆盖并直接退出，新手常在这里以为是脚本崩了。

4.2 从 report 文件追到 OTU：注释链路留存

report 文件的结构是“功能名 + 缩进的 OTU 明细”交替排列，格式规整但不利于直接读表，需要做一次清洗。

# 抓取每个 ASV 行及其上一行（即所属功能名），去掉分隔线 # 前缀 ASV_ 要与自己数据里的 OTU ID 前缀一致，用 Greengenes 风格命名就要换成对应的前缀 grep 'ASV_' -B 1 faprotax_report.txt | grep -v -P '^--$' > faprotax_report.clean # 用整理脚本把清洗后的报告拆成三类表 perl faprotax_report_sum.pl -i faprotax_report.clean -o faprotax_report

-B 1是关键，它把匹配行的前一行一起带出来，那一行正是功能名，否则你拿到的是一堆没有归属的 OTU 列表。grep -v -P '^--$'过滤掉报告里的横线分隔行，不然后面解析会错位。

这段清洗逻辑依赖两个前提，动手前务必确认：一是 OTU ID 前缀，如果你的表用的是OTU_1而不是ASV_1，grep 模式必须跟着改，否则清洗结果是空文件；二是 report 生成时的功能名行不含制表符缩进，如果版本差异导致格式变了，先head -30看一眼再改正则。

4.3 三类整理表与下游分析接口

整理脚本会吐出三张表，各自解决的问题不一样：

下游统计从faprotax.txt开始，通常是先转相对丰度再算差异，这里给一段 R 侧的处理骨架：

# 读入功能丰度矩阵，第一列是功能名 fap <- read.table("faprotax.txt", header = TRUE, row.names = 1, sep = "\t", check.names = FALSE) # 转相对丰度，避免测序深度带来的样本间偏差 fap_rel <- sweep(fap, 2, colSums(fap), "/") # 只保留在至少 20% 样本中出现过的功能，过滤掉偶发噪声 keep <- rowSums(fap_rel > 0) >= 0.2 * ncol(fap_rel) fap_filt <- fap_rel[keep, , drop = FALSE]

sweep那一行的含义是逐列除以该列总和，得到每个样本内各功能的相对占比；过滤阈值 20% 不是硬规定，样本量小的时候可以放宽到 10%，但一定要在方法里写清楚，别让审稿人猜。

出图部分，常用的两张是样本分组的功能丰度热图，以及环境因子与功能丰度的相关性泡泡图。后者记得做多重检验校正，Bonferroni 保守但稳妥，BH 更适合探索性分析；图上的星号必须说清对应的是哪种校正，这是审稿人必问的点。

5. 避坑与排查：注释率、格式与结果解读的常见翻车点

这一章是我自己踩过、也见过别人反复踩的地方。按“现象 → 原因 → 解决”整理，遇到问题可以直接对号入座。

5.1 注释率异常：从 leftovers 反推问题

现象一：Assigned 数字极低，leftovers 占九成以上。

原因通常是物种注释层级太浅。环境样本里大量未培养菌只注释到门，脚本沿层级往下退也碰不到数据库条目，全部落进 leftovers。

解决方式是先分层统计注释深度：算一下属级注释的 OTU 占总 OTU 的比例，如果这个比例本身就低，那不是 FAPROTAX 的问题，是上游注释的问题。可以换参考库重注释（SILVA 132 与新版数据库配合较好），或者接受现状，在方法里如实报告命中比例，把讨论限定在能注释到属的那部分菌群上。硬把三成命中率的结果讲成整体群落功能，答辩现场很容易被追穿。

现象二：换了参考库之后注释率反而下降。

原因多半是物种命名体系变了。不同库对同一属的写法不一致（有无g__前缀、是否有空格、属名是否重命名），而数据库里的物种名是固定写法，对不上就匹配失败。

解决办法是做一次名称规范化的比对：把数据库里的物种名列表提取出来，与注释表里的属名做交集，看交集大小。交集明显偏小就是命名体系不匹配，这时候换库比调参数有用。经验上 RDP 注释的结果与这套库匹配率最高，做预实验时值得花半天实测对比。

现象三：重复跑同一条命令，两次注释数不一致。

这种情况少见但有，通常是中间文件被--force覆盖时上一次的残留没清干净，或者输入 biom 被别的流程改过。做法是每次跑之前把输出文件删掉，或者干脆新建一个只读输入的目录，用绝对路径引用，避免手滑。

5.2 格式与脚本报错：八成出在 biom 上

现象四：脚本顺利跑完，但所有 OTU 都进了 leftovers，Assigned 为 0。

原因几乎可以锁定在--collapse_by_metadata的值与 biom 里的元数据列名不一致。脚本找不到那一列，自然一条都折叠不了，而且它不会明确报错。

排查方式是先验证 biom 里到底有什么：

# 查看 biom 的观测元数据表头，确认列名到底叫什么 biom table-ids -i otutab_rare_tax.biom --observations | head -3 biom show-install-info 2>/dev/null || true

更直接的办法是把 biom 转回文本看一眼：

# 转回 tsv 检查 taxonomy 列是否真的挂上去了 biom convert -i otutab_rare_tax.biom -o check.txt --to-tsv --header-key taxonomy head -3 check.txt

如果最后一列整列是空的，问题出在add-metadata那一步，回去核对--observation-header的列名与文件表头是否逐字一致，以及--observation-metadata-fp路径是否指对。

现象五：清洗报告时得到空文件。

原因就是 4.2 提到的前缀问题。grep 'ASV_'只在 OTU ID 以ASV_开头时有效，用OTU_1命名的表跑出来就是零行。

解决办法是先看两行原文再定模式：

# 看报告里 OTU 行的实际写法，再决定 grep 的前缀 head -10 faprotax_report.txt

一条head能省掉半小时的瞎猜，这类排查成本极低但最容易被跳过。

5.3 结果解读上的边界：别把预测当实测

现象六：把功能丰度当作绝对定量写进结论。

最常见的表述错误是“该样本固氮功能显著高于对照，说明固氮速率更高”。功能丰度反映的是具备该功能的类群相对占比，跟实际表达活性、速率是两回事，中间隔着转录和翻译两层。

写法上建议统一改成“固氮功能相关类群的相对丰度更高”，并在讨论里明确这是基于 16S 的预测结果，需要 qPCR 或宏基因组数据支撑时如实说明。这句话加上去不会削弱结论，反而显得专业。

现象七：把功能条目当成互斥分类。

一个 OTU 可能同时挂在多个功能下，因此功能矩阵的行与行之间并不独立，加起来也不等于 100%。做柱状堆叠图时尤其注意，堆起来超过总量不是数据错了，是数据结构本来如此。这类图我一般改成并排柱或热图，避免读者误读。

6. 验收与答辩：怎么证明这份功能预测站得住

功能预测类的工作在答辩现场最容易被打的点，不是结果不显著，而是“凭什么信你”。有三件事必须在提交前做完。

第一件，把输入侧的证据留全。抽平后的 OTU 表、注释文件、biom 的元数据截图、脚本的完整参数，全部归档，最好写成一个可复跑的 shell 脚本。答辩时被问到某个数字怎么来的，能当场打开文件，说服力远高于口头解释。

第二件，做一张验收自查表，逐项打勾再往下走：

第三件，准备一条能讲三层的证据链。第一层是功能层：哪个元素循环功能在组间显著差异，校正方法是什么。第二层是来源层：翻 report 说明这个功能主要由哪几个属贡献。第三层是文献层：这些属在已发表的可培养研究中被验证过什么代谢特征。三层讲完，一个“基于预测”的结果就变成了有据可查的推断，这是这类分析在答辩中最有分量的部分。

还有个小技巧：把faprotax_report.otu_func里几条典型 OTU 的功能注释单独拎出来做一张小表，附在附录。评审老师翻到附录看到你能追到单个 OTU 挂了几条功能、有没有明显矛盾的注释，对你数据处理细致程度的判断会完全不一样。

我早年做过一版没留 report 的结果，答辩时被问“这个甲烷氧化功能是哪些菌贡献的”，当场只能含糊过去，那场面至今记得。从那以后我每次跑 FAPROTAX 都强制走一遍“先自检输入表、再存 report、最后核对三类整理表”的流程，哪怕多花二十分钟。希望这套流程能帮到你。

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