简介这份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 上以及数据库文件。输出侧是两张表结构差异很大用途也完全不同。文件维度内容下游用途faprotax.txt功能 × 样本每个功能在各样本的丰度行是功能条目统计检验、热图、排序分析、与环境因子做相关faprotax_report.txt功能 → OTU 列表每个功能由哪些 OTU 贡献逐条列出追溯证据、讨论具体菌属、写结果解释功能矩阵的行名是功能描述不是 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 而不是 biomconda 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-typeOTU 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是两列制表符分隔的纯文本第一行是表头后面每行一个 OTUOTUID 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_1grep 模式必须跟着改否则清洗结果是空文件二是 report 生成时的功能名行不含制表符缩进如果版本差异导致格式变了先head -30看一眼再改正则。4.3 三类整理表与下游分析接口整理脚本会吐出三张表各自解决的问题不一样输出文件结构用途faprotax_report.matOTU × 功能 的 0/1 有无矩阵画 OTU 与功能的对应热图看哪些菌群在同一功能下聚集faprotax_report.otu_funcOTU 明细表逐条核对某 OTU 挂了哪些功能判断注释是否合理faprotax_report.func_otu功能 → OTU 汇总写结果时的证据来源直接引用到讨论里下游统计从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 都进了 leftoversAssigned 为 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 脚本。答辩时被问到某个数字怎么来的能当场打开文件说服力远高于口头解释。第二件做一张验收自查表逐项打勾再往下走检查项合格标准不合格时看哪里OTU 表是否抽平各样本序列数一致抽平日志与 OTU 表列和species 注释挂载biom 转回文本后 taxonomy 列非空add-metadata 的列名与路径注释命中率与自己数据同类型的文献水平相当异常低要解释report 里的 Assigned/leftovers功能矩阵维度行数为功能条目数列数与样本数一致collapse_table.py 的 -o 输出三类整理表mat / otu_func / func_otu 齐全且非空清洗步骤的前缀与分隔行可视化可复现换台机器用同一脚本能出同一张图相对丰度转换与过滤阈值第三件准备一条能讲三层的证据链。第一层是功能层哪个元素循环功能在组间显著差异校正方法是什么。第二层是来源层翻 report 说明这个功能主要由哪几个属贡献。第三层是文献层这些属在已发表的可培养研究中被验证过什么代谢特征。三层讲完一个“基于预测”的结果就变成了有据可查的推断这是这类分析在答辩中最有分量的部分。还有个小技巧把faprotax_report.otu_func里几条典型 OTU 的功能注释单独拎出来做一张小表附在附录。评审老师翻到附录看到你能追到单个 OTU 挂了几条功能、有没有明显矛盾的注释对你数据处理细致程度的判断会完全不一样。我早年做过一版没留 report 的结果答辩时被问“这个甲烷氧化功能是哪些菌贡献的”当场只能含糊过去那场面至今记得。从那以后我每次跑 FAPROTAX 都强制走一遍“先自检输入表、再存 report、最后核对三类整理表”的流程哪怕多花二十分钟。希望这套流程能帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/25 13:26:57 · 阅读 1.2k 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 上以及数据库文件。输出侧是两张表结构差异很大用途也完全不同。文件维度内容下游用途faprotax.txt功能 × 样本每个功能在各样本的丰度行是功能条目统计检验、热图、排序分析、与环境因子做相关faprotax_report.txt功能 → OTU 列表每个功能由哪些 OTU 贡献逐条列出追溯证据、讨论具体菌属、写结果解释功能矩阵的行名是功能描述不是 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 而不是 biomconda 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-typeOTU 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是两列制表符分隔的纯文本第一行是表头后面每行一个 OTUOTUID 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_1grep 模式必须跟着改否则清洗结果是空文件二是 report 生成时的功能名行不含制表符缩进如果版本差异导致格式变了先head -30看一眼再改正则。4.3 三类整理表与下游分析接口整理脚本会吐出三张表各自解决的问题不一样输出文件结构用途faprotax_report.matOTU × 功能 的 0/1 有无矩阵画 OTU 与功能的对应热图看哪些菌群在同一功能下聚集faprotax_report.otu_funcOTU 明细表逐条核对某 OTU 挂了哪些功能判断注释是否合理faprotax_report.func_otu功能 → OTU 汇总写结果时的证据来源直接引用到讨论里下游统计从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 都进了 leftoversAssigned 为 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 脚本。答辩时被问到某个数字怎么来的能当场打开文件说服力远高于口头解释。第二件做一张验收自查表逐项打勾再往下走检查项合格标准不合格时看哪里OTU 表是否抽平各样本序列数一致抽平日志与 OTU 表列和species 注释挂载biom 转回文本后 taxonomy 列非空add-metadata 的列名与路径注释命中率与自己数据同类型的文献水平相当异常低要解释report 里的 Assigned/leftovers功能矩阵维度行数为功能条目数列数与样本数一致collapse_table.py 的 -o 输出三类整理表mat / otu_func / func_otu 齐全且非空清洗步骤的前缀与分隔行可视化可复现换台机器用同一脚本能出同一张图相对丰度转换与过滤阈值第三件准备一条能讲三层的证据链。第一层是功能层哪个元素循环功能在组间显著差异校正方法是什么。第二层是来源层翻 report 说明这个功能主要由哪几个属贡献。第三层是文献层这些属在已发表的可培养研究中被验证过什么代谢特征。三层讲完一个“基于预测”的结果就变成了有据可查的推断这是这类分析在答辩中最有分量的部分。还有个小技巧把faprotax_report.otu_func里几条典型 OTU 的功能注释单独拎出来做一张小表附在附录。评审老师翻到附录看到你能追到单个 OTU 挂了几条功能、有没有明显矛盾的注释对你数据处理细致程度的判断会完全不一样。我早年做过一版没留 report 的结果答辩时被问“这个甲烷氧化功能是哪些菌贡献的”当场只能含糊过去那场面至今记得。从那以后我每次跑 FAPROTAX 都强制走一遍“先自检输入表、再存 report、最后核对三类整理表”的流程哪怕多花二十分钟。希望这套流程能帮到你。本文还有配套的精品资源点击获取