转录组分析流程（RNA-seq Analysis Pipeline）¶
一句话说明¶
转录组分析（RNA-seq）是把细胞在某个时刻"正在表达的基因"全部测出来，通过 质控 → 比对 → 定量 → 差异分析 → 功能富集 的流程，找到不同条件下哪些基因"开得更多"或"关得更多"，从而揭示疾病机制或生物学过程。你项目中用 featureCounts 做定量、DESeq2 做差异分析，就是这个流程的核心步骤。
目录¶
1. 核心概念（白话版）¶
1.1 什么是转录组（Transcriptome）¶
- 定义：转录组就是一个细胞/组织在某个时间点，所有被"转录"出来的 RNA 的总和。转录（Transcription）就是 DNA 上的基因被"抄写"成 mRNA 的过程。
- 白话比方：如果基因组（Genome）是一本完整的"菜谱大全"（包含所有菜的做法），那转录组就是"今天厨房实际在做的菜的清单"——虽然菜谱里有 2 万多道菜，但今天只做了其中几千道，每道菜做的份数还不一样。
- 为什么要研究转录组：
- 基因组是"静态"的（每个细胞的 DNA 序列基本一样）
- 转录组是"动态"的（不同组织、不同时间、健康 vs 疾病，表达的基因完全不同）
- 通过比较疾病组 vs 健康组的转录组差异，就能找到和疾病相关的关键基因
1.2 RNA-seq 实验流程概述¶
RNA-seq（RNA Sequencing，RNA 测序）是研究转录组最主流的技术：
| 步骤 | 做什么 | 白话解释 |
|---|---|---|
| 提取 RNA | 从细胞/组织中提取总 RNA | 把"菜单"从厨房拿出来 |
| mRNA 富集 | 两种方式二选一：① polyA 选择（用 oligo-dT 磁珠抓取有 polyA 尾巴的 mRNA，最常用）；② rRNA 去除（用试剂盒去掉占 80-90% 的核糖体 RNA）。rRNA：细胞里造蛋白质的机器上自带的 RNA，数量巨大但不是研究目标 | 只要 mRNA，不要其他杂 RNA |
| 逆转录 | 把 RNA 反转录成 cDNA（互补 DNA；cDNA：用酶把 RNA 反向抄回 DNA，因为测序仪只认 DNA） | 测序仪只能读 DNA，所以要把 RNA "翻译"回 DNA |
| 建库 | 打断 cDNA、加接头（adapter：像给每个片段贴上条形码和把手，测序仪靠它识别和固定片段） | 给 DNA 片段两头加上"条形码"，让测序仪识别 |
| 上机测序 | Illumina 测序仪进行边合成边测序 | 一个碱基一个碱基地读出序列 |
| 产出数据 | 得到 FASTQ 格式的原始数据 | 每条 read 约 150bp，几千万到上亿条 |
1.3 转录组 vs 基因组 vs 宏基因组的区别¶
| 维度 | 基因组（Genomics） | 转录组（Transcriptomics） | 宏基因组（Metagenomics） |
|---|---|---|---|
| 研究对象 | 一个物种的全部 DNA | 一个样本正在表达的 RNA | 一个环境中所有微生物的 DNA |
| 白话解释 | 整本菜谱 | 今天做了哪些菜、做了多少 | 整个美食城所有餐厅的菜谱合集 |
| 测什么 | DNA（静态） | mRNA（动态，反映实时状态） | 混合 DNA（多物种） |
| 关注点 | 突变、结构变异 | 基因表达量变化 | 物种组成 + 群落功能 |
| 比对参考 | 同一物种参考基因组 | 同一物种参考基因组 + 注释 | 微生物数据库 |
| 你项目中 | — | 用 featureCounts + DESeq2 | Kraken2 + HUMAnN3 |
面试陷阱：宏基因组测的是 DNA 不是 RNA，它反映的是微生物"有什么基因"，而不是"正在表达什么基因"。如果要研究微生物实际在干什么（正在表达的功能），需要做"宏转录组"（Metatranscriptomics）。
2. 转录组分析完整流程¶
流程图（文字版）¶
原始数据 (FASTQ)
│
▼
┌──────────────────┐
│ Step 1: 质控 │ FastQC 评估 + fastp 清洗
│ (Quality Control)│ 去接头、去低质量reads
└──────────────────┘
│
▼
┌──────────────────┐
│ Step 2: 比对 │ HISAT2 或 STAR
│ (Alignment) │ 把reads映射到参考基因组
└──────────────────┘
│
▼
┌──────────────────┐
│ Step 3: 定量 │ featureCounts / HTSeq / Salmon
│ (Quantification) │ 统计每个基因被多少reads覆盖
└──────────────────┘
│
▼
┌──────────────────┐
│ Step 4: 差异分析 │ DESeq2 / edgeR
│ (DE Analysis) │ 找出两组间表达差异显著的基因
└──────────────────┘
│
▼
┌──────────────────┐
│ Step 5: 功能富集 │ clusterProfiler (GO/KEGG)
│ (Enrichment) │ 差异基因参与了什么生物学功能
└──────────────────┘
Step 1: 质控（Quality Control）¶
做什么：检查原始测序数据的质量，去掉低质量的 reads 和接头序列（adapter）。
白话解释：原始数据就像刚拍完的照片，里面可能有模糊的、有水印的，质控就是"筛选照片"，只留下清晰的。
工具： - FastQC：只负责"看"质量，生成报告，不修改数据 - fastp：既能看又能"修"，自动检测并去除接头、过滤低质量 reads
关键指标：
| 指标 | 含义 | 合格标准 |
|---|---|---|
| Q20 | 碱基错误率 ≤1% | >95% 的碱基达到 Q20 |
| Q30 | 碱基错误率 ≤0.1% | >85% 的碱基达到 Q30 |
| GC 含量 | 鸟嘌呤+胞嘧啶的比例（GC含量：DNA 四个字母 ATCG 中，G 和 C 的占比） | 接近物种理论值（人类 ~42%） |
| 接头残留 | adapter 序列的比例 | 应接近 0% |
Step 2: 比对（Alignment / Mapping）¶
做什么：把清洗后的 reads "放回"参考基因组上，找到每条 read 来自基因组的哪个位置。
白话解释：你手里有几千万张拼图碎片（reads），参考基因组就是拼图的"原图"。比对就是把碎片一片一片放回原图对应的位置。
为什么 RNA-seq 比对特殊：mRNA 经过了"剪接"（Splicing，就是把内含子 intron 切掉、只保留外显子 exon），所以一条 read 可能跨越两个外显子，中间隔了一段很长的内含子。普通 DNA 比对工具（如 BWA）处理不了这种"跳跃式比对"，需要专门的 splicing-aware aligner（剪接感知比对器；splicing-aware：知道 RNA 会跳着比，能处理 read 横跨两个外显子中间隔了内含子的情况），比如 HISAT2 和 STAR。
主流比对工具：
| 特性 | HISAT2 | STAR |
|---|---|---|
| 全称 | Hierarchical Indexing for Spliced Alignment of Transcripts 2 | Spliced Transcripts Alignment to a Reference |
| 内存需求 | ~8 GB（人类基因组） | ~30-40 GB（人类基因组） |
| 速度 | 快 | 非常快（但需要大内存加载索引） |
| 适合场景 | 内存有限的服务器 | 内存充足、追求速度 |
| 二次比对 | 两步比对（2-pass 需手动） | 内置 2-pass 模式，自动发现新剪接位点 |
| 输出格式 | SAM/BAM | SAM/BAM |
面试答题要点：如果面试官问"所用的什么比对器"，说 HISAT2 就行，然后补充"HISAT2 内存占用小，适合我们实验室服务器配置；如果追求更高灵敏度可以用 STAR 的 2-pass 模式"。
Step 3: 定量（Quantification）—— 重点¶
做什么：统计每个基因被多少条 reads 覆盖（count），这个 count 数就是基因的表达量的"原始度量"。
白话解释：比对完之后，每条 read 都有了"地址"。定量就是数一数每个基因的"地址"上有多少条 read 落在那里——落得越多，说明这个基因表达越活跃。
三种主流定量方式对比：
| 工具 | 定量策略 | 速度 | 需要比对吗 | 特点 |
|---|---|---|---|---|
| featureCounts | 基于比对的 read 计数 | 极快 | 需要（输入 BAM） | 简单直接，你项目在用 |
| HTSeq-count | 基于比对的 read 计数 | 慢 | 需要（输入 BAM） | Python 实现，老牌工具 |
| Salmon | 轻量级比对（selective alignment，v1.1.0 起为默认；v0.9.x 及更早版本默认为 quasi-mapping） | 极快 | 不需要（直接输入 FASTQ） | 跳过比对步骤，直接估算 transcript 表达量 |
featureCounts 详解（你项目用的）¶
featureCounts 是 Subread 软件包中的定量工具，由 Liao et al. (2014) 发表在 Bioinformatics 上。
工作原理： 1. 读入 BAM 文件（比对结果）和 GTF 注释文件（记录每个基因的外显子坐标） 2. 对每条已比对的 read，检查它落在哪个基因的外显子区域内 3. 计数：每个基因获得的 read 数就是该基因的 raw count
关键概念： - Feature（特征）：GTF 文件中的每一行，通常是一个外显子（exon） - Meta-feature（元特征）：多个 feature 的集合，通常是一个基因（gene）。同一个基因的所有外显子的 reads 加起来，就是该基因的 count - 多映射 reads（multi-mapping reads）：一条 read 比对到基因组多个位置。默认不计数（避免重复计算），可用 -M 参数开启计数 - 跨基因 reads（ambiguous reads）：一条 read 同时覆盖两个基因。默认不计数，可用 -O 参数允许
表达量标准化指标：RPKM / FPKM / TPM¶
Raw counts 受两个因素影响：测序深度（样本间 reads 总量不同）和基因长度（长基因天然能捕获更多 reads）。以下三种指标都同时校正了这两个偏差：
| 指标 | 全称 | 适用测序类型 | 公式 |
|---|---|---|---|
| RPKM | Reads Per Kilobase per Million mapped reads | 单端测序（SE） | RPKM = reads / (总reads数/10^6 x 基因长度/10^3) |
| FPKM | Fragments Per Kilobase per Million mapped fragments | 双端测序（PE） | 与 RPKM 相同，但以 fragment（片段）而非 read 为单位，避免一个片段的两条 reads 被重复计算 |
| TPM | Transcripts Per Million | SE 和 PE 均可 | 先除基因长度（得到 RPK），再除以所有 RPK 之和/10^6 |
TPM vs RPKM/FPKM 的关键区别： - RPKM/FPKM 先校正测序深度，再校正基因长度 - TPM 先校正基因长度，再校正测序深度 - 结果：TPM 在所有样本中的总和相同（都等于 10^6），因此可以直接跨样本比较比例 - RPKM/FPKM 不同样本的总和可能不同，直接比较会有偏差
当前推荐： - 跨样本比较表达量：优先用 TPM（Salmon 默认输出 TPM） - 差异表达分析：使用 raw counts + DESeq2/edgeR 自带的标准化方法（median of ratios / TMM），不要用 TPM/FPKM 做差异分析
Step 4: 差异表达分析（Differential Expression Analysis）—— 重点¶
做什么：比较两组样本（如疾病组 vs 健康组）的基因表达量，找出表达量有"统计学显著差异"的基因。
白话解释：你有"糖尿病组"和"健康组"各 3 个样本的基因 counts 数据。某个基因在糖尿病组平均 1000 reads，在健康组平均 100 reads，看起来差异很大。但这到底是真的差异，还是仅仅因为个体间的随机波动？差异分析就是用统计方法来判断这个差异是否"真实可信"。
DESeq2 详解（你项目用的）¶
DESeq2 是最主流的 RNA-seq 差异分析 R 包，由 Love, Huber & Anders (2014) 发表在 Genome Biology 上。
DESeq2 做了什么（3 步核心）：
| 步骤 | 做什么 | 白话解释 |
|---|---|---|
| 1. 标准化（Normalization） | 计算 size factors，消除文库大小差异 | 张三测了 1 亿条 reads，李四只测了 5000 万条。同一个基因在张三那里 count 多不代表表达高，只是因为张三总量多。标准化就是把"总量"这个因素去掉 |
| 2. 离散度估计（Dispersion Estimation） | 估计每个基因的变异程度 | 有些基因天生"波动大"（比如免疫基因），有些"很稳定"（比如管家基因）。DESeq2 用负二项分布（Negative Binomial Distribution；负二项分布：比泊松分布多了一个允许数据波动更大的参数，更适合 RNA-seq 天然波动大的 count 数据）来建模这种变异 |
| 3. 统计检验（Wald Test；Wald Test：比较两组均值差异是否足够大的统计方法） | 对每个基因做假设检验 | 零假设："这个基因在两组间没有差异"。如果 p 值很小（<0.05），就拒绝零假设，认为差异显著 |
DESeq2 关键输出：
| 列名 | 含义 | 面试怎么说 |
|---|---|---|
| baseMean | 所有样本的平均标准化表达量 | 这个基因总体表达水平有多高 |
| log2FoldChange（LFC） | log2(组A表达量 / 组B表达量) | LFC=2 表示 A 组是 B 组的 4 倍；LFC=-1 表示 A 组是 B 组的 0.5 倍 |
| lfcSE | LFC 的标准误 | LFC 估计的不确定性 |
| stat | Wald 检验统计量 | — |
| pvalue | 原始 p 值 | 不能直接用，因为做了上万次检验，会有大量假阳性 |
| padj | 校正后 p 值（BH 方法） | 这才是我们用来筛选差异基因的，通常要求 padj < 0.05 |
为什么要用 padj 而不是 pvalue：
假设你有 20000 个基因，每个都做一次检验（p < 0.05）。即使所有基因实际上没有差异，纯靠运气也会有 20000 × 0.05 = 1000 个基因"假阳性"通过检验。所以必须做多重检验校正（Multiple Testing Correction），DESeq2 默认使用 Benjamini-Hochberg (BH) 方法，把 pvalue 校正成 padj（adjusted p-value）。
差异基因筛选标准（常用）： - padj < 0.05（统计显著） - |log2FoldChange| > 1（表达差异至少 2 倍）
Step 5: 功能富集分析（Functional Enrichment Analysis）¶
做什么：拿到差异基因列表后，分析这些基因集中参与了哪些生物学功能或通路。
白话解释：你找到了 500 个差异基因。这 500 个基因不是随机的——它们可能大量集中在"炎症反应"或"糖代谢"通路上。富集分析就是找到这些"扎堆"的模式。
两大数据库：
| 数据库 | 全称 | 白话解释 | 内容 |
|---|---|---|---|
| GO | Gene Ontology | 基因的"职能说明书" | 分 3 大类：BP（生物学过程）、MF（分子功能）、CC（细胞组分） |
| KEGG | Kyoto Encyclopedia of Genes and Genomes | 基因的"代谢地图" | 提供代谢通路图，看基因在通路中的位置 |
R 包推荐：clusterProfiler（余光创团队开发，中国人做的，Google Scholar 引用超 28,000 次）
3. 关键工具对比¶
3.1 STAR vs HISAT2（比对器对比）¶
| 对比维度 | STAR | HISAT2 |
|---|---|---|
| 发表时间 | 2013（Dobin et al.） | 2015（Kim et al.） |
| 索引算法 | 后缀数组（Suffix Array） | HGFM（层次化图 FM 索引） |
| 人类基因组索引大小 | ~30 GB | ~8 GB |
| 运行内存 | 30-40 GB | 8-10 GB |
| 速度 | 极快（~25 min/样本） | 快（~30-40 min/样本） |
| 新剪接位点发现 | 内置 2-pass 模式（推荐） | 可手动 2-pass |
| 适用场景 | 大内存集群、ENCODE/TCGA 大项目 | 个人服务器、内存有限环境 |
| ENCODE 项目用的 | STAR（官方推荐） | — |
| 准确率 | 两者在大多数基准测试中表现相当，STAR 在检测新剪接位点上略优 | |
| 你该怎么选 | 服务器内存 ≥32 GB 优先选 STAR | 内存不够就选 HISAT2 |
2024-2025 最新趋势：STAR 在大型项目（如 GTEx、TCGA、ENCODE）中是事实标准。但 HISAT2 因为内存友好，在教学和小实验室仍然广泛使用。两者结果高度一致，面试时说"所用的是 HISAT2，了解 STAR"就够了。
3.2 featureCounts vs HTSeq vs Salmon（定量方式对比）¶
| 对比维度 | featureCounts | HTSeq-count | Salmon |
|---|---|---|---|
| 语言/包 | C（Subread 包） | Python（HTSeq 包） | C++（独立工具） |
| 输入 | BAM + GTF | BAM + GTF | FASTQ（无需比对！） |
| 速度 | 极快（几分钟处理上亿 reads） | 慢（可能要几小时） | 极快（几分钟） |
| 定量单位 | gene-level counts | gene-level counts | transcript-level TPM/counts |
| 多线程 | 支持（-T 参数） | 不支持 | 支持 |
| 使用难度 | 简单，一行命令 | 简单但慢 | 需理解 selective alignment 概念（v1.1.0 起为默认，比旧版 quasi-mapping 更准确，最新版为 v1.11.0+） |
| 发表年份 | 2014 | 2015 | 2017 |
| 与 DESeq2 配合 | 直接输入 count 矩阵 | 直接输入 count 矩阵 | 需通过 tximport 导入 |
| 你项目用的 | featureCounts | — | — |
2024-2025 最新趋势：DESeq2 官方推荐的上游流程是 Salmon + tximport（无需比对，可校正 transcript 长度偏差）。但 featureCounts 因为简单直接，仍是最广泛使用的基因水平定量工具。
3.3 DESeq2 vs edgeR（差异分析对比）¶
| 对比维度 | DESeq2 | edgeR |
|---|---|---|
| 开发团队 | Michael Love, Simon Anders, Wolfgang Huber | Mark Robinson, Davis McCarthy, Gordon Smyth |
| 统计模型 | 负二项分布 + 收缩估计（shrinkage；shrinkage：把不太可靠的极端值往平均值方向拉一拉） | 负二项分布 + 经验贝叶斯（经验贝叶斯：先从全部基因学一个规律，再用这个规律帮助估计每个基因的波动） |
| 标准化方法 | Median of ratios（中位数比率法） | TMM（Trimmed Mean of M-values；TMM：先去掉极端值，再算修正系数消除测序深度差异） |
| 离散度估计 | 数据驱动的先验分布 | 经验贝叶斯收缩 |
| 假阳性控制 | 保守（假阳性少，但可能漏掉一些真差异基因） | 相对宽松 |
| 小样本表现 | 好（3 vs 3 都能用） | 好 |
| 大样本表现 | 好 | 好（大样本推荐 edgeR 的 QLF 检验） |
| 学习曲线 | 简单（函数少，一个 DESeq() 搞定核心分析） | 稍复杂（需手动调用多个函数） |
| LFC 收缩 | 内置 lfcShrink()，可视化更友好 | 需额外处理 |
| 引用次数 | ~86,000+（截至 2026 年初） | ~35000+ |
| 你项目用的 | DESeq2 | — |
实际结论：两个工具结果高度重合（差异基因 ~80% 重叠）。选哪个取决于实验室习惯。面试时说"我用 DESeq2，了解 edgeR，两者都基于负二项分布，结果高度一致"就很好。
4. 实操代码¶
4.1 HISAT2 比对命令¶
# ============================================
# HISAT2 比对流程
# ============================================
# --- 第一步：构建参考基因组索引 ---
# -p 8: 使用 8 个 CPU 核心并行处理
# genome.fa: 参考基因组 FASTA 文件（比如人类 hg38）
# genome_index: 输出的索引文件前缀
hisat2-build -p 8 genome.fa genome_index
# --- 第二步：比对 ---
# -x genome_index: 指定索引文件前缀
# -1 sample_R1.fq.gz: 双端测序的正向 reads（Read 1）
# -2 sample_R2.fq.gz: 双端测序的反向 reads（Read 2）
# --dta: 专门为下游 StringTie 或 featureCounts 优化输出
# （会在 BAM 中添加 XS 标签标记剪接方向）
# -p 8: 使用 8 个线程
# --new-summary: 输出新格式的比对统计摘要
# --summary-file: 统计摘要保存到这个文件
# | samtools sort: 管道符，直接把比对结果传给 samtools 排序
# -@ 4: samtools 使用 4 个线程
# -o sample.sorted.bam: 输出排序后的 BAM 文件
hisat2 -x genome_index \
-1 sample_R1.fq.gz \
-2 sample_R2.fq.gz \
--dta \
-p 8 \
--new-summary \
--summary-file sample_hisat2.summary \
| samtools sort -@ 4 -o sample.sorted.bam
# --- 第三步：建立 BAM 索引 ---
# 排序后的 BAM 文件必须建索引才能被后续工具快速访问
# 生成 sample.sorted.bam.bai 索引文件
samtools index sample.sorted.bam
# --- 第四步：查看比对统计 ---
# flagstat 会显示总 reads 数、比对率、配对率等
samtools flagstat sample.sorted.bam
4.2 featureCounts 定量命令（重点）¶
# ============================================
# featureCounts 基因表达定量
# ============================================
# 单样本定量：
# -T 8: 使用 8 个线程加速
# -p: 输入是双端测序数据（paired-end）
# --countReadPairs: 按"片段"（fragment）计数，而不是按"read"计数
# （双端测序一个片段有 2 条 reads，应该只算 1 次）
# -t exon: 只统计落在外显子（exon）区域的 reads
# （因为 mRNA 只包含外显子，内含子已经被剪掉了）
# -g gene_id: 按 gene_id 汇总（把同一个基因的所有外显子 counts 加起来）
# -a annotation.gtf: GTF 格式的基因注释文件
# （告诉 featureCounts 每个基因的外显子在基因组哪个位置）
# -o counts.txt: 输出的 count 矩阵文件
# sample.sorted.bam: 输入的比对结果文件
#
# 注意：如果文库是链特异性的（strand-specific），需要加 -s 参数：
# -s 1: 正链文库（read 方向与转录方向一致）
# -s 2: 反链文库（dUTP 文库常用，read 方向与转录方向相反，目前最常见）
# -s 0: 非链特异性（默认值，不区分正反链）
featureCounts -T 8 \
-p --countReadPairs \
-t exon \
-g gene_id \
-a annotation.gtf \
-o counts.txt \
sample.sorted.bam
# 多样本一起定量（推荐，保证 count 矩阵一致性）：
# 把所有样本的 BAM 文件一起传入，输出一个统一的 count 矩阵
featureCounts -T 8 \
-p --countReadPairs \
-t exon \
-g gene_id \
-a annotation.gtf \
-o all_samples_counts.txt \
T2D_1.sorted.bam T2D_2.sorted.bam T2D_3.sorted.bam \
Control_1.sorted.bam Control_2.sorted.bam Control_3.sorted.bam
# --- featureCounts 输出文件说明 ---
# counts.txt: 主文件，包含 Geneid、Chr、Start、End、Strand、Length、count 列
# counts.txt.summary: 统计摘要，显示 Assigned（成功分配）和各种未分配原因的 reads 数
# - Assigned: 成功分配到基因的 reads 数（通常应 >60-70%）
# - Unassigned_NoFeatures: 落在基因间区域（intergenic）
# - Unassigned_Ambiguity: 同时覆盖多个基因
4.3 DESeq2 差异分析 R 代码（重点）¶
# ============================================
# DESeq2 差异表达分析完整流程
# ============================================
# --- 加载 R 包 ---
library(DESeq2) # 差异分析核心包
library(ggplot2) # 画图
library(pheatmap) # 热图
library(EnhancedVolcano) # 火山图（可选，更好看）
# ============================================
# 第一步：读入 featureCounts 的结果
# ============================================
# featureCounts 输出的 counts.txt 前 6 列是基因信息，第 7 列开始才是 count 数据
# 所以我们需要跳过前面的注释行（第 1 行是命令信息，用 skip=1 跳过）
count_data <- read.table("all_samples_counts.txt",
header = TRUE, # 第一行是列名
row.names = 1, # 第一列（Geneid）作为行名
skip = 1, # 跳过第一行（featureCounts 命令信息）
sep = "\t") # tab 分隔
# 只保留 count 列（去掉 Chr, Start, End, Strand, Length 这 5 列）
count_data <- count_data[, 6:ncol(count_data)]
# 看看数据长什么样（前 5 个基因、前 3 个样本）
head(count_data[, 1:3])
# ============================================
# 第二步：准备样本信息表
# ============================================
# condition: 每个样本属于哪个组（T2D 糖尿病组 vs Control 健康组）
# 注意：行名必须和 count_data 的列名完全一致！
sample_info <- data.frame(
condition = factor(c("T2D", "T2D", "T2D", # 前 3 个是糖尿病组
"Control", "Control", "Control"), # 后 3 个是对照组
levels = c("Control", "T2D")) # 设定 Control 为参考水平
)
# 行名设为样本名（要和 count_data 的列名一致）
rownames(sample_info) <- colnames(count_data)
# 检查样本名是否一致（必须返回 TRUE）
all(rownames(sample_info) == colnames(count_data))
# ============================================
# 第三步：构建 DESeqDataSet 对象
# ============================================
# design = ~ condition: 告诉 DESeq2 按 condition 分组比较
dds <- DESeqDataSetFromMatrix(
countData = count_data, # count 矩阵
colData = sample_info, # 样本信息表
design = ~ condition # 实验设计公式
)
# ============================================
# 第四步：预过滤低表达基因
# ============================================
# 去掉那些在所有样本中 count 都很低的基因（没有分析价值，还浪费计算时间）
# 规则：至少在 3 个样本中 count >= 10（3 是最小组的样本数）
smallestGroupSize <- 3
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep, ]
# 打印过滤后剩多少基因
cat("过滤后剩余基因数:", nrow(dds), "\n")
# ============================================
# 第五步：运行 DESeq2 差异分析（核心一行代码！）
# ============================================
# 这一步内部自动完成：估计 size factors → 估计离散度 → Wald 检验
dds <- DESeq(dds)
# ============================================
# 第六步：提取结果
# ============================================
# contrast: 指定比较方向 → T2D vs Control
# 结果中 log2FoldChange > 0 表示 T2D 组表达更高
res <- results(dds, contrast = c("condition", "T2D", "Control"))
# 用 lfcShrink 收缩 log2FoldChange（推荐，让可视化更合理）
# type="apeglm": 最新推荐的收缩方法
res_shrunk <- lfcShrink(dds,
coef = "condition_T2D_vs_Control",
type = "apeglm")
# 查看结果摘要
summary(res)
# ============================================
# 第七步：筛选差异基因
# ============================================
# 条件：padj < 0.05 且 |log2FoldChange| > 1
sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
cat("显著差异基因数:", nrow(sig_genes), "\n")
# 按 padj 排序，最显著的排前面
sig_genes <- sig_genes[order(sig_genes$padj), ]
# 查看前 10 个最显著的差异基因
head(as.data.frame(sig_genes), 10)
# 导出差异基因结果到 CSV 文件
write.csv(as.data.frame(sig_genes),
file = "DEG_T2D_vs_Control.csv",
row.names = TRUE)
# 上调基因（T2D 中表达升高的）
up_genes <- subset(sig_genes, log2FoldChange > 1)
cat("上调基因数:", nrow(up_genes), "\n")
# 下调基因（T2D 中表达降低的）
down_genes <- subset(sig_genes, log2FoldChange < -1)
cat("下调基因数:", nrow(down_genes), "\n")
4.4 火山图绘制代码¶
# ============================================
# 火山图（Volcano Plot）
# ============================================
# 火山图是展示差异分析结果最经典的图：
# X 轴：log2FoldChange（表达变化倍数）
# Y 轴：-log10(padj)（统计显著性，值越大越显著）
# 颜色：标记上调/下调/不显著的基因
# --- 方法一：用 ggplot2 手动画（推荐掌握） ---
# 把 DESeq2 结果转成 data.frame
res_df <- as.data.frame(res)
res_df <- na.omit(res_df) # 去掉含 NA 的行
# 添加一列标记：上调 / 下调 / 不显著
res_df$change <- "Not Significant" # 默认不显著
res_df$change[res_df$padj < 0.05 & res_df$log2FoldChange > 1] <- "Up" # 上调
res_df$change[res_df$padj < 0.05 & res_df$log2FoldChange < -1] <- "Down" # 下调
res_df$change <- factor(res_df$change,
levels = c("Up", "Down", "Not Significant"))
# 画图
ggplot(res_df, aes(x = log2FoldChange,
y = -log10(padj),
color = change)) +
geom_point(alpha = 0.6, size = 1.2) + # 散点图，透明度 0.6
scale_color_manual(values = c("Up" = "red", # 上调基因红色
"Down" = "blue", # 下调基因蓝色
"Not Significant" = "grey")) + # 不显著灰色
geom_vline(xintercept = c(-1, 1), # 画 LFC 阈值线
linetype = "dashed", color = "black") +
geom_hline(yintercept = -log10(0.05), # 画 padj 阈值线
linetype = "dashed", color = "black") +
labs(title = "Volcano Plot: T2D vs Control", # 标题
x = "log2 Fold Change", # X 轴标签
y = "-log10(adjusted p-value)", # Y 轴标签
color = "Change") + # 图例标题
theme_bw() + # 白色背景主题
theme(plot.title = element_text(hjust = 0.5)) # 标题居中
# 保存图片
ggsave("volcano_plot.pdf", width = 8, height = 6)
ggsave("volcano_plot.png", width = 8, height = 6, dpi = 300)
# --- 方法二：用 EnhancedVolcano 包（一行代码，效果更好看） ---
library(EnhancedVolcano)
EnhancedVolcano(res,
lab = rownames(res), # 基因名标签
x = 'log2FoldChange', # X 轴
y = 'padj', # Y 轴
title = 'T2D vs Control', # 标题
pCutoff = 0.05, # padj 阈值
FCcutoff = 1, # LFC 阈值
pointSize = 2.0, # 点大小
labSize = 3.0) # 标签字号
5. 与项目的关联¶
该 T2D（2 型糖尿病）肠道菌群项目虽然主要是宏基因组方向，但在分析流程中用到了转录组分析的核心工具：
5.1 featureCounts 在你项目中的用法¶
- 场景：该宏基因组项目中使用 featureCounts 对比对到参考基因组后的 reads 进行基因水平定量
- 输入：经过 Bowtie2/BWA 比对后的 BAM 文件 + 注释文件
- 输出：基因 count 矩阵（每个基因在每个样本中有多少 reads）
- 参数选择：
-p --countReadPairs：双端测序按片段计数-t exon -g gene_id：按外显子统计、按基因汇总-T 8：多线程加速
5.2 DESeq2 在你项目中的用法¶
- 目的：比较 T2D 患者和健康对照的肠道微生物基因表达差异
- 分组设计：
design = ~ condition（T2D vs Control） - 筛选标准：
padj < 0.05 & |log2FoldChange| > 1 - 下游分析：差异基因做 GO/KEGG 富集分析，找到与 T2D 相关的代谢通路
5.3 面试时怎么串联¶
"该项目是宏基因组方向，研究 2 型糖尿病患者肠道菌群的组成和功能差异。在定量环节，我用 featureCounts 统计每个基因在各样本中的 read counts，因为 featureCounts 速度快、支持多线程，很适合批量处理。得到 count 矩阵后，我用 DESeq2 做差异分析，以 padj < 0.05、|LFC| > 1 为阈值筛选差异基因，最后对差异基因做 KEGG 富集分析，发现 T2D 组在短链脂肪酸代谢等通路上有显著变化。"
6. 面试怎么答¶
Q1: 请简述 RNA-seq 分析的基本流程¶
答：RNA-seq 分析主要分 5 步。第一步质控，用 FastQC 看数据质量、用 fastp 去接头和低质量 reads。第二步比对，用 HISAT2 或 STAR 把 reads 比对到参考基因组上，注意要用 splicing-aware 的比对器，因为 mRNA 有剪接。第三步定量，用 featureCounts 统计每个基因被多少 reads 覆盖，得到 count 矩阵。第四步差异分析，用 DESeq2 比较两组间的表达差异，筛选 padj < 0.05 且 LFC 绝对值大于 1 的差异基因。第五步功能富集，用 clusterProfiler 做 GO 和 KEGG 分析，看差异基因富集在什么通路上。
Q2: DESeq2 和 edgeR 有什么区别？你为什么选 DESeq2？¶
答：两个包都基于负二项分布建模 RNA-seq 的 count 数据，但标准化方法不同——DESeq2 用 median of ratios，edgeR 用 TMM。DESeq2 相对保守，假阳性控制更严格；edgeR 在大样本时可以用 QLF 检验。实际上两者的差异基因结果大约 80% 重叠。我选 DESeq2 主要是因为：第一，使用简单，核心就一个
DESeq()函数；第二，内置 LFC 收缩功能，可视化效果好；第三，文献引用最广泛，审稿人接受度高。
Q3: 为什么 RNA-seq 不能用普通的 DNA 比对器（如 BWA）？¶
答：因为 mRNA 经过了剪接，内含子被切掉了。所以一条 read 可能跨越两个外显子，在基因组上是不连续的，中间可能隔了几千甚至几万碱基的内含子。BWA 这类比对器假设 read 在基因组上是连续的，处理不了这种"跳跃式"比对。所以需要 HISAT2、STAR 这种 splicing-aware aligner，它们能识别并正确处理跨越剪接位点的 reads。
Q4: featureCounts 和 Salmon 有什么区别？¶
答：最大的区别是 featureCounts 需要先比对再计数，而 Salmon 不需要比对，直接从 FASTQ 文件估算 transcript 表达量。featureCounts 统计的是 gene-level 的 raw counts，直接可以给 DESeq2 用。Salmon 统计的是 transcript-level 的 TPM 和 estimated counts，需要通过 tximport 包导入 DESeq2。Salmon 的优势是速度快、能校正基因长度偏差，DESeq2 官方现在也推荐 Salmon + tximport 的流程。但 featureCounts 最简单直观，一行命令搞定，适合入门和常规分析。
Q5: DESeq2 中为什么要用 padj 而不是 pvalue 来筛选差异基因？¶
答：因为 RNA-seq 同时检验了上万个基因，这是典型的多重检验问题。如果直接用 p < 0.05 的阈值，2 万个基因里纯靠运气就会有 1000 个假阳性。所以必须做多重检验校正。DESeq2 默认用 Benjamini-Hochberg（BH）方法校正 p 值，得到 padj（adjusted p-value），它控制的是 FDR（False Discovery Rate，假发现率），也就是说如果我设 padj < 0.05，意思是"我接受筛出来的差异基因中最多有 5% 是假阳性"。
Q6: 什么是 log2FoldChange？LFC=2 是什么意思？¶
答：log2FoldChange 就是两组表达量比值取 log2。LFC=2 意味着实验组是对照组的 2^2 = 4 倍表达。LFC=1 是 2 倍，LFC=-1 是 0.5 倍（也就是下调了一半）。用 log2 的好处是让上调和下调对称：上调 2 倍是 LFC=1，下调 2 倍是 LFC=-1，在火山图上看起来对称美观。
Q7: RPKM、FPKM、TPM 有什么区别？差异分析应该用哪个？¶
答：三者都是同时校正测序深度和基因长度的表达量指标。RPKM 用于单端测序，FPKM 用于双端测序（以 fragment 为单位避免重复计数），两者计算顺序一样：先除测序深度再除基因长度。TPM 的不同之处是先除基因长度再除测序深度，好处是所有样本的 TPM 总和都等于 100 万，可以直接跨样本比较基因占总表达量的比例。RPKM/FPKM 不同样本总和可能不同，直接比较有偏差。所以现在推荐用 TPM。但做差异分析时，不应该用 TPM/FPKM，而应该用 raw counts 配合 DESeq2 或 edgeR 自带的标准化方法，因为这些工具需要原始计数数据来正确建立统计模型。
Q8: 火山图怎么看？¶
答：火山图 X 轴是 log2FoldChange，Y 轴是 -log10(padj)。右上角的红点是显著上调基因——变化大且统计显著。左上角的蓝点是显著下调基因。中间灰色的是不显著的基因。两条竖虚线是 LFC 阈值（通常 ±1），一条横虚线是 padj 阈值（通常 0.05）。超过这些线的基因才被认为是差异表达基因。
7. 延伸阅读¶
| 资源 | 链接 / 说明 |
|---|---|
| DESeq2 官方教程 | Bioconductor DESeq2 vignette（最权威的使用手册） |
| featureCounts 官网 | Subread/featureCounts 官方文档 |
| RNA-seq 分析工作流 | Bioconductor RNA-seq workflow（从 FASTQ 到差异分析的完整 R 流程） |
| HISAT2 论文 | Kim et al. (2015) Nature Methods |
| STAR 论文 | Dobin et al. (2013) Bioinformatics |
| DESeq2 论文 | Love et al. (2014) Genome Biology（必读） |
| edgeR 论文 | Robinson et al. (2010) Bioinformatics |
| Salmon 论文 | Patro et al. (2017) Nature Methods |
| clusterProfiler 原始版 | Yu et al. (2012) OMICS: A Journal of Integrative Biology |
| clusterProfiler 4.0 | Wu et al. (2021) The Innovation（推荐引用此版本，支持更多物种和通路数据库） |
| Galaxy RNA-seq 教程 | Galaxy Project 提供的图形界面 RNA-seq 分析教程，适合零代码基础入门 |
| 统计学背景 | 本知识库 11_统计学基础.md（p 值、多重检验校正等概念） |
| 测序技术原理 | 本知识库 13_测序技术原理.md（Illumina 测序原理） |
最后提醒：面试时最重要的不是背工具名称，而是能讲清楚"为什么这么做"。比如：为什么要标准化？因为不同样本测序深度不同。为什么用负二项分布？因为 RNA-seq 的 count 数据有过度离散（overdispersion），泊松分布不够用。为什么要做 LFC 收缩？因为低表达基因的 LFC 估计噪声大，收缩后更可靠。能说出这些"为什么"，面试官会觉得你是真懂，而不是只会跑命令。