转录组分析流程（RNA-seq Analysis Pipeline）¶

一句话说明¶

转录组分析（RNA-seq）是把细胞在某个时刻"正在表达的基因"全部测出来，通过 质控 → 比对 → 定量 → 差异分析 → 功能富集 的流程，找到不同条件下哪些基因"开得更多"或"关得更多"，从而揭示疾病机制或生物学过程。你项目中用 featureCounts 做定量、DESeq2 做差异分析，就是这个流程的核心步骤。

目录¶

1. 核心概念（白话版）
2. 转录组分析完整流程
3. 关键工具对比
4. 实操代码
5. 和你项目的关联
6. 面试怎么答
7. 延伸阅读

1. 核心概念（白话版）¶

1.1 什么是转录组（Transcriptome）¶

定义：转录组就是一个细胞/组织在某个时间点，所有被"转录"出来的 RNA 的总和。转录（Transcription）就是 DNA 上的基因被"抄写"成 mRNA 的过程。
白话比方：如果基因组（Genome）是一本完整的"菜谱大全"（包含所有菜的做法），那转录组就是"今天厨房实际在做的菜的清单"——虽然菜谱里有 2 万多道菜，但今天只做了其中几千道，每道菜做的份数还不一样。
为什么要研究转录组：
基因组是"静态"的（每个细胞的 DNA 序列基本一样）
转录组是"动态"的（不同组织、不同时间、健康 vs 疾病，表达的基因完全不同）
通过比较疾病组 vs 健康组的转录组差异，就能找到和疾病相关的关键基因

1.2 RNA-seq 实验流程概述¶

RNA-seq（RNA Sequencing，RNA 测序）是研究转录组最主流的技术：

实验流程（湿实验部分）：
提取 RNA → mRNA 富集（polyA 选择或 rRNA 去除） → 逆转录成 cDNA → 建库 → 上机测序 → 得到 FASTQ 文件

步骤	做什么	白话解释
提取 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 / 早期版本为 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 参数允许

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（余光创团队开发，中国人做的，引用超 2 万次）

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` 参数）	不支持	支持
使用难度	简单，一行命令	简单但慢	需理解 quasi-mapping 概念（quasi-mapping：不做完整比对只做粗定位，像扫一眼封面就知道是哪本书）
发表年份	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()`，可视化更友好	需额外处理
引用次数	~45000+（截至 2025）	~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: 火山图怎么看？¶

答：火山图 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
Galaxy RNA-seq 教程	Galaxy Project 提供的图形界面 RNA-seq 分析教程，适合零代码基础入门
统计学背景	本知识库 `11_统计学基础.md`（p 值、多重检验校正等概念）
测序技术原理	本知识库 `13_测序技术原理.md`（Illumina 测序原理）

最后提醒：面试时最重要的不是背工具名称，而是能讲清楚"为什么这么做"。比如：为什么要标准化？因为不同样本测序深度不同。为什么用负二项分布？因为 RNA-seq 的 count 数据有过度离散（overdispersion），泊松分布不够用。为什么要做 LFC 收缩？因为低表达基因的 LFC 估计噪声大，收缩后更可靠。能说出这些"为什么"，面试官会觉得你是真懂，而不是只会跑命令。