宏基因组功能注释（HUMAnN3 / eggNOG-mapper）¶

输入: 质控后的 fastq reads
          │
          ▼
┌─────────────────────────────────┐
│ Step 1: 物种鉴定（MetaPhlAn4）  │  → 确定样本中有哪些物种
└─────────────────────────────────┘
          │
          ▼
┌─────────────────────────────────┐
│ Step 2: 核酸比对（Bowtie2）     │  → reads 比对到已知物种的基因组
│         针对性比对，只比对已检测到  │    （ChocoPhlAn 定制库）
│         的物种的泛基因组           │
└─────────────────────────────────┘
          │ 未比对上的 reads
          ▼
┌─────────────────────────────────┐
│ Step 3: 翻译搜索（Diamond）     │  → 蛋白级别比对 UniRef 数据库
│         捕获未知/远缘物种的功能   │    （找到新功能）
└─────────────────────────────────┘
          │
          ▼
输出: 基因家族丰度 + 通路丰度 + 通路覆盖度

# ========== 方法1: conda 安装（推荐） ==========
conda create -n humann_env python=3.9    # 创建独立环境
conda activate humann_env                # 激活环境

# 安装 HUMAnN（会自动安装 MetaPhlAn4 作为依赖）
pip install humann                       # 安装 HUMAnN3 主程序

# 验证安装
humann --version                         # 应显示 humann v3.9
metaphlan --version                      # 应显示 MetaPhlAn 4.x

# ========== 方法2: bioconda ==========
conda install -c bioconda humann         # bioconda 渠道安装

# ========== 下载 ChocoPhlAn 核酸数据库（~16 GB） ==========
# 包含已知物种的泛基因组序列，用于核酸级别比对
humann_databases --download chocophlan full /path/to/databases/  # 下载完整版

# ========== 下载 UniRef 蛋白数据库 ==========
# uniref90: 精度高但数据库大（~21 GB，推荐正式分析用）
humann_databases --download uniref uniref90_diamond /path/to/databases/

# uniref50: 速度快数据库小（~7 GB，适合测试或资源有限时）
# humann_databases --download uniref uniref50_diamond /path/to/databases/

# ========== 下载 MetaPhlAn4 marker 数据库 ==========
metaphlan --install --bowtie2db /path/to/databases/metaphlan_db/

# ========== 更新 HUMAnN 配置指向数据库位置 ==========
humann_config --update database_folders nucleotide /path/to/databases/chocophlan/
humann_config --update database_folders protein /path/to/databases/uniref/

# 查看当前配置
humann_config --print                    # 确认数据库路径正确

# ========== 单样本运行 ==========
humann \
  --input sample_kneaddata.fastq.gz \    # 输入：KneadData 质控后的 reads
  --output results/humann/ \             # 输出目录
  --threads 8 \                          # 使用 8 个线程
  --memory-use maximum                   # 内存模式（maximum 更快）

# ========== 批量运行（循环所有样本） ==========
for sample in data/kneaddata/*_kneaddata.fastq.gz; do
    sample_name=$(basename ${sample} _kneaddata.fastq.gz)  # 提取样本名
    humann \
      --input ${sample} \                # 输入文件
      --output results/humann/ \         # 输出目录
      --threads 8 \                      # 线程数
      --output-basename ${sample_name}   # 指定输出文件名前缀
done

# ========== 关键参数说明 ==========
# --nucleotide-database   : 指定 ChocoPhlAn 路径（如果没配置 config）
# --protein-database      : 指定 UniRef 路径
# --bypass-translated-search : 跳过翻译搜索（快但不全）
# --taxonomic-profile     : 提供已有的物种 profile，跳过 MetaPhlAn 步骤

# 文件: sample_genefamilies.tsv
# Gene Family                                             RPKs
UniRef90_A6L0N6: Conserved protein                        67.0
UniRef90_A6L0N6|g__Bacteroides.s__Bacteroides_fragilis    57.0   # 物种贡献
UniRef90_A6L0N6|unclassified                              1.0    # 未知物种贡献

# 文件: sample_pathabundance.tsv
# Pathway                                                 Abundance
PWY-5484: glycolysis II (from fructose-6P)                54.7
PWY-5484|g__Bacteroides.s__Bacteroides_caccae             16.7   # 物种贡献
PWY-5484|unclassified                                     6.0

# 文件: sample_pathcoverage.tsv
# Pathway                                                 Coverage
PWY-5484: glycolysis II (from fructose-6P)                0.95

# ========== 1. 归一化 ==========
# 将 RPK 转换为 CPM（copies per million）或 relab（相对丰度）
humann_renorm_table \
  --input sample_genefamilies.tsv \      # 输入原始基因家族表
  --output sample_genefamilies_cpm.tsv \ # 输出归一化结果
  --units cpm                            # cpm 或 relab

humann_renorm_table \
  --input sample_pathabundance.tsv \
  --output sample_pathabundance_relab.tsv \
  --units relab                          # 通路用相对丰度

# ========== 2. 合并多样本 ==========
humann_join_tables \
  --input results/humann/ \              # 包含所有样本结果的目录
  --output merged_genefamilies.tsv \     # 合并后的基因家族表
  --file_name genefamilies_cpm           # 匹配文件名关键词

humann_join_tables \
  --input results/humann/ \
  --output merged_pathabundance.tsv \
  --file_name pathabundance_relab

# ========== 3. 基因家族重分组（UniRef → KEGG/GO/EC） ==========
# 将 UniRef 基因家族映射到 KEGG KO
humann_regroup_table \
  --input merged_genefamilies.tsv \      # 输入：UniRef 级别的基因家族
  --output merged_ko.tsv \               # 输出：KO 级别
  --groups uniref90_ko                   # 映射关系：UniRef90 → KO

# 映射到 EC 编号
humann_regroup_table \
  --input merged_genefamilies.tsv \
  --output merged_ec.tsv \
  --groups uniref90_level4ec             # UniRef90 → EC number

# 映射到 GO terms
humann_regroup_table \
  --input merged_genefamilies.tsv \
  --output merged_go.tsv \
  --groups uniref90_go                   # UniRef90 → GO

# ========== 4. 添加功能名称 ==========
humann_rename_table \
  --input merged_ko.tsv \                # 输入：KO 编号表
  --output merged_ko_named.tsv \         # 输出：带名称的表
  --names kegg-orthology                 # 名称来源

# ========== 5. 分离分层结果 ==========
# 拆分为社区总量和物种贡献两张表
humann_split_stratified_table \
  --input merged_pathabundance.tsv \     # 输入合并的通路表
  --output results/stratified/           # 输出目录
# 产生: merged_pathabundance_unstratified.tsv（社区总量）
#       merged_pathabundance_stratified.tsv（物种贡献）

# ========== 安装 ==========
pip install metaphlan                    # 通常随 HUMAnN 自动安装

# ========== 下载数据库 ==========
metaphlan --install --bowtie2db /path/to/metaphlan_db/

# ========== 运行 ==========
metaphlan \
  sample_kneaddata.fastq.gz \            # 输入：质控后的 reads
  --input_type fastq \                   # 输入类型
  --bowtie2db /path/to/metaphlan_db/ \   # 数据库路径
  --output_file sample_profile.tsv \     # 物种组成表
  --nproc 8                              # 线程数

# ========== 合并多样本 ==========
merge_metaphlan_tables.py \
  results/metaphlan/*.tsv \              # 所有样本的 profile
  > merged_metaphlan.tsv                 # 合并为一张表

# ========== 物种多样性分析 ==========
# 计算 alpha 多样性
metaphlan \
  sample.fastq.gz \
  --input_type fastq \
  --output_file profile.tsv \
  -t rel_ab_w_read_stats                 # 输出带 reads 统计的相对丰度

#clade_name                               relative_abundance
k__Bacteria                               99.8
k__Bacteria|p__Bacteroidetes              45.2
k__Bacteria|p__Bacteroidetes|...|s__Bacteroides_vulgatus    12.3
k__Bacteria|p__Firmicutes                 40.1
k__Bacteria|p__Firmicutes|...|s__Faecalibacterium_prausnitzii    8.7

# ========== 最新版本：v2.1.13（2025年6月） ==========

# ========== Step 1: 基因预测（Prodigal meta 模式） ==========
prodigal \
  -i assembled_contigs.fasta \           # 输入：组装后的 contigs
  -a predicted_proteins.faa \            # 输出：蛋白序列
  -p meta \                              # 宏基因组模式（必须！）
  -f gff \                               # GFF 格式坐标
  -o genes.gff

# ========== Step 2: eggNOG-mapper 注释 ==========
emapper.py \
  -i predicted_proteins.faa \            # 输入：预测的蛋白序列
  --itype proteins \                     # 输入类型：蛋白质
  -o metagenome_eggnog \                 # 输出前缀
  --data_dir /path/to/eggnog_db/ \       # eggNOG 数据库路径
  --cpu 16 \                             # CPU 数
  -m diamond \                           # 搜索引擎（diamond 最快）
  --dmnd_iterate yes \                   # 迭代搜索提高灵敏度
  --override                             # 覆盖已有输出

# ========== Step 3: 提取 KEGG 注释 ==========
# eggNOG-mapper 输出文件: metagenome_eggnog.emapper.annotations
# 第12列 = KEGG_ko（如 ko:K00844,ko:K01810）
# 第13列 = KEGG_Pathway（如 ko00010,ko00051）
# 第14列 = KEGG_Module（如 M00001,M00002）

# 用 awk 提取 KO 注释（去掉没注释到的行）
awk -F'\t' '$12 != "-"' metagenome_eggnog.emapper.annotations \
  | cut -f1,12 > gene_to_ko.tsv         # 基因 → KO 映射表

KO（KEGG Orthology）    →  最小功能单元（一个酶/蛋白）
     ↓ 组合
Module（功能模块）       →  完成一个小任务的酶组合（如 M00001 = 糖酵解核心模块）
     ↓ 组合
Pathway（代谢通路）      →  完整代谢途径（如 ko00010 = 糖酵解/糖异生）

# ========== HUMAnN 默认输出 MetaCyc 通路 ==========
# 需要重分组才能得到 KEGG 编号

# 基因家族 → KO
humann_regroup_table \
  --input merged_genefamilies_cpm.tsv \
  --output merged_ko_cpm.tsv \
  --groups uniref90_ko                   # UniRef90 → KO 映射

# 为 KO 添加功能描述
humann_rename_table \
  --input merged_ko_cpm.tsv \
  --output merged_ko_named.tsv \
  --names kegg-orthology

# ========== KO → Module/Pathway 映射（需额外工具） ==========
# 方法1: 使用 MinPath 或 KEGGDecoder
# 方法2: 使用 R 包 clusterProfiler 做 KEGG 富集分析
# 方法3: 自定义 Python 脚本映射 KO → Pathway

# 示例: Python 提取 KEGG module 丰度
# 原理：一个 module 的丰度 = 其包含的 KO 的最小值（保守估计）

# ========== R: KEGG 通路富集分析 ==========
library(clusterProfiler)                 # KEGG 富集分析包

# 读取差异 KO（如 T2D vs 健康的差异 KO 列表）
diff_ko <- read.table("diff_ko_list.txt", header=FALSE)$V1  # KO 编号列表
# 去掉 "ko:" 前缀
diff_ko <- gsub("ko:", "", diff_ko)      # K00844 格式

# KEGG pathway 富集
enrich_result <- enrichKEGG(
  gene = diff_ko,                        # 差异 KO 列表
  organism = "ko",                       # "ko" 表示直接用 KO 编号
  keyType = "kegg",                      # 输入类型
  pvalueCutoff = 0.05,                   # p 值阈值
  qvalueCutoff = 0.1                     # FDR 校正后阈值
)

# KEGG module 富集
module_result <- enrichMKEGG(
  gene = diff_ko,                        # 差异 KO
  organism = "ko",                       # KO 编号
  keyType = "kegg",
  pvalueCutoff = 0.05
)

# 查看结果
head(enrich_result@result)               # 显著富集的通路

# ========== 方法1: 从 KO 表中提取 SCFA 相关 KO ==========
# 丁酸相关 KO 列表
echo -e "K00929\nK01034\nK01035\nK00248\nK01782\nK00634" > butyrate_kos.txt

# 从合并的 KO 表中提取
grep -f butyrate_kos.txt merged_ko_named.tsv > butyrate_ko_abundance.tsv

# ========== 方法2: 从 MetaCyc 通路中找 SCFA 相关通路 ==========
# HUMAnN 输出的 MetaCyc 通路中，SCFA 相关的：
# PWY-5022: 4-aminobutyrate degradation V（丁酸相关）
# CENTFERM-PWY: pyruvate fermentation to butanoate（丁酸发酵）
# P163-PWY: L-lysine fermentation to acetate and butyrate
# PWY-5100: pyruvate fermentation to acetate and lactate II

grep -i "butyr\|butan\|acetate\|propanoate" merged_pathabundance.tsv \
  > scfa_pathways.tsv                    # 提取 SCFA 相关通路

#!/usr/bin/env python3
"""提取 SCFA 通路丰度，比较 T2D vs 健康对照"""
import pandas as pd                      # 数据处理
from scipy import stats                  # 统计检验

# 读取 HUMAnN 通路丰度表（unstratified，社区总量）
pathway_df = pd.read_csv(
    "merged_pathabundance_unstratified.tsv",
    sep="\t",
    index_col=0                          # 第一列是通路名
)

# 定义 SCFA 相关 MetaCyc 通路
scfa_pathways = [
    "CENTFERM-PWY",                      # 丙酮酸发酵产丁酸
    "P163-PWY",                          # 赖氨酸发酵产乙酸和丁酸
    "PWY-5100",                          # 丙酮酸发酵产乙酸
    "PWY-5022",                          # 4-氨基丁酸降解
    "PWY-7111",                          # 丙酸降解
]

# 提取 SCFA 通路子集
scfa_df = pathway_df[
    pathway_df.index.str.contains("|".join(scfa_pathways))
]

# 读取分组信息
metadata = pd.read_csv("metadata.tsv", sep="\t")  # 含 sample_id, group 列
t2d_samples = metadata[metadata["group"] == "T2D"]["sample_id"].tolist()
ctrl_samples = metadata[metadata["group"] == "Control"]["sample_id"].tolist()

# Wilcoxon 检验比较 T2D vs Control
for pathway in scfa_df.index:
    t2d_values = scfa_df.loc[pathway, t2d_samples]    # T2D 组丰度
    ctrl_values = scfa_df.loc[pathway, ctrl_samples]  # 对照组丰度
    stat, pval = stats.mannwhitneyu(t2d_values, ctrl_values)
    print(f"{pathway}: p={pval:.4f}, T2D_mean={t2d_values.mean():.2f}, "
          f"Ctrl_mean={ctrl_values.mean():.2f}")

#!/usr/bin/env python3
"""通路丰度热图"""
import seaborn as sns                    # 热图绑定
import matplotlib.pyplot as plt          # 画图基础库
import pandas as pd

# 读取归一化后的通路丰度
df = pd.read_csv("merged_pathabundance_relab_unstratified.tsv",
                 sep="\t", index_col=0)

# 选取 Top 30 通路
top30 = df.mean(axis=1).nlargest(30).index  # 平均丰度 Top 30
plot_df = df.loc[top30]

# 画热图
plt.figure(figsize=(12, 10))
sns.clustermap(
    plot_df,                             # 数据
    cmap="YlOrRd",                       # 配色：黄-橙-红
    standard_scale=0,                    # 按行标准化（Z-score）
    figsize=(14, 10),
    yticklabels=True,                    # 显示通路名称
    col_cluster=True,                    # 样本聚类
    row_cluster=True                     # 通路聚类
)
plt.savefig("pathway_heatmap.pdf", dpi=300, bbox_inches="tight")

# ========== HUMAnN 自带 barplot 工具 ==========
humann_barplot \
  --input merged_pathabundance.tsv \     # 输入：分层的通路丰度表
  --focal-feature "PWY-5484" \           # 指定关注的通路
  --focal-metadata group \               # 按组着色
  --last-metadata group \                # 分组变量
  --output glycolysis_barplot.png        # 输出柱状图

# ========== LEfSe: 找组间差异功能特征 ==========
# LEfSe = Linear discriminant analysis Effect Size
# 用途：找出在 T2D 组和健康组之间有显著差异的通路/物种

# Step 1: 准备 LEfSe 输入格式
# 第一行 = 分组信息，后面每行 = 一个 feature
# group   T2D   T2D   T2D   Ctrl  Ctrl  Ctrl
# PWY-1   0.5   0.3   0.4   0.1   0.05  0.08
lefse_format_input.py \
  merged_pathabundance_lefse.tsv \       # 输入：转换好格式的丰度表
  formatted.in \                         # 中间格式
  -c 1                                   # 第1行是分组

# Step 2: 运行 LEfSe
lefse_run.py \
  formatted.in \                         # 输入
  results.res \                          # 输出结果
  -l 2.0                                 # LDA score 阈值（>2 为显著）

# Step 3: 可视化
lefse_plot_res.py \
  results.res \                          # 输入结果
  lefse_barplot.png \                    # 输出柱状图
  --dpi 300

lefse_plot_cladogram.py \
  results.res \                          # 输入
  lefse_cladogram.png \                  # 输出进化分支图
  --dpi 300