简介围绕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 1.2.4具体元素循环功能条目人工整理的可培养菌文献碳氮磷硫循环里谁在动覆盖度低注释率常在 20%~60%PICRUSt2KO / EC / 代谢通路基因组数据库外推通路层面的整体代谢潜力中间步骤多结果难落到具体过程Tax4Fun2KO 通路百分比16S 参考序列库快速通路比较准确率高度依赖参考库覆盖表型预测类工具革兰氏、氧气利用、致病性已有表型数据库群落表型的粗略画像表型条目有限粒度粗选型逻辑其实很朴素题目落在元素生物地球化学循环上样本来自海洋、湖泊、沉积物、湿地、农田这类环境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 python3.9 -y conda activate faprotax # 装依赖biom 在 conda 里的包名是 biom-format conda install -c conda-forge numpy biom-format -y # 确认解释器路径落在新建的环境下避免装完还是用系统 python which pythonconda 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。报cant 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-typeOTU 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 的参数逐个说清参数不多但每一个都直接决定结果对不对值得逐条过一遍参数含义常见取值与注意事项-i输入 biom必须是挂了 taxonomy 的那个版本-g数据库文件指向FAPROTAX.txt路径写错会直接报错或零注释-o功能丰度矩阵行是功能、列是样本-r明细报告每个功能下命中的 OTU 列表与统计行--collapse_by_metadata用哪一列做归并写taxonomy必须和上一步的列名一致-v打印过程信息建议开着能看到进度和警告--force覆盖已有输出重跑时必加否则遇到同名文件直接退出对应的完整命令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 数没有意义测序深度和抽平方式都会影响结果。第一步永远是转相对丰度把每个样本的功能丰度除以该样本的总和再乘 100import pandas as pd df pd.read_csv(faprotax.txt, sep\t, index_col0) # 行是功能列是样本 df df.div(df.sum(axis0), axis1) * 100 # 按样本列归一化成百分比 df df.loc[df.sum(axis1) 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 数、以及关键功能的属级归属这三样抄进分析日志达不到预期就回头查输入绝不带着可疑的结果往下走。希望帮到你。本文还有配套的精品资源点击获取 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 1.2.4具体元素循环功能条目人工整理的可培养菌文献碳氮磷硫循环里谁在动覆盖度低注释率常在 20%~60%PICRUSt2KO / EC / 代谢通路基因组数据库外推通路层面的整体代谢潜力中间步骤多结果难落到具体过程Tax4Fun2KO 通路百分比16S 参考序列库快速通路比较准确率高度依赖参考库覆盖表型预测类工具革兰氏、氧气利用、致病性已有表型数据库群落表型的粗略画像表型条目有限粒度粗选型逻辑其实很朴素题目落在元素生物地球化学循环上样本来自海洋、湖泊、沉积物、湿地、农田这类环境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 python3.9 -y conda activate faprotax # 装依赖biom 在 conda 里的包名是 biom-format conda install -c conda-forge numpy biom-format -y # 确认解释器路径落在新建的环境下避免装完还是用系统 python which pythonconda 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。报cant 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-typeOTU 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 的参数逐个说清参数不多但每一个都直接决定结果对不对值得逐条过一遍参数含义常见取值与注意事项-i输入 biom必须是挂了 taxonomy 的那个版本-g数据库文件指向FAPROTAX.txt路径写错会直接报错或零注释-o功能丰度矩阵行是功能、列是样本-r明细报告每个功能下命中的 OTU 列表与统计行--collapse_by_metadata用哪一列做归并写taxonomy必须和上一步的列名一致-v打印过程信息建议开着能看到进度和警告--force覆盖已有输出重跑时必加否则遇到同名文件直接退出对应的完整命令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 数没有意义测序深度和抽平方式都会影响结果。第一步永远是转相对丰度把每个样本的功能丰度除以该样本的总和再乘 100import pandas as pd df pd.read_csv(faprotax.txt, sep\t, index_col0) # 行是功能列是样本 df df.div(df.sum(axis0), axis1) * 100 # 按样本列归一化成百分比 df df.loc[df.sum(axis1) 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 数、以及关键功能的属级归属这三样抄进分析日志达不到预期就回头查输入绝不带着可疑的结果往下走。希望帮到你。本文还有配套的精品资源点击获取