Bulk RNA-seq全流程指南：从测序数据到差异表达可视化

我最早接触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 readsSTAR会同时找到两个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和转录组序列做定量速度比STARfeatureCounts的组合快一个量级。它适合什么场景一种是样本量特别大上百个样本比对到基因组再定量的计算成本过高另一种是你只关心转录本或基因表达水平不需要处理变异信息。但Salmon有个前提要注意它基于参考转录组进行定量参考转录组的完整度直接决定结果质量。如果某个基因的可变剪接异构体没有包含在参考里它的表达量就可能被低估或者错配。另外Salmon定量后再转成count矩阵做差异分析还需要用tximport来汇总到基因水平。流程上多了一步但并没有多复杂。从实操角度看如果项目只有几个样本、几十个样本STARfeatureCounts仍然是更稳妥的选择因为它保留了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与shrinkageresults()返回的结果里有几列必须搞清楚baseMean是基因在所有样本中的平均归一化countslog2FoldChange是处理组相对于对照组的倍数变化取log2lfcSE是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. 富集分析把基因列表变回生物学意义差异基因列表本身只是基因名真正的生物学故事要靠富集分析来讲。最常用的是GOGene 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轴是log2FCy轴是-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分析流程本身已经相当标准化了真正的差异往往来自对每个环节原理的理解和对细节的把握——从实验设计时的重复数选择到链特异性参数的确认到差异基因筛选阈值的依据。把这些问题在动手之前想清楚远比多跑几个工具更重要。