宏基因组全流程（Metagenomics Pipeline）¶

一句话概述¶

宏基因组分析从原始测序数据（raw sequencing data）出发，经过 数据校验 → 质量评估 → 质量控制 → 去宿主 → 物种注释 → 功能注释 → 多样性分析 的完整流程，最终实现对环境样本中微生物群落组成和功能的全面解析。

目录¶

1. 核心知识点（流程总览表）
2. 各步骤详解
2.1 md5校验
2.2 质量评估
2.3 质量控制/清洗
2.4 去宿主
2.5 物种注释
2.6 功能注释
2.7 多样性分析
3. 实战命令/代码（完整流程）
4. 面试常问点
5. 易错/易混淆点
6. 补充知识

1. 核心知识点（流程总览表）¶

序号	步骤	工具	目的	输入	输出	关键参数
1	md5校验	`md5sum`	验证数据传输完整性，防止文件损坏	原始 fastq.gz 文件 + md5.txt	校验通过/失败报告	`-c`（check模式）
2	质量评估	`FastQC` + `MultiQC`	评估原始数据质量，发现潜在问题	fastq.gz 文件	HTML质量报告	`-t`（线程数）、`-o`（输出目录）
3	质量控制	`fastp`	去接头序列、过滤低质量reads、去除短reads	原始 R1/R2 fastq.gz	清洗后 clean R1/R2 fastq.gz	`--qualified_quality_phred 20`、`--length_required 50`、`--detect_adapter_for_pe`
4	去宿主	`Bowtie2` + `samtools`	去除人源序列污染，保留微生物reads	clean fastq.gz + 宿主基因组索引	非宿主 reads fastq.gz	`--very-sensitive`、`-p`（线程）、`-f 12`（保留双端未比对）
5	物种注释	`Kraken2` + `Bracken`	对reads进行物种分类和丰度估计	去宿主后 fastq.gz + Kraken2数据库	物种分类报告 + 丰度表	`--confidence 0.05`、`-r`（read length）、`-l`（分类层级）
6	功能注释	`HUMAnN3` / `eggNOG-mapper`	注释微生物群落的代谢通路和基因功能	去宿主后 fastq.gz 或组装后的基因序列	通路丰度表、基因家族表	`--threads`、数据库路径
7	多样性分析	`R`（vegan包）/ `Python`（scikit-bio）	评估群落组成差异和样本间关系	物种丰度表（OTU/species table）	Alpha/Beta多样性指数、PCoA图等	距离度量方法（Bray-Curtis等）

2. 各步骤详解¶

2.1 md5校验（Data Integrity Check）¶

什么是md5？

md5（Message Digest Algorithm 5）是一种哈希算法，可以为任意大小的文件生成一个固定长度的128位（32个十六进制字符）"指纹"。即使文件只改变1个字节，md5值也会完全不同。

为什么要做md5校验？

测序数据通常由测序公司通过网络传输，过程中可能出现数据损坏
测序数据文件往往很大（几GB到几十GB），传输中断或磁盘错误的概率不低
如果数据损坏却不知道，后续分析结果都不可靠

实际操作：

测序公司交付数据时会同时提供一个 md5.txt 文件，里面记录了每个fastq.gz文件对应的md5值。我们下载数据后，重新计算md5值并与之比较，一致则说明数据完整。

2.2 质量评估（Quality Assessment）¶

FastQC 报告的核心模块：

模块名称	含义	正常标准	异常信号
Per base sequence quality	每个位置的碱基质量分布	Q20以上占大多数，Q30>80%	3'端质量急剧下降
Per sequence quality scores	每条read的平均质量分布	峰值在Q30以上	双峰分布
Per base sequence content	每个位置的ATCG比例	四条线平行、接近25%	开头几个碱基波动（正常，来自random primer）
Per sequence GC content	GC含量分布	接近理论正态分布	多峰（可能有污染）
Sequence length distribution	reads长度分布	集中在固定长度（如150bp）	长度不均一
Adapter content	接头序列残留比例	接近0%	明显上升曲线（需要去接头）
Overrepresented sequences	高频重复序列	无或极少	大量重复（可能是接头或rRNA污染）

MultiQC 的作用：

当样本很多时（比如实习中做几十个样本的宏基因组），逐一看FastQC报告太麻烦。MultiQC 能把所有样本的FastQC报告汇总成一个交互式网页，方便横向比较。

关键质量指标阈值：

指标	好	一般	差
Q20 占比	>95%	90-95%	<90%
Q30 占比	>85%	80-85%	<80%
Adapter含量	<5%	5-15%	>15%
重复率	<30%	30-50%	>50%
GC分布	单峰正态	轻微偏移	多峰/严重偏移

2.3 质量控制/清洗（Quality Control）¶

fastp 是什么？

fastp 是一个高效的fastq文件预处理工具，可以同时完成质量过滤、接头去除、低质量碱基修剪等操作。相比老牌工具 Trimmomatic，fastp 速度更快（C++实现）、使用更简单，且自带HTML质量报告。

fastp 核心参数详解：

参数	含义	常用值	说明
`-i` / `-I`	输入R1/R2文件	fastq.gz	大写I对应R2
`-o` / `-O`	输出R1/R2文件	clean.fq.gz	大写O对应R2
`--detect_adapter_for_pe`	自动检测PE接头	开启	PE数据推荐开启
`--adapter_sequence`	指定R1接头序列	Illumina通用接头	知道接头序列时手动指定更准确
`--adapter_sequence_r2`	指定R2接头序列	Illumina通用接头	同上
`-q` / `--qualified_quality_phred`	碱基质量阈值	20	低于此值的碱基被认为是"不合格"
`-u` / `--unqualified_percent_limit`	不合格碱基最大比例	40	一条read中不合格碱基超过此比例则过滤
`-l` / `--length_required`	最短read长度	50	短于此值的read被丢弃
`--cut_front`	从5'端滑窗切除低质量	开启	类似Trimmomatic的LEADING
`--cut_tail`	从3'端滑窗切除低质量	开启	类似Trimmomatic的TRAILING
`--cut_window_size`	滑窗大小	4	滑窗内平均质量低于阈值则切除
`--cut_mean_quality`	滑窗内平均质量阈值	20	与滑窗配合使用
`-w` / `--thread`	线程数	8-16	根据服务器配置设置
`-j` / `-h`	JSON/HTML报告输出	指定路径	fastp自带质量报告，非常好用
`--n_base_limit`	允许的最大N碱基数	5	含N过多的read被丢弃
`--correction`	碱基校正（PE overlap区域）	开启	利用PE reads重叠区域纠错

常见Illumina接头序列：

TruSeq Adapter (Read 1): AGATCGGAAGAGCACACGTCTGAACTCCAGTCA
TruSeq Adapter (Read 2): AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT
Nextera Adapter:          CTGTCTCTTATACACATCT

reads保留率参考标准：

阶段	正常保留率	需要注意
fastp质控后	>90%	<85% 说明原始数据质量较差
去宿主后（人肠道样本）	70-99%（视样本类型）	<50% 可能宿主污染严重
去宿主后（土壤/水样本）	>95%	理论上宿主含量极低

2.4 去宿主（Host Removal）¶

为什么要去宿主？

宏基因组测序的对象是环境中的微生物DNA，但在样本采集和提取过程中，不可避免地会混入宿主（如人体肠道样本中的人源DNA）的序列。

不去宿主的后果： 1. 浪费计算资源：人源reads可能占10-50%甚至更高，这些reads参与后续分析纯属浪费 2. 干扰物种注释：人源序列可能被错误注释为微生物，产生假阳性 3. 影响丰度估计：分母（总reads数）被人源序列稀释，微生物丰度被低估 4. 功能注释污染：人类基因功能被混入微生物功能分析中

使用什么参考基因组？

人源宿主：GRCh38（hg38），即人类参考基因组第38版
下载地址：NCBI（GCA_000001405.15）或 UCSC（hg38.fa.gz）
小鼠样本：GRCm39（mm39）
其他动物样本：对应物种的参考基因组

Bowtie2 去宿主的原理：

先用 bowtie2-build 对宿主参考基因组建立索引（只需做一次）
将质控后的reads比对到宿主基因组
提取未比对上的reads（即非宿主、属于微生物的reads）
用 samtools 从BAM文件中提取未比对的reads

关键概念：samtools flag 过滤

flag	含义
4	该read未比对上（unmapped）
8	配对的mate未比对上
12	该read和mate都未比对上（4+8=12）

对于PE数据，我们用 samtools view -b -f 12 提取双端都没比对上宿主基因组的read对。

Bowtie2 比对模式：

模式	参数	速度	灵敏度	适用场景
very-fast	`--very-fast`	最快	最低	粗筛
fast	`--fast`	快	较低	快速分析
sensitive	`--sensitive`	中等	中等	默认模式
very-sensitive	`--very-sensitive`	慢	最高	去宿主推荐

去宿主时推荐 --very-sensitive，因为我们希望尽可能多地识别和去除宿主序列，宁可多去一些也不要漏掉。

2.5 物种注释（Taxonomic Profiling）¶

Kraken2¶

原理：

Kraken2 使用精确的 k-mer 匹配算法进行物种分类。它将每条read切割成连续的 k-mer（默认k=35），然后在预建的数据库中查找每个k-mer属于哪个物种。最终通过LCA（Lowest Common Ancestor，最低公共祖先）算法确定read的分类。

数据库类型：

数据库	大小	覆盖范围	适用场景
Standard	~70GB	RefSeq完整库（细菌、古菌、病毒、人类）	常规分析、精度要求高
MiniKraken	~8GB	标准库的缩减版	内存有限时使用
PlusPF	~70GB+	Standard + 真菌 + 原生动物	需要分析真核微生物时
自建库	可变	自定义	特殊需求

confidence 参数详解：

--confidence 参数控制分类的严格程度，取值范围 0-1。

含义：一条read被分类到某个物种时，支持该分类的k-mer比例必须 >= confidence值
默认值：0（不过滤，所有能分类的read都保留分类结果）

confidence值	效果	适用场景
0	最宽松，分类最多reads，但假阳性最高	不推荐用于最终分析
0.05	常用推荐值，轻微过滤，平衡灵敏度和精确度	常规宏基因组分析
0.1	中等严格	需要较高可信度时
0.2	较严格，只保留高置信分类，减少假阳性但可能丢失低丰度物种	要求高精度、不关心低丰度物种
0.5	非常严格，大量reads变为"未分类"	特殊验证用途

实际选择建议： - 常规分析：--confidence 0.05 或 --confidence 0.1 - 0.05 保留更多信息，适合探索性分析 - 0.2 更保守，适合需要高精度的发表级分析 - 具体选择还取决于数据库大小和样本类型

Bracken¶

为什么需要Bracken？

Kraken2 的分类结果存在一个问题：很多reads只能被分类到较高的分类层级（如属level），而不是种level。这是因为不同物种之间共享大量k-mer，Kraken2通过LCA算法会把这些reads推到它们的共同祖先节点。

Bracken（Bayesian Reestimation of Abundance with KrakEN）通过贝叶斯算法，将这些高层级的reads重新分配到种level，从而得到更准确的物种水平丰度估计。

Bracken 关键参数：

参数	含义	常用值
`-r`	read长度	150（与实际测序长度一致）
`-l`	分类层级	S（种）、G（属）、P（门）
`-t`	最小reads数阈值	10（少于此数的物种被过滤）

分类层级代码：

代码	层级	英文
D	域	Domain
P	门	Phylum
C	纲	Class
O	目	Order
F	科	Family
G	属	Genus
S	种	Species

2.6 功能注释（Functional Annotation）¶

为什么要做功能注释？

物种注释告诉我们"有哪些微生物"，功能注释则告诉我们"这些微生物在做什么"。功能注释将基因序列映射到已知的功能数据库，揭示微生物群落的代谢潜力。

HUMAnN3¶

HUMAnN3（The HMP Unified Metabolic Analysis Network 3）是宏基因组功能注释的主流工具。

工作流程： 1. 先用MetaPhlAn进行物种鉴定 2. 根据鉴定结果选择对应的泛基因组数据库 3. 用Bowtie2将reads比对到泛基因组 4. 未比对上的reads用DIAMOND比对到UniRef蛋白数据库 5. 整合结果，输出通路丰度和基因家族丰度

主要输出：

输出文件	内容	用途
genefamilies.tsv	基因家族丰度（UniRef90/UniRef50）	基因水平功能分析
pathabundance.tsv	代谢通路丰度（MetaCyc）	通路水平功能分析
pathcoverage.tsv	通路覆盖度	判断通路是否完整存在

eggNOG-mapper¶

eggNOG-mapper 通过将基因序列比对到eggNOG数据库，提供多种功能注释： - COG/KOG 分类：蛋白功能正交分类 - GO (Gene Ontology)：基因本体注释 - KEGG pathway：代谢通路注释 - CAZy：碳水化合物活性酶

HUMAnN3 vs eggNOG-mapper：

特点	HUMAnN3	eggNOG-mapper
输入	原始reads或clean reads	基因序列（通常需要先组装+预测）
速度	较慢	较快
数据库	MetaCyc + UniRef	eggNOG
注释类型	通路为主	多种注释体系（COG/GO/KEGG）
常用场景	基于reads的功能分析	基于基因的功能分析

2.7 多样性分析（Diversity Analysis）¶

Alpha多样性（Alpha Diversity）¶

Alpha多样性衡量的是单个样本内部的物种多样性。

常用指标：

指数	英文	含义	特点
Observed Species	观测到的物种数	样本中检测到的物种总数	最直观，但受测序深度影响大
Shannon	Shannon-Wiener Index	综合考虑物种丰富度和均匀度	最常用，值越大多样性越高
Simpson	Simpson Index	随机抽取两个个体属于不同物种的概率	对优势物种更敏感
Chao1	Chao1 Estimator	根据低丰度物种估计总物种数	常用于16S，宏基因组中较少用
Pielou	Pielou's Evenness	物种均匀度（0-1）	1表示完全均匀
ACE	ACE Estimator	基于丰度的覆盖度估计	与Chao1类似

通俗解释Shannon指数：

想象一个花园： - 花园A：100朵花，分别是50朵玫瑰、30朵百合、20朵菊花 → Shannon较高 - 花园B：100朵花，95朵玫瑰、3朵百合、2朵菊花 → Shannon较低 - Shannon指数同时考虑了"有多少种花"（丰富度）和"每种花的数量是否均衡"（均匀度）

Beta多样性（Beta Diversity）¶

Beta多样性衡量的是不同样本之间的物种组成差异。

常用距离/相异度指标：

度量	英文	特点	适用场景
Bray-Curtis	Bray-Curtis Dissimilarity	考虑丰度信息，范围0-1	最常用，适合丰度数据
Jaccard	Jaccard Distance	只考虑有无，不考虑丰度	关注物种组成差异
UniFrac (weighted)	Weighted UniFrac	考虑丰度+进化关系	16S数据（需要进化树）
UniFrac (unweighted)	Unweighted UniFrac	只考虑有无+进化关系	16S数据
Euclidean	Euclidean Distance	常见距离度量	一般用于标准化后的数据

常用可视化方法：

方法	英文	用途
PCoA	Principal Coordinates Analysis	展示样本间距离关系（最常用）
NMDS	Non-metric Multidimensional Scaling	展示样本排列，适合非线性关系
PCA	Principal Component Analysis	降维展示，适合标准化数据
热图	Heatmap	展示物种丰度矩阵

统计检验：

方法	英文	用途
PERMANOVA	Permutational MANOVA	检验组间Beta多样性是否有显著差异
ANOSIM	Analysis of Similarities	类似PERMANOVA，非参数检验
Wilcoxon	Wilcoxon Rank Sum Test	Alpha多样性组间比较
Kruskal-Wallis	Kruskal-Wallis H Test	多组Alpha多样性比较

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 0.23.4
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被丢弃
# --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 表示组间差异显著

4. 面试常问点¶

★ 宏基因组分析的基本流程是什么？¶

参考答案（结合实习经历）：

在百迈客实习期间，我参与了多个宏基因组项目的分析。基本流程如下：

首先是数据校验，用md5sum验证下机数据的完整性；然后用FastQC和MultiQC做质量评估，查看碱基质量、GC含量、接头残留等指标；接着用fastp做质量控制，去除接头序列和低质量reads；之后用Bowtie2将clean reads比对到人类参考基因组GRCh38进行去宿主，用samtools提取未比对上的微生物reads；再用Kraken2做物种分类注释，配合Bracken进行种水平丰度重估计；如果需要功能注释，会用HUMAnN3分析代谢通路；最后是多样性分析，包括Alpha多样性（Shannon、Simpson指数等）和Beta多样性（Bray-Curtis距离、PCoA等）。

在我的个人项目中，我还独立搭建了一个肠道宏基因组分析流程，从FastQC到多样性分析全流程自动化，通过shell脚本串联各个步骤。

★ 为什么要去宿主？不去会怎样？¶

参考答案：

去宿主是宏基因组分析中非常关键的一步。宏基因组测序的目标是微生物DNA，但样本采集时不可避免会混入宿主DNA，比如人肠道样本中可能含有10-50%的人源序列。

如果不去宿主，会造成几个问题： 1. 浪费计算资源：大量人源reads参与后续分析，既耗时又无意义 2. 干扰物种注释：人源序列可能被Kraken2错误地分类为微生物，产生假阳性结果 3. 影响丰度估计：分母（总reads数）被人源序列稀释，导致微生物真实丰度被低估 4. 功能注释污染：人类基因的功能信息混入微生物功能分析

我们一般用Bowtie2的 --very-sensitive 模式将reads比对到人类参考基因组GRCh38，然后用samtools提取双端都没比对上的reads对（flag 12），这些就是微生物reads。选择very-sensitive模式是因为去宿主时我们宁可多去一些也不要漏掉。

★ Kraken2的confidence参数怎么选？0.05和0.2有什么区别？¶

参考答案：

Kraken2的 --confidence 参数控制物种分类的严格程度，取值范围0到1。它的含义是：一条read被分类到某个物种时，支持该分类的k-mer比例必须达到这个阈值。

0.05（常用推荐值）：比较宽松，只要有5%以上的k-mer支持就接受分类。这样做的好处是灵敏度高，能检测到更多低丰度物种，但可能有一些假阳性。适合探索性分析和常规宏基因组分析。

0.2（较严格）：需要20%以上的k-mer支持，大幅减少假阳性，但会丢失一些低丰度物种或与近缘物种k-mer共享较多的species。适合需要高精度的发表级分析。

实际选择时要根据项目需求：如果是初步探索、想尽可能全面地了解群落组成，用0.05；如果是准备发文章、需要高可信度，可以用0.1或0.2。在百迈客的项目中，我们一般用0.05作为默认值，如果客户对假阳性敏感会适当提高。

需要注意的是，调高confidence不会改善Kraken2数据库本身的局限性，比如数据库中没有的物种无论confidence设多少都检测不到。

★ fastp的主要质控参数有哪些？你怎么设置的？¶

参考答案：

fastp是我们做宏基因组质控的主要工具，主要参数和我的常用设置如下：

接头去除：开启 --detect_adapter_for_pe 自动检测PE数据接头，也可以用 --adapter_sequence 手动指定Illumina TruSeq接头序列

碱基质量过滤：--qualified_quality_phred 20 表示质量值低于Q20的碱基被视为不合格（Q20对应1%的错误率）；--unqualified_percent_limit 40 表示一条read中不合格碱基超过40%就整条丢弃

长度过滤：--length_required 50，处理后短于50bp的read丢弃，因为太短的read比对效率低且容易产生假阳性

滑窗修剪：--cut_front 和 --cut_tail 开启首尾滑窗修剪，配合 --cut_window_size 4 和 --cut_mean_quality 20，在读取的5'端和3'端用4bp滑窗检测，平均质量低于20就切掉

N碱基过滤：--n_base_limit 5，含N碱基超过5个的read丢弃

碱基纠错：--correction 利用PE reads重叠区域进行碱基校正

一般经过这些参数处理后，reads保留率在90%以上就算正常。fastp还自带HTML报告，可以直观看到质控前后的效果对比。

★ 如何判断测序数据质量好不好？（FastQC报告解读）¶

参考答案：

拿到FastQC报告后，我会重点关注以下几个模块：

Per base sequence quality（最重要）：看每个位置的碱基质量分布箱线图。正常情况下Q值应该在20以上，Q30以上的碱基占比应该超过85%。如果3'端质量急剧下降，说明需要做trimming

Adapter content：看接头残留情况。如果曲线明显上升，说明有接头污染，质控时需要去接头

Per sequence GC content：看GC含量分布是否接近正态分布。如果出现多峰或严重偏移，可能有其他物种DNA污染

Sequence length distribution：正常的原始数据应该集中在固定长度（如150bp）。如果分布不均一可能有问题

Overrepresented sequences：如果有大量高频重复序列，可能是接头残留或rRNA污染

Per base sequence content：前几个位置ATCG比例有波动是正常的（random primer导致），但后面应该趋于稳定

总体判断标准：Q20>95%、Q30>85%、接头含量<5%、GC分布正常、无大量重复序列，就认为数据质量是好的。如果是多个样本，我会用MultiQC汇总查看，方便横向比较发现异常样本。

★ reads保留率多少算正常？¶

参考答案：

reads保留率需要分两个阶段看：

fastp质控后：正常保留率应该在90%以上。如果低于85%，说明原始数据质量较差，可能测序出了问题，需要跟测序公司沟通。保留率过低也可能是质控参数设置过严导致的

去宿主后：这个取决于样本类型。人肠道样本一般保留70-99%，因为采样部位和方法不同，人源DNA含量差异很大（粪便样本人源含量通常较低，活检组织样本人源含量高）。土壤、水体等环境样本去宿主后保留率通常>95%，因为基本没有宿主DNA

在实际工作中，如果fastp后保留率低于85%或去宿主后保留率异常低（如人肠道样本低于50%），我会回头检查原始数据质量或参数设置，必要时联系测序公司确认是否需要重测。

★ 宏基因组和16S扩增子的区别？¶

参考答案：

这是两种不同的微生物群落研究方法，主要区别如下：

对比维度宏基因组（Metagenomics） 16S扩增子（16S Amplicon）

测序对象 样本中所有DNA（全基因组）仅16S rRNA基因的特定区域（如V3-V4）

信息量 基因组+功能信息仅物种分类信息

分辨率 可达种甚至菌株水平通常到属水平，种水平不太可靠

功能分析 可以做（HUMAnN3等）只能做预测（PICRUSt2，不够准确）

成本较高（需要大量测序数据）较低（测序量小）

分析复杂度 较高（计算资源需求大）较低（分析流程相对简单）

覆盖范围 细菌、古菌、真菌、病毒等主要是细菌和古菌（16S基因特异）

偏好性 偏好小（直接测全基因组）受引物偏好影响（不同引物扩增效率不同）

常用工具 Kraken2, HUMAnN3, MEGAHIT等 QIIME2, DADA2, mothur等

常见数据量 每样本5-20G 每样本5-10万条reads

简单总结：16S像是"人口普查"——知道有哪些微生物；宏基因组像是"全面调查"——不仅知道有谁，还知道它们在做什么。

在百迈客实习时，两种方法的项目我都做过。16S项目主要是用QIIME2做OTU聚类和多样性分析；宏基因组项目流程更复杂，但获得的信息更丰富。客户预算有限时推荐16S，需要功能信息或高分辨率分类时推荐宏基因组。

★ Kraken2和MetaPhlAn的区别？（加分题）¶

参考答案：

对比 Kraken2 MetaPhlAn

原理 k-mer精确匹配 clade-specific marker gene比对

数据库 全基因组（RefSeq）标记基因

速度非常快较慢

内存需求 较大（标准库约70GB）较小（几GB）

灵敏度 高（但假阳性较多）较低（但精确度高）

丰度估计 需配合Bracken 自带丰度估计

适用场景 快速筛查、大规模样本精确定量、小样本深入分析

★ 组装（Assembly）在宏基因组中的作用是什么？（加分题）¶

参考答案：

宏基因组分析有两条路线： 1. 基于reads的分析（read-based）：直接用reads做物种注释和功能注释，如Kraken2、HUMAnN3 2. 基于组装的分析（assembly-based）：先将reads拼接成更长的contigs，再做基因预测和注释

组装常用工具如 MEGAHIT 或 metaSPAdes，组装后可以用 Prodigal 做基因预测，再用 eggNOG-mapper 或 DIAMOND 做功能注释。

组装的优势：得到更完整的基因序列，功能注释更准确；可以做binning（分箱）获得MAGs（Metagenome-Assembled Genomes），即直接从宏基因组数据中恢复单个微生物的基因组。

组装的劣势：低丰度物种可能组装不出来；嵌合体contigs（chimeric contigs）问题；计算资源需求大。

5. 易错/易混淆点¶

⚠️ 1. Kraken2和Bracken的区别¶

易混淆点： 初学者常把两者混为一谈，或者只用Kraken2不用Bracken。

正确理解： - Kraken2 是分类器（classifier）：将每条read分配到一个分类单元，但很多reads只能到属甚至更高层级 - Bracken 是丰度估计器（abundance estimator）：在Kraken2结果基础上，用贝叶斯算法把高层级reads重新分配到种水平 - 两者是配合使用的关系：先Kraken2分类，再Bracken估计丰度 - 如果只用Kraken2，种水平的丰度会严重低估（因为大量reads停留在属及以上层级）

⚠️ 2. 去宿主前后reads数变化的解读¶

易错点： 看到去宿主后reads大幅减少就认为"数据不好"。

正确理解： - 人肠道样本的宿主比例变异很大（5-50%），粪便样本通常人源含量较低，但活检或灌洗液样本可能很高 - 去宿主后reads减少不一定是坏事，反而说明成功去除了大量人源污染 - 需要关注的是去宿主后剩余的微生物reads数量是否足够做后续分析（通常需要至少几百万条） - 如果去宿主后保留率极低（如<30%），才需要考虑是否采样有问题

⚠️ 3. fastp的--qualified_quality_phred含义¶

易混淆点： 误以为 --qualified_quality_phred 20 表示"所有低于Q20的碱基都直接去掉"。

正确理解： - 这个参数定义的是"什么叫不合格碱基"的阈值，低于Q20的碱基被标记为不合格 - 但仅仅标记为不合格不会直接丢弃，还要看 --unqualified_percent_limit - 只有当一条read中不合格碱基的比例超过了 --unqualified_percent_limit（默认40%），这条read才会被丢弃 - 真正负责"切掉低质量碱基"的是 --cut_front、--cut_tail 等滑窗参数 - 理解三者关系：qualified_quality_phred 定义"什么是低质量"；unqualified_percent_limit 决定"低质量碱基太多时丢弃整条read"；cut_front/cut_tail 负责"修剪掉末端的低质量区域"

⚠️ 4. 数据库版本不一致的影响¶

易错点： 不同样本或不同时间的分析使用了不同版本的数据库（如Kraken2标准库），然后把结果直接合并比较。

正确理解： - 数据库更新会增加或修改物种分类信息，不同版本数据库的注释结果不具有直接可比性 - 同一批项目的所有样本必须使用完全相同版本的数据库 - 建议在项目开始时记录数据库版本和下载日期 - 发表文章时必须注明使用的数据库版本（如 "Kraken2 standard database, 2023-06-05 release"） - 如果需要重分析历史数据，应该用新数据库对所有样本重新分析，而不是只分析新样本

⚠️ 5. SE和PE数据处理差异¶

易混淆点： 用PE数据的参数和命令直接处理SE数据，或者反过来。

正确理解： - PE（Paired-End）：一条DNA片段从两端分别测序，产生R1和R2两个文件 - SE（Single-End）：只从一端测序，产生一个文件 - fastp处理SE数据不需要 -I/-O 参数和 --detect_adapter_for_pe - Bowtie2处理SE数据用 -U 参数而不是 -1/-2 - samtools提取未比对reads时，SE数据用 -f 4 而不是 -f 12 - Kraken2处理SE数据不需要 --paired 参数 - 宏基因组分析绝大多数情况是PE数据，SE在宏基因组中很少用

⚠️ 6. 比对率 ≠ 分类率¶

易混淆点： 混淆 Bowtie2 的比对率和 Kraken2 的分类率。

正确理解： - Bowtie2比对率：reads比对到宿主参考基因组的比例 → 这个越高，说明宿主污染越严重 - Kraken2分类率：reads被成功分类到物种的比例 → 这个越高，说明数据库覆盖越好 - 两者完全不同的概念，不要混淆 - Kraken2分类率通常在50-80%，有20-50%的reads无法分类是正常的（数据库不可能覆盖所有微生物）

⚠️ 7. --confidence 不等于准确率¶

易混淆点： 以为设置 --confidence 0.95 就能保证95%的准确率。

正确理解： - confidence 是指支持某分类的k-mer占该read所有k-mer的比例阈值 - 它是一个过滤标准，不是准确率的保证 - 设太高（如0.5以上）会导致大量reads变为"未分类"，可能丢失有价值的信息 - 准确率还受数据库质量、物种相似度等因素影响 - 选择合适的confidence需要在灵敏度和精确度之间权衡

⚠️ 8. 多样性分析中丰度表的标准化¶

易混淆点： 直接用原始reads计数做多样性分析，不做标准化。

正确理解： - 不同样本的测序深度（总reads数）可能差异很大 - 不做标准化，测序深度高的样本会天然显示更高的物种丰富度 - Alpha多样性分析前通常需要做稀释（rarefaction）或标准化 - Beta多样性分析中，Bray-Curtis距离自带相对丰度标准化，但最好还是先将计数转为相对丰度 - 不同的标准化方法（rarefaction、TSS、CSS、TMM等）可能导致不同结果，需要根据具体情况选择

6. 补充知识¶

6.1 常见文件格式¶

格式	后缀	内容	用途
FASTQ	.fq / .fastq	reads序列 + 质量值	原始/清洗后的测序数据
FASTA	.fa / .fasta	序列（无质量值）	参考基因组、组装结果
SAM	.sam	序列比对结果（文本）	比对信息，可读
BAM	.bam	序列比对结果（二进制）	SAM的压缩版，常用
GFF/GTF	.gff / .gtf	基因注释信息	基因坐标、功能注释
TSV/CSV	.tsv / .csv	表格数据	丰度表、统计结果
BIOM	.biom	生物观测矩阵	QIIME2等工具的标准格式

6.2 质量值（Phred Quality Score）换算¶

Phred值	错误率	准确率	含义
Q10	10%	90%	每10个碱基有1个错误
Q20	1%	99%	每100个碱基有1个错误
Q30	0.1%	99.9%	每1000个碱基有1个错误
Q40	0.01%	99.99%	每10000个碱基有1个错误

计算公式：Q = -10 * log10(P)，其中P是错误概率

6.3 常用Linux命令速查¶

# 查看fastq文件的reads数（每4行一条read）
zcat sample.fq.gz | wc -l | awk '{print $1/4}'

# 查看fastq文件前几条reads
zcat sample.fq.gz | head -12   # 查看前3条reads（每条4行）

# 统计文件大小
du -sh *.fq.gz

# 后台运行并记录日志
nohup bash pipeline.sh > pipeline.log 2>&1 &

# 查看后台任务
jobs -l

# 查看磁盘空间（分析前确认空间足够）
df -h

# 统计某个物种在Kraken2 report中的丰度
grep "Escherichia coli" sample_kraken2.report

6.4 分析流程时间和资源参考¶

步骤	单样本时间（PE150，10G数据）	内存需求	磁盘需求
FastQC	5-10分钟	1GB	很小
fastp	5-15分钟	2-4GB	与输入相当
Bowtie2建索引	1-2小时（一次性）	8GB+	~10GB
Bowtie2比对去宿主	30-60分钟	4-8GB	中间文件较大
Kraken2	10-30分钟	50-70GB（标准库）	较小
Bracken	1-2分钟	很小	很小
HUMAnN3	2-6小时	16GB+	较大

注意：Kraken2标准数据库需要大量内存（约70GB），这是宏基因组分析中对服务器配置要求最高的步骤之一。内存不足时可以使用MiniKraken数据库（约8GB）。

快速复习卡片¶

流程口诀：
校（md5校验）→ 评（FastQC评估）→ 控（fastp质控）→ 去（Bowtie2去宿主）→ 分（Kraken2分类）→ 功（功能注释）→ 多（多样性分析）

关键数字：
- Q20 > 95%, Q30 > 85%
- fastp保留率 > 90%
- Kraken2 confidence: 0.05（常规）/ 0.2（严格）
- 去宿主用 --very-sensitive
- samtools -f 12（PE双端未比对）
- Bracken -r 150 -l S -t 10

🔄 最新版本动态（2026年4月更新）¶

工具	最新版本	关键更新
Kraken2	v2.17.0	新增PrackenDB数据库（2026.01），支持多NCBI数据库同时分类，FASTA/Q解析器换用kseq大幅加速
Kraken2数据库	PrackenDB (2025.10)	包含所有NCBI参考组装（GenBank+RefSeq）的细菌/古菌/原生生物/真菌，加人基因组和病毒
Bracken	配合Kraken2	Standard-8/PlusPF-8/PlusPFP-8预构建数据库仅需~8GB内存
fastp	v0.24+	持续维护中，仍是最主流的质控工具

📌 面试加分：提到Kraken2最新版v2.17.0和PrackenDB，说明你关注工具更新。Standard-8数据库只需8GB内存，适合内存有限的服务器。

最后更新：2026-04-27 适用岗位：生物信息分析工程师关联知识库：02_微生物多样性分析.md、03_功能注释数据库.md

对比维度	宏基因组（Metagenomics）	16S扩增子（16S Amplicon）
测序对象	样本中所有DNA（全基因组）	仅16S rRNA基因的特定区域（如V3-V4）
信息量	基因组+功能信息	仅物种分类信息
分辨率	可达种甚至菌株水平	通常到属水平，种水平不太可靠
功能分析	可以做（HUMAnN3等）	只能做预测（PICRUSt2，不够准确）
成本	较高（需要大量测序数据）	较低（测序量小）
分析复杂度	较高（计算资源需求大）	较低（分析流程相对简单）
覆盖范围	细菌、古菌、真菌、病毒等	主要是细菌和古菌（16S基因特异）
偏好性	偏好小（直接测全基因组）	受引物偏好影响（不同引物扩增效率不同）
常用工具	Kraken2, HUMAnN3, MEGAHIT等	QIIME2, DADA2, mothur等
常见数据量	每样本5-20G	每样本5-10万条reads

对比	Kraken2	MetaPhlAn
原理	k-mer精确匹配	clade-specific marker gene比对
数据库	全基因组（RefSeq）	标记基因
速度	非常快	较慢
内存需求	较大（标准库约70GB）	较小（几GB）
灵敏度	高（但假阳性较多）	较低（但精确度高）
丰度估计	需配合Bracken	自带丰度估计
适用场景	快速筛查、大规模样本	精确定量、小样本深入分析