16S扩增子FAPROTAX元素循环功能预测：流程与避坑

简介：围绕16S扩增子功能预测与元素循环解读的PDF讲义，聚焦FAPROTAX将OTU分类表与人工整理的功能注释数据库关联，用于海洋、湖泊等环境样本中碳、氢、氮、磷、硫循环的功能预测。内容适合具备一定微生物组分析基础的学生与研究者，用于理解FAPROTAX的数据库来源、预测边界、结果矩阵形式，以及功能丰度热图、环境因子相关性泡泡图等呈现方式。包内共1个PDF文件，约2.47MB，篇幅集中，便于作为专题课配套资料查阅。资源还结合哺乳动物与人类皮肤微生物、水稻氮利用等案例，讨论功能预测在差异挖掘与基因/代谢通路关联中的用法，并提及ImageGP一键分析与Linux本地分析两条路径的适用场景及依赖、输入表等常见问题。已有1148人学习，适合需要快速建立FAPROTAX知识框架、补充元素循环分析思路的读者。

1. 16S 数据做完物种注释之后：为什么我先补一层元素循环功能预测

扩增子跑完，OTU 表和物种注释摆在面前，紧接着一定会有人问一句：这些菌到底在干什么，碳氮磷硫循环里谁在干活。PICRUSt2 给的是一堆 KO 编号，Tax4Fun 给的是通路百分比，而 FAPROTAX 直接把 OTU 翻译成「硝化」「反硝化」「甲烷氧化」「硫酸盐还原」这种能直接写进结论的功能条目。它源自 2016 年那篇把海洋微生物组功能与分类解耦的工作，数据库靠人工从可培养菌文献整理，全部代码就一个collapse_table.py。覆盖度确实不高，但每条注释都能追到文献，做环境样本、农田土壤、湖泊沉积物的元素循环题目时，比一堆 KO 好解释得多。

2. 拆开 FAPROTAX 的黑匣子：一张二值表加一个 python 脚本

动手之前先把这套东西的内部结构看清楚，后面报错时才知道该往哪儿找。FAPROTAX 没有模型、没有隐层、没有任何机器学习成分，它就是一次查表加一次求和，把所有不确定性都压在「数据库里到底有没有你这个属」这一件事上。

2.1 数据库长什么样：菌与功能的二值映射

FAPROTAX.txt是纯文本、制表符分隔。第一列是分类单元名，绝大多数落在属一级，少量到种或科；后面每一列是一个功能条目，交叉位置用 1/0 标记该分类单元是否被文献支持具备这项功能。整张表本质是一个二部图：一个属可以同时挂上硝化、亚硝酸盐氧化、好氧氨氧化好几个标签，一个功能条目下面也拴着几十上百个属。脚本做的工作就是把这个二部图和你的 OTU 表做乘法——OTU 按注释归到属，属映射到功能，同一功能下所有 OTU 的丰度加起来。

这里有个必须提前知道的细节：功能条目的粒度差别很大。nitrification是粗线条的汇总条目，下面还挂着aerobic_ammonia_oxidation、aerobic_nitrite_oxidation这些细分条目，跑完的结果里两类会同时出现。如果你直接拿去画热图，同一个 OTU 会被数两遍，图上呈现的就是自相关，审稿人一眼能看出来。我一般会在出图前把汇总级别的条目单独拎出来，只留细分条目做展示。

还有一类条目要特别小心，像fermentation、chemoheterotrophy、aerobic_chemoheterotrophy这种范围极宽的，几乎每个样本都有、丰度还都不低，做组间差异分析永远不显著。它们不是错，只是信息量太低，属于「知道就行、别往结论里写」的那一档。

2.2 和 KO 通路、表型预测工具的分工

同一个 OTU 表可以喂给好几个预测工具，选哪个取决于你要回答什么问题。下面这张表是我自己在项目里做选型时用的对照，参数和输出形态都按实际跑出来的结果填：

选型逻辑其实很朴素：题目落在元素生物地球化学循环上，样本来自海洋、湖泊、沉积物、湿地、农田这类环境，FAPROTAX 是对口的；题目落在宿主相关样本上也能用，但覆盖度会明显掉一截，因为它的底子是可培养的环境菌。有研究把它用在哺乳动物和人类皮肤菌群的比较上，也有人把它和籼稻根系氮利用的表型关联起来，这两类工作能发出来，靠的不是覆盖度，而是「命中的那部分功能确实对得上已知机制」。

2.3 输入输出的形状：四个文件之间的关系

输入侧只有两样东西：一份不带注释的 OTU 丰度表（转成 biom 格式），和一份把 OTUID 对上 taxonomy 字符串的元数据文件。脚本默认取注释字符串的最后一级去和数据库匹配，通常就是属名，所以如果你的注释只到科，匹配率会掉得很难看。

输出侧是两个文件。faprotax.txt是功能 × 样本的丰度矩阵，行是功能条目，列是样本，数值是落在该功能下的所有 OTU 丰度之和。faprotax_report.txt是明细报告，每个功能条目下列出命中了哪些 OTU，末尾还有一行统计，格式类似Assigned 979 records to groups, 1924 records were leftovers。别小看这个 report，后面写讨论、做图注、解释「为什么某个功能突然冒起来」，全靠它把功能追回到具体序列。

中间还有一个容易被忽略的产物：biom 里挂上去的 observation metadata。它既不是输入也不是输出，但整条链路成不成就看它，第 5 章里一半的坑都出在这一步。

3. 本地环境搭起来：python3、numpy、biom 与 1.2.4 版脚本自检

这套脚本没有网络版依赖，纯本地跑，好处是数据不出内网、参数完全可控，坏处是环境得自己收拾。1.2.4 版是 2020 年 11 月底更新的，数据库做了扩充，脚本本身也从 python2 迁到了 python3，这个迁移是很多人翻车的源头。

3.1 下载与解包：先认清目录里有什么

压缩包解压后是一个以版本号命名的目录，里面至少要有三样东西：collapse_table.py主脚本、FAPROTAX.txt数据库、以及一份说明文件。三者的相对位置不要动，脚本按相对路径去找数据库是常见写法，挪了位置就得在命令行里显式指定。

# 假设压缩包已经下载到当前目录，先解压 unzip FAPROTAX_1.2.4.zip -d . # 看一眼目录结构，确认关键文件都在 ls -lh FAPROTAX_1.2.4/ # 预期看到 collapse_table.py、FAPROTAX.txt、README 之类

操作本身没什么难度，容易出问题的是「版本混用」：网上 2019 年以前流传的教程，命令里带的是 python2 的写法，直接抄到 python3 环境里会报语法错。看到报错信息里出现print语句相关的SyntaxError，基本可以断定你抄了老教程，或者手里拿的是老版本脚本。判断依据很简单，python --version输出 3.x，脚本也必须是 1.2.4 这一版，两者对齐才谈得上后面的事。

3.2 conda 环境与依赖安装：为什么不建议装在 base 里

依赖只有两个大头：numpy和biom。听起来简单，但这两者在不同渠道的包名和版本约束不一样，装在 base 环境里很容易和你已有的其他生信工具打架。我一般会单开一个环境，用完就扔，出问题直接删掉重建，比修依赖快得多。

# 新建一个干净的 python3 环境，名字随意 conda create -n faprotax python=3.9 -y conda activate faprotax # 装依赖，biom 在 conda 里的包名是 biom-format conda install -c conda-forge numpy biom-format -y # 确认解释器路径落在新建的环境下，避免装完还是用系统 python which python

conda create里指定 python 版本这一步别省，直接继承系统 python 有时候会带进一堆历史包袱。biom-format这个包名值得记一下，很多教程里写的是conda install biom，能不能装上取决于频道里有没有同名包，用-c conda-forge加全名是最稳的写法。

如果你机器上已经有一个装好 QIIME2 的环境，里面自带 numpy 和 biom-format，也可以直接激活它来跑，省一次安装。代价是 QIIME2 的依赖树比较深，偶尔会有版本锁死的情况，真出问题时排查成本比自己建环境高。选哪条路取决于你手头有没有现成的环境，没有就从零建。

3.3 自检：一行命令弹出帮助才算装好

装完不要急着上真实数据，先用不带参数的方式跑一次脚本。这类命令行工具通常在不给参数时打印帮助信息然后退出，能正常弹出帮助，就说明解释器、依赖、脚本三样都通了。

# 不带任何参数运行，预期打印 usage 帮助信息 python FAPROTAX_1.2.4/collapse_table.py

如果这一步报ModuleNotFoundError: No module named 'numpy'，说明当前激活的环境里没有 numpy，回到上一步执行conda install -c conda-forge numpy -y。报No module named 'biom'同理，装biom-format。报can't open file则是路径写错了，用ls确认脚本的实际位置。

还有一种情况是命令能跑但立刻退出、什么也不打印，这多半是脚本读参数时提前 return 了，把-h显式加上再试一次。自检这一步花三十秒，能省掉后面半小时的无效排查，属于典型的后悔药。

4. 从 otutab.txt 到 faprotax.txt：三步转换加一次预测

整条链路只有三段：文本表转 biom、把物种注释挂上去、跑 collapse。每一段都有它的格式约束，任何一环不合规，最坏的结果不是报错，而是静默地给你一张全零矩阵。

4.1 txt 转 biom json：两个格式约束先满足

原始 OTU 表通常是制表符分隔的纯文本，第一行是#OTU ID加一串样本名，下面是数字。转换前先做两件事：把最后一列如果有的话——通常是 taxonomy 字符串——单独拆出去存成另一个文件；确认所有单元格是纯数字，不带千分位逗号、不带空格、不带NA。

# 文本 OTU 表转 biom 的 json 格式 biom convert -i otutab_rare.txt -o otutab_rare.biom \ --table-type="OTU table" --to-json

--table-type这个参数必须写，写错了后面的工具链会按错误的类型解析；--to-json生成的是 json 版本而不是 HDF5，好处是纯文本、能直接head看内容、出问题时肉眼可查，对调试非常友好。输入表里的样本顺序和分组表顺序不要求一致，biom 靠列名匹配，但列名里的空格、中文、特殊符号建议提前清掉，这些在后续各种工具里都是雷。

转换完成后建议立刻验证一遍：用biom summarize-table -i otutab_rare.biom看每个样本的序列数、总 OTU 数是否和原表对得上。这一步不做，等到预测结果全零时再回头找，就得从头查起了。

4.2 add-metadata 把 taxonomy 挂上去：表头和分隔符是关键

这一步是整个流程里最容易翻车的地方。需要准备一个两列的元数据文件，第一列是 OTUID，必须和 biom 里的 observation ID 逐字符一致；第二列是 taxonomy 字符串，各级之间用分号分隔。

# 把物种注释作为 observation metadata 挂到 biom 上 biom add-metadata -i otutab_rare.biom \ --observation-metadata-fp taxonomy2.txt \ -o otutab_rare_tax.biom \ --sc-separated taxonomy \ --observation-header OTUID,taxonomy

--sc-separated taxonomy告诉程序这一列内部是按分号切分的多级结构，切完会在 biom 里存成列表；--observation-header OTUID,taxonomy是元数据文件的表头声明，这两个词必须一字不差。写成OTU_ID、#OTUID、OTUID,taxonomy_name都不会报错，但挂上去的 taxonomy 是空的，直到你跑完预测看到一张全零的矩阵才会发现。我见过不止一次有人在这上面耗掉一整天。

挂完之后同样要验：biom table -i otutab_rare_tax.biom --observations随便看几条，确认每个 OTU 后面都跟着一段注释字符串。另外注意注释字符串里最后一级是不是属名，中间有没有混进g__、unclassified、norank这类前缀或占位词，这些在匹配阶段都是纯噪声。

4.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

跑完之后终端会打印一行统计，形如Assigned 979 records to groups, 1924 records were leftovers。前半句是成功映射到功能的记录数，后半句是没匹配上的。这个比例就是你的注释率，它决定了后面所有分析的可信度，看到数字的第一反应应该是判断它合不合理，而不是直接进入画图环节。环境样本的注释率能到 40%~60% 算正常，宿主样本掉到 20% 以下也不奇怪，如果环境样本只有个位数百分比，那一定是前面某一步出了问题。

4.4 读 report 判断质量：先看统计行，再看命中列表

faprotax_report.txt的价值被严重低估。它按功能条目分块，每个块下列出命中的 OTU 及其归属，末尾带统计。拿它做三件事：判断注释率是否合理、看某个关键功能是被哪些 OTU 撑起来的、以及在讨论部分解释「这个功能的升高来自哪几个属」。

如果要把 report 整理成可用的矩阵，常见做法是先把每个功能块的首行和它下面的 ASV 行筛出来，去掉分隔线，再按功能分组展开。下面的脚本做的是同一件事，用 python 解析比记一堆 grep 参数更省心：

import collections # faprotax_report.clean 是筛掉分隔线后的报告：功能名行 + 缩进的 OTU 行交替出现 func2otu = collections.OrderedDict() cur = None with open("faprotax_report.clean") as fh: for line in fh: line = line.rstrip("\n") if not line.strip(): continue if line.startswith("ASV_") or line.startswith("OTU_"): if cur: # 当前功能下追加一条命中记录 func2otu[cur].append(line.split("\t")[0]) else: cur = line.strip() # 遇到功能名行，切换当前功能 func2otu.setdefault(cur, []) # 输出「功能 - 命中OTU数」的速查表，用来判断哪些功能是靠单个OTU撑起来的 for f, otus in func2otu.items(): print(f, len(otus))

这段脚本的输出可以直接当质控表用。命中 OTU 数只有一两个的功能条目，在组间差异里极不稳定，换个抽平深度可能就消失了，写结论时要留余地。反过来，命中几十个 OTU 的功能更稳健，但往往也就是那些宽泛条目，差异不显著。两头都要具体分析，不能一刀切。

5. 五条踩坑记录：注释率、报错与结果对不上号的排查路径

下面这些是我和周围人真实踩过、并且在别人数据上也复现过的问题。每条按「现象→原因→解决」写，遇到对得上的直接照方抓药。

5.1 现象：Assigned 只有一百多条，leftovers 两三千

原因通常有三个。最常见的是 taxonomy 只注到科，脚本取最后一级拿到的是科名，数据库里根本没有这一级；其次是属名带了一堆前缀和占位符，比如g__、unclassified、norank、uncultured，这些字符串拿去匹配当然一条都中不了；第三是生境本身特殊，数据库的底子是可培养的环境菌，你的样本来自极端环境或宿主肠道，库里确实没有对应条目。

解决顺序建议从便宜到贵：先写个脚本把所有注释字符串的最后一级抽出来，统计有多少是unclassified之类的占位词；再把g__前缀统一去掉、属名首字母大写规范化；如果注释只到科，换一套能到属的参考库重新注释，经验上 132 版的 SILVA 和这套数据库配合得比较好；都做完注释率还是低，那就是生境问题，接受它并在方法里如实写明覆盖度。

5.2 现象：跑脚本报 ModuleNotFoundError，明明刚装过

原因几乎都是解释器和环境不匹配。装依赖时装在 A 环境，跑脚本时激活的是 B 环境，或者干脆用了系统自带的 python。conda 环境切换不彻底是很常见的操作失误。

解决：运行which python和which conda，确认前者路径里带envs/你的环境名。如果不对，conda deactivate之后重新conda activate，再不行就直接用环境里的绝对路径调用解释器，类似~/miniconda3/envs/faprotax/bin/python collapse_table.py。这条路径写死在命令行里，虽然难看，但绝对不会跑错环境。

5.3 现象：预测跑完了，faprotax.txt 全是 0

原因八成在上一步的 add-metadata。要么--observation-header写错导致 taxonomy 列没挂上，要么元数据文件里的 OTUID 和 biom 里的不一致——多了个#、少了ASV_前缀、大小写不同、行尾带\r，任何一种都能让匹配全军覆没，而且整个过程不报错。

解决：先用biom table -i otutab_rare_tax.biom --observations | head看前几条，确认 taxonomy 字段真的有内容。如果为空，回去查表头；如果有内容但预测结果还是全零，就对比两边的 ID，我一般会cut -f1 taxonomy2.txt | sort > a.txt、把 biom 的 observation ID 也导出来排序，diff一下，差异一目了然。Windows 上生成的文件记得dos2unix。

5.4 现象：换一套物种注释重跑，两次结果相关性只有 0.6

原因不是脚本不稳定，而是输入变了。不同参考库对同一个 OTU 的属级归属本来就可能不同，能命中数据库的 OTU 集合随之改变，功能丰度自然整体漂移。如果两次流程里还有一次做了抽平、另一次没做，差异会更大。

解决：同一批数据从头到尾锁死一套注释流程，中途不换库；确实需要比较时，先把两边都转成相对丰度再算相关；更严谨的做法是把两次 report 里命中的 OTU 列表取交集，只在交集上比较功能丰度，把「注释差异」这个变量控制住。

5.5 现象：某个功能在所有样本里丰度几乎一模一样

原因有两种。一种是这个功能命中的 OTU 极少且丰度很低，抽平之后被随机性主导，看起来就像常数；另一种是命中的是那几个无处不在的宽泛条目，本来就没什么组间变化。

解决：先用 4.4 那段脚本看该功能命中了几个 OTU，只有一两个的直接从图里删掉，别硬画；如果是宽泛条目，保留在表里但在正文里不展开。低丰度 OTU 在抽平前先做一次过滤，比如保留在至少 10% 样本中出现的 OTU，能明显减少这类噪声。

6. 拿到功能矩阵之后：归一化、追到 OTU、用实测指标兜底

faprotax.txt里是绝对丰度加和，直接跨样本比 reads 数没有意义，测序深度和抽平方式都会影响结果。第一步永远是转相对丰度，把每个样本的功能丰度除以该样本的总和，再乘 100：

import pandas as pd df = pd.read_csv("faprotax.txt", sep="\t", index_col=0) # 行是功能，列是样本 df = df.div(df.sum(axis=0), axis=1) * 100 # 按样本列归一化成百分比 df = df.loc[df.sum(axis=1) > 0.01] # 剔掉几乎不出现的功能 df.to_csv("faprotax_rel.txt", sep="\t")

归一化之后才谈得上热图、气泡图、组间差异。做环境因子相关的泡泡图时，功能这一侧建议只保留与碳氮磷硫直接相关的细分条目，把宽泛条目挡在外面，否则图上会挤满一堆不显著的点。

第二步是把功能追回到 OTU。这一步很多人跳过，结果讨论部分只能写「硝化功能在组间显著升高」，说不清是谁升高的。用第 4.4 节的脚本把 report 拆成「功能—OTU」对应表，再和 OTU 表的物种注释一join，就能落到具体的属。我做氮循环相关题目时，这一步几乎必做，因为它能帮你判断观察到的升高是真信号还是某个低丰度 OTU 的随机波动。

第三步，也是最能提升说服力的一步：用实测指标做锚点。功能预测终究是预测，如果条件允许，同一批样本补测几个关键基因的 qPCR，硝化补 amoA、反硝化补 nirK/nirS、甲烷氧化补 pmoA，把预测丰度和实测拷贝数做个相关。方向一致，后面的机制讨论才站得住；方向相反，那就是个明确的信号，说明这套数据库在你的生境里不适用，趁早换工具比硬写强。没有 qPCR 条件的，退一步用环境因子的实测值做相关也行，铵态氮、硝态氮、有机碳这些理化指标很多项目本来就测了。

血泪教训是：我曾经拿着一份注释率只有 12% 的结果写了一整节机制讨论，投稿后被审稿人一句「你这条通路只由一个未培养属的单个 OTU 支撑」直接打回。从那以后我每跑完一次collapse_table.py，都强制先把 Assigned 比例、每个功能的命中 OTU 数、以及关键功能的属级归属这三样抄进分析日志，达不到预期就回头查输入，绝不带着可疑的结果往下走。希望帮到你。

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