简介一份聚焦微生物群落功能预测的专题 PDF围绕 FAPROTAX 在 16S 扩增子分析中的应用展开面向微生物生态、环境科学方向的研究生与科研人员帮助理解环境样本中碳、氢、氮、磷、硫等元素循环的功能注释思路也适合作为组会汇报与课题设计时的概念参考。内容从数据库的文献构建依据讲起说明 OTU 分类表与功能注释的关联方式并延伸到功能—样本矩阵、分组功能丰度热图、环境因子相关性泡泡图等结果呈现同时对比 ImageGP 一键分析与 Linux 本地运行两种实现路径也提示了预测准确度与覆盖度的适用边界。压缩包仅含 1 个 PDF 文件大小约 2.47MB轻量便于查阅。已有 1148 人学习读者可据此建立 FAPROTAX 预测流程的整体框架并通过哺乳动物皮肤微生物、水稻氮利用等案例把握功能预测在差异挖掘与机制关联中的用法。1. 拿到 16S 结果之后功能预测这一步怎么落地扩增子测序做到最后ASV 表和物种注释都齐了合作方或者审稿人往往不会满足于“这批样本里拟杆菌门多、变形菌门少”他们要的是——这些菌在干哪些活氮循环、硫循环、碳降解到底有没有差异。FAPROTAXFunctional 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 python3.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 转 biomJSON 格式注意表类型写 OTU table biom convert -i otutab_rare.txt -o otutab_rare.biom \ --table-typeOTU 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,taxonomybiom 不会报错但会把 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进分析之前要做标准化。因为功能组丰度受测序深度影响样本间直接比大小没有意义。常见做法是转成相对丰度每列除以该列总和或者用 TSStotal 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 NR3{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表附上审稿人问“这个功能由谁贡献”一查就有。下面这张表是我常用的自检清单检查项看什么合格线经验值匹配覆盖度Assigned / 总记录数环境样本 30%~60%低于 15% 需复查注释抽平一致性三次抽平后功能组排序相关性Spearman 0.9功能组冗余每个功能组对应的 ASV 数多数功能组应由多个 ASV 贡献单一 ASV 主导需说明差异结果稳健性换分组检验方法Wilcoxon / t 检验主要结论方向一致外部验证与理化指标或 qPCR 的相关方向关键功能组方向不矛盾还有一个具体技巧值得单独说做组间比较前先把faprotax.txt里在所有样本中相对丰度都低于 0.01% 的功能组过滤掉。这些低丰度功能组本身由一两条 ASV 撑着稍一抽平就消失留着只会制造假阳性让多重检验校正后的结果更难看。过滤阈值别拍脑袋按你数据的总丰度分布来定写进方法里就行。最后提一句和实测数据的交叉验证。如果项目里同时测了总氮、铵态氮、pH把硝化、反硝化、固氮这几个功能组单独拎出来做散点或相关分析方向对得上故事就闭环了方向对不上先别急着删数据回头查抽平和注释库很多时候问题就出在这两处。从那以后我每次跑 FAPROTAX都强制先看 Assigned 比例、再抽平三次自检、最后反查一遍 report 里的贡献 OTU这三步不走完后面的图一张都不出。希望帮到你。本文还有配套的精品资源点击获取