宏基因组全流程 — 实战命令与代码¶

3. 实战命令/代码（完整流程）¶

以下命令基于 conda 环境，假设已安装所需工具。数据为 Illumina PE150 测序，人肠道宏基因组样本。

3.0 环境准备¶

# ============================================================
# 步骤0：创建并激活conda环境
# ============================================================

# 创建宏基因组分析专用环境（建议Python 3.9+）
conda create -n metagenome python=3.9 -y

# 激活环境
conda activate metagenome

# 安装所需工具（通过bioconda渠道）
conda install -c bioconda -c conda-forge \
    fastqc multiqc fastp bowtie2 samtools \
    kraken2 bracken -y

# 验证安装
fastqc --version    # 应输出类似 FastQC v0.12.1
fastp --version     # 应输出类似 fastp 1.3.3
bowtie2 --version   # 应输出版本信息
kraken2 --version   # 应输出版本信息
samtools --version   # 应输出版本信息

3.1 项目目录结构¶

# ============================================================
# 步骤0.5：创建标准化的项目目录结构
# ============================================================

# 设置项目根目录（根据实际修改）
PROJECT_DIR="/path/to/project"

# 创建标准目录结构
mkdir -p ${PROJECT_DIR}/{00_rawdata,01_fastqc,02_clean,03_hostfree,04_taxonomy,05_function,06_diversity}

# 目录说明：
# 00_rawdata/   - 原始下机数据（fastq.gz + md5.txt）
# 01_fastqc/    - FastQC 和 MultiQC 质量评估报告
# 02_clean/     - fastp 质控后的 clean data
# 03_hostfree/  - 去宿主后的微生物 reads
# 04_taxonomy/  - Kraken2/Bracken 物种注释结果
# 05_function/  - 功能注释结果
# 06_diversity/ - 多样性分析结果

3.2 md5校验¶

# ============================================================
# 步骤1：md5校验 - 验证数据完整性
# ============================================================

# 进入原始数据目录
cd ${PROJECT_DIR}/00_rawdata

# 查看md5文件内容（通常由测序公司提供）
# 格式示例：d41d8cd98f00b204e9800998ecf8427e  sample1_R1.fq.gz
cat md5.txt

# 方法1：使用 -c 参数自动比对（推荐）
# -c 表示 check 模式，自动读取md5.txt中的校验值并逐一比对
md5sum -c md5.txt

# 输出示例：
# sample1_R1.fq.gz: OK
# sample1_R2.fq.gz: OK
# sample2_R1.fq.gz: OK
# sample2_R2.fq.gz: OK
# 如果出现 FAILED，说明该文件损坏，需要重新下载

# 方法2：手动计算并比对（不推荐，了解原理即可）
md5sum sample1_R1.fq.gz
# 输出的md5值与md5.txt中对应行比较

# 统计校验结果
md5sum -c md5.txt 2>&1 | grep -c "OK"      # 通过的文件数
md5sum -c md5.txt 2>&1 | grep -c "FAILED"   # 失败的文件数

3.3 质量评估（FastQC + MultiQC）¶

# ============================================================
# 步骤2：质量评估 - 评估原始数据质量
# ============================================================

# --- 2a. 对原始数据运行 FastQC ---

# -t 8：使用8个线程并行处理（根据服务器配置调整）
# -o：指定输出目录
# 00_rawdata/*.fq.gz：对所有fastq.gz文件运行
fastqc -t 8 \
    -o ${PROJECT_DIR}/01_fastqc/ \
    ${PROJECT_DIR}/00_rawdata/*.fq.gz

# FastQC 会为每个文件生成两个结果文件：
# sample1_R1_fastqc.html  - 可视化报告（浏览器打开）
# sample1_R1_fastqc.zip   - 原始数据（供MultiQC使用）

# --- 2b. 用 MultiQC 汇总所有 FastQC 报告 ---

# 将所有样本的FastQC结果整合为一个交互式网页
multiqc ${PROJECT_DIR}/01_fastqc/ \
    -o ${PROJECT_DIR}/01_fastqc/multiqc_report \
    -n raw_data_multiqc       # 报告名称

# 打开 raw_data_multiqc.html 查看汇总结果
# 重点关注：
# 1. 各样本的 Q20/Q30 比例
# 2. 接头含量（Adapter Content）
# 3. GC 分布是否正常
# 4. 是否有过多的重复序列

3.4 质量控制（fastp）¶

# ============================================================
# 步骤3：质量控制 - fastp 清洗数据
# ============================================================

# --- 单个样本处理示例 ---

SAMPLE="sample1"

fastp \
    -i ${PROJECT_DIR}/00_rawdata/${SAMPLE}_R1.fq.gz \
    -I ${PROJECT_DIR}/00_rawdata/${SAMPLE}_R2.fq.gz \
    -o ${PROJECT_DIR}/02_clean/${SAMPLE}_clean_R1.fq.gz \
    -O ${PROJECT_DIR}/02_clean/${SAMPLE}_clean_R2.fq.gz \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --unqualified_percent_limit 40 \
    --length_required 50 \
    --n_base_limit 5 \
    --cut_front \
    --cut_tail \
    --cut_window_size 4 \
    --cut_mean_quality 20 \
    --correction \
    --thread 8 \
    --json ${PROJECT_DIR}/02_clean/${SAMPLE}_fastp.json \
    --html ${PROJECT_DIR}/02_clean/${SAMPLE}_fastp.html

# 参数详解：
# -i / -I                     : 输入的R1和R2文件（大写I对应R2）
# -o / -O                     : 输出的R1和R2文件（大写O对应R2）
# --detect_adapter_for_pe     : 自动检测PE数据的接头序列并去除
# --qualified_quality_phred 20: 质量值低于20（即错误率>1%）的碱基被视为"不合格"
# --unqualified_percent_limit 40: 一条read中不合格碱基超过40%就丢弃该read
# --length_required 50        : 处理后长度<50bp的read被丢弃
#   # 纯质控场景50可接受；若下游使用 MetaPhlAn/Bracken，建议设为 70 以上
# --n_base_limit 5            : 含N碱基超过5个的read被丢弃
# --cut_front                 : 从5'端开始滑窗，遇到低质量区域则切除
# --cut_tail                  : 从3'端开始滑窗，遇到低质量区域则切除
# --cut_window_size 4         : 滑窗大小为4个碱基
# --cut_mean_quality 20       : 滑窗内平均质量低于20则切除该窗口
# --correction                : 对PE reads重叠区域进行碱基校正
# --thread 8                  : 使用8个线程
# --json / --html             : 输出fastp自带的质量报告

# --- 批量处理所有样本（推荐） ---

# 假设样本名列表文件 samples.txt，每行一个样本名
# sample1
# sample2
# sample3

for SAMPLE in $(cat ${PROJECT_DIR}/samples.txt); do
    echo "Processing: ${SAMPLE}"
    fastp \
        -i ${PROJECT_DIR}/00_rawdata/${SAMPLE}_R1.fq.gz \
        -I ${PROJECT_DIR}/00_rawdata/${SAMPLE}_R2.fq.gz \
        -o ${PROJECT_DIR}/02_clean/${SAMPLE}_clean_R1.fq.gz \
        -O ${PROJECT_DIR}/02_clean/${SAMPLE}_clean_R2.fq.gz \
        --detect_adapter_for_pe \
        --qualified_quality_phred 20 \
        --unqualified_percent_limit 40 \
        --length_required 50 \
        --n_base_limit 5 \
        --cut_front \
        --cut_tail \
        --cut_window_size 4 \
        --cut_mean_quality 20 \
        --correction \
        --thread 8 \
        --json ${PROJECT_DIR}/02_clean/${SAMPLE}_fastp.json \
        --html ${PROJECT_DIR}/02_clean/${SAMPLE}_fastp.html
    echo "${SAMPLE} done."
done

# --- 质控后再次运行 FastQC 比较效果 ---
mkdir -p ${PROJECT_DIR}/01_fastqc/clean
fastqc -t 8 \
    -o ${PROJECT_DIR}/01_fastqc/clean/ \
    ${PROJECT_DIR}/02_clean/*_clean_R*.fq.gz

multiqc ${PROJECT_DIR}/01_fastqc/clean/ \
    -o ${PROJECT_DIR}/01_fastqc/clean/multiqc_report \
    -n clean_data_multiqc

3.5 去宿主（Bowtie2）¶

# ============================================================
# 步骤4：去宿主 - 去除人源序列
# ============================================================

# --- 4a. 下载人类参考基因组（只需做一次） ---

# 下载 GRCh38（hg38）参考基因组
# 方法1：从NCBI下载
# wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

# 存放位置（根据实际修改）
HOST_DB="/path/to/host_db"
mkdir -p ${HOST_DB}

# 解压
# gunzip GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz

# --- 4b. 建立 Bowtie2 索引（只需做一次，耗时约1-2小时） ---

# bowtie2-build 参数：
# 第一个参数：参考基因组fasta文件
# 第二个参数：索引前缀名
# --threads：线程数
bowtie2-build \
    --threads 16 \
    ${HOST_DB}/GRCh38.fna \
    ${HOST_DB}/GRCh38

# 建索引后会生成6个文件：
# GRCh38.1.bt2, GRCh38.2.bt2, GRCh38.3.bt2,
# GRCh38.4.bt2, GRCh38.rev.1.bt2, GRCh38.rev.2.bt2

# --- 4c. 比对去宿主 ---

SAMPLE="sample1"

# 第一步：将 clean reads 比对到人类基因组
# --very-sensitive：最灵敏模式，尽可能找到所有人源reads
# -p 16：使用16个线程
# -x：Bowtie2索引前缀
# -1 / -2：PE数据的R1和R2
# -S：输出SAM文件（后续转为BAM）
bowtie2 \
    --very-sensitive \
    -p 16 \
    -x ${HOST_DB}/GRCh38 \
    -1 ${PROJECT_DIR}/02_clean/${SAMPLE}_clean_R1.fq.gz \
    -2 ${PROJECT_DIR}/02_clean/${SAMPLE}_clean_R2.fq.gz \
    -S ${PROJECT_DIR}/03_hostfree/${SAMPLE}_mapped.sam \
    2> ${PROJECT_DIR}/03_hostfree/${SAMPLE}_bowtie2.log

# Bowtie2运行日志示例：
# 10000000 reads; of these:
#   10000000 (100.00%) were paired; of these:
#     2000000 (20.00%) aligned concordantly 0 times  ← 未比对上（微生物reads）
#     7500000 (75.00%) aligned concordantly exactly 1 time  ← 比对上宿主
#     500000 (5.00%) aligned concordantly >1 times
# 80.00% overall alignment rate  ← 宿主比例

# 第二步：SAM 转 BAM（压缩存储）
# -b：输出BAM格式
# -S：输入是SAM格式（较新版samtools可省略）
# -@：线程数
samtools view -bS \
    -@ 8 \
    ${PROJECT_DIR}/03_hostfree/${SAMPLE}_mapped.sam \
    -o ${PROJECT_DIR}/03_hostfree/${SAMPLE}_mapped.bam

# 第三步：提取未比对上宿主的reads（即微生物reads）
# -b：输出BAM格式
# -f 12：只保留 flag 包含 12 的 reads（read本身和mate都未比对上）
# -F 256：排除次要比对（secondary alignment）
samtools view -b \
    -f 12 \
    -F 256 \
    -@ 8 \
    ${PROJECT_DIR}/03_hostfree/${SAMPLE}_mapped.bam \
    -o ${PROJECT_DIR}/03_hostfree/${SAMPLE}_unmapped.bam

# 第四步：BAM排序（按read名排序，便于后续提取PE reads）
# -n：按read名排序（不是按坐标）
samtools sort -n \
    -@ 8 \
    ${PROJECT_DIR}/03_hostfree/${SAMPLE}_unmapped.bam \
    -o ${PROJECT_DIR}/03_hostfree/${SAMPLE}_unmapped_sorted.bam

# 第五步：将BAM转回FASTQ（分别输出R1和R2）
# -1：输出R1 reads
# -2：输出R2 reads
# -0：输出单端reads（通常很少）
# -s：输出singleton reads
# -n：保持read名中的 /1 /2 标识
samtools fastq \
    -@ 8 \
    -1 ${PROJECT_DIR}/03_hostfree/${SAMPLE}_hostfree_R1.fq.gz \
    -2 ${PROJECT_DIR}/03_hostfree/${SAMPLE}_hostfree_R2.fq.gz \
    -0 /dev/null \
    -s /dev/null \
    -n \
    ${PROJECT_DIR}/03_hostfree/${SAMPLE}_unmapped_sorted.bam

# 第六步：清理中间文件（节省磁盘空间）
rm -f ${PROJECT_DIR}/03_hostfree/${SAMPLE}_mapped.sam
rm -f ${PROJECT_DIR}/03_hostfree/${SAMPLE}_mapped.bam
rm -f ${PROJECT_DIR}/03_hostfree/${SAMPLE}_unmapped.bam
rm -f ${PROJECT_DIR}/03_hostfree/${SAMPLE}_unmapped_sorted.bam

# --- 4d. 统计去宿主效果 ---

# 统计 clean reads 数
echo "Clean reads (R1):"
zcat ${PROJECT_DIR}/02_clean/${SAMPLE}_clean_R1.fq.gz | wc -l | awk '{print $1/4}'

# 统计去宿主后 reads 数
echo "Host-free reads (R1):"
zcat ${PROJECT_DIR}/03_hostfree/${SAMPLE}_hostfree_R1.fq.gz | wc -l | awk '{print $1/4}'

# 计算保留率（人肠道样本通常保留70-99%）

3.6 物种注释（Kraken2 + Bracken）¶

# ============================================================
# 步骤5：物种注释 - Kraken2 + Bracken
# ============================================================

# --- 5a. 准备 Kraken2 数据库（只需做一次） ---

# 数据库存放路径（根据实际修改）
KRAKEN2_DB="/path/to/kraken2_db/standard"

# 方法1：下载预建数据库（推荐，节省时间）
# 常用下载源：https://benlangmead.github.io/aws-indexes/k2
# wget https://genome-idx.s3.amazonaws.com/kraken/k2_standard_20230605.tar.gz
# tar -xzf k2_standard_20230605.tar.gz -C ${KRAKEN2_DB}

# 方法2：自己构建（耗时较长，需要大内存）
# kraken2-build --standard --db ${KRAKEN2_DB} --threads 16

# --- 5b. 运行 Kraken2 ---

SAMPLE="sample1"

kraken2 \
    --db ${KRAKEN2_DB} \
    --threads 16 \
    --confidence 0.05 \
    --paired \
    --gzip-compressed \
    --output ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_kraken2.output \
    --report ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_kraken2.report \
    ${PROJECT_DIR}/03_hostfree/${SAMPLE}_hostfree_R1.fq.gz \
    ${PROJECT_DIR}/03_hostfree/${SAMPLE}_hostfree_R2.fq.gz

# 参数详解：
# --db             : Kraken2数据库路径
# --threads 16     : 使用16个线程
# --confidence 0.05: 置信度阈值，支持该分类的k-mer比例需>=5%
# --paired         : 输入是双端测序数据
# --gzip-compressed: 输入文件是gzip压缩的
# --output         : 逐条read的分类结果（很大，通常不用保留）
# --report         : 汇总分类报告（重要！后续分析用这个）

# --- Kraken2 report 格式说明 ---
# 列1: 该分类单元下reads占总reads的百分比
# 列2: 该分类单元下reads数（含子节点）
# 列3: 直接分到该分类单元的reads数
# 列4: 分类等级代码（D/P/C/O/F/G/S等）
# 列5: NCBI Taxonomy ID
# 列6: 物种学名

# --- 5c. 运行 Bracken（丰度重估计） ---

# Bracken 在 Kraken2 结果基础上，将高层级reads重分配到种水平
bracken \
    -d ${KRAKEN2_DB} \
    -i ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_kraken2.report \
    -o ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_bracken_S.tsv \
    -w ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_bracken_S.report \
    -r 150 \
    -l S \
    -t 10

# 参数详解：
# -d : Kraken2数据库路径（Bracken使用同一个数据库）
# -i : Kraken2的report文件
# -o : Bracken输出的丰度表（TSV格式）
# -w : Bracken生成的新report文件
# -r 150 : read长度（必须与实际一致，PE150就写150）
# -l S   : 分析层级，S表示种（Species）
# -t 10  : 最小reads阈值，低于10条reads的物种被过滤

# 也可以生成属水平丰度
bracken \
    -d ${KRAKEN2_DB} \
    -i ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_kraken2.report \
    -o ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_bracken_G.tsv \
    -w ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_bracken_G.report \
    -r 150 \
    -l G \
    -t 10

# --- 5d. 查看结果 ---

# 查看Bracken种水平丰度表
# 列：name, taxonomy_id, taxonomy_lvl, kraken_assigned_reads, added_reads, new_est_reads, fraction_total_reads
head -20 ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_bracken_S.tsv

# 查看分类统计
echo "=== Classification Summary ==="
grep -c "C" ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_kraken2.output   # 已分类reads数
grep -c "U" ${PROJECT_DIR}/04_taxonomy/${SAMPLE}_kraken2.output   # 未分类reads数

3.7 功能注释（HUMAnN3 简要）¶

# ============================================================
# 步骤6：功能注释 - HUMAnN3（简要示例）
# ============================================================

# --- 安装 HUMAnN3（在单独的conda环境中，避免依赖冲突） ---
# conda create -n humann3 python=3.9 -y
# conda activate humann3
# pip install humann

# --- 下载数据库（首次使用时下载） ---
# humann_databases --download chocophlan full /path/to/humann_db --update-config yes
# humann_databases --download uniref uniref90_diamond /path/to/humann_db --update-config yes

# --- 运行 HUMAnN3 ---

SAMPLE="sample1"

# HUMAnN3 接受合并的PE reads文件
# 先合并R1和R2（cat拼接，HUMAnN3内部处理）
cat ${PROJECT_DIR}/03_hostfree/${SAMPLE}_hostfree_R1.fq.gz \
    ${PROJECT_DIR}/03_hostfree/${SAMPLE}_hostfree_R2.fq.gz \
    > ${PROJECT_DIR}/05_function/${SAMPLE}_merged.fq.gz

# 运行HUMAnN3
humann \
    --input ${PROJECT_DIR}/05_function/${SAMPLE}_merged.fq.gz \
    --output ${PROJECT_DIR}/05_function/${SAMPLE}_humann3/ \
    --threads 16

# 主要输出文件：
# ${SAMPLE}_genefamilies.tsv    - 基因家族丰度（RPK）
# ${SAMPLE}_pathabundance.tsv   - 代谢通路丰度
# ${SAMPLE}_pathcoverage.tsv    - 通路覆盖度

# --- 结果标准化 ---

# 将丰度标准化为 CPM（copies per million）
humann_renorm_table \
    --input ${PROJECT_DIR}/05_function/${SAMPLE}_humann3/${SAMPLE}_merged_genefamilies.tsv \
    --output ${PROJECT_DIR}/05_function/${SAMPLE}_genefamilies_cpm.tsv \
    --units cpm

humann_renorm_table \
    --input ${PROJECT_DIR}/05_function/${SAMPLE}_humann3/${SAMPLE}_merged_pathabundance.tsv \
    --output ${PROJECT_DIR}/05_function/${SAMPLE}_pathabundance_cpm.tsv \
    --units cpm

# --- 合并多个样本的结果 ---
# humann_join_tables -i ${PROJECT_DIR}/05_function/ -o merged_genefamilies.tsv --file_name genefamilies_cpm

3.8 多样性分析（R脚本示例）¶

# ============================================================
# 步骤7：多样性分析 - 基于R（示例代码）
# ============================================================

# 在R环境中运行以下脚本
# 需要安装：vegan, ggplot2, dplyr, ape 等包

# ============================================================
# R脚本：宏基因组多样性分析
# ============================================================

# --- 加载R包 ---
library(vegan)     # 生态学多样性分析核心包
library(ggplot2)   # 画图
library(dplyr)     # 数据处理
library(ape)       # 进化分析（PCoA用）

# --- 读取丰度表 ---
# 假设已将所有样本的Bracken种水平丰度合并为一个矩阵
# 行=物种，列=样本，值=reads数
abundance <- read.csv("species_abundance_matrix.csv", row.names = 1)

# 转置为样本在行、物种在列的格式（vegan要求）
abundance_t <- t(abundance)

# --- Alpha多样性 ---

# 计算Shannon指数
shannon <- diversity(abundance_t, index = "shannon")

# 计算Simpson指数
simpson <- diversity(abundance_t, index = "simpson")

# 计算观测到的物种数（richness）
richness <- specnumber(abundance_t)

# 计算Pielou均匀度（Shannon / ln(物种数)）
pielou <- shannon / log(richness)

# 汇总Alpha多样性结果
alpha_div <- data.frame(
    Sample = rownames(abundance_t),
    Shannon = round(shannon, 3),
    Simpson = round(simpson, 3),
    Richness = richness,
    Pielou = round(pielou, 3)
)
print(alpha_div)

# 写出结果
write.csv(alpha_div, "alpha_diversity.csv", row.names = FALSE)

# --- Alpha多样性箱线图 ---
# 假设有分组信息
# group_info <- read.csv("group_info.csv")  # 列：Sample, Group
# alpha_div <- merge(alpha_div, group_info, by = "Sample")
#
# ggplot(alpha_div, aes(x = Group, y = Shannon, fill = Group)) +
#     geom_boxplot() +
#     theme_minimal() +
#     labs(title = "Shannon Diversity", y = "Shannon Index")

# --- Beta多样性 ---

# 计算Bray-Curtis距离矩阵
bray_dist <- vegdist(abundance_t, method = "bray")

# PCoA分析（主坐标分析）
pcoa_result <- pcoa(bray_dist)

# 提取前两个主坐标轴
pcoa_df <- data.frame(
    Sample = rownames(pcoa_result$vectors),
    PC1 = pcoa_result$vectors[, 1],
    PC2 = pcoa_result$vectors[, 2]
)

# 计算各轴解释的变异比例
var_explained <- pcoa_result$values$Relative_eig * 100

# PCoA可视化
# pcoa_df <- merge(pcoa_df, group_info, by = "Sample")
#
# ggplot(pcoa_df, aes(x = PC1, y = PC2, color = Group)) +
#     geom_point(size = 3) +
#     stat_ellipse(level = 0.95) +
#     theme_minimal() +
#     labs(
#         title = "PCoA (Bray-Curtis)",
#         x = paste0("PC1 (", round(var_explained[1], 1), "%)"),
#         y = paste0("PC2 (", round(var_explained[2], 1), "%)")
#     )

# --- PERMANOVA 统计检验 ---
# 检验不同分组间的Beta多样性是否有显著差异
# group_info <- read.csv("group_info.csv")
# permanova_result <- adonis2(bray_dist ~ Group, data = group_info, permutations = 999)
# print(permanova_result)
# p值 < 0.05 表示组间差异显著