MSBS7101 Lecture 5 · NGS 序列比对（BWA 与 IGV）

一、本讲信息

课程：MSBS7101 Essential Skills in High-throughput Sequencing Data Analysis（高通量测序数据分析核心技能，与 SBMS7201 合授）

讲次：2026-09-29 · NGS Sequence Alignment

教师：Prof. Jason Wong（黄），School of Biomedical Sciences, HKU；邮箱 jwhwong@hku.hk

形式：讲授 + Google Colab 命令行实操（NGS_sequence_alignment_command_line.ipynb）

课件：NGS_alignment_2026.pdf（47 页）

二、学习目标

依据课件第 2 页的 Lecture objectives：

能够描述 FASTQ 格式（FASTQ format）的特征；

理解短读段比对（short read alignment）面临的挑战；

能够使用 BWA 将短读段比对到参考基因组；

能够在整合基因组浏览器（Integrative Genomics Viewer, IGV）中可视化已比对的短读段。

三、核心内容

1. 回顾：NGS 测序如何产生数据

下一代测序（next-generation sequencing, NGS）包含多种技术，本课程以 Illumina 测序仪为例。Illumina 采用边合成边测序（sequencing by synthesis）的策略：DNA 复制时每掺入一个带荧光标记的碱基，仪器就记录一次荧光信号，序列信息由荧光信号推断而来（课件 p3–4）。由此得到的不是完整的基因组序列，而是大量较短的测序读段（short reads）——本讲练习中的数据为 36 bp（ChIP-seq 单端）与 100 bp（WGS 双端）。

测序得到的原始短读段本身只是零散碎片，必须通过比对（alignment）确定每条读段在参考基因组上的位置，才能进行后续分析（课件 p5）。

2. FASTQ 格式：测序数据的存储形式

原始读段通常存为 FASTQ 格式（.fastq 或 .fq，常压缩为 .fq.gz）。每条读段固定占 4 行 （课件 p11）：

行	内容
1	读段标识符（read identifier），以 `@` 开头
2	碱基序列（sequence）
3	分隔行（spacer），以 `+` 开头
4	质量分数（quality score），与第 2 行等长

读段标识符的字段 （课件 p12），例如 HWI-ST216_307:3:1101:1704:1834#0/1：

HWI-ST216_307：仪器名称（instrument name）

3：流动槽泳道编号（flowcell lane）

1101：tile 上的 x 坐标

1704：tile 上的 y 坐标

1834：多路复用标签（multiplex tag，即 index/barcode）

/1 或 /2：读段配对标签（read pair label），/1 为双端中的第一条，/2 为第二条

质量分数（quality score） 采用 Phred 标度，表示每个碱基被正确判读的概率。Phred 分数 Q 与错误概率 P 的关系为 P = 10^(−Q/10)：Q = 30 对应错误率 0.1%，即 99.9% 的概率判读正确；Illumina 序列的最高分为 40（课件 p13）。

ASCII 编码 （课件 p14）：为节省空间，FASTQ 用一个 ASCII 字符存储一个质量分：分值 = 字符的 ASCII 码 − 33。例如字符 C（ASCII 67）表示 67 − 33 = 34 分；字符 I（ASCII 73）表示 73 − 33 = 40 分。

实操提示（课件 p9）：文件名中的 ChIP-seq_H3K27ac_example 里包含小写字母 l 而不是数字 1，输入命令时须注意区分。

3. 质量控制（Quality Control）

正式分析前需用 FastQC 程序检查读段质量：运行 fastqc <文件名>.fq.gz，它会生成 HTML 报告，下载后可在浏览器中查看（课件 p17–19）。FastQC 报告包含每条读段的碱基质量分布、GC 含量、接头（adapter）污染等指标，用于判断数据是否可用、是否需要修剪。

4. 短读段比对的挑战与 Burrows–Wheeler 变换

挑战（课件 p20）：如何把数百万条短序列匹配到人类基因组上？困难在于——序列短、数量巨大、且必须容许非精确匹配（测序有错误，个体基因组与参考基因组也存在真实差异）。

Burrows–Wheeler 变换（Burrows–Wheeler transform, BWT） （课件 p21）：

基因组太大，无法对每条读段做标准双序列比对（即动态规划，dynamic programming）——计算量不可行；

思路类似 BLAST：先对参考基因组建索引（index），用于快速定位种子（seeding）；

但搜索算法使用 Burrows–Wheeler 变换。BWT 本是一种压缩算法，其搜索经过优化，特别适合短序列的近似精确匹配 （near-exact matches）；

因此对于短读段，BWA（Burrows–Wheeler Aligner）比 BLAST 更快、更准确。

5. 用 BWA 进行比对

BWA 是常用的短读段比对软件，核心子命令为 bwa mem。本课程参考基因组使用人类基因组 hg38 版本；为节省练习时间，数据库只包含 2 号染色体（课件 p24）：

chr2.fa：参考基因组序列（仅 chr2）

chr2.fai：基因组序列索引

chr2.fa.*：由 bwa index 生成的 BWA 索引文件，运行比对所必需

比对主命令（默认参数，课件 p25）：

6. SAM/BAM 格式与 samtools

SAM 格式（Sequence Alignment Map format） （课件 p27）：bwa mem 的输出。文件开头是头部（header），以 @ 开头，记录所用的参考序列（@SQ）与比对命令（@PG）等信息；随后每一行是一条读段的比对结果，包括比对到的染色体、位置、质量、是否配对等。

BAM 格式 ：SAM 是纯文本，体积大，通常压缩为 BAM（二进制格式）存储；大多数下游工具直接使用 BAM 而非 SAM。samtools 是处理 SAM/BAM 最常用的工具集（课件 p28）。

常用操作流程（课件 p30–31）：

本次 ChIP-seq 数据的 flagstat 结果：共 9,151 条读段，其中 9,057 条成功比对（mapped）。基因组浏览器加载的是排序并建好索引 的 BAM 文件。

7. IGV 基因组浏览器可视化

IGV 是运行在本机（client side）的基因组浏览器，功能与 UCSC Genome Browser 类似，但更适合查看 NGS 比对结果而非基因组 track（课件 p33）。在 Colab 中通过 igv-notebook 加载 track，指定基因组版本 hg38 与位置（locus）即可。

界面要点（课件 p35–38）：

覆盖度（coverage）track ：显示每个位置上方叠放的读段数量；

比对读段（aligned reads） ：以短横条显示每条读段；

基因注释（gene annotation） ：细条为内含子（intron），粗条为外显子（exon）；进一步放大后可见外显子粗细分级——稍粗的条为非编码外显子（如 5′UTR），最粗的条为编码氨基酸的外显子；

输入坐标（如 chr2:47783034-47783438）或用缩放按钮可放大到单个碱基；放大后彩色竖条标记错配（mismatch） 位置；在 settings 中选择 "Show all bases" 可显示读段的实际碱基。

生物学示例（课件 p36）：H3K27ac 组蛋白修饰标记基因组中的活跃启动子/增强子，因此 ChIP-seq 读段在基因启动子区域富集，形成上一讲所说的「峰（peak）」。

8. 单端与双端测序，WGS 双端数据比对

单端测序（single-end sequencing） 只从 DNA 片段一端测序：

文库制备成本更低；

同样费用可测更多 DNA 片段；

适合 ChIP-seq 等应用。

双端测序（paired-end sequencing） 从同一 DNA 片段两端各测一条读段（课件 p39）：

定位更准确（两条读段的位置互相约束）；

能够发现参考序列中的断裂（即结构变异，structural variants）或 RNA-seq 中由剪接（splicing）产生的连接；

WGS 数据集即双端数据，故有 WGS_example_1.fq.gz（配对中的第一条读段）与 WGS_example_2.fq.gz（第二条读段）两个文件。

双端比对用管道（pipe）串联可省去中间文件（课件 p41）：

WGS 数据统计（课件 p42）：共 80,792 条读段，80,589 条比对成功；注意部分读段会比对到多个位置（supplementary，补充比对），因此配对读段计数可高于比对读段数。在 IGV 中加载 WGS track 后输入 MSH6 或 chr2:47781678-47809690 即可查看该基因；通过 settings 中的 "View as pairs" 成对显示读段、"Color alignments by > read strand" 按链着色，右键还可按不同选项对读段排序（课件 p45）。

9. 自主练习（Self-exercise）

Datasets/ 文件夹中还有两组数据 WXS_example_*.fq.gz 与 WXS_example_cancer_*.fq.gz，分别是一位结直肠癌（colorectal cancer）患者血液与癌组织的测序数据，下周将用于单核苷酸变异（single nucleotide variant, SNV）检测；本周的目标只是把它们比对好。题目见 Moodle 与 NGS_exercise.ipynb。

四、关键图解

课件 p4：Illumina 边合成边测序原理示意——DNA 合成过程中由荧光信号推断碱基序列，产出短测序读段。

课件 p5：比对的概念——将零散的短读段定位到参考基因组上，得到堆积在基因区域的已比对读段（aligned reads）。

课件 p11：FASTQ 格式中每条读段占 4 行——读段标识符、碱基序列、分隔行与质量分数。

课件 p14：质量分数的 ASCII 编码——字符 ASCII 码减 33 即 Phred 分值（C=34，I=40）。

课件 p27：SAM 文件结构——@ 开头的头部（含 @SQ 参考序列、@PG 比对命令）之后逐行记录每条读段的比对信息。

课件 p36：IGV 中 H3K27ac ChIP-seq 读段在基因启动子区域富集形成「峰」；基因注释中细条为内含子、粗条为外显子。

课件 p37：IGV 放大后可见单条读段，彩色竖条标记错配碱基；外显子按粗细区分非编码区（如 5′UTR）与编码区。

五、课前思考 / 自测

以下题目取自课件与 NGS_exercise.ipynb（答案参照 NGS_exercise_answers_2026.ipynb）：

为什么每个样本有两个文件（如 WXS_example_1.fq.gz 与 WXS_example_2.fq.gz）？——因为这是双端测序数据，两个文件分别对应 DNA 片段两端测出的读段。

FASTQ 文件中每条读段占几行？第 4 行质量分数 I（ASCII 73）对应的 Phred 分值与正确率是多少？——4 行；73 − 33 = 40 分，即最高分，错误率 10^(−40/10)。

ChIP-seq 文件共有 9,151 条读段，若已知总行数，如何用一条 Linux 命令统计读段数？——每条读段 4 行，zcat < 文件 | wc -l 后除以 4。

samtools flagstat 显示 WGS 数据 80,792 条读段中 80,589 条比对成功，为什么配对（paired）计数会高于比对读段数？——因为部分读段产生补充（supplementary）比对，一条读段可记录多条比对结果。

WXS 数据在 IGV 中查看时，读段主要覆盖基因的哪些区域？为什么？——主要覆盖外显子；WXS 是全外显子组测序（whole exome sequencing），只测编码区。

六、参考材料

主课件：Materials/NGS_alignment_2026.pdf（Prof. Jason Wong，2026）

课堂实操：Materials/NGS_sequence_alignment_command_line.ipynb

练习与答案：Materials/NGS_exercise.ipynb、Materials/NGS_exercise_answers_2026.ipynb