我最早接触Bulk RNA-seq的时候，脑子里其实只有个很模糊的概念：把细胞里的RNA提出来测一测，然后看看哪些基因表达变高了、哪些变低了。真到自己手里要跑完一整个项目时才发现，从实验设计到最终拿到一张可以放进文章里的热图，中间隔着一大串需要做判断的环节——参考基因组选哪个版本、比对用STAR还是Hisat2、定量该相信featureCounts还是Salmon、差异分析为什么DESeq2给出来的padj跟别人给的p-value对不上。任何一个环节理解不到位，最后跑出来的结果就可能在组会上被问得支支吾吾。
这篇内容就是围绕小鼠Bulk RNA-seq的全流程展开的，适合刚入手转录组分析、或者已经跑通了流程但对每一步的原理和选型逻辑还不太清楚的人。我会按照实际操作顺序，把实验设计阶段就要想清楚的事、原始数据下来之后的质控、比对定量的原理与命令、差异表达分析的核心逻辑、富集分析和可视化这几大块一次讲透，然后把那些文档里不会写、但实际做起来非常容易踩的坑单独拎出来说。
1. 开始之前：搞清楚你手上的数据到底意味着什么
很多人容易把单细胞和Bulk混在一起想，实际上两者的分析逻辑差异非常大。Bulk RNA-seq测的是某一群细胞的平均转录状态，就像把一锅汤搅匀了再尝一口，你能知道整锅平均是什么味道，但你不知道里面哪块肉咸了、哪根菜淡了。单细胞则是把每个细胞当成独立样品分别去测，信息粒度完全不同。这意味着Bulk RNA-seq解决的核心问题是"两组之间的平均表达差异"，而不是"哪个细胞亚群在响应刺激"。
搞清楚这个前提之后，再来看整个流程，会发现所有步骤其实都指向同一个目标：得到一个可靠的基因表达矩阵，然后再从这个矩阵里找出有统计学依据的差异基因。
1.1 生物学重复：最容易省、也最不能省的一环
实验设计阶段最大的坑，是把技术重复当成生物学重复来用。我见过有人把同一个RNA样本分开建两次库、跑两次测序，以为这就是两个重复了。这不是重复，这只是把同样一份材料测了两遍。技术重复解决的是操作误差，而生物学重复解决的是个体差异——同一种处理下不同小鼠之间的天然波动。
在小鼠模型里，生物学重复数量最好不低于3个，条件允许就做5个。重复数直接决定差异分析能不能出显著基因。三个重复勉强够用，但组间本身表达差异不大的话，会发现padj很难达到0.05；五个重复的时候，统计效力就会明显改善。原因倒不复杂：DESeq2这类工具本质是在估计生物学变异的大小，重复数太少，变异估计不准，离散度参数就会飘，最后检验统计量的可靠性自然就打了折扣。
1.2 批次效应要从实验源头就控制
做小鼠实验时，不同批次的动物、不同时间点取样、不同人操作，都会引入批次效应。转录组分析的下游虽然可以用ComBat-seq这类工具去矫正批次，但这是最后一道防线，而不是一开始就拿来兜底的。最好的策略是在实验设计时把处理组和对照组在同一天、用同一批小鼠、随机分组来完成。如果实在做不到，至少把批次信息记录清楚，后面在DESeq2的设计公式里用~ batch + condition这样的方式纳入模型，比事后矫正要自然得多。
1.3 建库策略与链特异性：现在就要做的决定
RNA提取之后建库，核心要选的是两件事：富集方式和链特异性。
- polyA富集：适合mRNA分析，可以去掉rRNA占比，数据利用率高，是目前转录组研究的主流方案。缺点是检测不了没有polyA尾巴的RNA（比如部分lncRNA和circRNA）。
- rRNA去除：保留更全的转录本种类，适合需要看lncRNA或者降解程度较高的FFPE样本。缺点是数据里rRNA残留比例可能较高，浪费一部分测序量。
链特异性也是一个看似不起眼、实则影响很大的选项。链特异性文库（strand-specific）能告诉我们转录本来自哪条DNA链，这个信息在基因注释冲突区、反义转录本分析时非常关键。如果建库时做了链特异性，但下游比对定量时没有指定对应参数，会产生大量错误定位的reads，表现为基因注释率很低，或者反义链上也出现莫名其妙的信号。
2. 数据落地后的第一关：原始reads质控
测序公司交付的数据通常是fastq.gz格式。拿到数据之后，第一件事不是急着比对，而是做质控。这个环节的目标有两个：一是确认测序质量到底行不行，二是决定哪些reads需要修剪。
2.1 FastQC报告怎么看，别只盯着绿色勾
FastQC会生成一份HTML报告，里面有十来个模块。刚入门的同学最容易犯的错是只看Per base sequence quality是不是绿色，然后就觉得自己数据没问题了。真正需要认真看的是下面几项：
- Per base sequence content：正常情况下每个位点四种碱基的比例应该比较接近，如果前几个碱基出现明显偏移，多半是接头残留或者引物污染。
- Adapter Content：直接反映接头残留比例，这个值如果超过5%，就有必要修剪。
- Overrepresented sequences：如果某条序列出现超高比例，有可能是接头二聚体、rRNA残留或者建库时的PCR重复。
- Sequence Duplication Levels：高重复率可能来自PCR扩增过度，也可能是本身丰度就高的RNA种类（比如rRNA），要结合Overrepresented sequences一起判断。
用MultiQC把所有样本的FastQC报告汇总成一份总报告，比一个个打开HTML文件效率高得多。运行方式没什么门槛：
multiqc ./它会自动扫描当前目录下的fastqc_html和fastqc_zip，输出一个汇总页面。
2.2 修剪与过滤的标准流程
质控确认有问题之后，用cutadapt或者Trimmomatic做修剪。这里以cutadapt为例，最常用的参数长这样：
cutadapt \ -a AGATCGGAAGAGC \ -j 8 \ -q 20 \ -m 20 \ -o clean.fastq.gz \ input.fastq.gz解释一下每一步的含义：-a指定要切除的接头序列，这个序列要跟建库试剂盒对应，通常Illumina TruSeq的接头是AGATCGGAAGAGC，但具体还是要看建库报告；-j是线程数；-q 20表示从reads末端开始切除质量低于20的碱基；-m 20表示修剪后长度不足20bp的reads直接丢弃。
修剪之后需要重新跑一次FastQC，确认接头残留问题已经解决了再往下走。这个回头看的动作很关键，很多人修剪完直接比对，结果发现比对率上不去，又要回头查是不是修剪参数出了问题。
2.3 别急着跑比对：先核对基因组版本与注释文件
这一步是整个流程里最容易被忽视、却也最容易导致"全盘皆输"的地方。小鼠参考基因组现在常见的两个版本是mm10/GRCm38和mm39/GRCm39。mm39是2020年发布的，修复了不少mm10的错误组装和gap区域。如果你手里的GTF注释文件是mm10的，参考基因组用的却是mm39，那STAR在比对时就会报一堆warning，甚至大量reads比对不上。
我自己的习惯是先在项目目录下建一个清晰的目录结构，把参考基因组、注释文件、原始数据分开存放，同时写一个README记录版本信息：
data/ raw_fastq/ clean_fastq/ reference/ genome_mm39.fa gencode.vM32.annotation.gtf scripts/ results/小鼠的注释文件我推荐用GENCODE的版本，小鼠对应gencode.vM系列，人类对应gencode.vH系列。选版本的时候要注意GTF和参考基因组要来自同一套版本体系，不要混搭。
3. 比对与定量：从fastq到表达矩阵
这个环节是整个Bulk RNA-seq分析的核心。目标是两件事：把reads比对到参考基因组上，然后统计每个基因覆盖了多少reads，最终得到表达矩阵。
3.1 STAR的种子-延伸策略：为什么它能对付剪接事件
STAR是目前Bulk RNA-seq比对的主流选择。它之所以快，核心在于两步走的种子-延伸策略。第一步，它会在参考基因组上找一段足够长的完全匹配片段（seed），基因组的索引以suffix array方式组织，这一步可以非常快地定位reads的大致位置；第二步，遇到reads跨过外显子边界的情况（也就是junction reads），STAR会同时找到两个exon端的匹配，然后用动态规划算法把中间的剪接事件连起来。
这个设计的直接好处是：RNA-seq数据里很大比例的reads是跨剪接位点的，传统的线性比对工具（比如BWA）是为DNA测序设计的，遇到这种reads会直接比对失败，导致基因定量严重偏低。STAR专门绕开了这个问题。
建立索引和比对的标准命令如下：
# 建立索引（只用做一次） STAR --runMode genomeGenerate \ --genomeDir ./reference/star_index/ \ --genomeFastaFiles ./reference/genome_mm39.fa \ --sjdbGTFfile ./reference/gencode.vM32.annotation.gtf \ --sjdbOverhang 149 \ --runThreadN 16 # 比对 STAR --runMode alignReads \ --genomeDir ./reference/star_index/ \ --readFilesIn clean_1.fastq.gz clean_2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outFileNamePrefix sample1_ \ --runThreadN 16--sjdbOverhang这个参数一般设为读长减1。比如150bp双端测序，就填149。这个值是给剪接位点数据库预留的长度匹配边界用的，填错了不会报错，但会轻微影响junction reads的比对准确度。
3.2 featureCounts做定量：GTF和链式参数是两大雷区
得到BAM文件之后，下一步就是统计每个基因的reads数。featureCounts因为速度快、内存占用低，是绝大多数人的第一选择。
featureCounts \ -a gencode.vM32.annotation.gtf \ -o counts.txt \ -s 2 \ -T 16 \ -p \ --countReadPairs \ sample1_sorted.bam sample2_sorted.bam sample3_sorted.bam这里有两个参数非常关键。第一个是-s，它指定链式信息：0代表非链特异性，1代表链特异性文库且reads来自和注释链相同的方向，2代表反向。如果你建库的时候用的是链特异性试剂盒（常见的是dUTP法），那-s 2通常是正确的选择。如果填反了，你会看到一个诡异的结果：exon区的reads计数只剩下一半，大量reads被算到基因间区去了。
第二个关键是GTF注释版本必须和STAR索引时的一致。用mm39的基因组搭配vM32的注释，没问题；但如果你不小心用了vM25（那是mm10时代的注释），那featureCounts会忽略掉大量新注释的转录本，最后得到的基因数量和总reads数都会缩水。
3.3 要不要用Salmon走alignment-free路线
Salmon这类alignment-free工具不生成BAM文件，直接基于k-mer和转录组序列做定量，速度比STAR+featureCounts的组合快一个量级。它适合什么场景？一种是样本量特别大（上百个样本），比对到基因组再定量的计算成本过高；另一种是你只关心转录本或基因表达水平，不需要处理变异信息。
但Salmon有个前提要注意：它基于参考转录组进行定量，参考转录组的完整度直接决定结果质量。如果某个基因的可变剪接异构体没有包含在参考里，它的表达量就可能被低估或者错配。另外，Salmon定量后再转成count矩阵做差异分析，还需要用tximport来汇总到基因水平。流程上多了一步，但并没有多复杂。
从实操角度看，如果项目只有几个样本、几十个样本，STAR+featureCounts仍然是更稳妥的选择，因为它保留了BAM文件，后续想看某个基因的比对情况、检查可疑区域，都可以直接回到BAM里查。Salmon的定量结果要回看比对情况就不那么直观了。
3.4 为什么表达量最好用TPM，别再用RPKM/FPKM
拿到counts之后，你会发现不同基因的reads数直接比较是不公平的——基因越长，理论上能比对上的reads就越多；测序深度越深，所有基因的reads数都会成比例增加。RPKM/FPKM先把reads数除以基因长度，再做总reads数归一化。问题在于，它用样本的总reads数做分母，而在总reads数固定的情况下，少数极高表达的基因会把"分母"撑大，导致其他基因的RPKM被系统性压低。
TPM的计算顺序恰好反过来：先按基因长度归一化每个基因的reads数（得到per-kilobase值），再按基因之间的比例去归一化到"每百万条转录本"。这样做的好处是，不同样本的TPM总和恒定为1e6，样本之间可以直接比较。所以如果你要画表达量箱线图、做样本间相关性、或者展示单个基因的表达，用TPM。但如果要做差异分析，就别自己"归一化"，把原始counts直接交给DESeq2去处理。
4. 差异表达分析：DESeq2的工作逻辑与实操
拿到counts矩阵，进入真正的核心环节——差异表达分析。现在主流的工具有DESeq2、edgeR、limma-voom三个，它们的统计模型和处理逻辑各有侧重，但对常规的小鼠Bulk RNA-seq来说，DESeq2是使用最广泛、结果最稳定、社区资料最丰富的一个。
4.1 为什么不能只看倍数变化，还要有离散度模型
最朴素的差异分析思路是：算每个基因在两组之间的平均表达倍数变化，然后用t检验看p值。但转录组数据有个特点：低表达基因的变异往往比高表达基因大得多，而且counts数据是离散的，不是连续的正态分布。如果直接对counts做t检验，会得到大量假阳性。
DESeq2的核心逻辑是，先把每个基因的counts建模为负二项分布，这个分布有两个关键参数：均值表达量和离散度。离散度描述了同一个基因在不同样本之间的波动程度。在样本量少的情况下，单独估计每个基因的离散度不准，DESeq2的做法是把所有基因的离散度-均值关系拟合一条曲线，然后把每个基因的离散度往曲线上收缩（shrinkage），让低表达、低重复数基因的离散度估计更稳定。这一步直接决定了最终的p值和padj质量。
4.2 在R里完整跑一遍DESeq2
下面是一段可以直接照用的代码：
library(DESeq2) library(tximport) # 读入featureCounts结果 countdata <- read.table("counts.txt", header = TRUE, row.names = 1, sep = "\t") countdata <- countdata[, 6:ncol(countdata)] # 构建样本信息表 coldata <- data.frame( condition = factor(c("ctrl", "ctrl", "ctrl", "treat", "treat", "treat"), levels = c("ctrl", "treat")) ) rownames(coldata) <- colnames(countdata) # 构建DESeq2对象 dds <- DESeqDataSetFromMatrix( countData = round(countdata), colData = coldata, design = ~ condition ) # 过滤低表达基因（至少10个reads总和的基因） keep <- rowSums(counts(dds)) >= 10 dds <- dds[keep, ] # 核心运行 dds <- DESeq(dds) # 提取结果，规定比较方向：treat相对于ctrl res <- results(dds, contrast = c("condition", "treat", "ctrl"))你可能会注意到我用了round(countdata)。这是因为DESeq2要求输入整数counts，而featureCounts输出的通常是整数，但如果你之前对counts做了任何算术操作，就可能出现小数，必须舍入。
4.3 结果筛选：padj、log2FoldChange与shrinkage
results()返回的结果里，有几列必须搞清楚：baseMean是基因在所有样本中的平均归一化counts；log2FoldChange是处理组相对于对照组的倍数变化取log2；lfcSE是log2FoldChange的标准误；pvalue和padj分别是原始p值和多重检验校正后的p值（Benjamini-Hochberg方法）。
筛差异基因的标准组合通常是：
res_df <- as.data.frame(res) res_df$gene <- rownames(res_df) sig <- res_df[which(res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1), ]padj < 0.05控制假阳性率，|log2FC| > 1对应两倍以上变化，是生物学意义和统计学意义的双重门槛。如果你发现筛出来的基因太少，可以尝试放宽到|log2FC| > 0.58，也就是1.5倍变化，但必须在文章中明确说明阈值。
还有一个很容易被忽略的细节：低表达基因的log2FoldChange非常不稳定，两个样本间一个reads的差异就能产生很大的倍数变化。这时候需要看shrunken estimates（收缩后的倍数变化）。用lfcShrink函数：
res_shrink <- lfcShrink(dds, contrast = c("condition", "treat", "ctrl"), type = "apeglm")收缩后的log2FoldChange更保守，画火山图时那些低表达基因不会因为极端的fold change跑到图的边缘去干扰视线。
5. 富集分析：把基因列表变回生物学意义
差异基因列表本身只是基因名，真正的生物学故事要靠富集分析来讲。最常用的是GO（Gene Ontology）和KEGG通路富集。用R里的clusterProfiler包，配以小鼠的org.Mm.eg.db注释库，即可完成整个分析。
5.1 GO与KEGG：两个数据库的使用边界要求
GO注释分为三个本体：生物过程（Biological Process）、分子功能（Molecular Function）和细胞组分（Cellular Component）。做转录组富集时，最常看的是BP，因为差异基因是"启动了某些生物学程序"的反映。
KEGG提供的是代谢和信号通路层面的注释，直接对应"通路"这个概念。它跟GO最大的区别在于，KEGG通路之间存在层级关系，比如某个基因既属于"PI3K-Akt signaling pathway"又属于"FoxO signaling pathway"，这是正常的，因为通路本身是互相交错的网络。
使用上有一个常见坑：做KEGG富集时，clusterProfiler需要联网下载KEGG数据，如果网络环境受限，会一直卡住。解决方法是先下载好KEGG db文件，后面离线也可以用。
5.2 clusterProfiler实操，从基因名到富集结果
library(clusterProfiler) library(org.Mm.eg.db) # 把差异基因的symbol转换为ENTREZID gene_symbols <- sig$gene gene_entrez <- bitr(gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Mm.eg.db) # GO富集 ego <- enrichGO(gene = gene_entrez$ENTREZID, OrgDb = org.Mm.eg.db, keyType = "ENTREZID", ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05) # KEGG富集 ekegg <- enrichKEGG(gene = gene_entrez$ENTREZID, organism = "mmu", pvalueCutoff = 0.05)KEGG富集时organism = "mmu"对应小鼠，如果用了人类对应的"hsa"，全部基因映射不上，结果为空，这种错误经常出现。
5.3 富集结果的常见误读方式
富集分析给出来的是"统计上显著富集的条目"，不等于"生物学上最重要的通路"。p值只告诉你这个通路的基因比例比随机抽样更高，它不管通路的上下游关系。真正的解读，需要你把富集到的通路放回实验背景里看逻辑链：比如炎症模型中富集到TNF signaling pathway和NF-kappa B signaling pathway，这两个通路本身就有明确的上下游关系，这时候能讲出一个完整的故事。如果几个富集通路之间毫无关联，那更可能是筛选出的基因集合本身特异性不够，或者富集分析用错了基因列表。
另一个误读是把富集条的p值当成效应量。GO条目的p值受通路基因数量的影响较大，大通路的p值天然更容易显著。多个通路之间比较时，不要只凭p值判断哪个通路"变化最大"，要看GeneRatio（差异基因中属于该通路的比例）。
6. 可视化：用来投稿和汇报的几张图
分析做完，最终要落实到图片上。组会汇报、论文投稿、以及补充材料里最常用的几张图有：PCA/样本相关性热图、火山图、差异基因热图、富集分析气泡图。每一张的用法和画法都值得单独说。
6.1 PCA与样本相关性：先确认整体数据质量
任何下游分析之前，我都会先画PCA和样本相关性的图。PCA的主要用途是看样本的分组趋势：同一处理组的样本应该在PCA空间里聚在一起，如果对照组和处理组完全没有分开，那后续差异分析即使出了结果，也要谨慎对待。更严重的情况是PCA按批次而不是按处理分开，说明批次效应非常强。
library(ggplot2) library(DESeq2) # 使用vst变换后的数据做PCA，避免高表达基因主导 vsd <- vst(dds, blind = TRUE) pca_data <- plotPCA(vsd, intgroup = "condition", returnData = TRUE) percentVar <- round(100 * attr(pca_data, "percentVar")) ggplot(pca_data, aes(PC1, PC2, color = condition)) + geom_point(size = 3) + xlab(paste0("PC1: ", percentVar[1], "% variance")) + ylab(paste0("PC2: ", percentVar[2], "% variance")) + theme_classic()样本相关性热图可以用全部基因的TPM值来算Spearman相关系数，组内样本相关性一般要高于0.95。如果出现某个样本跟其它重复的相关性特别低，要考虑是不是样本标签搞混了、文库构建失败或者测序深度严重偏低。
6.2 火山图和热图的操作细节
火山图展示的是所有基因的p值和倍数变化关系，用ggplot2画，x轴是log2FC，y轴是-padj取log10，把显著基因标色。
热图展示的是关键基因在不同样本里的表达模式。请特别注意以下几点：第一个，热图里展示的应当是归一化后的表达量（vst或TPM），直接用原始counts会导致文库大小差异反映成表达差异。第二个，画热图前对每行做z-score归一化，把基因的表达模式变成"相对自身均值的偏离"，否则高表达和低表达基因放在同一色阶里，低表达基因的变化根本看不出来。第三个，小鼠的基因symbol首字母大写、其余小写，跟人类的写法不同（Gapdh vs GAPDH），投稿时基因名的写法要符合物种命名规则。
library(pheatmap) # 提取vst表达矩阵并对显著差异基因取行子集 mat <- assay(vsd)[sig$gene, ] mat <- t(scale(t(mat))) # 标注样本的分组信息 annotation_col <- data.frame(condition = coldata$condition) rownames(annotation_col) <- colnames(mat) pheatmap(mat, scale = "none", annotation_col = annotation_col, show_rownames = TRUE, show_colnames = TRUE, color = colorRampPalette(c("navy", "white", "firebrick3"))(50))这个配色方案是在白色背景下把低表达映射到深蓝、高表达映射到红色，是转录组热图里最常用的配色之一。需要注意scale = "none"，因为我们已经手动做了行标准化，这里不要再重复缩放。
6.3 从"跑完"到"可复现交付"：参数记录的重要性
最后一个建议可能看起来不那么"技术"，但在我眼里是最有价值的：从拿到fastq数据的那一刻起，就用笔记软件记录每一步的软件版本和关键参数。STAR是哪个版本、featureCounts的-s值是多少、DESeq2用的参考基因组是mm39还是mm10、筛差异基因的阈值是多少。原因很简单：审稿人或者师兄师姐随时可能问你要这些细节。转录组分析流程的复现性完全依赖这些"元数据"，而这些信息如果当时没记，事后靠回忆基本补不回来。
我自己吃过一次亏。某次项目复盘时想要重新生成结果，发现当时GTF注释文件下载的是哪个版本已经想不起来了，只能从比对率、基因数量这些指标倒推排错，折腾了好几天。后来我就在服务器上养成了一个习惯：每个项目的根目录放一个params.yaml，记录基因组版本、索引路径、比对软件版本、定量参数、差异分析阈值这些关键信息，随时可以追溯。
还有一个直接影响效率的小经验：如果同一个项目中多个样本是并行运行的，STAR比对时一定要让不同样本在独立的临时目录下运行，或者确保--outFileNamePrefix的样本名前缀唯一，否则会产生文件覆盖或者写冲突。跑几十个样本的时候，这种低级错误很容易在脚本里悄悄出现，最后比对结果文件的大小都不对，排查起来很费劲。
最后的最后，回到方法论层面说一句。Bulk RNA-seq分析流程本身已经相当标准化了，真正的差异往往来自对每个环节原理的理解和对细节的把握——从实验设计时的重复数选择，到链特异性参数的确认，到差异基因筛选阈值的依据。把这些问题在动手之前想清楚，远比多跑几个工具更重要。
