宏基因组对比：Kraken2/MetaPhlAn/mOTUs¶

# === 安装 ===
conda install -c bioconda kraken2 bracken   # Kraken2 + Bracken
conda install -c bioconda metaphlan         # MetaPhlAn4
conda install -c bioconda motus             # mOTUs3

# === 数据库下载 ===
# Kraken2 标准数据库（~70GB，含细菌/古菌/病毒/人）
kraken2-build --standard --db kraken2_db --threads 8
# 或下载预建数据库（推荐，更快）:
# https://benlangmead.github.io/aws-indexes/k2

# Bracken数据库（基于Kraken2数据库构建）
bracken-build -d kraken2_db -t 8 -k 35 -l 150  # 读长150bp

# MetaPhlAn4 数据库（自动下载，~1.5GB）
metaphlan --install --bowtie2db metaphlan_db    # 首次运行自动下载

# mOTUs3 数据库（自动下载）
motus downloadDB                                # 下载数据库

# === Kraken2 分类 ===
kraken2 \
  --db kraken2_db \                        # 数据库路径
  --threads 8 \                            # 线程数
  --paired \                               # 双端测序
  --output sample_kraken.out \             # 分类结果（每条读段）
  --report sample_kraken.report \          # ★汇总报告
  --confidence 0.05 \                      # 置信度阈值（减少假阳性）
  sample_R1.fq.gz sample_R2.fq.gz        # 输入文件

# Kraken2报告格式：
# %reads | #reads | #reads_direct | rank | taxid | name
# 0.95    5000     3000           S      562    Escherichia coli

# === Bracken 丰度重估计（★推荐配合使用） ===
bracken \
  -d kraken2_db \                          # 数据库
  -i sample_kraken.report \               # Kraken2报告
  -o sample_bracken.txt \                 # Bracken输出
  -r 150 \                                # 读长
  -l S \                                  # 分类层级: S=种 G=属 F=科
  -t 10                                   # 最低读段数阈值

# 合并多样本结果
combine_bracken_outputs.py \
  --files sample1_bracken.txt sample2_bracken.txt \
  --output merged_bracken.txt

# === MetaPhlAn4 分析（★研究首选） ===
metaphlan \
  sample_R1.fq.gz,sample_R2.fq.gz \     # 输入（逗号分隔双端）
  --input_type fastq \                    # 输入类型
  --bowtie2db metaphlan_db \             # 数据库路径
  --nproc 8 \                            # 线程数
  --output_file sample_profile.txt \     # ★输出物种丰度表
  --bowtie2out sample.bowtie2.bz2 \     # 比对文件（可保存重用）
  --unclassified_estimation               # 估计未分类reads比例

# 输出格式:
# #SampleID  Metaphlan_Analysis
# k__Bacteria|p__Firmicutes|...|s__Faecalibacterium_prausnitzii  12.5

# === 合并多样本 ===
merge_metaphlan_tables.py \
  sample1_profile.txt sample2_profile.txt sample3_profile.txt \
  > merged_abundance.txt

# === 提取特定层级 ===
# 只看种水平
grep -E "(^#)|(s__)" merged_abundance.txt | grep -v "t__" > species_table.txt

# 只看属水平
grep -E "(^#)|(g__)" merged_abundance.txt | grep -v "s__" > genus_table.txt

# === MetaPhlAn4 可视化 ===
# GraPhlAn分类树图
export2graphlan.py --skip_rows 1 \
  -i merged_abundance.txt \
  --tree merged.tree.txt \
  --annotation merged.annot.txt

graphlan_annotate.py --annot merged.annot.txt merged.tree.txt merged.xml
graphlan.py merged.xml graphlan.png --dpi 300

# === mOTUs3 分析 ===
motus profile \
  -f sample_R1.fq.gz \                    # 正向读段
  -r sample_R2.fq.gz \                    # 反向读段
  -t 8 \                                  # 线程数
  -o sample_motus.txt \                   # 输出文件
  -n sample1 \                            # 样本名
  -k species \                            # 分类层级
  -p                                      # 输出相对丰度

# 合并多样本
motus merge \
  -i sample1_motus.txt,sample2_motus.txt \
  -o merged_motus.txt

# mOTUs特色：能检测ref_mOTU（已知物种）和meta_mOTU（★未知物种）
# meta_mOTU = 在数据库中没有参考基因组，但通过保守基因确认是新物种

# === R语言对比分析 ===
library(tidyverse)

# 读取各工具结果
metaphlan <- read.table("species_table.txt", header=TRUE, sep="\t")
bracken <- read.table("merged_bracken.txt", header=TRUE, sep="\t")
motus <- read.table("merged_motus.txt", header=TRUE, sep="\t", comment.char="#")

# 比较检测到的物种数
cat("MetaPhlAn4检测物种:", nrow(metaphlan), "\n")
cat("Kraken2+Bracken检测物种:", nrow(bracken), "\n")
cat("mOTUs3检测物种:", nrow(motus), "\n")
# 通常: Kraken2 >> MetaPhlAn4 ≥ mOTUs3（Kraken2有更多假阳性）

# 推荐：取MetaPhlAn4和至少一种其他方法的交集

需求是什么？
├── 快速筛查/临床检测 → Kraken2+Bracken（速度优先）
├── 研究发表/精准定量 → MetaPhlAn4（★首选）
├── 探索未知物种 → mOTUs3（meta_mOTU优势）
├── 需要功能注释 → HUMAnN3（基于MetaPhlAn4）
└── 最佳实践 → MetaPhlAn4 + Kraken2取交集

# === 物种分类速查 ===
# Kraken2+Bracken（快速）
kraken2 --db db --paired --report report R1.fq.gz R2.fq.gz
bracken -d db -i report -o bracken.txt -r 150 -l S

# MetaPhlAn4（★精准，研究首选）
metaphlan R1.fq.gz,R2.fq.gz --input_type fastq --nproc 8 -o profile.txt

# mOTUs3（可检测未知物种）
motus profile -f R1.fq.gz -r R2.fq.gz -o motus.txt

# 合并: merge_metaphlan_tables.py / combine_bracken_outputs.py
# 选择: 临床→Kraken2 | 研究→MetaPhlAn4 | 未知物种→mOTUs3