简介这是一份微生物群落功能预测专题资料聚焦FAPROTAX在环境样本元素循环分析中的应用适合从事16S扩增子测序、宏基因组功能注释的科研人员与生信学习者。资料以PDF形式系统讲解FAPROTAX的数据库原理、OTU表转换流程、结果矩阵解读并结合Science、PNAS、Nature Biotechnology等经典案例展示其在海洋、湖泊、皮肤微生物及水稻氮利用研究中的具体用法同时涉及ImageGP一键分析平台与Linux本地命令行操作对比帮助读者按需选择分析路径。压缩包为单份PDF大小2.47MB轻量易得已有1147人浏览学习。对希望从物种分类快速推断碳、氢、氮、磷、硫等元素循环功能的用户这份材料能提供从数据库下载、软件安装到结果可视化的完整思路是微生物生态学功能预测入门与进阶的实用参考。1. 为什么16S功能预测要单独给FAPROTAX留一章手里握着一张16S扩增子的OTU表想知道样本里的微生物到底在做什么、谁在参与氮循环、谁在做硫氧化很多人的第一反应是打开PICRUSt把KO通路跑一遍。但PICRUSt跑出来的是“推断基因组”的结果靠的是参考基因组和学习模型相当于把黑匣子当骰子用。FAPROTAX是另一条路它不推断只查表把OTU的物种注释直接映射到手工整理的文献功能数据库上命中就报没命中就不报。准确度有文献兜底代价是覆盖度低。这篇笔记把FAPROTAX的下载、环境配置、输入格式、核心命令collapse_table.py的每个参数以及我实际跑数据踩过的坑拆开讲一遍适合手上有OTU表、关心碳氢氮磷硫元素循环、想把功能预测从玄学变成可追溯的人。2. FAPROTAX原理与定位手工文献库如何做到“准而窄”2.1 它不是预测是逐级查表FAPROTAX的全称是Functional Annotation of Prokaryotic Taxa2016年发表在Science上作者团队做的事情很朴素把已发表的、经过培养验证的微生物文献里的物种-功能关系一条条手工整理成数据库。每一行记录是一个功能分组比如“nitrate_reduction”“sulfate_respiration”“methanotrophy”后面跟着一组能执行这个功能的物种分类名。运行的时候collapse_table.py把你输入的OTU表逐行读出来取出每个OTU的物种注释从域到属、到种逐级和数据库比对命中了就把这个OTU的丰度累加进对应的功能分组。所以FAPROTAX本质上是查表不是统计推断。这带来两个直接结果第一数据库里没有的功能它不会给你编一个出来这就保证了“报出来的结果都有文献出处”第二它的覆盖度完全取决于数据库收录了多少物种环境样本里大量未培养的类群文献里没写过自然就匹配不上。我最早用的时候以为是自己OTU表格式错了后来才想明白覆盖度低就是它的设计边界不是bug。这个逻辑和PICRUSt截然不同。PICRUSt的思路是用16S序列位置去比对参考基因组然后推断整个基因组的KO通路。看起来什么都能预测但每一条通路都是“推断”出来的里面有多少是真实存在的需要通过宏基因组或者qPCR去验证。FAPROTAX是“只报有把握的”宁可少报报出来的东西基本靠谱。2.2 四种常用16S功能预测工具怎么选做16S功能预测市面上高频出现的是PICRUSt/PICRUSt2、FAPROTAX、Tax4Fun和bugbase这四类经常有人纠结选哪个。我把它们放在一张表里对比结论会更直观工具原理输入要求输出内容优势局限PICRUSt/PICRUSt216S序列比对参考基因组推断OTU/ASV表推荐用Greengenes或SILVA注释KO通路、EC酶丰度覆盖全基因功能层面推断误差难以评估宿主相关样本偏差大FAPROTAX物种注释映射手工文献库带物种注释的biom格式OTU表元素循环等功能分组丰度矩阵准确度高可追溯文献覆盖度低只覆盖可培养类群Tax4Fun16S序列映射到SILVA参考基因组OTU/ASV表SILVA注释KO、KEGG通路百分比基于SILVA兼容性好输出百分比绝对值意义弱bugbase基于Greengenes的表型推断OTU表Greengenes注释革兰氏、氧利用、移动元件等表型结果直观、易解释表型特征有限覆盖面窄我的选择习惯是如果问题聚焦在元素循环比如“水稻根际的氮循环功能是否受基因型影响”“湖泊沉积物中硫氧化功能在空间上怎么分布”直接用FAPROTAX因为它输出的功能分组恰好就是“nitrogen_fixation”“sulfate_respiration”这类生态过程名词不需要再去KO表里一个个查通路名。如果问题是“整体代谢潜力怎么样”那就PICRUSt2打底FAPROTAX做交叉验证。2.3 FAPROTAX的适用边界环境样本是主场FAPROTAX最典型的使用场景是海洋、湖泊、土壤、沉积物这类环境样本。原因是这些生态位里的元素循环过程研究得早、文献积累厚数据库里的功能分组基本都是围绕碳、氢、氮、磷、硫循环建的。做水稻根际微生物研究时有人发现籼稻中富集的9个高丰度通路里有6个和氮相关用的就是FAPROTAX加上差异分析这种“功能分组-环境因子-样本分组”的联动是它的强项。但如果你做的是人的肠道菌群、皮肤菌群这类宿主相关样本要谨慎。宿主相关微生物里大量物种是未培养的或者说没有功能验证文献FAPROTAX的匹配率会明显下降。一篇关于哺乳动物和人类皮肤微生物的研究里用过它但那是跨物种比较功能组成差异用的也是数据库中真正命中的部分不是拿全矩阵硬解释。数据库版本方面2020年11月27日更新到了1.2.4版脚本运行环境更新为Python3。我个人的经验是1.2.4版注释比例比老版本高不少注释后建议优先配SILVA132两者兼容性更好。但也有人说RDP注释出来的匹配率最高这点我在这章先不展开第4章避坑部分细说。3. 本地跑通FAPROTAX从OTU表到功能矩阵的完整命令链3.1 环境准备conda装numpy和biom一个都不能少FAPROTAX本体是一个Python脚本加一个数据库文件不需要编译但依赖numpy和biom这两个Python包。我习惯先建一个独立环境不和qiime2主环境混在一起避免依赖冲突把环境搞坏。# 下载FAPROTAX 1.2.4压缩包并解压 wget -c https://pages.uoregon.edu/slouca/LoucaLab/archive/FAPROTAX/SECTION_Download/MODULE_Downloads/CLASS_Latest%20release/UNIT_FAPROTAX_1.2.4/FAPROTAX_1.2.4.zip unzip FAPROTAX_1.2.4.zip # 创建独立Python3环境并激活 conda create -n faprotax python3.7 -y conda activate faprotax # 安装依赖numpy和biom-format注意包名是biom-format不是biom conda install -c conda-forge numpy biom-format -y # 测试脚本是否可运行弹出帮助信息即正常 python FAPROTAX_1.2.4/collapse_table.py这里有个容易翻车的地方biom包在conda里的名字是biom-format如果直接conda install biom会报包找不到或者装到同名但是完全无关的包。我最初就是因为这个卡了半小时后来发现装完后命令行工具叫biom但包名必须写biom-format。另外如果你已经有qiime2环境也可以直接在里面跑conda activate qiime2-2020.6之后不需要额外装东西因为qiime2自带numpy和biom。但我不建议在qiime2环境里手动conda install numpy版本很容易被覆盖导致整个环境出问题。3.2 输入文件制作txt转biom再挂物种注释collapse_table.py只认biom格式的OTU表所以第一步是把常规的txt格式OTU表转成biom json格式然后给每个OTU挂上物种注释。# 第一步txt格式OTU表转biom json格式 biom convert -i otutab_rare.txt -o otutab_rare.biom --table-typeOTU table --to-json # 第二步为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第一条命令里otutab_rare.txt是稀有化后的OTU表行是OTU/ASV列是样本第一列是OTU ID。--table-typeOTU table是biom转换时的必填参数确保输出的biom文件被识别为OTU丰度表--to-json指定输出为json格式。第二条命令是给biom表加观测元数据也就是物种注释。注意taxonomy2.txt这个文件必须是两列第一列OTU ID第二列物种注释表头是OTUID和taxonomy列之间用Tab分隔。--sc-separated taxonomy的意思是我告诉biom物种注释这一列是用分号分隔的多个分类级别有的版本写法是--sc-separated ;如果你遇到参数识别报错优先检查biom版本然后改成带分隔符的写法。最关键的点物种注释的层级顺序必须从域开始格式类似“Bacteria;Proteobacteria;Gammaproteobacteria;Burkholderiales;Burkholderiaceae;Cupriavidus”。FAPROTAX是逐级匹配的从宽泛分类往精细分类搜注释不完整或者顺序乱了会直接影响匹配结果。3.3 核心转换collapse_table.py逐参数拆解做好了带物种注释的biom表就可以运行核心脚本了。这一步的输出就是功能丰度矩阵。# 运行FAPROTAX功能预测 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格式OTU表必须带物种注释-g指定FAPROTAX数据库文件位置注意路径里一定包含FAPROTAX.txt这个文件名--collapse_by_metadata taxonomy指定按biom观测元数据里的哪一列进行聚合这里填的就是3.2里加注释时的列名taxonomy如果列名不一致这里会直接报错。-v是verbose模式运行时会打印详细信息比如“Assigned 979 records to groups, 1924 records were leftovers”。这句话特别关键Assigned数量表示有多少个OTU被成功分配到了功能分组leftovers表示没匹配上的OTU数量。--force表示输出文件存在时直接覆盖不加这个参数而输出文件已存在的话脚本会拒绝执行。-o指定输出的功能丰度矩阵文件名-r指定report文件名。report文件记录的是每个功能分组具体由哪些OTU贡献的调试时必看。运行完看一眼终端输出如果Assigned比例过低比如只有几百个说明你的物种注释格式或数据库版本有问题不要急着往下走先回头检查注释文件。3.4 输出文件解读faprotax.txt和faprotax_report.txtfaprotax.txt是主结果格式和OTU表一样行是功能分组列是样本每个单元格是对应功能在该样本中的丰度值。可以直接拿去做热图、CCA、PERMANOVA或者差异分析。功能分组名字本身就带着生态过程语义比如“nitrogen_fixation”“nitrate_reduction”“methanotrophy”读起来非常直观。faprotax_report.txt是明细展开看能发现每个功能分组下面挂了哪些OTU。比如“nitrate_reduction”这个功能分组列出了所有贡献该功能的OTU的物种注释和ID。这个文件常用于两件事一是追溯某个功能变化是由哪些微生物驱动的二是做第5章说的OTU-功能有无矩阵。实际使用中faprotax.txt可以直接和样本分组信息组合做差异分析比如用STAMP或者R里的edgeR/limma处理。report文件通常要再做一步清洗把OTU级别的对应关系拆出来这一步放在第5章讲。4. FAPROTAX避坑指南覆盖度低、匹配为零等5个常见坑4.1 运行成功但输出功能表只有一行标题没有任何数据现象脚本正常跑完终端没有报错但是faprotax.txt里只有表头和样本名一行功能数据都没有。原因最常见的两个一是biom表里的物种注释列名和--collapse_by_metadata指定的名字不一致二是注释文件里的分类层级分隔符不是分号FAPROTAX逐级匹配时整个字符串拿去比比不上就是零命中。解决先biom head -i otutab_rare_tax.biom看注解的列名是什么再head -5 taxonomy2.txt看注释格式。确认列名是taxonomy且分隔符是分号。我自己的习惯是在运行collapse_table.py之前先用biom summarize-tables -i otutab_rare_tax.biom看一眼观测数确认OTU数量还在排除输入文件本身为空的情况。4.2 Assigned数量极低leftovers占绝大多数现象终端打印“Assigned 300 records to groups, 4000 records were leftovers”匹配率不到10%。原因这种情况通常不是格式错误而是物种注释精度不够或者样本类型超出了数据库覆盖范围。FAPROTAX数据库收录的是有功能文献证据的属和种你注释到科甚至门级别匹配率自然会崩。另外宿主相关样本里大量未培养类群数据库里根本没有对应条目。解决优先换RDP训练集重新注释OTU这是经验上匹配率最高的方案改用SILVA时尽量选SILVA132配合FAPROTAX 1.2.4。如果重新注释后匹配率还是低不要强行解释低丰度功能分组只关注那些有足够OTU支持的功能分组。做环境样本时匹配率30%左右是可以正常用的样本间比较的是相对丰度趋势不是绝对定量。4.3 conda安装biom-format后命令行还是找不到biom现象conda install -c conda-forge biom-format顺利装完但在终端输入biom提示command not found。原因环境没激活或者把包装到了base环境但当前在别的环境里执行命令。我见过更隐蔽的情况是conda的bin目录没有正确加到PATH。解决先conda activate faprotax再执行which biom确认路径。如果还是找不到用python -m biom的方式调用biom的命令行入口也可以这样执行。另外检查一下是不是真的在目标环境里conda list | grep biom看输出。4.4 --force不加重跑直接失败加了好结果被覆盖现象第一次跑成功修改参数后想重跑脚本报错拒绝运行加了--force之后重跑发现上一次的产出被覆盖想找回来看不到了。原因collapse_table.py默认输出文件若已存在就拒绝执行这是防止误覆盖的设计--force虽然能强制覆盖但不会帮你保留旧文件。解决不用固定输出名每次带上日期或版本后缀。比如-o faprotax_20250112.txt -r faprotax_report_20250112.txt养成保留过程文件的习惯。这个教训我是付过时间成本的重跑一次要等很久结果旧数据没了新数据感觉也不太对只能重新来一轮。4.5 换了SILVA注释后匹配率反而比RDP低现象同时用RDP和SILVA注释同一批OTU分别跑FAPROTAX发现SILVA注释的匹配率明显低于RDP。原因FAPROTAX数据库里的物种名以经典的培养物种拉丁名为主。SILVA注释里经常出现“Candidatus”“uncultured”这类候选分类或未培养标记这些名字在文献库里根本没有对应条目所以匹配时全被跳过。解决在SILVA和RDP之间做FAPROTAX我会优先选RDP。如果必须用SILVA先对注释做清洗把“uncultured”“Candidatus”等前缀去掉或者把注释截断到属水平能提高一部分匹配率。但这里也有个度把注释强行截到属水平会丢失种水平的有用信息FAPROTAX里不少功能细分是种级别才有的截太狠反而会漏掉这部分。4.6 快速排查清单如果跑出来的结果怎么看怎么不对劲按下面顺序排查能省很多时间排查事项检查方法正常状态脚本环境运行python FAPROTAX_1.2.4/collapse_table.py弹出帮助信息输入文件格式biom summarize-tables -i otutab_rare_tax.biomOTU数量与原始表一致物种注释列名biom head -i otutab_rare_tax.biom列名包含taxonomy数据库路径ls FAPROTAX_1.2.4/FAPROTAX.txt文件存在且非空匹配情况查看终端-v输出及report文件Assigned占OTU总数比例15%5. 进阶用法让功能矩阵经得起审稿人追问的三件事5.1 从report文件拆出OTU功能有无矩阵faprotax_report.txt除了用于追溯还能拆成一个OTU级别的功能有无矩阵用来做网络图或者功能富集排序。做法是先清洗report文件把每个OTU对应的功能分组关系提取出来。# 筛选ASV/OTU注释行及其前面的功能分组标题行 grep ASV_ -B 1 faprotax_report.txt | grep -v -P ^--$ faprotax_report.clean # 用脚本整理为三类表格有无矩阵、OTU功能列表、功能OTU列表 perl faprotax_report_sum.pl -i faprotax_report.clean -o faprotax_report生成的faprotax_report.mat就是OTU×功能的有无矩阵1表示该OTU参与了对应功能。这个矩阵里能直接看出哪些OTU是多功能型的、哪些功能组只由少数OTU支撑做共现网络之前先用它过滤低支持度的功能分组可以避免无意义的连线。5.2 画热图和气泡图之前先查三个数字功能丰度矩阵做可视化之前我先查三个数字一是总丰度是否做了归一化输入OTU表使用的是稀有化后的表跑出来的功能丰度才有样本间可比性二是每个功能分组对应的OTU数量支持OTU数量太少的功能组波动大画热图时容易被极端值带偏建议过滤掉只有1-2个OTU支撑的分组三是功能分组丰度的分布如果某个分组丰度占比超过80%说明群落功能高度单一这时候更适合做分组间差异分析而不是做聚类热图。环境因子与功能丰度的相关性泡泡图也类似先对功能丰度做标准化再算Spearman相关避免高丰度功能主导相关系数。5.3 和PICRUSt做交叉验证的取舍FAPROTAX的命中是基于文献证据但覆盖度低PICRUSt覆盖广但每一条都是推断。我自己的习惯是同一个数据集同时跑FAPROTAX和PICRUSt2当两个工具在氮循环或者碳循环相关指标上方向一致时这个结论我才敢写进文章里。比如FAPROTAX里nitrate_reduction丰度在某种处理下显著上升PICRUSt里对应的narG、napA相关KO通路也同步上升两个证据互为印证。如果出现冲突优先采信FAPROTAX的方向然后从report里找出驱动变化的具体OTU用这些OTU的丰度和功能关系做进一步解释。从那以后我每次做16S功能预测都强制走一遍“FAPROTAX出方向、PICRUSt做验证、report文件追溯OTU”的完整流程不再让一个工具单独说了算。希望帮到你。本文还有配套的精品资源点击获取