MSBS7101 Lecture 8 · RNA 测序数据分析（DESeq2 与 GSEA）

一、本讲信息

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

日期：2026 年 10 月 27 日

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

材料：讲义 RNA-seq_2026.pdf（32 页，本讲主课件）；Colab 实操笔记 RNA_Seq_Pipeline.ipynb（比对与定量）；R 实操讲义 Differential_expression.Rmd（差异表达与 GSEA）

本讲定位 ：前几周课程围绕 DNA 层面的变异检测（WGS/WXS），本讲转向 RNA 层面，完整走一遍 RNA 测序（RNA-seq）基因表达分析的标准流程。课程明确只聚焦最基本的基因表达分析，不展开 isoform、非编码 RNA、RNA 编辑等进阶应用。

二、学习目标

依据讲义第 2 页的学习目标并结合本讲内容：

能够按顺序描述 RNA-seq 基因表达分析的完整步骤：质控 → 比对 → 定量 → 差异表达 → 可视化与通路分析。

理解 RNA-seq 比对与 DNA-seq 比对的关键区别（剪接感知），以及 STAR 比对器的工作原理。

理解 featureCounts 如何借助 GTF 注释把 reads 汇总为基因水平计数矩阵，并关注其关键参数。

掌握 DESeq2 差异表达的统计思想：负二项分布、离散度估计、广义线性模型与多重检验校正。

会用 PCA、火山图（volcano plot）和基因集富集分析（GSEA）解读差异表达结果，并能用 R 与 Google Colab 完成上述分析。

三、核心内容

3.1 为什么测 RNA：应用与定位

同一个体内不同细胞（如肌肉细胞与血细胞）的 DNA 几乎相同，细胞类型与功能的差异主要由转录产生的 RNA 决定，因此测量 RNA 可以直接反映"细胞当前在做什么"。讲义第 3 页列出 RNA-seq 的主要应用：

基因表达谱（gene expression profiling） 与样本间差异表达（differential expression） 分析——本讲主题；

判断细胞类型身份（包括单细胞转录组）；

与表观遗传学、基因调控数据整合；

发现新基因，尤其是人类的非编码 RNA（non-coding RNA）——不编码蛋白质的 RNA 分子；

可变剪接（alternative splicing） 与 RNA 编辑/突变检测。

课堂强调：比较肿瘤与正常组织在 RNA 层面的差异，可以寻找肿瘤的分子机制、弱点与治疗靶点；讲师提到 RNA-seq 是 NGS 中仅次于变异检测的常用应用（课堂录音）。

3.2 从样本到计数矩阵：总体流程

RNA-seq 建库流程与 DNA-seq 基本一致，区别仅在最初需分离 RNA 并通过逆转录（reverse transcription） 把 RNA 转成更稳定的互补 DNA（cDNA）再进行测序。测序仪通过荧光信号读取碱基序列，最终产出数百万条 DNA 片段（reads）。

本讲实操的标准流程（讲义第 5 页流程图 + notebook）：

课堂数据为 6 号染色体部分区域的肿瘤-正常配对样本，已下采样（downsampled）以加快运行；比对用参考基因组索引由 hg38 FASTA 与 MANE GTF 预先构建。

3.3 质控：FastQC

FastQC 是对 FASTQ 文件做质量控制的工具，输出 HTML 报告。课件数据为 150 bp 双端（paired-end）reads——即每条 DNA 片段从两端各测一次——且已完成** 修剪（trimming）**，即去除了接头（adapter）与低质量碱基。课堂展示的报告要点（课堂录音）：

reads 长度集中在约 50–55 bp，质量分数良好；

GC 含量分布偏离理论曲线、出现重复与过表达序列，这在 RNA-seq 中并不意外：少数高表达基因会产生大量相同或相近的 reads；

由于原始等长 reads 已被修剪，报告没有 adapter 警告，课堂因此跳过修剪步骤。

3.4 比对：为什么必须用剪接感知比对器

比对（alignment） 是把每条 read 定位到参考基因组上的过程。RNA-seq 比对的核心难点在于RNA 剪接（splicing） ：基因在基因组上由外显子（exon，保留部分） 与内含子（intron，被切除部分） 组成，成熟的 mRNA 只含外显子；逆转录得到的 cDNA 因此是"外显子相连"的序列。若用 BWA 等常规 DNA 比对器把 RNA reads 比到基因组，被删去的内含子会被误判为大片段 gap，比对效果差（讲义第 8 页；课堂录音）。

两种分析策略（讲义第 8–10 页）：

比对到参考基因组 ：需要 STAR、HISAT2、TopHat2 等剪接感知（splice-aware） 比对器，能够跨过剪接位点比对；适用于大多数场景，包括发现新转录本。

比对到转录组（cDNA） ：cDNA 已无内含子，可用常规比对器甚至更快的直接计数算法；若只关心已知基因的表达量，此法更快。

比对器选择还涉及运行时间与内存的权衡、以及精度需求：例如检测新型变异或 RNA 编辑需要比差异表达更高的比对精度。

3.5 STAR 比对器的工作原理

STAR 是本讲使用的剪接感知比对器（讲义第 11 页）：

对参考基因组直接构建后缀数组（suffix array） 索引并常驻内存——速度快但内存消耗大；

对每条 read 寻找基因组中最长的精确匹配片段（种子 seed ），而非逐碱基动态规划；

把位置相近的种子聚类并"拼接（stitch）"起来；

若拼接中出现数百至数千碱基的大缺口，STAR 将其解释为剪接位点（splice junction） ，即外显子–内含子边界；

按打分方案选择最佳比对：匹配加分，错配、gap、软修剪（soft clipping）扣分；剪接位点需有经典剪接信号 GT–AG 支持才被允许。

课堂日志显示，约 1 万条 reads 中约六七千条跨剪接位点，符合"外显子较短、reads 常横跨外显子连接处"的预期；绝大多数剪接位点为 GT–AG 经典信号（课堂录音）。注意 STAR 需要专用索引，不能复用 BWA 索引；实际分析中应检查众多参数的默认值是否适合自己的数据（课堂录音）。

输出（讲义第 12–13 页）：主输出为按基因组坐标排序的 BAM 文件（2-0019A_Aligned.sortedByCoord.out.bam）——BAM 是比对结果的二进制压缩格式，供下游计数与可视化使用；其余为日志与 reads 分布统计文件。

3.6 定量：featureCounts 与 GTF

featureCounts （subread 软件包内）把 BAM 中的 reads 与 GTF 注释中的外显子区域取交集，再把同一基因各外显子的计数按 gene_id 汇总，得到基因水平的** 计数矩阵（count matrix）**——即"每个基因在每个样本中被数到多少条 reads"的整数表格（讲义第 14 页）。

GTF 是基因注释文件，充当"基因在基因组上的地图"（本讲用 MANE 注释），包含 seqname、source、feature 类型、起止坐标、链方向等列；定量时只用到其中的 exon 记录。关键参数（讲义第 14 页；课堂录音）：

-p：双端数据按片段（fragment，即一对 reads） 计数，避免同一对 reads 被重复计算；

-s 0：链特异性设置；本数据为 0 表示文库非链特异性，正反链 reads 都计数；

-g gene_id：按基因聚合外显子计数。

featureCounts 能处理外显子重叠、reads 跨剪接位点、同时命中多个基因等复杂情况，并可设置此类歧义 reads 的分配或忽略规则，这是它比简单的区间取交工具更适合定量的原因（课堂录音）。

输出（讲义第 16 页）：主文件 featureCounts.txt 包含 Geneid、染色体、起止、链、长度及每个 BAM 一列整数计数，可直接作为 DESeq2/edgeR/limma 的输入；*.summary 文件给出质控统计：Assigned（成功分配到基因）、Unassigned_Unmapped（未比对上）、Unassigned_NoFeatures（落在无注释区域）、Unassigned_Ambiguity（歧义）等。课堂 summary 显示相当一部分 reads 未被分配，真正用于定量的只是一部分（课堂录音）。

3.7 可视化验证：IGV

IGV（Integrative Genomics Viewer） 用于在基因组浏览器中检查 BAM。课堂观察（课堂录音）：覆盖度峰主要集中在外显子区域；跨剪接位点的 reads 显示为从上一个外显子末端直接跳到下一个外显子；按链着色时正链为红色、反链为蓝色。这一步用于人工核对比对与定量是否合理。

3.8 差异表达分析：DESeq2

流程从 Google Colab 切换到 R/RStudio（讲义第 18–19 页）。DESeq2 是最常用的差异表达工具之一（另一为 edgeR），集标准化、方差稳定与统计检验于一体。

为什么不能用 t 检验 （讲义第 21 页；课堂录音）：t 检验假设数据服从正态分布，而 RNA-seq 计数含有大量零值、高表达端迅速下降，服从负二项分布（negative binomial distribution） 。DESeq2 的做法：

为每个基因估计离散度（dispersion）——衡量该基因计数方差超过泊松预期的程度，不同基因的生物学变异程度不同，需逐基因估计；

对每个基因拟合广义线性模型（GLM） 并比较条件（如 Tumour vs Normal）；设计公式可纳入配对信息与年龄、性别等协变量，而不仅限于两组比较；

输出 log2 倍数变化（log2 fold change） 与经多重检验校正的 FDR 校正 P 值（adjusted p-value）——同时检验约两万个基因时必须校正，否则仅凭偶然也会出现"显著"基因。

文库标准化（library normalization） （讲义第 22 页）：使计数在不同样本间可比，需解决两个问题：测序深度不同（文库大小）与文库组成不同（如不同组织类型中高表达基因的构成差异）。DESeq2 内置按文库大小的标准化并在过程中抑制离群基因的影响。

实操数据 （讲义第 23–24 页；Differential_expression.Rmd）：乳腺癌数据集，6 对肿瘤/配对正常样本共 12 列；原始 33 646 个基因去除 Ensembl ID 的非编码基因后剩 25 281 个。步骤：构建样本信息表（Conditions）→ DESeqDataSetFromMatrix 构建 DESeq2 数据集（设计公式 ~ Type + Sample，其中 Sample 项实现配对检验）→ vst() 做方差稳定化变换（variance stabilising transformation） ：原始计数中高表达基因方差远大于低表达基因，vst 使方差近似与均值无关（同方差性），供 PCA、聚类热图使用 → DESeq() 拟合模型并做 Wald 检验 → results() 提取结果并按显著性筛选排序。

3.9 结果解读：PCA、火山图与 GSEA

讲义第 25 页指出：仅有一份差异基因清单难以解释生物学意义，需要三类互补分析。

PCA（主成分分析，Principal Component Analysis） （讲义第 26–27 页）：把每个样本视为约 2.5 万维空间（每维一个基因）中的一个点，PCA 将其压缩为少数主成分（PC）——PC1 捕捉最大方差来源，PC2 次之。作用是在做差异检验之前检查样本是否按预期分组、标签是否正确、有无离群样本或隐藏混杂因素。课堂结果：PC1 解释 42% 方差、PC2 解释 11%，肿瘤与正常样本沿 PC1 清晰分离（仅一个正常样本例外），正常样本组内部异质性更大——分组分离良好预示会有较多显著差异基因（课堂录音）。

火山图（volcano plot） （讲义第 28 页）：横轴为 log2 fold change（效应大小与方向），纵轴为统计显著性，每点一个基因，可快速识别既统计显著又变化幅度大的基因；观察点的整体分布形态（聚集在中下部或向两侧升起）有助于判断效应与置信度的平衡。

基因集富集分析（GSEA, Gene Set Enrichment Analysis） （讲义第 29–31 页）：上千个差异基因难以逐一解释，GSEA 回答"某些通路或功能是否整体上调或下调"。做法：

用 DESeq2 结果给全部基因排序（按 Wald 统计量——带方向的检验统计量，正值越大表示在肿瘤中越上调；无生物学重复时可直接用 fold change 排序，见 Rmd）；

选择预定义基因集数据库（Hallmark、KEGG、MSigDB 等；Hallmark 约含 50 组通路，每组最多约 200 个基因）；

检验各基因集的成员是否集中在排序列表顶端（上调）或底端（下调），保留显著通路。

结果解读三要素：NES（Normalized Enrichment Score，归一化富集分数）——方向与强度，正值为顶部富集（上调），绝对值越大信号越强；** 校正 P 值（FDR q 值）**——主要阈值（如 < 0.05）；leading edge——驱动该通路富集信号的核心基因子集。课堂结果：E2F 靶基因、G2M 检查点等通路在肿瘤中正向富集；血管生成（angiogenesis）通路整体不显著——即使个别成员基因变化很大，若基因集成员分散在整个排序列表中，整体仍判定为不显著，说明要看整组基因的分布而非单个基因（课堂录音）。

差异基因数量 ：讲义第 25 页为 padj < 0.05 时 2 813 个显著差异基因；课堂演示用更严格阈值（校正 P < 0.01）得到约 1 297 个（课堂录音）；Rmd 讲义示例则用 padj < 0.001 且 |log2FC| > 2 筛选。数量随阈值变化，属正常现象。

3.10 作业二（Assignment 2）

截止 2026 年 11 月 13 日（周五）23:59（讲义第 32 页）：从癌症样本的全外显子组测序（WXS）数据中识别致癌突变；在 Google Colab 上用模板 .ipynb 完成，需用到此前课堂生成的 WXS_germline.vcf 与 WXS_example_cancer_sorted.bam（缺失时运行初始代码块会自动下载）。课堂补充（录音）：在 Moodle 下载作业 notebook，填入姓名与 UId，补全代码并写入结果；第 7 题需在 UCSC Genome Browser 生成并保存 session 文件；最终将 .ipynb 与 session 文件一并提交到 Moodle。

四、关键图解

讲义第 5 页：RNA-seq 分析总流程图，从原始数据经质控、比对、定量到差异表达与下游分析的完整环节，是本讲的结构主线。

讲义第 8 页：cDNA（转录组）与基因组 DNA 比对的对比。mRNA 已切除内含子，reads 比对到基因组时必须跨越剪接位点，说明 RNA-seq 不能直接使用常规 DNA 比对器。

讲义第 11 页：STAR 比对算法——先找最长精确匹配种子，聚类拼接，大缺口被解释为剪接位点，并按 GT–AG 经典信号打分选择最佳比对。

讲义第 14 页：featureCounts 将 BAM 中 reads 与 GTF 外显子区域取交集并按 gene_id 汇总为基因水平计数矩阵的过程。

讲义第 21 页：DESeq2 以负二项分布估计离散度、用广义线性模型比较条件，产出 log2 fold change 与 FDR 校正 P 值。

讲义第 26 页：PCA 把数万维基因表达压缩为少数主成分，用于检查样本分组、重复一致性与离群样本。

讲义第 28 页：火山图横轴为 fold change、纵轴为显著性，可同时评估效应大小与置信度，快速定位关键基因。

讲义第 29 页：GSEA 将基因按差异统计量排序，检验预定义基因集是否集中在列表一端，从而判断通路的整体上调或下调。

五、课前思考 / 自测

RNA-seq 的 reads 来自已剪接的 mRNA（逆转录为 cDNA）。若坚持用 BWA 等常规 DNA 比对器把 reads 比到参考基因组，比对结果会出现什么问题？为什么剪接感知比对器能避免？

featureCounts 的 -p（按片段计数）与 -s 0（非链特异性）两个参数分别解决什么问题？如果双端数据误按单条 read 计数，计数矩阵会怎样被影响？

为什么 RNA-seq 计数不适合用 t 检验做差异分析？请从分布形态（负二项分布）、逐基因离散度、多重检验校正三个角度说明 DESeq2 的对应处理。

做差异检验之前先画 PCA 图的目的是什么？如果 PCA 中肿瘤与正常样本严重重叠，你对接下来的差异分析结果应有什么预期？

一个基因集中有个别基因变化幅度非常大，但 GSEA 判定该基因集不显著。这说明 GSEA 的判断依据是什么？

六、参考材料

主课件：Materials/RNA-seq_2026.pdf（RNA-seq, SBMS7201/MSBS7101, Prof. Jason Wong, 32 页）

实操材料：Materials/RNA_Seq_Pipeline.ipynb（Colab：FastQC、STAR、featureCounts、IGV）；Materials/Differential_expression.Rmd（R：DESeq2、PCA、火山图、fgsea/clusterProfiler）

课堂录音：Recordings/MSBS7101-1027-summary.md、Recordings/MSBS7101-1027-transcript.md

讲义中引用的公开来源：StatQuest RNA-seq 视频（p4、p22）；Conesa et al., Genome Biology 2016（p9）；HBC Training 教程（p11、p14、p21）；BioRender 火山图解读模板（p28）

七、课堂要点补充（据课堂录音整理）

以下为讲义未展开、老师课堂口头强调的内容：

应用背景 ：人类基因组约 2 万个蛋白编码基因，同时存在大量非编码 RNA；此前作业关注 DNA 变异，本讲关注 RNA 差异，两者互补。

工具生态 ：剪接感知比对器除 STAR 外还有 TopHat2、HISAT2；RNA-seq 文件通常比 DNA-seq 小，整体分析更快。

参数意识 ：STAR 与 featureCounts 参数众多，实际分析务必检查默认值是否适合自己的数据；课堂数据已预先 trimmed，故跳过修剪步骤——实际研究中是否修剪需视数据而定。

定量细节 ：GTF 中虽有 CDS、起始密码子、UTR 等注释，定量只用 exon；基因长度信息在后续标准化中会用到。

统计直觉 ：DESeq2 的 GLM 设计公式可纳入配对样本（~ Type + Sample）及年龄、性别等协变量；这是它比两组直接比较更灵活之处。

GSEA 案例解释 ：E2F 靶基因与 G2M 检查点通路上调与肿瘤增殖活跃一致；angiogenesis 通路虽与肿瘤相关，但课堂数据中其基因集整体不显著，再次说明要看整组基因在排序列表中的分布。

课程安排 ：本讲为讲师本学期最后一讲，随后四周为 BGI 讲座；课程录像将上传并把链接放到课程平台。