16S扩增子PICRUSt2功能预测¶

一句话概述¶

PICRUSt2能根据16S rRNA基因测序数据"猜测"微生物群落有什么功能基因——不需要做昂贵的宏基因组测序，用便宜的16S数据就能预测功能，虽然是"猜"但准确度还不错。

核心知识点表格¶

知识点	说明
PICRUSt2	基于系统发育关系预测微生物群落功能组成的工具
原理	利用已知基因组的功能注释，通过16S系统发育推断未知菌的功能
输入	16S ASV/OTU序列 + 丰度表
输出	基因家族丰度、EC号丰度、MetaCyc通路丰度
最新版本	v2.6.0+（官方，内置PICRUSt2-SC数据库），PICRUSt2-SC使用GTDB r214数据库（2025年发布）
局限性	只是功能"预测"，不是真实测量，需谨慎解释

白话解释¶

想象你参加一个聚会，只知道每个人的姓氏（= 16S序列），不知道每个人的工作（= 功能基因）。

PICRUSt2的做法： 1. 先查"族谱"（参考基因组数据库），找出和聚会上每个人最像的亲戚 2. 这些亲戚的工作（功能基因）是已知的 3. 根据"亲戚做什么工作"来推测聚会上每个人可能做什么工作 4. 最后统计整个聚会的"职业分布"

关键限制： 这终究是"推测"——如果聚会上来了一个和族谱里所有人都不像的陌生人，PICRUSt2就猜不准了。

安装与配置¶

# 方法1：用conda安装（推荐）
conda create -n picrust2 -c bioconda -c conda-forge picrust2  # 创建独立环境
conda activate picrust2  # 激活环境

# 验证安装
picrust2_pipeline.py --help  # 查看帮助信息

# 方法2：用pip安装
pip install picrust2  # pip安装

# 检查版本
picrust2_pipeline.py --version  # 查看版本号

# 系统要求：
# - Linux/MacOS（Windows需要WSL）
# - 至少16GB内存
# - Python 3.8+

输入数据准备¶

从QIIME2获取输入¶

# PICRUSt2需要两个输入文件：
# 1. study_seqs.fna — ASV代表序列（FASTA格式）
# 2. study_seqs.biom — ASV丰度表（BIOM格式）

# 如果你是用QIIME2分析的16S数据：

# 导出ASV代表序列
qiime tools export \
    --input-path rep-seqs.qza \       # QIIME2的代表序列文件
    --output-path exported_seqs/      # 导出目录
# 得到 exported_seqs/dna-sequences.fasta

# 导出feature表（BIOM格式）
qiime tools export \
    --input-path table.qza \          # QIIME2的feature表
    --output-path exported_table/     # 导出目录
# 得到 exported_table/feature-table.biom

# 重命名文件
mv exported_seqs/dna-sequences.fasta study_seqs.fna    # 重命名序列文件
mv exported_table/feature-table.biom study_seqs.biom    # 重命名表文件

数据预处理注意事项¶

# 重要：PICRUSt2对输入数据有要求

# 1. 必须是去噪后的ASV或聚类的OTU，不能是原始reads
# 2. 建议去除单例（只出现一次的ASV），减少噪音
# 3. 不要做稀释（rarefaction）后再输入PICRUSt2
# 4. 序列不需要trim到统一长度

# 在QIIME2中过滤低频ASV
qiime feature-table filter-features \
    --i-table table.qza \             # 输入feature表
    --p-min-frequency 10 \            # 最低总频率10
    --p-min-samples 2 \               # 至少出现在2个样本中
    --o-filtered-table filtered_table.qza  # 输出过滤后的表

运行PICRUSt2¶

一键运行（推荐新手）¶

# 最简单的运行方式：一个命令搞定一切
picrust2_pipeline.py \
    -s study_seqs.fna \        # 输入ASV序列
    -i study_seqs.biom \       # 输入丰度表
    -o picrust2_output \       # 输出目录
    -p 16 \                    # 线程数
    --stratified \             # 输出按物种分层的结果
    --per_sequence_contrib     # 输出每个ASV对通路的贡献

# 输出文件说明：
# picrust2_output/
# ├── EC_metagenome_out/        — EC号（酶）预测结果
# │   └── pred_metagenome_unstrat.tsv.gz  — 非分层丰度表
# ├── KO_metagenome_out/        — KEGG KO预测结果
# │   └── pred_metagenome_unstrat.tsv.gz  — 非分层丰度表
# ├── pathways_out/             — MetaCyc通路预测结果
# │   ├── path_abun_unstrat.tsv.gz  — 通路丰度表
# │   └── path_abun_strat.tsv.gz   — 按物种分层的通路丰度
# ├── out.tre                   — 放置ASV到参考树的结果
# └── intermediate/            — 中间文件

分步运行（理解每一步）¶

# 第一步：序列放置（将ASV放到参考系统发育树上）
place_seqs.py \
    -s study_seqs.fna \        # 输入ASV序列
    -o out.tre \               # 输出放置后的树
    -p 16 \                    # 线程数
    --intermediate intermediate/placement  # 中间文件目录

# 第二步：隐藏状态预测（预测每个ASV的基因家族）
hsp.py \
    -i 16S \                   # 输入类型：16S
    -t out.tre \               # 输入放置后的树
    -o EC_predicted.tsv.gz \   # 输出EC号预测
    -n \                       # 使用最近邻插值
    -p 16                      # 线程数

hsp.py \
    -i KO \                    # 输入类型：KO
    -t out.tre \
    -o KO_predicted.tsv.gz \
    -n -p 16

# 第三步：宏基因组预测（结合丰度表计算群落功能）
metagenome_pipeline.py \
    -i study_seqs.biom \       # 输入丰度表
    -m EC_predicted.tsv.gz \   # 输入EC预测结果
    -o EC_metagenome_out \     # 输出目录
    --strat_out                # 输出分层结果

metagenome_pipeline.py \
    -i study_seqs.biom \
    -m KO_predicted.tsv.gz \
    -o KO_metagenome_out \
    --strat_out

# 第四步：通路推断（用MinPath推断通路丰度）
pathway_pipeline.py \
    -i EC_metagenome_out/pred_metagenome_contrib.tsv.gz \  # 输入EC丰度
    -o pathways_out \          # 输出通路目录
    -p 16                      # 线程数

结果分析与可视化¶

# 解压结果文件
gunzip picrust2_output/KO_metagenome_out/pred_metagenome_unstrat.tsv.gz
gunzip picrust2_output/pathways_out/path_abun_unstrat.tsv.gz

# 用R分析PICRUSt2结果
Rscript << 'REOF'
# 加载必要的R包
library(ggplot2)  # 绑定ggplot2包

# 读取KO丰度表
ko_table <- read.delim(
    "picrust2_output/KO_metagenome_out/pred_metagenome_unstrat.tsv",
    row.names = 1  # 第一列作为行名
)

# 读取通路丰度表
pathway_table <- read.delim(
    "picrust2_output/pathways_out/path_abun_unstrat.tsv",
    row.names = 1
)

cat("预测到的KO数量:", nrow(ko_table), "\n")      # 打印KO数量
cat("预测到的通路数量:", nrow(pathway_table), "\n")  # 打印通路数量
cat("样本数量:", ncol(ko_table), "\n")              # 打印样本数量

# 差异分析示例（用Wilcoxon检验比较两组）
# 假设metadata有Group列（T2D/Control）
metadata <- read.delim("metadata.tsv", row.names = 1)  # 读取样本信息

# 对每条通路做差异检验
pvals <- apply(pathway_table, 1, function(x) {
    group1 <- x[metadata$Group == "T2D"]      # T2D组的丰度
    group2 <- x[metadata$Group == "Control"]  # 对照组的丰度
    wilcox.test(group1, group2)$p.value        # Wilcoxon检验
})

# FDR校正
padj <- p.adjust(pvals, method = "BH")  # BH方法校正p值

# 筛选显著差异的通路
sig_pathways <- names(padj[padj < 0.05])  # 校正后p<0.05
cat("显著差异的通路数量:", length(sig_pathways), "\n")  # 打印数量
REOF

将KO号转换为KEGG通路描述¶

# 添加功能描述
add_descriptions.py \
    -i picrust2_output/KO_metagenome_out/pred_metagenome_unstrat.tsv.gz \  # 输入KO表
    -m KO \                   # 描述类型：KO
    -o KO_with_descriptions.tsv  # 输出带描述的表

# 添加EC号描述
add_descriptions.py \
    -i picrust2_output/EC_metagenome_out/pred_metagenome_unstrat.tsv.gz \
    -m EC \
    -o EC_with_descriptions.tsv

# 添加MetaCyc通路描述
add_descriptions.py \
    -i picrust2_output/pathways_out/path_abun_unstrat.tsv.gz \
    -m METACYC \
    -o pathways_with_descriptions.tsv

在QIIME2中使用PICRUSt2¶

# QIIME2插件版本的PICRUSt2
# 安装q2-picrust2插件
conda install -c bioconda q2-picrust2  # 安装插件

# 在QIIME2中运行PICRUSt2
qiime picrust2 full-pipeline \
    --i-table table.qza \              # 输入feature表
    --i-seq rep-seqs.qza \             # 输入代表序列
    --output-dir q2-picrust2-output \  # 输出目录
    --p-threads 16 \                   # 线程数
    --p-hsp-method mp \                # 预测方法：最大简约
    --p-max-nsti 2                     # 最大NSTI值（排除无参考的ASV）

# 输出的QIIME2 artifact可以直接用于下游分析
# 如差异分析、PCoA可视化等

PICRUSt2-SC：2025年最新数据库更新¶

# PICRUSt2-SC (Sugar-Coated) 使用GTDB r214数据库
# 包含27,870个参考基因组（比原始数据库更多更新）

# 安装PICRUSt2-SC（2025年4月发表）
# 需要单独下载更新的参考树和基因数据
# 详见：https://academic.oup.com/bioinformatics/article/41/5/btaf269

# 使用PICRUSt2-SC的自定义数据库运行
picrust2_pipeline.py \
    -s study_seqs.fna \
    -i study_seqs.biom \
    -o picrust2sc_output \
    --custom_trait_tables /path/to/picrust2sc/ko_counts.tsv.gz \  # 自定义KO表
    --marker_gene_table /path/to/picrust2sc/16S_counts.tsv.gz \   # 自定义16S表
    --ref_dir /path/to/picrust2sc/ \  # 自定义参考目录
    -p 16

常见报错与解决¶

报错	原因	解决方案
`Sequences not in FASTA format`	序列文件格式不对	确保是标准FASTA格式，不是FASTQ
`No sequences overlap with reference`	ASV序列和参考树不匹配	检查序列是否是16S，是否去了引物
`NSTI value too high`	ASV与参考基因组差异太大	用 `--max_nsti 2` 过滤掉NSTI>2的ASV
内存不足（OOM）	ASV太多，内存不够	过滤低频ASV，减少输入数量
`biom format error`	BIOM文件版本不兼容	`biom convert` 转换版本
运行非常慢	线程数设太低或ASV太多	增加 `-p` 线程数，过滤稀有ASV

速查表¶

# PICRUSt2完整命令
picrust2_pipeline.py -s seqs.fna -i table.biom -o output -p 16

# 输出文件对应关系
KO_metagenome_out → KEGG直系同源基因丰度
EC_metagenome_out → 酶分类号丰度
pathways_out → MetaCyc代谢通路丰度

# 支持的功能数据库
KO — KEGG直系同源基因
EC — 酶分类号
COG — 直系同源聚类
PFAM — 蛋白质家族
TIGRFAM — TIGR蛋白家族

# NSTI值解读（Nearest Sequenced Taxon Index）
NSTI < 0.03 → 参考基因组很近，预测可靠
NSTI 0.03-0.15 → 一般可靠
NSTI 0.15-2.0 → 预测不太可靠
NSTI > 2.0 → 预测不可靠，建议排除

# 安装验证
picrust2_pipeline.py --version  # 查看版本
place_seqs.py --help           # 查看帮助

面试高频问题¶

Q1：PICRUSt2的原理是什么？

A：PICRUSt2利用已测序基因组的功能注释和16S系统发育关系来预测功能。具体步骤：①将ASV序列放置到参考系统发育树上；②利用祖先状态重建（隐藏状态预测）推断每个ASV对应的基因家族组成；③结合ASV丰度表计算群落整体的功能丰度；④用MinPath推断代谢通路。

Q2：PICRUSt2的结果可靠吗？有什么局限性？

A：PICRUSt2预测的是功能"潜力"而非实际表达。局限性：①依赖参考数据库覆盖度，新物种预测不准；②不能反映基因表达水平（转录调控）；③假设同一物种的不同菌株功能相似（忽略株水平差异）。NSTI值可以评估预测可靠性。

Q3：PICRUSt2和宏基因组测序的功能分析有什么区别？

A：PICRUSt2基于16S数据"预测"功能，便宜快速但是推断性质的；宏基因组直接测序所有DNA来"测量"功能基因，准确但贵。PICRUSt2适合生成假设，宏基因组适合验证假设。

Q4：什么是NSTI？NSTI值高说明什么？

A：NSTI（Nearest Sequenced Taxon Index）表示ASV与最近已测序参考基因组的进化距离。NSTI越高说明参考越远，预测越不可靠。通常NSTI>2的ASV建议排除。

Q5：在你的T2D项目中用过PICRUSt2吗？

A：（结合自身项目回答）16S数据可以用PICRUSt2做初步的功能预测分析，关注糖代谢、短链脂肪酸合成等与T2D相关的通路。但重要发现建议用宏基因组数据验证。