简介：这份PDF讲义面向微生物生态与16S扩增子分析学习者，聚焦FAPROTAX在环境样本元素循环功能预测中的应用。内容围绕将物种分类与功能注释关联的思路展开：FAPROTAX依据已发表、经过验证的可培养菌文献手动构建，预测准确度较好但覆盖度有限；讲义说明如何由OTU分类表获得功能描述与样本对应的预测矩阵，并重点讨论碳、氢、氮、磷、硫等生物地球化学循环的注释结果解读。讲义还涉及功能丰度热图、环境因子与功能相关性泡泡图等呈现方式，以及ImageGP一键分析与Linux本地命令行两条实践路径，并借哺乳动物皮肤微生物、水稻氮利用等案例展示差异挖掘与基因/通路关联。包内为1个PDF文件，压缩包约2.47MB，篇幅集中、便于按主题查阅。已有1148人学习，适合希望理解FAPROTAX原理、结果形式与适用边界的读者用作专题参考。

1. 从一份 16S 表到功能矩阵：FAPROTAX 1.2.4 本地安装的难点在哪

第一次把这套东西装到本地，是因为一批沉积物样品的 16S 数据要讲碳氮硫循环的故事，在线一键版出图确实快，但默认跑的是 1.1 版数据库，注释比例一直上不去，审稿意见问的就是“这些功能你凭什么说是氮循环”。FAPROTAX 的定位很窄也很准：拿 OTU/ASV 的物种注释去查一张手工整理的文献表，预测微生物群落参与的元素循环功能，碳、氢、氮、磷、硫这些过程是它的主场。窄意味着覆盖度低，准意味着好解释。

装不好时的症状往往不是报错，而是跑完一看，几千条记录全是 leftovers，功能矩阵空得没法做热图，还得回头怀疑自己的注释库。这篇记录的是 1.2.4 版在 Linux 本地的完整落地过程：Python3 环境怎么搭、biom 依赖怎么装、带物种注释的输入表怎么造、collapse_table.py 的参数怎么给，以及跑挂在哪几个地方。

2. 装前把环境列清：Python3、numpy 与 biom 的版本账

FAPROTAX 本身没有二进制安装包，也没有 pip/conda 包，它就是一个 Python 脚本加一个纯文本数据库。这意味着环境问题全部落在解释器和两个依赖库上，装错版本的代价比装错库更隐蔽——脚本能启动，跑出来的结果是空的。

2.1 为什么这版脚本必须锁 Python3

1.2.4 是 2020 年 11 月 27 日更新的版本，数据库同步做了修订，脚本运行环境从 Python2 迁到了 Python3。这个时间点很关键，因为网上流传的大量教程还停留在 1.1 时代，命令长得几乎一样，但直接把老脚本丢进 Python3 解释器，第一行 print 就会抛SyntaxError: Missing parentheses in call to 'print'。更麻烦的是有些旧环境里同时躺着 1.1 和 1.2.4，路径写错一个字母，跑的是旧脚本配新库，日志不报错，功能矩阵却少一大半。

我一般会在下载目录里单独放一个FAPROTAX_1.2.4文件夹，不和其他版本混放，后面所有命令都写相对这个目录的路径，避免“以为跑的是新的”这种问题。版本这件事上没什么可折中的，要么全套 1.2.4，要么别用本地跑，去用在线版。

2.2 用 conda 起一个只装三件套的环境

不推荐往 base 环境里装。原因很实际：biom-format 对 numpy 版本有一定要求，而做扩增子的机器上往往已经装了 QIIME2 或者别的生信套件，版本冲突排查起来非常费时间。另起一个干净环境，装错了直接删掉重来，这就是后悔药。

# 新建一个 Python3 环境，名字随意，这里叫 faprotax conda create -n faprotax python=3.8 -y conda activate faprotax # 两个依赖：numpy 做矩阵运算，biom-format 提供 biom convert / add-metadata conda install -c conda-forge -c bioconda numpy biom-format -y # 确认当前解释器确实来自刚建的环境，而不是系统自带的 which python python -V

which python这一步别省。conda 激活失败、shell 里PYTHONPATH被别的东西污染、或者前一个环境没退干净，都会让python指向另一个解释器，表现出来就是“明明装了 numpy 还是 ModuleNotFoundError”。老教程里写的依赖是conda install numpy和conda install biom，biom这个包名在 bioconda 通道下现在通常对应biom-format，直接写biom有时能解析成别名、有时报 PackagesNotFoundError，统一写biom-format更稳。

如果机器上已经有一套 QIIME2 的 2020.6 环境，也可以直接conda activate qiime2-2020.6借用，它自带 numpy 和 biom，省一次安装。代价是这个环境比较重，一旦在里面动了包，可能影响后续的 QIIME2 分析，我个人更偏向单独建环境。

2.3 下载与首次自检

软件包官方只提供 zip 压缩包，下载地址里带%20这种空格转义，wget 遇到会截断 URL，所以整条地址必须加引号。

mkdir -p ~/biosoft && cd ~/biosoft # 地址整体加引号，否则 %20 会被 shell 和 wget 各自解释一遍 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 -q FAPROTAX_1.2.4.zip ls FAPROTAX_1.2.4

解压后应该能看到collapse_table.py和FAPROTAX.txt这两个核心文件，前者是转换脚本，后者就是那张文献整理的数据库。-c参数是断点续传，网络抖动时重跑不会从零开始。

自检只需要一条命令，不带任何参数，正常情况会弹出 usage 帮助：

python FAPROTAX_1.2.4/collapse_table.py

弹出用法说明就说明解释器和脚本本身没问题。这里如果抛的是 ImportError，那就是依赖没装进当前环境；抛的是 SyntaxError，那基本可以确定手里这个脚本不是 1.2.4 版。两种情况在第 5 章都有对应处理。

3. 输入表准备：把 txt 丰度表转成带物种注释的 biom

FAPROTAX 的输入不是普通的 txt 丰度表，而是带 observation metadata 的 biom 格式。这个设计是有道理的：biom 能把“每个 OTU 属于哪个物种”和“每个 OTU 在各样本里有多少”绑在同一个文件里，脚本解析时不会串行。所以中间必然多出两步转换，这两步也是最容易出格式问题的地方。

3.1 biom convert 之前，先把 OTU 表整理干净

输入的 txt 表要求是制表符分隔，第一行是#OTU ID加样本名，之后每行是一个 OTU/ASV，后面的数值是丰度。建议直接用抽平之后的 rare 表，因为 FAPROTAX 输出的是功能相对丰度矩阵，样本间测序深度差太多会让下游热图失真。

先做几项检查，别急着转：

# 看前几行的行尾和分隔符，^I 是制表符，出现 ^M 说明是 Windows 换行 head -3 otutab_rare.txt | cut -f1-6 | cat -A # 表头有多少个样本列 awk -F'\t' 'NR==1{print "样本列数:", NF-1}' otutab_rare.txt # 检查除表头外每行列数是否一致，输出不是 1 就说明有脏行 awk -F'\t' 'NR>1{print NF}' otutab_rare.txt | sort -u | wc -l

第三条命令的返回值是 1 才算干净。如果返回多个值，通常是文件末尾有空行、某一行的注释列多写了一个字段，或者混进了从 Excel 里复制出来的多余制表符。这类脏行会让 biom convert 直接失败，报错信息还很不具体，先在这里清掉比后面猜要省事。

物种注释文件是另一个独立的两列文本，第一列是 OTU/ASV 编号，第二列是注释字符串，两级之间用分号分隔：

head -2 taxonomy2.txt # ASV_1 d__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;...;g__Pseudomonas

注释库用哪个不强求，SILVA132 和 1.2.4 版的数据库配合得比较顺，而实操经验里 RDP 注释的结果匹配率往往最高。这两者不冲突，容易匹配是因为 RDP 在属级命名上更贴近 FAPROTAX 数据库里的写法。如果手上只有 Greengenes 的注释，也能跑，只是 leftovers 会多一些，心里要有预期。

3.2 添加物种注释：add-metadata 的两个必填项

转换分两条命令，第一条把 txt 转成 biom，第二条把注释挂上去。

# 1) txt 丰度表转 biom json 格式 biom convert -i otutab_rare.txt -o otutab_rare.biom \ --table-type="OTU table" --to-json # 2) 把物种注释作为 observation metadata 挂到每条记录上 biom add-metadata -i otutab_rare.biom \ --observation-metadata-fp taxonomy2.txt \ -o otutab_rare_tax.biom \ --sc-separated taxonomy \ --observation-header OTUID,taxonomy

参数里有两个地方写错就前功尽弃。--table-type="OTU table"是 biom 的枚举值，必须一字不差，写成OTU_table或者otu table都会报 Invalid table type。--sc-separated taxonomy告诉 biom 这一列是以分号分隔的多级结构，不加这个参数，整串注释会被当成一个不可分割的字符串，FAPROTAX 拿不到属名，匹配率直接掉到底。--observation-header OTUID,taxonomy声明的是注释文件本身的表头列名，必须和taxonomy2.txt第一行的实际写法一致，文件里叫OTUID命令里就不能写OTU_ID。

--to-json的作用是输出可读的 JSON 文本而不是 HDF5 二进制，好处是出问题时能直接 grep 看内容，代价是文件体积大一些。对几十个样本的扩增子数据来说完全无所谓，排查方便更重要。

3.3 转换结果怎么自查

转换完不要直接往下跑，先确认注释真的挂上去了。最轻量的办法是用 biom 自带的摘要命令看一眼：

biom summarize-table -i otutab_rare_tax.biom | head -20

它能告诉你样本数、每个样本的最小/最大深度。但摘要不显示 metadata，要确认注释挂载情况，用几行 Python 更直接：

import biom # 读入带注释的 biom 表 t = biom.load_table("otutab_rare_tax.biom") md = t.metadata(axis="observation") # 统计缺注释的条目 ids = t.ids(axis="observation") missing = [i for i in ids if not md.get(i, {}).get("taxonomy")] print("观测条目数:", len(ids)) print("缺注释条目数:", len(missing)) print("示例注释:", md[ids[0]].get("taxonomy"))

metadata(axis="observation")取的就是每条 OTU 上挂的属性字典，如果missing的数量等于总条目数，说明第二步 add-metadata 根本没生效，多半是--observation-header和文件表头不一致。如果 missing 是 0 但示例注释是一整串没有分号的乱码，说明--sc-separated没加或者列名写错。这两个检查点花不了两分钟，能挡掉后面大半的无效运行。

4. 跑通 collapse_table.py：参数逐条拆解与三张结果表

真正执行功能预测只有一条命令，但这条命令的参数没有默认值兜底，写错一个就是空结果。跑完之后还要把 report 文件拆成能进下游作图的三张表，这一步官方没给现成脚本的场合比较多，得自己补。

4.1 一条命令的七个参数

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

--collapse_by_metadata 'taxonomy'里的taxonomy必须和 3.2 步--sc-separated声明的列名完全一致，这是两条命令之间唯一的隐性契约。-r那个 report 文件很多人嫌它大就不生成，等到审稿人问“你这个氮固定功能是哪几个 ASV 贡献的”时，就只剩一条路——全部重跑。

执行完屏幕上会打印一行统计，形如：

Assigned 979 records to groups, 1924 records were leftovers

这行是整套流程里信息量最大的输出。前半句是成功匹配到至少一个功能组的记录数，后半句是没匹配上的。leftovers 占比高不一定是你操作错了，FAPROTAX 的数据库基于已发表的可培养菌文献构建，覆盖度天然受限，环境样本里大量未培养类群本来就查不到。但如果 leftovers 超过九成，或者 Assigned 直接是 0，那就是流程问题而不是覆盖度问题。

4.2 看懂 faprotax.txt 与 faprotax_report.txt

faprotax.txt是主结果，结构跟 OTU 表、KO 表一个样子：第一行是功能名加样本名，之后每行是一个功能分组，数值是该功能在各样本中的汇总丰度。这个表可以直接进 R 做热图，也可以按实验设计分组做差异比较。它的行数远少于 OTU 表，因为一个功能组会吸收多个物种的贡献。

faprotax_report.txt是明细，格式上是分块的：每一块以功能组名开头，下面跟着所有被归入该组的 OTU/ASV 行，每行带上完整的物种注释，块与块之间用--分隔。这个文件的价值在于可解释性。做出来一张“氮循环功能在组间显著富集”的图之后，能立刻回到 report 里查是哪些菌在撑这个信号，写讨论时这就是现成素材。我一般会把它和主结果一起归档，不单独删。

还有一点值得提前知道：中间文件里单个 OTU 的功能注释结果，除了做统计，还可以拿来做展示。比如挑几个关键功能，把贡献它们的物种列成表放进补充材料，比干巴巴一张热图有说服力得多。

4.3 从 report 拆出三张表

report 是给人看的，进下游还得转成矩阵。如果手头有配套的汇总脚本，一条命令就够：

# 先清洗：只保留“功能组名 + 其下 ASV 行”，去掉分隔线 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

grep 'ASV_' -B 1的意思是命中所有 ASV 行，同时带出它上面一行，也就是所属的功能组名；grep -v -P '^--$'把纯分隔线剔掉。这里有个硬性前提：你的 OTU 编号前缀必须是ASV_。如果表里用的是OTU_1这种命名，第一条 grep 会一条都匹配不到，清洗文件是空的，后面所有步骤都白做，所以执行完先wc -l看一下行数。

产物是三张：faprotax_report.mat是 OTU 层面的功能有无矩阵，适合画 Presence/Absence 热图；faprotax_report.otu_func是 OTU 到功能的对应列表；faprotax_report.func_otu反过来，按功能列出包含的 OTU。

手上没有这个 perl 脚本时，用 awk 加 R 也能自己造，逻辑不复杂。先生成两张表之间的关系：

awk -F'\t' 'BEGIN{OFS="\t"} /^ASV_/ {print F, $1; next} # 命中 ASV 行：输出"当前功能 + 该 ASV" /^--$/ {next} # 跳过组分隔线 {F = $0} # 其余行是功能组名，记到变量 F 里 ' faprotax_report.clean | sort -u > otu_func.tsv

F = $0是关键，它把功能组名暂存下来，供紧接着的 ASV 行使用；sort -u去掉同一个 ASV 在同一功能里重复出现的情况。拿到otu_func.tsv后，用 R 把同功能的 ASV 丰度加起来：

# otu_func.tsv 两列：功能名、ASV 编号 otu <- read.delim("otutab_rare.txt", row.names = 1, check.names = FALSE) map <- read.delim("otu_func.tsv", header = FALSE, col.names = c("func", "asv")) # 只保留原始丰度表里真实存在的 ASV，防止注释表里有孤儿条目 map <- map[map$asv %in% rownames(otu), ] # 按功能分组，把同组 ASV 的丰度逐样本相加 mat <- rowsum(otu[map$asv, , drop = FALSE], group = map$func) write.table(mat, "faprotax_mat.tsv", sep = "\t", quote = FALSE, col.names = NA)

rowsum的语义就是按 group 分组求和，正好对应“把一个功能组下所有物种的丰度合并”这个动作。drop = FALSE不能省，否则当某个功能只有一个 ASV 时，R 会自动降维成向量，rowsum直接报错。col.names = NA是为了让写出的表首行首列留空，符合丰度表在 R 里的常规读法。

5. 排错记录：FAPROTAX 本地跑挂的五个典型现场

这套东西的报错分两派：一派是环境派，报错信息明确但指向的包名不直观；另一派是数据派，不报错，结果全是空。后者更耗时间，下面每条都按现象、原因、解决三段记。

5.1 脚本一跑就 SyntaxError

现象：python FAPROTAX_1.2.4/collapse_table.py直接抛SyntaxError: Missing parentheses in call to 'print'，或者提示某一行invalid syntax。

原因：手里的脚本不是 1.2.4 版。1.1 及更早的脚本是 Python2 语法，print 还是语句而不是函数，Python3 解释器第一步就过不去。

解决：重新从官方地址下载 1.2.4 压缩包，别用别处转存的旧包。如果因为某些原因必须用旧数据库，那就老实建一个 Python2.7 环境跑旧脚本，但新旧混合的坑更多，不推荐。确认版本最省事的办法是解压后看一眼脚本开头，有没有from __future__之类的兼容声明，以及 print 有没有带括号。

5.2 依赖报错与 conda 包名

现象：ModuleNotFoundError: No module named 'numpy'，或者No module named 'biom'，但明明记得装过。

原因有三层。第一层是装进了别的环境，base 里装了但当前不是 base；第二层是PYTHONPATH或PATH里有系统 Python 抢在前面；第三层是包名写错，conda install biom有时候能解析成别名，有时候直接 PackagesNotFoundError。

解决：先which python和which pip确认解释器来源，再conda list | grep -i biom看当前环境里到底有没有。装的时候统一用conda install -c conda-forge -c bioconda numpy biom-format，包名写全，通道显式指定，比让 conda 自己猜要稳。

5.3 leftovers 高得离谱或 Assigned 为 0

现象：跑完日志是Assigned 0 records to groups, 2841 records were leftovers，或者 leftovers 占了九成以上，功能矩阵几乎没有内容。

原因排序：最常见的是--observation-header与注释文件表头不一致，导致 taxonomy 列压根没挂到 biom 上，这时 FAPROTAX 看到的每条记录都没有物种信息；其次是--collapse_by_metadata里的列名写错，跟 biom 里实际的 metadata 列对不上；再次是注释字符串里只有门纲水平、没有属名，FAPROTAX 靠属级名匹配，没有属名就查不到；最后才是注释库本身的覆盖度问题。

解决：按 3.3 的 Python 脚本先确认注释挂载率和内容形态，能定位到第一步。如果注释挂得没问题，那就换注释库重跑一次，SILVA132 和 1.2.4 版数据库配合较好，RDP 的匹配率在实操里最高。这一步建议在正式分析前用一小批样本先试，别等跑了全量才发现。

5.4 覆盖失败与输出路径问题

现象：重跑时提示输出文件已存在，脚本退出；或者-o指到一个不存在的子目录，直接报错。

原因：collapse_table.py默认不覆盖已有结果，这是防止误删的保护；而-o只接受文件路径，不会自动创建父目录。

解决：日常重跑统一加--force，批量分析时用带日期的输出目录，比如-o result_20220109/faprotax.txt之前先mkdir -p。另外--force会无提示覆盖，历史结果要留就提前改名，别指望它做备份。

5.5 biom 转换时行列错位

现象：biom convert报Invalid table type，或者转换成功但biom summarize-table显示的样本数比原表少。

原因：--table-type的值必须严格写成OTU table，大小写和空格都不能变；样本数变少通常是原 txt 里有某几行列数不一致，biom 解析时把异常行丢掉了。

解决：回到 3.1 的列数检查，把awk -F'\t' 'NR>1{print NF}' | sort -u的返回值压到只剩 1 个数字再转。另外注意文件编码，从 Excel 另存的 txt 有时带 BOM，第一列列名会变成\ufeff#OTU ID，grep 表头就匹配不上，用dos2unix或sed -i '1s/^\xEF\xBB\xBF//'清一下。

6. 装完怎么验证跑通：三条自检与一次结果核对

装完不等于跑通。脚本一旦返回 0 退出，界面安静得像成功了，但结果可能是空表。我习惯在交付前固定走三条检查，花不了五分钟，能挡掉绝大多数“看着跑完了但不能用”的情况。

第一条查总量守恒。功能矩阵的丰度总和应当小于等于输入 OTU 表的总和，因为只有被注释到的部分才会计入，但两者不应该差出一个数量级。

# 功能表总丰度：跳过表头行 awk -F'\t' 'NR>1{for(i=2;i<=NF;i++) s+=$i} END{printf "function total: %.2f\n", s}' faprotax.txt # OTU 表总丰度：跳过 # 开头的表头 awk -F'\t' '!/^#/{for(i=2;i<=NF;i++) s+=$i} END{printf "otu total: %.2f\n", s}' otutab_rare.txt

两个数字如果完全相等，反而是可疑的，说明有可能把未注释的记录也算了进去；功能总量占 OTU 总量的三到七成是比较常见的区间。具体表头格式以实际输出为准，必要时把NR>1换成!/^#/。

第二条查 leftover 比例，直接看运行时那行日志，或者回读 report 的规模。这个比例反映的是注释库与样本的契合度，不是脚本质量问题，但它是判断“结果能不能讲故事”的直接依据。

第三条做一次反向抽查，挑一个你关心的功能，看看撑起这个功能的物种在生物学上说不说得通。

# 提取某个功能组下的所有 ASV 明细，块与块之间以 -- 分隔 awk '/^nitrogen_fixation/{f=1; next} /^--$/{f=0} f' faprotax_report.txt | head -20

抽出来的物种如果全是固氮相关的已知类群，这条流程就是可信的；如果冒出来一堆完全不相干的属，那就得回头查注释文件是不是和 OTU 编号对错了行。这一步做过一次之后，后面再跑同类项目会踏实很多。

验证通过后再往下游走。faprotax.txt可以直接导进 STAMP 做组间差异，也可以按分组信息做功能丰度热图，或者拿去和环境因子做相关性分析。第一次跑通时我建议先用一小批样本、用已知注释库的公开数据试一遍，把这条链路走顺，再上自己的项目数据。

血泪经验是有一回赶进度，report 文件没生成就交了结果，后面想补一个“这些功能由哪些菌贡献”的说明，只能把整条流程从头重跑一遍，中间隔了两周，前面用的哪个版本注释库都记不清了。从那以后我每次跑 FAPROTAX，都强制走一遍“总量核对、leftover 记录、随机反查一个功能”这三步，并且把命令、版本号、注释库来源原样写进项目目录的 README，谁接手都不用猜。希望帮到你。

本文还有配套的精品资源，点击获取