简介FAPROTAXFunctional Annotation of Prokaryotic Taxa功能预测课件聚焦微生物群落16S扩增子数据的功能注释面向环境微生物组、海洋与湖泊生态等研究方向的科研人员。该预测方法基于已发表、可培养菌文献构建物种-功能映射数据库通过Python脚本将OTU分类表转换为功能预测矩阵重点注释碳、氢、氮、磷、硫等元素循环通路。包内为1个PDF文件压缩包大小2.47MB内容系统讲解FAPROTAX的数据库原理、软件安装、biom格式OTU表生成与转换、collapse_table.py运行流程并结合Science、PNAS和Nature Biotechnology等期刊案例展示功能丰度热图、环境因子相关性泡泡图以及差异菌与氮循环功能挖掘等典型应用。已有1147人学习适合希望快速上手功能预测的读者参考既能理解从OTU表到功能注释的完整分析思路也能对比ImageGP在线平台与本地Linux命令行两种实现方式为后续统计分析和论文图表绘制提供方法支撑。1. FAPROTAX1 功能预测从 OTU 表到元素循环的一张功能名片手头刚跑完一批某湿地土壤的 16S 扩增子数据分类学注释表躺在文件夹里下一步最想知道的不是「谁在」而是「它们到底在帮土壤转什么元素」。碳、氮、磷、硫这些元素循环功能靠测序注释本身看不出来得做功能预测。标题里的「31功能预测2-元素循环FAPROTAX1」就是分析流程里常见的一环用 FAPROTAX 把 OTU/ASV 表映射成一张元素循环功能表回答「这个群落里谁在固氮、谁在反硝化、谁在产甲烷」。它解决的是物种注释到功能注释的最后一公里适合正做到功能预测环节、想快速拿到可发表热图或差异检验结果的从业者。FAPROTAX 的核心是查表比 PICRUSt2 那种基因推断方案轻得多代价是只覆盖「有文献依据」的原核类群理解这一点后面所有的坑都好解释了。2. FAPROTAX 的映射逻辑与选型为什么元素循环这块它比 PICRUSt2 顺手2.1 一条 OTU 是怎么变成「硝酸盐还原」功能的FAPROTAX 全称 Functional Annotation of Prokaryotic Taxa名字里就写明白了它不是预测基因而是给原核生物分类单元做功能注释。它的工作方式非常直白社区维护者把已发表培养实验、富集培养和同位素示踪文献里验证过的「分类单元—功能表型」关系整理成一张映射表常见文件名 FAPROTAX.txt每一行是某个属或种对应的生态功能比如某属能够进行 nitrate_reduction某科参与 methanogenesis。运行脚本时它把你 OTU 表里每个 OTU 的物种谱系拿出来逐级去映射表里匹配命中一个功能就把这个 OTU 的丰度累加到对应功能列上。这个机制决定了两个重要特性。第一它是「功能表型」注释而不是「基因潜力」注释。硝化、反硝化、硫酸盐呼吸这类过程往往由多基因、多步骤协作完成单看一两个功能基因并不等于过程真的在发生而 FAPROTAX 直接落到「文献证明这类群能干这件事」所以对元素循环过程类功能的注释比单基因推断更贴近实测通量。第二它是透明的查表不是黑匣子。每个 OTU 被分配了什么功能可以一路回溯到映射表里那一条文献证据这对审稿人来说非常友好也是我在流程里坚持用它兜底的原因。映射的颗粒度基本卡在属和种两级。数据库里绝大多数条目写的是「属名功能」或「种名功能」少数覆盖到科。这也解释了为什么注释只到门纲目水平的 OTU 大概率命中不了映射表里根本没有门纲目的 key。做元素循环方向的样本建议上游尽力把物种注释做到属级以上否则 FAPROTAX 这步基本白跑。注释级别这一条是整个流程里性价比最高的投入点。2.2 同赛道工具对比PICRUSt2、Tax4Fun、MicFunction 谁更适合元素循环功能预测的常用工具不止 FAPROTAX 一个选型决定了后面结果解读的口径。我把几个主流方案放在一张表里对比结论会清晰很多。工具推断方式输入要求输出颗粒度元素循环适配度主要坑FAPROTAX文献查表映射OTU表物种注释功能表型固氮、反硝化等高C/N/S 循环是强项覆盖窄依赖注释分辨率PICRUSt216S拷贝数参考基因组推断OTU表系统发育信息KEGG/EC/MetaCyc 基因丰度中需要自己筛通路低丰度 OTU 预测噪声大Tax4Fun16S 与参考数据库最近邻匹配OTU表参考数据库KEGG KO中参考数据库版本极度敏感MicFunction数据库注释映射宏基因组或扩增子多组功能类别中参数多上手慢我的习惯是PICRUSt2 和 FAPROTAX 都跑前者看基因层面的潜力有没有 narG、nirK 这些基因后者看生态过程层面的结论这个样本有没有反硝化功能群。两者不是替代关系而是互相校验。PICRUSt2 强在能推断未培养类群的基因弱在低丰度 OTU 的预测容易飘FAPROTAX 强在结论可靠、可回溯弱在覆盖范围完全取决于数据库里收录了多少类群。标题里的「元素循环FAPROTAX1」通常就是指流程里专门用 FAPROTAX 做元素循环功能注释的这个模块负责出生态过程结论。选型上还有一条经验如果你的样本来自极端环境或者大量未培养类群主导的系统比如某些深海沉积物、特殊工业废水FAPROTAX 的命中率会明显下降这时候要以 PICRUSt2 为主、FAPROTAX 为辅并在文章里说明注释覆盖率的局限。反过来如果样本来土壤、湿地、活性污泥这类研究充分的生态系统FAPROTAX 的命中率通常能到一半以上元素循环结论用它就够硬了。判断依据很简单跑完看被注释上的 OTU 丰度占总量多少低于三成就得换主次策略。3. 跑通 FAPROTAX1 最小流程输入转换、脚本调用与三个关键参数3.1 输入格式整理与转换BIOM 表、TSV 表和注释文件怎么对齐FAPROTAX 的输入就三样东西丰度表、物种注释表、映射数据库文件。丰度表常见来源是 QIIME2 的 feature-table.biom注释表常见来源是 taxonomy.qza 导出的 TSV。这两个文件里的 OTU ID 必须完全一致一个字符都不能差这是第一步对齐要解决的核心问题。我一般先把 BIOM 转成 TSV再单独整理注释文件让两个文件都只有两列核心信息OTU ID 和对应的值。# 从 QIIME2 导出丰度表feature-table.biom 换成你的实际文件名 qiime tools export --input-path feature-table.biom --output-path exported # 把 BIOM 转成 TSVbiom 2.x 语法 biom convert -i exported/feature-table.biom -o otu_table.txt --to-tsv # 查看转换结果第一列是 OTU ID其余列是样本名和丰度 head -n 5 otu_table.txt # 物种注释单独整理成两列OTU_ID 完整谱系 # 常见 taxonomy.tsv 是 QIIME2 导出的格式保留前两列即可 cut -f1,2 taxonomy.tsv tax_for_faprotax.txt head -n 5 tax_for_faprotax.txt这里有一个容易踩的格式细节QIIME2 导出的 BIOM 转成 TSV 后文件开头会带一行# Constructed from biom file之类的注释行FAPROTAX 脚本读取时对这类注释行的容忍度不同。我见过某开发者因为没删这行导致整个输入被判定为空表。稳妥做法是在跑脚本前用tail -n 2 otu_table.txt otu_table_clean.txt去掉首行或者直接手写一个小循环把第一列是#的行过滤掉。别小看这步它排掉了输入格式里最隐蔽的一个坑。注释文件的格式也有讲究。FAPROTAX 映射表里存的是类似Nitrosomonas europaea这样的双名法名称你的注释谱系如果写的是d__Bacteria;p__Proteobacteria;c__Betaproteobacteria;o__Nitrosomonadales;f__Nitrosomonadaceae;g__Nitrosomonas;s__europaea脚本需要按分隔符拆开之后逐级匹配。不同流程导出的注释分隔符不一样常见的是;分号加空格或;只有分号。拿到文件先统一成分号分隔别在匹配阶段让分隔符不一致把命中断了。这一步直接用批量替换搞定五分钟能省后面一小时的排查。3.2 跑通 FAPROTAX 脚本与三个关键参数输入文件就绪后接下来是调用脚本。常见做法是从 FAPROTAX 发布页下载脚本和数据库文件脚本通常是run_FAPROTAX.py映射数据库是FAPROTAX.txt。不同发行版的参数名略有差异我第一次跑的时候吃过参数不匹配的亏所以现在固定先python run_FAPROTAX.py --help看一眼再填参数。下面这个调用是常见形态以你下载版本的帮助输出为准。# 常见调用形态-i 输入丰度表-o 输出前缀-t 注释文件 # 具体参数名以 python run_FAPROTAX.py --help 为准 python run_FAPROTAX.py \ -i otu_table_clean.txt \ -o faprotax_output \ -t tax_for_faprotax.txt # 跑完检查输出文件列表 ls -lh faprotax_output*运行结束后会生成两个核心结果文件一个带 grouped 标记一个带 collapsed 标记。grouped 文件保留的是「每个 OTU × 每个功能」的明细行同一个功能可能因为多个 OTU 都命中而出现多行collapsed 文件把同一功能的所有 OTU 丰度合并成一行是下游分析真正要用的表。我第一次跑完直接拿 grouped 当结果表用加和功能丰度时怎么算都不对后来才知道要用 collapsed。后面第 4 章会细说这两个文件的结构差异。脚本有三个参数值得单独说。第一个是列过滤阈值控制某功能至少要出现在多少个样本里才保留调高能滤掉只在个别样本里出现一次的功能避免结果表列数爆炸第二个是行丰度阈值低于该丰度的 OTU 不参与映射适合理掉极低丰度的噪声 OTU第三个是 verbose 输出打开后可以看每个 OTU 的匹配过程排查「为什么这个 OTU 没被注释上」时非常有用。这三个参数没有通用最优值我的习惯是先不开任何过滤跑一遍看注释率再根据列数和噪声决定阈值。如果官方脚本的环境依赖装起来费劲或者你想彻底搞清楚匹配逻辑可以写一个最小匹配器做验证。下面这段 Python 脚本不做正式注释只用来验证「我的 OTU 表和注释文件能不能被正确读取、匹配逻辑有没有问题」跑通了再用官方脚本出正式结果。# 最小 FAPROTAX 匹配验证脚本确认输入文件可读、注释能对上 import pandas as pd # 读入丰度表第一列是 OTU ID abun pd.read_csv(otu_table_clean.txt, sep\t, index_col0) # 读入注释表两列otu, lineage tax pd.read_csv(tax_for_faprotax.txt, sep\t, headerNone, names[otu, lineage]) tax_dict dict(zip(tax[otu], tax[lineage].fillna())) # 用几个典型元素循环功能做关键词验证 func_keys { nitrogen_fixation: [nitrogen_fixation, nif], nitrate_reduction: [nitrate_reduction, nar], sulfate_respiration: [sulfate_respiration, dsr], } for func, kws in func_keys.items(): hit [oid for oid in abun.index if oid in tax_dict and any(k in tax_dict[oid] for k in kws)] print(func, 命中 OTU 数:, len(hit), 总丰度:, abun.loc[hit].sum().sum() if hit else 0)这段脚本的逻辑是把每个 OTU 的谱系字符串拿出来检查里面是否包含某个功能的典型关键词统计命中的 OTU 数量和总丰度。注意它只是验证文件通路的占位脚本不是 FAPROTAX 的正式匹配逻辑正式结果一定以官方脚本输出为准。它的价值在于把「输入文件格式错误」和「OTU 本身没被数据库收录」这两类问题分开——如果这里都读不到数据那官方脚本大概率也会静默失败。参数说明就一条func_keys 里的关键词按你自己关心的循环功能改验证完这步再放心跑正式脚本。4. 输出表拆解哪些列是碳氮硫磷循环哪些是噪声4.1 读懂 grouped 与 collapsed一行功能数据的来龙去脉跑完 FAPROTAX 之后结果表里会出现一堆似懂非懂的功能名光看列名就能劝退一半人。grouped 文件的每一行代表「某个 OTU 被分配到了某个功能」所以同一功能会重复出现多行每行对应不同的 OTU 或分类单元。collapsed 文件则是把这些行按功能名合并一个功能只剩一行各样本的丰度是该功能下所有 OTU 丰度之和。我做下游分析只用 collapsedgrouped 只在排查单 OTU 去向时才打开看。collapsed 表里常见的前几列是功能名和分类学层级信息后面每一列对应一个样本。不同版本输出列的排列略有差别但核心信息一致。下面列几个我在湿地和活性污泥样本里最常见的元素循环功能列名方便你对号入座元素循环典型功能列名生态意义氮循环nitrogen_fixation生物固氮把 N2 转成氨氮循环nitrate_reduction硝酸盐还原第一步还原氮循环nitrification硝化氨氧化到硝酸盐氮循环denitrification 系列反硝化N2O/N2 释放碳循环methanogenesis产甲烷严格厌氧发酵末端碳循环methanotrophy甲烷氧化甲烷汇碳循环cellulolysis纤维素分解有机碳降解硫循环sulfate_respiration硫酸盐呼吸厌氧硫还原碳循环fermentation发酵兼性厌氧有机物降解解读这些列之前必须记住一个规则功能列之间不是互斥的。同一个 OTU 可能同时被注释为 nitrate_reduction 和 fermentation因为它在不同条件下确实能干这两件事映射表按文献证据把两个功能都挂上了。所以千万不要拿功能表做整表归一化然后声称各功能占比加起来是 100%——它们是一张多标签表不是单标签分类。我见过有人直接把功能矩阵丢进 PCA 画图第一轴全是总丰度差异后来才发现是没处理多标签和包含关系。另外要留意那些「万金油」类型的功能比如 aerobic_chemoheterotrophy、fermentation。这类功能覆盖面极广几乎每个异养菌都能挂上在表里往往占掉一大半丰度。它们不是没用但在做元素循环特异性分析时会把真正的循环信号稀释掉。我的习惯是先看总表再手动筛出氮、硫、甲烷相关的功能子集做重点分析把好氧异养这类大而全的功能单独放在一边不掺进循环强度计算里。4.2 按元素循环归类功能并做丰度归一化把 collapsed 表里的功能列按元素循环归类这一步决定了下游所有统计的可靠性。我一般写一小段 Python用关键词把功能名匹配到氮、硫、碳三个大类顺便把磷循环的情况摸清楚——注意 FAPROTAX 对磷循环的覆盖明显弱于碳氮硫因为数据库的条目主要围绕生物地球化学循环里的 C/N/S 和甲烷展开磷相关的条目寥寥无几硬要用它做磷循环结论不够稳这时候得换宏基因组方案补。下面这段是常用的归类与归一化代码。# collapsed 表按元素循环归类并做样本内相对丰度归一化 import pandas as pd import re collapsed pd.read_csv(faprotax_output.collapsed, sep\t, index_col0) # 只保留数值矩阵忽略功能描述类列 num collapsed.select_dtypes(includenumber) cycles { nitrogen: rnitr|denitri|fixation|ammon|anammox, sulfur: rsulfur|sulfate|sulfite|sulfide|thiosulf, carbon_methane: rmethan|cellul|ferment|hydrocarbon|methylotroph, } cycle_abund {} for name, pattern in cycles.items(): cols [c for c in num.columns if re.search(pattern, c, re.I)] if cols: cycle_abund[name] num[cols].sum(axis1) cycle_df pd.DataFrame(cycle_abund, indexnum.index) # 相对丰度每个样本内三类循环功能占该样本总捕获功能的比例 rel cycle_df.div(cycle_df.sum(axis1), axis0) * 100 rel.to_csv(cycle_relative_abundance.csv) print(rel.round(1))代码逻辑分三步先用正则把功能列名分到氮、硫、碳甲烷三个循环桶再把桶内所有功能丰度按行加和得到每个样本的循环功能总丰度最后除以样本内三类循环之和做归一化。这样得到的是一张三列的样本×循环强度表直接就能画箱线图做分组比较。参数说明里最值得留意的是正则的覆盖范围我写的 pattern 尽量覆盖了常见功能名后缀但不同 FAPROTAX 数据库版本的功能命名可能有差异比如有的版本写 nitrous_oxide_denitrification有的简写成 nitrite_denitrification跑完先打印num.columns人工核对一遍缺了就补进 pattern。归一化这一步还有个容易被忽略的点要不要除以总丰度取决于你的科学问题。如果样本间测序深度差异大除以每个样本的全部功能丰度能消除深度影响但代价是「高活性样本的整体功能丰度」信息丢失。如果问题是「氮循环和硫循环谁更强」用我上面的相对方案合适如果问题是「处理组氮循环功能绝对丰度是否显著更高」那应该用原始丰度加协变量做 DESeq2 或 edgeR。两种口径在 FAPROTAX 输出上都成立只是一定要在方法部分写清楚用的是哪一种审稿人最烦看到相对丰度结论却用绝对丰度的逻辑解释。5. FAPROTAX 实战避坑五条能让元素循环结果翻车的细节5.1 现象所有 OTU 都没有被注释上输出表几乎是空的这是我见过最多的情况新手一跑发现结果表只有表头立刻以为工具坏了。真正原因大概率在注释分辨率FAPROTAX 映射表里的 key 到属和种而你的注释只到门纲目匹配自然全军覆没。另一个常见原因是 OTU ID 对不上比如丰度表里是OTU_1注释表里是OTU001肉眼看着像同一个程序眼里完全不同。解决方法是先head检查注释文件里有没有属级名称再用grep -c g__ tax_for_faprotax.txt统计带属注释的行数占比低于一半就先补注释再做功能预测。ID 一致性用comm -12 (cut -f1 otu_table_clean.txt) (cut -f1 tax_for_faprotax.txt)直接比对能快速算出两个文件共有多少 OTU。5.2 现象同一份数据两次跑结果功能数量不一样有次我在两个环境里各跑了一遍同样的 OTU 表输出功能列数差了十多个一度以为脚本有随机性。后来比对发现两个环境里放的 FAPROTAX.txt 不是同一个版本数据库条目更新过映射结果自然不同。原因很简单FAPROTAX 是持续维护的数据库不同时间下载的文件内容不一样。解决方法是锁定数据库版本记录下载日期、文件大小和校验值最好把 FAPROTAX.txt 固定放在项目目录里随分析一起归档而不是每次现下。这样结果才能复现文章方法部分也能写清楚用的是哪个版本的映射库。5.3 现象一大半丰度堆在 aerobic_chemoheterotrophy元素循环信号被淹没土壤样本经常出现这种局面表里最大的列永远是 aerobic_chemoheterotrophy氮循环、硫循环相关的列占比都很小画出来的热图几乎只有一种颜色。原因是这类「好氧化能异养」功能被数据库挂在大量常见类群名下覆盖面广命中率高。这不是 bug是数据库设计使然。解决方法是不要拿全表做热图或排序先按第 4.2 节的方法筛出氮、硫、甲烷等功能子集再分析解释结论时把范围限定在特定循环过程不要用「群落功能多样性」这种试图覆盖一切的表述。我在文章里通常写「在本研究中关注的氮循环功能类群中」明确的边界反而是加分项。5.4 现象grouped 表手工加总的结果跟 collapsed 表对不上有同学拿着 grouped 表自己按功能名加总发现怎么都对不上 collapsed 的值怀疑脚本算错了。原因有两层一是 grouped 表里同一 OTU 可能同时命中同一功能的不同映射条目数据库里冗余条目合并方式不同二是手工加总时容易把描述列算进去或者某些版本会对功能行做额外归一化。解决方法是不要在 grouped 上做任何统计汇总collapsed 就是官方帮你合并好的最终结果下游分析一律以它为准。如果非要自己合并只对数值列求和处理先把select_dtypes(includenumber)之外的所有列丢掉。5.5 现象FAPROTAX 和 PICRUSt2 对同一个样本的氮循环结论对不上这是工具联用时最常出现的困惑PICRUSt2 里硝酸盐还原基因丰度很高FAPROTAX 的 nitrate_reduction 却几乎没有于是怀疑某个工具出错。实际两者口径完全不同——PICRUSt2 推的是基因潜力哪怕细菌只是携带基因没表达也算FAPROTAX 推的是文献支持的功能表型只有数据库里明确收录了这类群才给注释。如果样本由大量未培养类群主导PICRUSt2 靠系统发育推断能给出基因预测FAPROTAX 却只能静默跳过差异就出来了。解决方法是不要用两个工具互相证伪而是定义清楚各自回答的问题基因潜力用 PICRUSt2生态过程用 FAPROTAX差异过大时回到物种注释看是否以未培养类群为主并在文章里如实报告 FAPROTAX 注释覆盖率。6. 把功能预测接进下游元素循环差异检验与热图可视化的一个顺手做法功能表跑出来只是中间产物真正的结论要靠差异检验和可视化收口。我一般先把功能表聚合到元素循环大类再做两件事一是按分组做 perMANOVA 看群落功能结构是否有显著差异二是画热图看具体哪些循环功能在组间分化。下面这段 R 代码是我在模拟项目X里反复用的一套流程改动最少、出图最快。# 读入 collapsed 表和分组信息 library(pheatmap) library(vegan) dat - read.delim(faprotax_output.collapsed, row.names 1, check.names FALSE) meta - read.delim(metadata.txt, row.names 1, check.names FALSE) # 1) 筛选氮循环相关功能子集避免被好氧异养淹没 n_funcs - grep(nitr|denitri|fixation|ammon, colnames(dat), value TRUE, ignore.case TRUE) dat_n - dat[, n_funcs, drop FALSE] # 2) hellinger 转换后再算 Bray-Curtis降低零丰度主导 dat_hel - decostand(dat_n, hellinger) dist_bc - vegdist(dat_hel, method bray) # 3) 分组差异检验 adonis2(dist_bc ~ Group, data meta, permutations 999) # 4) 相对丰度热图 dat_rel - as.data.frame(t(t(dat_n) / rowSums(dat_n) * 100)) pheatmap(dat_rel, scale column, clustering_method ward.D2, annotation_col meta[Group], fontsize_row 8)这段代码里最值得说的是两步处理。第一步筛功能子集只留氮循环相关的列这比全表直接热图干净得多第二步做 hellinger 转换再算 Bray-Curtis是因为 FAPROTAX 输出有大量零值和偏态分布原始丰度算距离会被高丰度功能主导hellinger 转换后再算距离能稳定很多。热图里我加了scalecolumn让每个功能自己的组间差异显现出来而不是让高丰度功能霸占整个色阶。perMANOVA 的置换次数 999 在正式分析里可以提到 9999结论更稳。一个让我印象很深的教训有次拿功能表直接怼进 PCA第一主成分全是总丰度差异组间差异怎么画都不显著后来才意识到是没先聚合大类、没做归一化白白折腾了两天。从那以后我养成了固定习惯——功能表到手先聚合到循环大类再决定归一化口径最后才碰可视化。做元素循环功能预测真正值钱的部分不是你跑通了多少脚本而是你能不能把一张多标签的复杂表压缩成审稿人能直接看懂的循环强度结论。希望这篇对你有用后面跑数据遇到玄学问题先回头查注释分辨率和数据库版本八成问题都在那里。本文还有配套的精品资源点击获取