14. 比对与组装工具

一句话说明

比对（Alignment）是把测序reads"贴回"参考基因组找位置，组装（Assembly）是没有参考时把reads像拼图一样拼成长序列——二者是宏基因组分析的基础操作。

---

核心概念（白话版）

什么是序列比对（Alignment）

白话比方： 想象你有一本完整的书（参考基因组），现在有人把书的某些句子撕下来打乱了（测序reads）。比对就是把这些碎纸条一张张找回原书中对应的位置。

技术定义： 将测序产生的短序列（reads）与已知的参考基因组进行匹配，确定每条read来自基因组的哪个位置。

什么是基因组组装（Assembly）

白话比方： 你没有原书，只有一堆碎纸条。你要根据纸条之间的重叠部分，把它们拼回完整的文章。这就是从头组装（de novo assembly）。

技术定义： 不依赖参考基因组，利用reads之间的序列重叠关系，将短序列拼接成较长的连续序列（contigs）。

比对 vs 组装：什么时候用哪个

场景	用比对	用组装
有参考基因组	是	-
没有参考基因组	-	是
宏基因组（混合菌群）	去宿主/定量时用比对	拼基因组时用组装
变异检测（SNP/InDel）	是	-
发现新基因/新物种	-	是
你的T2D项目	Bowtie2比对去宿主	MEGAHIT组装菌群基因组

参考基因组比对 vs de novo 组装

• 参考比对：速度快，适合已有参考的物种（如人、大肠杆菌）
• de novo 组装：计算量大，但能发现全新序列，宏基因组必须用这个来获得MAGs（宏基因组组装基因组，就是从混合菌群数据里拼出来的"单个细菌的近完整基因组"）

---

常用比对工具详解

1. BWA / BWA-MEM2（全基因组比对的"黄金标准"）

原理白话： BWA用一种叫BWT（Burrows-Wheeler Transform，相当于把书里所有句子按字母重新排列，查找起来快得多）的算法把参考基因组压缩建索引，就像给书做了一个超级目录。查找时通过这个目录快速定位。BWA-MEM2是BWA的加速版，利用CPU的SIMD（CPU的一种加速技巧，一次能同时处理多条数据）指令集，速度提升1.3-3.1倍，结果完全一样。

适用场景： - 全基因组重测序（WGS） - 变异检测（SNP calling） - 外显子组测序（WES，基因组中编码蛋白质的那一小部分，约1-2%）

优缺点： | 优点 | 缺点 | |------|------| | 准确度高，是GATK（基因组分析工具包，变异检测的行业标准软件）流程标配 | 不支持剪接比对（不适合RNA-seq） | | BWA-MEM2速度提升显著 | 索引文件较大（人基因组~10GB） | | 支持PE（双端测序）reads | 对超短reads（<70bp）表现一般 |

命令示例：

# ============ BWA-MEM2 比对流程 ============

# 第一步：建索引（把参考基因组做成"目录"，只需做一次）
# -p prefix 指定输出前缀名
bwa-mem2 index -p human_ref human_g1k_v37.fasta

# 第二步：比对（把reads贴回参考基因组）
# -t 8 表示用8个CPU线程加速
# -R 添加Read Group信息（GATK必须要）
# human_ref 是索引前缀
# sample_R1.fq.gz 和 sample_R2.fq.gz 是双端reads
bwa-mem2 mem -t 8 \
  -R '@RG\tID:sample1\tSM:sample1\tPL:ILLUMINA' \
  human_ref \
  sample_R1.fq.gz sample_R2.fq.gz \
  > aligned.sam

---

2. Bowtie2（短读比对，你项目的核心工具）

原理白话： Bowtie2也是基于BWT/FM-index（一种基于BWT的快速查找索引），但它的特点是索引文件特别小、内存占用低。它支持局部比对（local alignment），就像允许碎纸条只有一部分能贴上原书也算匹配成功。

为什么你的T2D项目选了Bowtie2： - 宏基因组分析中，"去宿主"这一步需要把人的reads过滤掉 - Bowtie2内存占用小（人基因组索引仅~3.2GB），适合普通服务器 - 速度快，支持 --un-conc 参数直接输出未比对的reads（即微生物reads） - 在宏基因组领域使用率极高，很多pipeline默认用Bowtie2去宿主

适用场景： - 宏基因组去宿主污染（你的项目就是这个） - ChIP-seq（看蛋白质在基因组哪些位置结合的技术）比对 - 16S/ITS扩增子比对 - 一般短读DNA比对

优缺点： | 优点 | 缺点 | |------|------| | 索引文件小（人基因组~3.2GB） | 对长reads（>500bp）不如minimap2 | | 内存友好，普通电脑就能跑 | 不支持剪接比对 | | --un-conc可直接输出未比对reads | 速度略慢于BWA-MEM2 | | 支持local模式处理部分比对 | - | | 宏基因组领域标准工具 | - |

命令示例（完整去宿主流程）：

# ============ Bowtie2 去宿主比对流程（你T2D项目的核心步骤） ============

# 第一步：建索引（给人参考基因组建"目录"）
# bowtie2-build <参考基因组fasta> <输出索引前缀>
# 只需要运行一次，之后可以重复使用
bowtie2-build human_genome.fa human_index

# 第二步：比对 + 去宿主
# --very-sensitive 使用最灵敏模式（慢一点但更准，少漏人源reads）
# --threads 8 用8个线程
# -x human_index 指定索引前缀
# -1 和 -2 指定双端reads输入
# --un-conc-gz 把"没比上"的reads（即微生物reads）输出为gz压缩文件
# -S 输出SAM格式比对结果
bowtie2 --very-sensitive \
  --threads 8 \
  -x human_index \
  -1 clean_R1.fq.gz \
  -2 clean_R2.fq.gz \
  --un-conc-gz microbe_reads_%.fq.gz \
  -S host_aligned.sam

# 说明：
# microbe_reads_1.fq.gz 和 microbe_reads_2.fq.gz 就是去掉人源后的微生物reads
# 这些reads接下来送去做宏基因组组装或物种注释

# 第三步：看比对统计（多少reads比上了人基因组？）
# 比对率越高说明宿主污染越严重
# 一般粪便样本宿主比对率在 1-30% 不等

# ---------- 其他常用参数 ----------

# end-to-end模式（默认，要求整条read都比上）
bowtie2 -x index -1 r1.fq -2 r2.fq -S out.sam

# local模式（允许reads两端软裁剪（read两头有碱基没比上参考，但保留在数据里不删），更宽松）
bowtie2 --local -x index -1 r1.fq -2 r2.fq -S out.sam

# 快速模式（牺牲灵敏度换速度）
bowtie2 --very-fast -x index -1 r1.fq -2 r2.fq -S out.sam

---

3. HISAT2（转录组比对专用）

原理白话： HISAT2专门为RNA-seq设计。RNA经过剪接（splicing），一条mRNA可能来自基因的多个外显子拼接。HISAT2能"跳着比"，识别剪接位点。就像碎纸条上的文字在原书中不是连续的几行，而是跨了好几页拼起来的。

适用场景： - RNA-seq转录组分析 - 差异表达基因分析 - 可变剪接分析

优缺点： | 优点 | 缺点 | |------|------| | 支持剪接比对（splice-aware） | 只适合RNA-seq，不适合DNA | | 索引极小（人基因组~8GB，含SNP和transcript信息） | 对长reads支持不如minimap2 | | 速度快于STAR，内存低于STAR | 转录组专用，范围窄 |

命令示例：

# ============ HISAT2 转录组比对 ============

# 下载或构建索引（HISAT2官网提供预建索引）
# hisat2-build <参考基因组> <索引前缀>
hisat2-build genome.fa genome_index

# 比对RNA-seq数据
# --dta 为下游StringTie等转录本组装工具优化输出
# --rna-strandness RF 指定链特异性文库类型
hisat2 --dta \
  --rna-strandness RF \
  -p 8 \
  -x genome_index \
  -1 rnaseq_R1.fq.gz \
  -2 rnaseq_R2.fq.gz \
  -S transcriptome.sam

---

4. minimap2（三代长读比对的"瑞士军刀"）

原理白话： minimap2用minimizer（最小化子，从序列里挑出有代表性的小片段作为关键词，先粗定位再精确比对）做种子，然后做链式延伸。它就像一个万能工具，既能处理错误率高达15%的三代长reads（几千到几万bp），也能比对短reads，甚至做RNA剪接比对。一个工具覆盖多个场景。

适用场景： - PacBio CLR/HiFi reads比对 - Oxford Nanopore reads比对 - 长读RNA-seq（Iso-Seq） - 基因组间比对（assembly-to-assembly） - 也支持Illumina短读（-ax sr模式）

优缺点： | 优点 | 缺点 | |------|------| | 极快（比BWA-MEM快几十倍处理长reads） | 对短reads准确度不如专用工具 | | 支持多种数据类型（预设模式丰富） | 不擅长检测小外显子 | | 容错率高，15%错误率也能比 | 索引参数不同数据要分开建 | | 2200+ GitHub stars，社区活跃 | - |

命令示例：

# ============ minimap2 长读比对 ============

# PacBio HiFi reads比对（错误率<1%的高质量长reads）
minimap2 -ax map-hifi ref.fa pacbio_hifi.fq.gz > aln_hifi.sam

# Oxford Nanopore reads比对（错误率较高的长reads）
minimap2 -ax map-ont ref.fa nanopore.fq.gz > aln_ont.sam

# 短reads比对（Illumina，但一般还是推荐BWA/Bowtie2）
minimap2 -ax sr ref.fa read1.fq read2.fq > aln_sr.sam

# 长读RNA-seq（剪接比对）
minimap2 -ax splice ref.fa iso_seq.fq > aln_rna.sam

# 基因组间比对（两个组装结果互相比）
minimap2 -ax asm5 ref.fa assembly.fa > asm_aln.sam

---

工具选择对比表

工具	适用数据	读长	内存(人基因组)	速度	核心场景
BWA-MEM2	Illumina DNA	70-300bp	~10GB	快	WGS变异检测
Bowtie2	Illumina DNA	50-300bp	~3.2GB	中	宏基因组去宿主、ChIP-seq
HISAT2	Illumina RNA	50-300bp	~8GB	快	RNA-seq转录组
minimap2	PacBio/ONT/all	1kb-100kb+	~6GB	极快(长读)	三代测序、全场景

选择口诀： - DNA短读 + 变异检测 → BWA-MEM2 - DNA短读 + 宏基因组去宿主 → Bowtie2（你的项目） - RNA短读 + 转录组 → HISAT2 - 三代长读（PacBio/Nanopore） → minimap2

---

常用组装工具详解

1. MEGAHIT（宏基因组组装首选）

原理白话： MEGAHIT用"简洁de Bruijn图"（succinct de Bruijn graph）来组装。把reads切成k-mer（长度为k的小片段），然后像搭积木一样，通过k-mer之间的重叠关系拼成长序列。MEGAHIT的特别之处是用多个k值（从小到大），先用小k捕获低丰度物种，再用大k提高高丰度物种的连续性。

适用场景： - 宏基因组组装（你的T2D项目核心步骤） - 大规模数据集（内存极省） - 土壤/肠道等复杂样本

优缺点： | 优点 | 缺点 | |------|------| | 内存极低（比SPAdes省10倍以上） | 连续性略低于metaSPAdes | | 速度快 | 对低丰度物种可能丢失 | | 专为宏基因组优化 | 不输出scaffold graph | | 支持超大数据集 | - |

命令示例：

# ============ MEGAHIT 宏基因组组装 ============

# 基本用法：双端reads输入
# -1 和 -2 分别是正向和反向reads
# -o 输出目录
# -t 线程数
# --min-contig-len 500 只保留≥500bp的contig（太短的没意义）
megahit \
  -1 microbe_R1.fq.gz \
  -2 microbe_R2.fq.gz \
  -o megahit_output \
  -t 16 \
  --min-contig-len 500

# 复杂宏基因组（如土壤）使用预设
megahit --presets meta-large \
  -1 soil_R1.fq.gz -2 soil_R2.fq.gz \
  -o soil_assembly

# 多个样本共组装（co-assembly，增加深度）
megahit \
  -1 s1_R1.fq.gz,s2_R1.fq.gz,s3_R1.fq.gz \
  -2 s1_R2.fq.gz,s2_R2.fq.gz,s3_R2.fq.gz \
  -o co_assembly

# 输出文件：megahit_output/final.contigs.fa 就是组装结果

2. SPAdes / metaSPAdes

原理白话： SPAdes也用de Bruijn图，但它比MEGAHIT更"精细"——会做错误纠正、用paired-end信息做scaffold（把多个contig按顺序串起来，中间允许有空缺）。metaSPAdes是它的宏基因组版本，针对不均匀覆盖度做了优化。

适用场景： - 单个细菌基因组组装（SPAdes） - 小规模宏基因组样本（metaSPAdes） - 需要高连续性的场景

优缺点： | 优点 | 缺点 | |------|------| | 组装连续性高（N50更大，想象积木按大到小排列，堆到总重量一半时那块积木的重量就是N50） | 内存需求巨大（几百GB） | | 内置错误纠正 | 速度慢 | | scaffold质量好 | 不适合超大数据集 |

命令示例：

# ============ metaSPAdes 宏基因组组装 ============

# 基本用法
# metaspades.py 本身就是宏基因组模式，不需要加 --meta
# （另一种等价写法：spades.py --meta，两者效果完全一样）
# --only-assembler 跳过错误纠正（已用fastp质控过就不需要）
metaspades.py \
  --only-assembler \
  -1 microbe_R1.fq.gz \
  -2 microbe_R2.fq.gz \
  -t 16 \
  -m 200 \
  -o spades_output

# 输出文件：spades_output/scaffolds.fasta

MEGAHIT vs metaSPAdes 选择： | 条件 | 选MEGAHIT | 选metaSPAdes | |------|-----------|--------------| | 内存 < 64GB | 是 | 否 | | 数据量 > 50GB | 是 | 否 | | 要求最高连续性 | 否 | 是 | | 多样本共组装 | 是 | 否 | | 你的T2D项目 | 是 | 否 |

3. 组装质量评估

# ============ QUAST：评估contig统计指标 ============
# -o 输出目录
# -t 线程
# --min-contig 500 只统计≥500bp的contig
quast.py \
  megahit_output/final.contigs.fa \
  -o quast_report \
  -t 8 \
  --min-contig 500
# 关注指标：N50（越大越好）、总长度、contig数量

# ============ CheckM：评估MAGs（宏基因组组装基因组）质量 ============
# 看MAGs的完整度（completeness）和污染度（contamination）
checkm lineage_wf \
  bins_dir/ \
  checkm_output/ \
  -t 16 \
  -x fa
# 高质量MAG标准：完整度>90%，污染<5%

---

SAM/BAM 格式详解

SAM（Sequence Alignment/Map）是比对结果的标准文本格式，BAM是它的二进制压缩版。

SAM文件核心字段（11个必须字段）

列号	字段名	含义	白话解释
1	QNAME	读序名称	这条read叫什么名字
2	FLAG	比对标志位	用数字编码的状态信息（是否配对、是否反向互补等）
3	RNAME	参考序列名	贴在了哪条染色体上
4	POS	位置（1-based）	贴在染色体的第几个碱基位置
5	MAPQ	比对质量	这次比对有多可靠（0-60，越高越好）
6	CIGAR	比对操作串	描述如何对齐（M匹配、I插入、D缺失、S软裁剪（read两头有碱基没比上参考，但保留在数据里不删））
7	RNEXT	配对read的染色体	另一端在哪
8	PNEXT	配对read的位置	另一端在第几个碱基
9	TLEN	模板长度	两端reads之间的距离（插入片段大小）
10	SEQ	序列	read本身的碱基序列
11	QUAL	质量值	每个碱基的测序质量

常用FLAG值速查

FLAG	含义
4	read未比对上（unmapped）
16	read比对到反向互补链
99	正常配对，正向read是第一条
147	正常配对，反向read是第二条

常用 samtools 命令

# ============ samtools 常用操作 ============

# SAM转BAM（压缩，节省空间）
# -b 输出BAM格式
# -S 输入是SAM（新版本可省略）
# -@ 8 用8个线程
samtools view -bS -@ 8 aligned.sam > aligned.bam

# 排序（按染色体位置排序，后续分析必须）
samtools sort -@ 8 aligned.bam -o sorted.bam

# 建索引（很多工具要求BAM必须有索引）
samtools index sorted.bam

# 查看比对统计
samtools flagstat sorted.bam
# 输出示例：
# 10000000 + 0 in total (所有reads数)
# 9500000 + 0 mapped (95.00% : mapped) （比对上的百分比）

# 只提取未比对的reads（用于提取微生物reads）
# -f 4 表示FLAG包含4（未比对）
samtools view -b -f 4 sorted.bam > unmapped.bam

# 只提取已比对的reads
# -F 4 表示FLAG不包含4（即已比对）
samtools view -b -F 4 sorted.bam > mapped.bam

# 查看某个区域的比对情况
samtools view sorted.bam chr1:1000000-1001000

# 计算覆盖深度
samtools depth sorted.bam > depth.txt

---

实操代码：完整比对流程（建索引→比对→排序→统计）

#!/bin/bash
# ============================================================
# 完整比对流程演示（以Bowtie2去宿主为例，模拟你的T2D项目）
# 输入：质控后的reads（已用fastp处理）
# 输出：去除人源后的微生物reads
# ============================================================

# ---------- 0. 设置变量 ----------
REF="human_genome_hg38.fa"          # 人参考基因组
INDEX="bt2_human_index"             # Bowtie2索引前缀
R1="clean_R1.fq.gz"                 # fastp质控后的正向reads
R2="clean_R2.fq.gz"                 # fastp质控后的反向reads
THREADS=8                           # CPU线程数
OUTDIR="alignment_results"          # 输出目录

mkdir -p ${OUTDIR}

# ---------- 1. 建索引（只需运行一次） ----------
echo "Step 1: Building Bowtie2 index..."
bowtie2-build --threads ${THREADS} ${REF} ${INDEX}
# 输出：human_index.1.bt2, .2.bt2, .3.bt2, .4.bt2, .rev.1.bt2, .rev.2.bt2

# ---------- 2. 比对（核心步骤） ----------
echo "Step 2: Aligning reads to human genome..."
bowtie2 --very-sensitive \
  --threads ${THREADS} \
  -x ${INDEX} \
  -1 ${R1} \
  -2 ${R2} \
  --un-conc-gz ${OUTDIR}/microbe_%.fq.gz \
  -S ${OUTDIR}/host_aligned.sam \
  2> ${OUTDIR}/bowtie2_log.txt
# --un-conc-gz 输出未比对reads = 微生物reads
# 2> 把比对统计信息保存到log文件

# ---------- 3. SAM转BAM + 排序 ----------
echo "Step 3: Converting and sorting..."
samtools view -bS -@ ${THREADS} ${OUTDIR}/host_aligned.sam | \
  samtools sort -@ ${THREADS} -o ${OUTDIR}/sorted.bam
# 管道操作：SAM → BAM → 排序，一步到位

# ---------- 4. 建BAM索引 ----------
samtools index ${OUTDIR}/sorted.bam

# ---------- 5. 统计比对结果 ----------
echo "Step 5: Generating alignment statistics..."
samtools flagstat ${OUTDIR}/sorted.bam > ${OUTDIR}/flagstat.txt

# 查看结果
echo "=== 比对统计 ==="
cat ${OUTDIR}/flagstat.txt
echo ""
echo "=== Bowtie2日志 ==="
cat ${OUTDIR}/bowtie2_log.txt

# ---------- 6. 清理中间文件（节省磁盘） ----------
rm ${OUTDIR}/host_aligned.sam  # SAM文件很大，已转BAM可删除

echo "Done! 微生物reads在: ${OUTDIR}/microbe_1.fq.gz 和 microbe_2.fq.gz"

---

和你项目的关联：T2D 项目中 Bowtie2 做了什么

在你的2型糖尿病肠道菌群宏基因组项目中，Bowtie2 的角色是去宿主污染：

1. 为什么要去宿主？ 粪便样本测序后，有一部分reads来自人的肠道上皮细胞DNA，不是细菌的。如果不去掉，后续分析会把人的序列误认为是微生物。

2. 流程位置： fastp质控 → Bowtie2去宿主 → MEGAHIT组装 / Kraken2物种注释

3. 你用的参数： --very-sensitive 模式确保尽量把人源reads都找出来，宁可多去也不能漏。

4. 结果解读： 如果你的样本宿主比对率是15%，说明85%的reads是微生物来源的，数据质量不错。

5. 面试加分点： 能说出"为什么选Bowtie2而不是BWA"——因为Bowtie2的 --un-conc 参数可以直接导出未比对reads，操作方便；内存小适合批量处理多个样本。

---

面试怎么答

Q1：比对（alignment）和组装（assembly）有什么区别？

比对是"有答案对答案"——你有参考基因组，把reads找到对应位置。组装是"没答案自己拼"——没有参考，靠reads之间的重叠关系拼出连续序列。我的T2D项目两个都用了：Bowtie2做比对去宿主，MEGAHIT做组装拼菌群基因组。

Q2：你为什么选Bowtie2做去宿主而不是BWA？

三个原因：第一，Bowtie2有 --un-conc-gz 参数可以直接输出未比对的reads，操作一步到位；第二，Bowtie2索引只有3.2GB，内存友好，我们有几十个样本要批量处理；第三，在宏基因组去宿主这个场景下，Bowtie2是领域内最常用的工具，很多发表的pipeline都用它。

Q3：SAM和BAM有什么区别？MAPQ是什么？

SAM是文本格式，人能直接看但文件大；BAM是二进制压缩版，小很多但需要samtools解析。MAPQ是比对质量分数，0到60分，表示这条read比对位置的可靠程度。MAPQ=0表示可能比到了多个位置（多比对），MAPQ=60表示位置唯一确定很可靠。一般过滤时用 MAPQ >= 20。

Q4：MEGAHIT和SPAdes怎么选？

看数据量和内存。MEGAHIT内存省、速度快，专为宏基因组设计，几十GB数据也能跑。SPAdes/metaSPAdes组装连续性更好（N50更高），但内存要求高达几百GB，适合小数据集或单基因组。我的T2D项目数据量大、样本多，选MEGAHIT是标准做法。

Q5：minimap2是什么？和BWA/Bowtie2有什么区别？

minimap2是专为三代长读（PacBio、Nanopore）设计的比对工具，能处理几千到几万bp的reads，容忍15%的错误率。BWA和Bowtie2是二代短读工具（100-300bp），要求错误率低。minimap2比BWA快几十倍处理长reads。如果未来宏基因组用Nanopore测序，就要换成minimap2。

Q6：什么是N50？怎么评价组装质量？

N50就是"把所有contig按长度从大到小排列，累加到总长度一半时那条contig的长度"。白话说就是"中位水平的contig有多长"。N50越大说明组装越连续。评价宏基因组组装还看：总长度、contig数量、最长contig。对MAGs还要用CheckM看完整度（>90%好）和污染度（<5%好）。

Q7：Bowtie2的end-to-end模式和local模式有什么区别？

end-to-end要求整条read从头到尾都比上参考序列，比较严格。local模式允许read两端被"软裁剪"（soft-clip），只要中间部分比上就算成功。去宿主时用end-to-end + very-sensitive，确保找到所有人源reads；如果做其他分析需要更宽松匹配则用local。

---

延伸阅读

1. BWA-MEM2 官方仓库：https://github.com/bwa-mem2/bwa-mem2 （v2.3, 2025年发布，速度比BWA快1.3-3.1x）
2. Bowtie2 官方文档：https://bowtie-bio.sourceforge.net/bowtie2/manual.shtml
3. HISAT2 官方：http://daehwankimlab.github.io/hisat2/
4. minimap2 论文：Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics, 34:3094-3100
5. MEGAHIT 论文：Li et al. (2015) MEGAHIT: An ultra-fast single-node solution for large and complex metagenomics assembly. Bioinformatics
6. SAM格式规范：https://samtools.github.io/hts-specs/SAMv1.pdf
7. samtools 文档：http://www.htslib.org/doc/samtools.html
8. CheckM 文档：https://github.com/Ecogenomics/CheckM

---

速记卡片

比对 = 有参考，把reads贴回去找位置
组装 = 没参考，把reads拼成长序列
BWA-MEM2 = WGS变异检测标配，快
Bowtie2 = 宏基因组去宿主首选，内存省
HISAT2 = RNA-seq专用，能跳着比（剪接比对）
minimap2 = 三代长读万能刀，极快
MEGAHIT = 宏基因组组装首选，内存极省
SAM/BAM = 比对结果标准格式，BAM是压缩版
N50 = 组装质量核心指标，越大越好