简介这份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 为什么这版脚本必须锁 Python31.2.4 是 2020 年 11 月 27 日更新的版本数据库同步做了修订脚本运行环境从 Python2 迁到了 Python3。这个时间点很关键因为网上流传的大量教程还停留在 1.1 时代命令长得几乎一样但直接把老脚本丢进 Python3 解释器第一行 print 就会抛SyntaxError: Missing parentheses in call to print。更麻烦的是有些旧环境里同时躺着 1.1 和 1.2.4路径写错一个字母跑的是旧脚本配新库日志不报错功能矩阵却少一大半。版本脚本运行环境数据库状态备注1.1 及更早Python2旧在线一键版长期默认此版本1.2.4Python32020-11-27 修订本地分析应固定使用此版本我一般会在下载目录里单独放一个FAPROTAX_1.2.4文件夹不和其他版本混放后面所有命令都写相对这个目录的路径避免“以为跑的是新的”这种问题。版本这件事上没什么可折中的要么全套 1.2.4要么别用本地跑去用在线版。2.2 用 conda 起一个只装三件套的环境不推荐往 base 环境里装。原因很实际biom-format 对 numpy 版本有一定要求而做扩增子的机器上往往已经装了 QIIME2 或者别的生信套件版本冲突排查起来非常费时间。另起一个干净环境装错了直接删掉重来这就是后悔药。# 新建一个 Python3 环境名字随意这里叫 faprotax conda create -n faprotax python3.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 -Vwhich python这一步别省。conda 激活失败、shell 里PYTHONPATH被别的东西污染、或者前一个环境没退干净都会让python指向另一个解释器表现出来就是“明明装了 numpy 还是 ModuleNotFoundError”。老教程里写的依赖是conda install numpy和conda install biombiom这个包名在 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 丰度表转成带物种注释的 biomFAPROTAX 的输入不是普通的 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 NR1{print 样本列数:, NF-1} otutab_rare.txt # 检查除表头外每行列数是否一致输出不是 1 就说明有脏行 awk -F\t NR1{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-typeOTU 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-typeOTU 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(axisobservation) # 统计缺注释的条目 ids t.ids(axisobservation) 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(axisobservation)取的就是每条 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参数含义写错会怎样-i带物种注释的 biom 输入表指向无注释的 biom全部 leftovers-gFAPROTAX.txt 数据库路径路径错直接 FileNotFoundError--collapse_by_metadata指定按哪个 metadata 列做归并列名写错则无可用注释-v输出详细过程信息不加则看不到匹配统计--force强制覆盖已有输出重跑时卡在文件已存在-o功能 × 样本丰度矩阵不指定则不会生成主结果-r每个功能组包含哪些 OTU 的明细不指定则无法回溯--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.txtfaprotax.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_reportgrep 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.tsvF $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 NR1{print NF} | sort -u的返回值压到只剩 1 个数字再转。另外注意文件编码从 Excel 另存的 txt 有时带 BOM第一列列名会变成\ufeff#OTU IDgrep 表头就匹配不上用dos2unix或sed -i 1s/^\xEF\xBB\xBF//清一下。6. 装完怎么验证跑通三条自检与一次结果核对装完不等于跑通。脚本一旦返回 0 退出界面安静得像成功了但结果可能是空表。我习惯在交付前固定走三条检查花不了五分钟能挡掉绝大多数“看着跑完了但不能用”的情况。第一条查总量守恒。功能矩阵的丰度总和应当小于等于输入 OTU 表的总和因为只有被注释到的部分才会计入但两者不应该差出一个数量级。# 功能表总丰度跳过表头行 awk -F\t NR1{for(i2;iNF;i) s$i} END{printf function total: %.2f\n, s} faprotax.txt # OTU 表总丰度跳过 # 开头的表头 awk -F\t !/^#/{for(i2;iNF;i) s$i} END{printf otu total: %.2f\n, s} otutab_rare.txt两个数字如果完全相等反而是可疑的说明有可能把未注释的记录也算了进去功能总量占 OTU 总量的三到七成是比较常见的区间。具体表头格式以实际输出为准必要时把NR1换成!/^#/。第二条查 leftover 比例直接看运行时那行日志或者回读 report 的规模。这个比例反映的是注释库与样本的契合度不是脚本质量问题但它是判断“结果能不能讲故事”的直接依据。第三条做一次反向抽查挑一个你关心的功能看看撑起这个功能的物种在生物学上说不说得通。# 提取某个功能组下的所有 ASV 明细块与块之间以 -- 分隔 awk /^nitrogen_fixation/{f1; next} /^--$/{f0} f faprotax_report.txt | head -20抽出来的物种如果全是固氮相关的已知类群这条流程就是可信的如果冒出来一堆完全不相干的属那就得回头查注释文件是不是和 OTU 编号对错了行。这一步做过一次之后后面再跑同类项目会踏实很多。验证通过后再往下游走。faprotax.txt可以直接导进 STAMP 做组间差异也可以按分组信息做功能丰度热图或者拿去和环境因子做相关性分析。第一次跑通时我建议先用一小批样本、用已知注释库的公开数据试一遍把这条链路走顺再上自己的项目数据。血泪经验是有一回赶进度report 文件没生成就交了结果后面想补一个“这些功能由哪些菌贡献”的说明只能把整条流程从头重跑一遍中间隔了两周前面用的哪个版本注释库都记不清了。从那以后我每次跑 FAPROTAX都强制走一遍“总量核对、leftover 记录、随机反查一个功能”这三步并且把命令、版本号、注释库来源原样写进项目目录的 README谁接手都不用猜。希望帮到你。本文还有配套的精品资源点击获取