小鼠Bulk RNA-seq全流程实操指南：从实验设计到差异表达分析

做小鼠的 Bulk RNA-seq最怕的不是不会跑 pipeline而是跑完了发现实验设计有问题或者中间某个环节埋了雷最后样本全废哭着回来补做。我自己最早入坑生信就是从小鼠转录组开始的那时候一边看教程一边手动敲命令质控、比对、定量、差异分析一步步走下来踩过的坑比吃过的盐还多。所以这篇指南我想直接对着完整流程讲一遍从拿到测序数据开始到最终拿到差异基因和富集结果每一步该做什么、为什么这么做、有哪些参数不能随便动全部按我实际操作的习惯来写。这篇内容适合几类人一是刚接触 Bulk RNA-seq想快速跑通全流程的同学二是湿实验出身、数据已经送测但不知道怎么分析的研究生三是已经跑过流程但结果总是不对劲想回来排查问题的朋友。我不会只丢给你一串命令而是把每一步背后的逻辑拆开讲清楚读完你至少能明白自己在做什么而不是机械地 copy 代码。1. 实验设计先行不给下游分析埋雷很多人拿到测序数据就开始跑流程结果到差异分析时发现样本分组混乱、重复数不足、批次效应严重根本没法补救。Bulk RNA-seq 的下游分析看起来是生信问题但绝大多数致命问题都出在实验设计阶段。所以第一步不是装软件而是先把设计捋清楚。1.1 分组设计中最容易被忽略的三件事第一件是生物学重复的数量。我见过不少同学为了省钱每个组只做 2 个重复甚至有人只做 1 个问就是“差异应该很大不需要重复”。这种想法会直接卡死在差异分析那一步。Bulk RNA-seq 的差异检验依赖组内方差估计重复数少于 3edgeR 和 DESeq2 根本算不出可靠的离散度结果就是 P 值大得离谱或者反过来因为过拟合给出完全不可信的显著基因。我的建议是常规处理组和对照组各 3 到 5 个生物学重复如果预期效应量小至少 5 个。这里的成本不是浪费是在买统计功效。第二件是批次效应的预分配。如果样本不是同一天提 RNA、不是同一次建库、不是同一张 flow cell 测序就会引入批次差异。最糟糕的情况是对照组全在第一批次、处理组全在第二批次那下游无论用什么算法都很难把处理效应和批次效应拆开。正确的做法是在实验设计时就把不同组的样本打散分配到不同批次里让批次与处理因素尽量正交。这样即使有批次差异也能在后续分析中用模型矫正。第三件是混杂因素的设计。比如处理组全是雄性、对照组全是雌性那差异基因到底是处理造成的还是性别造成的没人说得清。类似的问题还包括体重、年龄、饲养笼等。能用协变量在模型里矫正的前提是你在设计时记录下来了如果根本没记录那分析阶段很难补。小鼠实验还有一个容易被忽略的细节同笼饲养的小鼠有“笼效应”建议把同一巢的个体也尽可能分到不同组。1.2 建库策略与测序参数选择的底层逻辑Bulk RNA-seq 有两种主流建库方式链特异性文库stranded和非链特异性文库unstranded。链特异性文库保留了转录本的方向信息在判断基因反义转录、区分同源基因时非常重要。大部分公司默认给的就是链特异性建库但你在对比对参数和定量参数时一定要确认这一点。featureCounts 里有-s参数链特异性文库设为 1 或 2非链特异性文库设为 0这个参数如果设错定量结果会很奇怪。测序深度也是一个需要提前想清楚的问题。如果目标是找差异表达基因每个样本 10M 到 20M reads 基本够用如果要看低表达基因、可变剪接或者融合基因推荐 30M 以上。我个人的习惯是常规转录组测序设置为每个样本 20M pair-end reads读长为 PE150。这个配置在成本和结果稳定性之间比较平衡。还有一个很多人忽略的问题文库片段大小选择。常规转录组一般选择 250bp 左右的插入片段但如果你下游要做转录本异构体分析插入片段太短会导致剪接事件信息不足。这一点最好在建库前想清楚因为测序结束后无法补救。2. 上游分析从原始数据到表达矩阵拿到下机数据之后一般会是一堆.fastq.gz文件也可能公司已经帮你做了初步质控。上游分析的目标非常明确把每个样本的 reads 转化为一个基因表达矩阵。这中间包括质控过滤、参考基因组比对、基因水平定量三大步每一步都有不少细节需要处理。2.1 原始数据质控与过滤fastp 的参数取舍我习惯用 fastp 做一步式的 QC 过滤因为它会把接头识别、质量过滤、长度过滤、单端/双端合并处理都包在一个命令里速度比自己串 fastqc cutadapt trimmomatic 快很多。对于小鼠转录组的常规 PE150 数据我常用的命令大致长这样fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \ -o sample_R1.clean.fastq.gz -O sample_R2.clean.fastq.gz \ --detect_adapter_for_pe \ --cut_front --cut_tail \ --cut_front_mean_quality 20 \ --cut_tail_mean_quality 20 \ --qualified_quality_phred 20 \ --length_required 50 \ --thread 16 \ --json sample.fastp.json --html sample.fastp.html参数看起来多但逻辑很清晰--detect_adapter_for_pe让 fastp 自动识别接头序列并切除双端测序的接头识别一般比手动指定更稳--cut_front和--cut_tail分别切除两端质量偏低的位置--qualified_quality_phred 20表示 Q20 标准--length_required 50过滤掉长度小于 50bp 的 reads。这个标准不是特别严格适合转录组分析因为过度过滤反而可能扔掉有用的短转录本信息。如果你第一次接触 fastp建议跑完后打开.html报告检查几个核心指标Q20/Q30 比例是否达到 90% 以上、reads 长度分布是否均匀、duplication rate 是否过高、是否有残余接头。我遇到过一个问题某次建库质量不好R2 端整体质量明显低于 R1这时候就需要单独注意 R2 的质量过滤否则比对率会断崖式下降。质控后我还习惯把所有样本的清理报告汇总一下用 MultiQC 直接读取 fastp 的 json 输出生成一个总览页面。这一步虽然没有技术含量但在审稿时需要展示质控指标或者在组内讨论数据质量时非常方便。2.2 参考基因组与注释文件的选择小鼠参考基因组和注释文件看起来只是“下载一个文件”实际上这里踩坑的人非常多。最经典的问题是版本不匹配比如基因组用 Ensembl 的 GRCm38注释却用了 UCSC 的 mm10或者基因组用了 GRCm38注释用了 GRCm39比对和定量的时候程序不会直接报错但结果里会出现大量基因无法正确注释甚至 reads 定位到错误的基因区间。我目前推荐的做法是去 Ensembl 官网下载对应的 primary assembly 基因组文件和 GTF 注释文件。以小鼠为例Ensembl 当前主推的版本是 GRCm39下载时确认Mus_musculus.GRCm39.dna.primary_assembly.fa.gz和Mus_musculus.GRCm39.110.gtf.gz版本号一致。注意110是 Ensembl 的 release 号不同 release 对应的注释内容有差异尽量用较新的 release但一旦开始分析就不要中途更换版本否则所有样本需要重新比对。还有一个容易被忽略的坑基因组文件里的染色体命名方式。Ensembl 的染色体名称是1, 2, 3 ... X, Y, MT而 UCSC 是chr1, chr2 ...。如果比对索引和 GTF 文件来源不一致STAR 在生成索引时虽然不会报错但比对结果中来自线粒体基因组的 reads 可能无法正确分配到基因注释上。建议全程使用同一数据源并在比对前用grep -c 基因组文件检查一下染色体名称的格式。2.3 比对工具的选择STAR、Hisat2 与转录本定量的差别上游比对工具有很多选择目前在 Bulk RNA-seq 领域用得最多的三个是 STAR、Hisat2 和直接做转录本定量的 Salmon。选哪个工具取决于你的下游需求。如果你要做基因水平差异表达分析加 Spliced reads 的准确定位很关键那我推荐 STAR。它的速度和准确性都很优秀尤其是处理长 reads 和复杂剪接事件时表现稳定。STAR 的标准流程分两步第一步生成基因组索引第二步比对。生成索引的命令大致如下STAR --runMode genomeGenerate \ --genomeDir /path/to/star_index \ --genomeFastaFiles /path/to/Mus_musculus.GRCm39.dna.primary_assembly.fa \ --sjdbGTFfile /path/to/Mus_musculus.GRCm39.110.gtf \ --sjdbOverhang 149 \ --runThreadN 16这里--sjdbOverhang设置为读长减 1PE150 就填 149这是 STAR 手册里明确建议的参数直接决定剪接位点检测的灵敏度。比对样本的命令如下STAR --genomeDir /path/to/star_index \ --readFilesIn sample_R1.clean.fastq.gz sample_R2.clean.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outSAMstrandField intronMotif \ --outFileNamePrefix sample_ \ --runThreadN 16--outSAMtype BAM SortedByCoordinate直接输出按坐标排序的 BAM 文件省去后续排序步骤。--outSAMstrandField intronMotif会为每条记录加一个链方向标签虽然主流定量工具不一定需要这个字段但如果你之后想用 IGV 查看或者跑一些需要链信息的工具这个参数能省不少事。Hisat2 的优势是内存占用更低、速度更快适合在配置一般的服务器上跑但它的定位准确性还是比 STAR 稍弱一些。而 Salmon 这类工具干脆不做基因组比对而是把 reads 直接映射到转录本序列上速度极快、内存占用小在估算转录本丰度和 TPM 值时非常合适。但如果你需要拿到 BAM 文件做可视化、需要识别新剪接事件那还是要走 STAR。我的建议是常规差异表达分析直接用 STAR 比对 featureCounts 定量这条路径最稳妥、最容易被审稿人接受。如果你只想快速算个表达量、不想做剪接分析那 Salmon 也可以但建议先把比对这关走扎实。2.4 featureCounts 基因定量参数细节与常见坑拿到排序后的 BAM 文件下一步就是基因水平定量。我最常用的工具是 featureCounts来自 Subread 包。核心命令featureCounts \ -a /path/to/Mus_musculus.GRCm39.110.gtf \ -o counts.txt \ -T 8 \ -p \ --countReadPairs \ -s 1 \ -t exon \ -g gene_id \ sample1.bam sample2.bam ... sampleN.bam我重点说几个容易设置错误的参数。-p表示双端测序早期版本中 featureCounts 对双端数据的默认行为是只统计 read pairs 中的一条也就是说最终计数约等于片段数。设置--countReadPairs后它会统计整个片段但要注意不同版本的默认行为可能有差异跑完看一下counts.txt头部的 Summary 行确认统计值是不是接近预期。-s表示链特异性1 表示正向链特异性文库2 表示反向链特异性文库0 表示非链特异性。这个值一定要和建库方式对应上否则基因计数会丢掉一大半。-t exon -g gene_id的意思是统计落在 exon 特征上的 reads并按gene_id汇总到基因水平。这套组合是最标准的操作因为基因表达定量本质上是把所有外显子区域的 reads 数加起来。这里有个常见误区如果 GTF 注释版本和比对使用的索引版本不一致featureCounts 通常不会报错但结果中会有很多基因计数为 0或者某些 reads 明明比对上却没有被分配到任何基因这时候优先检查版本一致性。featureCounts 输出的counts.txt包含多列信息真正需要的是各样本的计数列。我一般直接用 R 读取这个文件去掉前 6 列注释信息后用DESeq2或edgeR构建矩阵。这里还有个技巧在跑 featureCounts 之前用samtools flagstat看一眼 BAM 文件的比对率和 rRNA 比例如果 rRNA 比例超过 5%说明建库时 rRNA 去除不彻底下游分析前要考虑额外过滤。3. 下游差异表达与功能富集表达矩阵拿到手后好戏才刚开始。从这里开始涉及统计模型和生物学解释信息密度也一下子高起来。差异表达分析不只是一个函数调用背后的归一化逻辑、过滤逻辑和多重检验矫正逻辑都需要心里有数。3.1 差异基因分析edgeR 和 DESeq2 怎么选目前比较主流的小鼠 Bulk RNA-seq 差异分析工具是 edgeR 和 DESeq2两个都是基于负二项分布模型。它们的核心思想类似用所有基因的计数数据估计每个基因的离散度再对基因表达量做假设检验。但它们的归一化方法和对小样本的处理方式不同所以结果会有细微差异。如果是常规 3 对 3 的实验设计我一般首选 edgeR TMM 归一化。TMM 归一化对文库组成差异的矫正比较稳健计算量也小。核心流程如下library(edgeR) counts - read.delim(counts.txt, row.names1, headerTRUE) group - factor(c(ctrl,ctrl,ctrl,treat,treat,treat)) y - DGEList(countscounts, groupgroup) keep - filterByExpr(y, groupgroup) y - y[keep,,keep.lib.sizesFALSE] y - calcNormFactors(y) design - model.matrix(~group) y - estimateDisp(y, design) fit - glmQLFit(y, design) res - glmQLFTest(fit, coef2) topTags(res, nInf)filterByExpr这一步可以自动过滤掉低表达基因保留的基因数通常是 1.5 万到 2 万。这个过滤非常重要一来能减少多重检验矫正的负担二来能提高离散度的估计稳定性。很多人上来就直接跑差异检验不做过滤结果就是大量低表达基因给出了非常不可靠的高显著 P 值。DESeq2 的优势在于它使用 shrinkage 方法对离散度进行估计对样本量小、组内方差大的情况更稳。如果你只有 2 个重复或者明显有离群样本DESeq2 通常比 edgeR 更稳。DESeq2 的典型流程是library(DESeq2) countData - read.delim(counts.txt, row.names1, headerTRUE) colData - data.frame(condition factor(c(ctrl,ctrl,ctrl,treat,treat,treat))) dds - DESeqDataSetFromMatrix(countDatacountData, colDatacolData, design~condition) dds - DESeq(dds) res - results(dds, contrastc(condition,treat,ctrl))这里有个小坑DESeq2 要求 count 数据是整数如果 featureCounts 输出的矩阵是整数的那就没问题。另外DESeq2 会做独立过滤所以results()返回的 padj 已经自动调整了较低表达基因的过滤不需要手动再走一遍 filterByExpr。我个人实验里两个工具跑出来的显著基因集合通常有 70% 以上的重叠差异主要出现在低表达基因和边缘显著的基因上。审稿时两个工具的结果可以互为验证但最终报告里选一个作为主要结果就行另一个放补充材料。3.2 功能富集分析GO、KEGG 与 GSEA 的思路差异拿到差异基因列表后下一步就是功能富集分析。这一步看似简单但很多人把 GO/KEGG 富集和 GSEA 混为一谈导致结果解读出现偏差。GO 和 KEGG 富集分析属于“过表达分析”输入是一组显著差异基因检验这些基因是否在某些功能通路中富集。常用的工具是 R 包 clusterProfiler。对于小鼠数据需要先加载org.Mm.eg.db包。我的大致流程是library(clusterProfiler) library(org.Mm.eg.db) deg - rownames(subset(res, padj 0.05 abs(log2FoldChange) 1)) ego - enrichGO(gene deg, OrgDb org.Mm.eg.db, keyType ENTREZID, ont BP, pAdjustMethod BH, qvalueCutoff 0.05) ekegg - enrichKEGG(gene deg, organism mmu, keyType kegg, pvalueCutoff 0.05)这里有两个非常重要的小细节。第一如果差异基因列表用的是 Entrez IDKEGG 的organism参数在小鼠里要写成mmu不是mouse写错了会报错。第二GO 富集的ont参数有三个选项BP生物学过程、CC细胞组分、MF分子功能。很多人只跑 BP但其实 CC 和 MF 有时能提供非常重要的线索比如差异基因集中定位在某个细胞器中这本身就暗示了功能方向。但过表达分析有一个天然缺陷它只关注显著差异基因把表达量变化幅度不大但很稳定的基因全部忽略掉了。GSEA 则是用所有基因的表达水平排序然后看某一个基因集合在排序中是否显著富集在顶部或底部。这个思路更接近真实生物学很多调控基因表达变化可能不大但它们所在通路的整体活性已经发生偏移。我常用 fgsea 包跑 GSEA。核心思路是先用所有基因的 log2FC 做一个排序列表然后对每个 KEGG 或 Hallmark 基因集做富集分析。这里注意两点一是基因排序指标用 signal-to-noise 或 log2FC 都行但如果是小样本实验log2FC 更容易受个别极端值影响可以配合shrinkage后的 log2FC 使用二是 GSEA 需要的是基因的数值型排序向量而不是显著的二进制列表这是它和 GO/KEGG 富集最大的区别。3.3 可视化要点的实操建议差异分析和富集分析的最终呈现通常依赖几张图PCA 图、火山图、热图和 GSEA 富集图。PCA 图用来展示样本的整体关系一般用所有基因的表达量做 PCA看组间是否分离、组内是否聚集。如果组内样本分散严重就要警惕批次效应。火山图最常用EnhancedVolcano包横轴是 log2FC纵轴是-log10(padj)两条虚线分别代表 fold change 阈值和显著性阈值。热图我一般用pheatmap对显著差异基因做 z-score 标准化后展示。如果基因数太多可以先按方差筛选前 1000 到 2000 个基因再画否则热图密密麻麻根本没法看。GSEA 结果图形一般用enrichplot::gseaplot2画能同时展示富集分数曲线、基因排序位置和核心基因的热图。这里要注意的是GSEA 图里的核心基因集leading edge才是这个通路里真正贡献富集信号的基因写文章时把这部分基因提出来重点讨论会更有说服力。4. 常见问题排查与避坑实录跑 Bulk RNA-seq 流程最容易被卡住的地方其实不是命令本身而是结果异常时无从下手。我把自己这几年遇到的高频问题整理成一份排查清单基本覆盖了从原始数据到下游差异分析的大部分坑。4.1 比对率低或 rRNA 污染严重怎么处理正常小鼠转录组数据的比对率应该在 85% 到 95% 之间。如果比对率明显低于这个范围先看 MultiQC 或 fastp 报告里的 GC 含量和 duplication 情况再检查是否有过多 rRNA reads。如果你的比对率在 70% 左右可以先用 Bowtie 2 把 reads 比对到 rRNA 序列上看看是不是建库时 rRNA 去除不彻底。如果是 rRNA 污染建议找公司补做 rRNA depletion 或用商品化的 rRNA 去除试剂盒重新建库单纯靠生信过滤很难彻底挽救。如果只是小幅偏低比如 80% 到 85%可以先检查 STAR 比对时是否因为基因组注释文件版本不一致导致大量 reads 落在“未注释区域”这种情况可以通过samtools view看 BAM 文件中比对到线粒体基因组的比例来辅助判断。线粒体基因比例高其实是正常的通常占 5% 到 20%太高则说明建库时线粒体 rRNA 没有去除干净。4.2 批次效应什么时候用 ComBat-seq什么时候不能用很多时候你跑完 PCA 发现样本没有按处理分组聚在一起而是按测序批次聚在一起这就是典型的批次效应。如果实验设计阶段已经把组别均匀分配在不同批次里这时候可以用sva包的 ComBat-seq 函数对 count 数据做批次矫正。ComBat-seq 的优势是可以直接作用于 count 矩阵不破坏整数属性且比 ComBat 更适合转录组数据。但 ComBat-seq 不是万能药。如果批次和处理因素完全混杂比如对照组全在第一批次、处理组全在第二批次这时候任何批次矫正算法都很难有效甚至会把真实的处理效应也一并矫正掉。这种情况最好的策略是回到实验设计阶段重新补样本或者尽可能增加对照组和处理组在批次间的交叉。如果已经不可补救那只能在文章中如实说明批次混杂的限制并谨慎解释差异分析结果。4.3 生物学重复不足时的补救策略如果重复数真的不够比如每组只有 2 个样本不是完全不能分析但要非常谨慎。edgeR 和 DESeq2 在重复数为 2 时都能跑但检验功效极低很容易漏掉真实差异。我有一招应对这种情况先跑 DESeq2用betaPriorTRUE或默认的收缩估计获取 log2FC把显著性阈值放宽一些优先关注那些在两个重复中方向一致且表达变化幅度大的基因。这种基因即使 P 值不显著也值得在后续 qPCR 或 Western blot 中验证。另外一个思路是合并处理效应相近的样本组。比如你有一个低剂量处理组和一个高剂量处理组如果各只有 2 个重复可以考虑先用趋势检验如contrast设计为线性趋势检验剂量效应而不是两两比较。这样能在不增加样本的情况下用上所有数据的信息。4.4 内存与计算资源不够怎么办跑 STAR 生成基因组索引是最吃内存的一步GRCm39 全基因组索引大约需要 30GB 内存如果你的服务器内存不够有两个办法一是在生成索引时加--genomeSAindexNbases 13或 12降低索引占用二是直接用 Salmon 这类不依赖基因组长索引的比对定量工具速度更快、内存占用小很多但代价是拿不到 BAM 文件做剪接分析。对于常规差异分析Salmon 的结果完全够用。如果机器只有 8GB 内存我建议走 salmon tximport 这个组合时间大概只占 STAR 路径的三分之一。用 tximport 把转录本水平定量汇总到基因水平下游的 edgeR/DESeq2 流程完全不受影响。5. 全流程的完整命令串联建议这里我把自己的日常流程做一个完整串联方便你照着跑通一遍。假设测序数据已经在raw_data/目录下文件名是S1_R1.fastq.gz和S1_R2.fastq.gz。第一步对于每个样本运行 fastpfor s in S1 S2 S3 C1 C2 C3; do fastp -i raw_data/${s}_R1.fastq.gz -I raw_data/${s}_R2.fastq.gz \ -o clean_data/${s}_R1.fastq.gz -O clean_data/${s}_R2.fastq.gz \ --detect_adapter_for_pe --cut_front --cut_tail \ --qualified_quality_phred 20 --length_required 50 --thread 16 \ --json clean_data/${s}.json --html clean_data/${s}.html done第二步生成 STAR 索引STAR --runMode genomeGenerate \ --genomeDir star_index \ --genomeFastaFiles ref/Mus_musculus.GRCm39.dna.primary_assembly.fa \ --sjdbGTFfile ref/Mus_musculus.GRCm39.110.gtf \ --sjdbOverhang 149 --runThreadN 16第三步比对所有样本for s in S1 S2 S3 C1 C2 C3; do STAR --genomeDir star_index \ --readFilesIn clean_data/${s}_R1.fastq.gz clean_data/${s}_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outSAMstrandField intronMotif \ --outFileNamePrefix bam/${s}_ --runThreadN 16 done第四步featureCounts 定量featureCounts -a ref/Mus_musculus.GRCm39.110.gtf \ -o counts.txt -T 8 -p --countReadPairs -s 1 \ -t exon -g gene_id \ bam/S1_Aligned.sortedByCoord.out.bam bam/S2_Aligned.sortedByCoord.out.bam \ bam/S3_Aligned.sortedByCoord.out.bam bam/C1_Aligned.sortedByCoord.out.bam \ bam/C2_Aligned.sortedByCoord.out.bam bam/C3_Aligned.sortedByCoord.out.bam第五步把counts.txt导入 R跑 edgeR 或 DESeq2。具体代码前面已经写过这里不再重复。持续跑到这一步其实整个流程已经从原始数据到了差异表达列表后面的富集分析只是调用函数的问题。我个人的习惯是每一个关键步骤结束后都记录一下输出文件的 checksum、版本信息和参数配置。生物信息学分析不是跑完就结束很多结果需要回溯版本和参数信息对后续排查和写方法部分非常重要。即使没有专门做 project management也至少把命令存成一个 shell 脚本放进项目目录方便之后复现。从零开始跑完小鼠 Bulk RNA-seq 全流程其实花不了太长时间真正花时间的是排查那些看起来玄学的异常。我第一次跑通整个流程用了整整两周其中一半时间都花在版本不匹配和参数设置上。后来踩的坑多了才慢慢总结出这套相对稳健的路线。希望这份指南能帮你少走一些弯路如果真的遇到了我上面写过的坑那你直接对照着排查就行。