187_微生物组分类器评估¶

一句话概述¶

微生物组分类器评估是对16S/ITS/宏基因组数据中物种分类工具（如QIIME2分类器、Kraken2、MetaPhlAn4等）进行准确性、灵敏度、特异性和计算效率对比的系统性方法，使用模拟数据集和标准参考社区进行benchmarking。

核心知识点表格¶

知识点	说明
分类器类型	基于比对（BLAST）、基于k-mer（Kraken2）、基于标记基因（MetaPhlAn）
评估指标	精确度(Precision)、召回率(Recall)、F1-score、假阳性率
模拟数据	使用已知组成的人工数据集评估
Mock community	ZymoBIOMICS等标准参考社区
QIIME2分类器	sklearn Naive Bayes、VSEARCH、BLAST+
Kraken2	K-mer精确匹配，速度极快
MetaPhlAn4	物种特异标记基因比对
Bracken	Kraken2的丰度重估计工具

步骤详解¶

第一步：准备评估数据集¶

白话解释：要评价分类器好不好，需要一个"答案已知"的数据集。有两种方式：(1)用电脑模拟已知物种组成的测序数据；(2)使用商业标准菌群（如ZymoBIOMICS Mock Community）。

# 方法1：使用ART模拟16S测序数据
# 首先准备已知物种的16S序列
cat << 'EOF' > mock_composition.tsv
Species Relative_Abundance
Escherichia_coli    0.12
Bacillus_subtilis   0.12
Staphylococcus_aureus   0.12
Lactobacillus_fermentum 0.12
Pseudomonas_aeruginosa  0.12
Salmonella_enterica 0.12
Enterococcus_faecalis   0.12
Listeria_monocytogenes  0.12
Saccharomyces_cerevisiae    0.02
Cryptococcus_neoformans 0.02
EOF

# 使用InSilicoSeq模拟宏基因组数据
pip install insilicoseq
iss generate \
    --genomes mock_genomes/ \
    --abundance mock_abundances.txt \
    --model HiSeq \
    --output simulated \
    --n_reads 1000000 \
    --cpus 8

# 方法2：使用CAMI评估数据集
# https://data.cami-challenge.org/
wget https://data.cami-challenge.org/participate  # CAMI2数据

第二步：运行多个分类器¶

# 1. Kraken2 + Bracken
kraken2 --db kraken2_db \
    --paired simulated_R1.fq simulated_R2.fq \
    --output kraken2_output.txt \
    --report kraken2_report.txt \
    --threads 16

bracken -d kraken2_db \
    -i kraken2_report.txt \
    -o bracken_output.txt \
    -r 150 \
    -l S \
    -t 10

# 2. MetaPhlAn4
metaphlan simulated_R1.fq,simulated_R2.fq \
    --input_type fastq \
    --bowtie2db metaphlan_db/ \
    --output_file metaphlan_output.txt \
    --nproc 16

# 3. Kaiju（蛋白质级别分类）
kaiju -t nodes.dmp -f kaiju_db.fmi \
    -i simulated_R1.fq -j simulated_R2.fq \
    -o kaiju_output.txt \
    -z 16

# 4. QIIME2 classify-sklearn（16S数据）
qiime feature-classifier classify-sklearn \
    --i-classifier silva-138-99-nb-classifier.qza \
    --i-reads rep-seqs.qza \
    --o-classification taxonomy.qza \
    --p-n-jobs 16

# 5. Centrifuge
centrifuge -x centrifuge_db \
    -1 simulated_R1.fq -2 simulated_R2.fq \
    -S centrifuge_output.txt \
    --report-file centrifuge_report.txt \
    -p 16

第三步：计算评估指标¶

import pandas as pd
import numpy as np
from sklearn.metrics import precision_score, recall_score, f1_score
import matplotlib.pyplot as plt

# 读取真实组成
true_comp = pd.read_csv("mock_composition.tsv", sep="\t")
true_species = set(true_comp['Species'].str.replace('_', ' '))
true_abundance = dict(zip(
    true_comp['Species'].str.replace('_', ' '),
    true_comp['Relative_Abundance']
))

def evaluate_classifier(predicted_df, true_species, true_abundance, tool_name):
    """评估分类器性能"""
    pred_species = set(predicted_df['Species'])

    # 真阳性、假阳性、假阴性
    tp = true_species & pred_species
    fp = pred_species - true_species
    fn = true_species - pred_species

    precision = len(tp) / (len(tp) + len(fp)) if (len(tp) + len(fp)) > 0 else 0
    recall = len(tp) / (len(tp) + len(fn)) if (len(tp) + len(fn)) > 0 else 0
    f1 = 2 * precision * recall / (precision + recall) if (precision + recall) > 0 else 0

    # 丰度误差（L1距离）
    abundance_error = 0
    for sp in true_species | pred_species:
        true_ab = true_abundance.get(sp, 0)
        pred_ab = predicted_df[predicted_df['Species'] == sp]['Abundance'].values
        pred_ab = pred_ab[0] if len(pred_ab) > 0 else 0
        abundance_error += abs(true_ab - pred_ab)

    return {
        'Tool': tool_name,
        'True_Positives': len(tp),
        'False_Positives': len(fp),
        'False_Negatives': len(fn),
        'Precision': precision,
        'Recall': recall,
        'F1': f1,
        'L1_Error': abundance_error,
        'Total_Predicted': len(pred_species)
    }

# 读取各工具结果（需根据实际格式调整）
# 示例：解析Bracken输出
bracken = pd.read_csv("bracken_output.txt", sep="\t")
bracken_eval = evaluate_classifier(
    pd.DataFrame({'Species': bracken['name'], 'Abundance': bracken['fraction_total_reads']}),
    true_species, true_abundance, 'Kraken2+Bracken'
)

# 汇总所有工具结果
results = pd.DataFrame([bracken_eval])  # 添加其他工具的结果
print(results.to_string(index=False))

第四步：可视化对比¶

import matplotlib.pyplot as plt
import seaborn as sns

# 模拟的评估结果
eval_results = pd.DataFrame({
    'Tool': ['Kraken2', 'MetaPhlAn4', 'Centrifuge', 'Kaiju', 'QIIME2-NB'],
    'Precision': [0.85, 0.95, 0.80, 0.75, 0.90],
    'Recall': [0.90, 0.80, 0.85, 0.88, 0.82],
    'F1': [0.87, 0.87, 0.82, 0.81, 0.86],
    'Speed_reads_per_sec': [500000, 50000, 300000, 100000, 10000],
    'Memory_GB': [40, 8, 30, 15, 4]
})

fig, axes = plt.subplots(1, 3, figsize=(18, 5))

# 1. F1-score对比
axes[0].barh(eval_results['Tool'], eval_results['F1'], color='steelblue')
axes[0].set_xlabel('F1 Score')
axes[0].set_title('Classification Accuracy')
axes[0].set_xlim(0.5, 1.0)

# 2. Precision vs Recall
axes[1].scatter(eval_results['Recall'], eval_results['Precision'], s=100, c='coral')
for i, tool in enumerate(eval_results['Tool']):
    axes[1].annotate(tool, (eval_results['Recall'][i]+0.01, eval_results['Precision'][i]))
axes[1].set_xlabel('Recall')
axes[1].set_ylabel('Precision')
axes[1].set_title('Precision vs Recall')
axes[1].set_xlim(0.5, 1.0)
axes[1].set_ylim(0.5, 1.0)

# 3. 速度对比
axes[2].barh(eval_results['Tool'], eval_results['Speed_reads_per_sec'], color='green', alpha=0.7)
axes[2].set_xlabel('Reads/second')
axes[2].set_title('Processing Speed')
axes[2].set_xscale('log')

plt.tight_layout()
plt.savefig('classifier_comparison.png', dpi=300)

第五步：不同分类学层级的评估¶

# 在不同分类层级评估
levels = ['Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species']

level_results = []
for level in levels:
    # 在每个层级聚合真实和预测结果
    # 通常随着层级深入，分类准确度下降
    for tool in ['Kraken2', 'MetaPhlAn4', 'Centrifuge']:
        # 模拟的结果
        f1 = max(0.5, 0.95 - levels.index(level) * 0.05 + np.random.normal(0, 0.02))
        level_results.append({'Level': level, 'Tool': tool, 'F1': f1})

level_df = pd.DataFrame(level_results)

plt.figure(figsize=(10, 6))
for tool in ['Kraken2', 'MetaPhlAn4', 'Centrifuge']:
    subset = level_df[level_df['Tool'] == tool]
    plt.plot(subset['Level'], subset['F1'], 'o-', label=tool, linewidth=2)
plt.xlabel('Taxonomic Level')
plt.ylabel('F1 Score')
plt.title('Classification Accuracy by Taxonomic Level')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('accuracy_by_level.png', dpi=300)

实战命令速查¶

# 模拟数据生成
iss generate --genomes refs/ --abundance abs.txt --model HiSeq -o sim --n_reads 1M

# Kraken2+Bracken快速评估
kraken2 --db db --paired R1.fq R2.fq --report report.txt -t 16
bracken -d db -i report.txt -o bracken.txt -l S

# MetaPhlAn4快速评估
metaphlan R1.fq,R2.fq --input_type fastq --output_file out.txt --nproc 16

# OPAL评估框架（CAMI标准化评估工具）
opal.py -g gold_standard.profile -o opal_output/ tool1.profile tool2.profile

面试常问点¶

Q1: Kraken2和MetaPhlAn4的核心区别是什么？ A: Kraken2基于k-mer精确匹配整个数据库，覆盖所有已测序物种，速度极快但需要大内存（标准数据库~60GB）。MetaPhlAn4只使用物种特异性标记基因（约1.1M markers），内存需求小（~15GB），特异性更高但可能漏掉新物种。Kraken2适合全面筛查，MetaPhlAn4适合精确定量。

Q2: 如何选择适合自己数据的分类器？ A: 考虑因素：(1)数据类型：16S用QIIME2分类器，宏基因组用Kraken2/MetaPhlAn4；(2)精度需求：高精度用MetaPhlAn4，高灵敏度用Kraken2；(3)计算资源：资源有限用MetaPhlAn4；(4)数据库覆盖：非模式环境考虑Kraken2+自定义数据库。

Q3: 为什么需要用模拟数据评估分类器？ A: 真实数据的"真实组成"是未知的，无法计算precision和recall。模拟数据提供了已知的ground truth，可以精确量化分类器的各项性能指标。但模拟数据不能完全代替真实数据，应结合mock community实验数据。

Q4: Bracken对Kraken2结果做了什么改进？ A: Kraken2倾向于将reads分类到最近公共祖先（LCA），导致很多reads被分到较高的分类层级。Bracken利用reads在物种水平的分布概率，将这些"高层级"reads重新分配到具体物种，提供更准确的物种级丰度估计。

Q5: 分类器的假阳性和假阴性哪个更严重？ A: 取决于应用。临床诊断中假阳性更危险（可能导致不必要的治疗）。生态研究中假阴性更糟（遗漏重要物种）。一般建议高精度（低假阳性）为优先，因为假阳性比假阴性更容易误导。

易错点¶

评估时用训练数据库中的序列：应使用独立于分类器训练集的评估数据
只在物种水平评估：应在多个分类层级评估，某些工具在属级别很好但种级别差
忽略丰度准确性：只看是否检测到物种，不评估丰度估计的准确性
模拟数据不代表真实情况：模拟数据通常比真实数据简单，性能会被高估
数据库版本不一致：不同工具使用不同版本的数据库，影响公平比较

补充知识¶

主流分类器对比¶

工具	方法	速度	内存	精度	适用数据
Kraken2	K-mer匹配	极快	大(40-60GB)	中高	WGS/WES
MetaPhlAn4	标记基因	中	小(15GB)	高	WGS
Centrifuge	FM-index	快	中(8-16GB)	中	WGS/WES
Kaiju	蛋白质匹配	中	大(40GB)	中	WGS
QIIME2-NB	朴素贝叶斯	慢	小(4GB)	中	16S/ITS
mOTUs3	标记基因	中	小	高	WGS