★ FEATURED ARTICLE RNA-seq大概是高通量测序数据分析里最“亲切”的入口了。很多人一听到“高通量测序数据分析”这几个字第一反应是命令行、服务器、成堆的FASTQ文件感觉门槛高得吓人。但真正跑过一轮RNA-seq的标准分析之后你会发现它的流程其实非常固定从原始测序数据到表达矩阵再到差异基因和功能富集每一步都有成熟的工具和约定俗成的参数。与其说它是高深的算法研究不如说是一套需要细心和耐心的“数据加工流水线”。我早期带团队的时候经常有新人拿着一个RNA-seq项目来问第一步该干什么。这个问题听起来简单但如果没人给你梳理全貌很容易一头扎进细节里出不来——比如纠结某个比对工具的某个参数到底该设成多少结果跑到最后才发现参考基因组版本和注释文件根本对不上。所以这篇文章我不打算只讲某个单一工具的操作而是把我平时处理RNA-seq项目的一套完整流程、每个环节的工具选型理由、需要避开的坑、以及各类报错的排查思路整理出来。不管你是刚接触高通量测序数据分析的学生还是已经跑过一两个项目但想系统梳理一遍的从业者希望这篇内容能帮你把整个RNA-seq分析的框架搭起来。1. 先建立RNA-seq分析的整体框架1.1 为什么RNA-seq的流程如此“标准化”RNA-seq的核心目标是检测转录组层面的基因表达情况。简单来说就是把细胞或组织里的RNA反转录成cDNA然后用高通量测序仪读出大量短序列reads再通过生物信息学方法把这些序列比对回参考基因组或转录组统计每个基因对应了多少reads最终得到一张“基因表达量表”。这听起来很直接但实际操作里每一步都有讲究。比如reads里混入了测序接头序列怎么办测序质量差的碱基会不会影响比对结果比对时reads跨越了剪接位点该怎么处理不同样本测序深度不一样怎么公平比较……这些问题都有对应的工具和策略于是也就形成了一套相对固定的标准流程。这套标准流程的骨架大致是原始数据质控 → 比对或定量 → 表达定量 → 差异表达分析 → 功能富集分析。无论是人的组织样本、肿瘤样本还是模式动物的转录组数据前四步几乎都是一样的只有到功能富集环节才会根据研究目的产生分化。这也是为什么我总建议新手先把这个骨架背下来——一旦你清楚自己正处于哪个环节每个环节该用什么工具、会得到什么文件整个项目就不会乱。1.2 两条主流技术路线的选择目前RNA-seq的数据分析有两条主流路线。第一条是传统路线先将reads比对到参考基因组上再基于比对结果统计每个基因的reads数。代表工具是STAR和HISAT2定量环节一般用featureCounts或HTSeq-count。第二条是伪比对路线不把reads逐一比对到基因组上而是直接把reads比对到转录组序列集合上通过快速算法估算转录本丰度代表工具是Salmon、kallisto和RSEM。这两条路线各有优劣。传统路线的优势是结果直观、可控性强后续断点识别、融合基因检测等都可以复用比对结果缺点是对计算资源要求较高跑到STAR比对这一步时内存经常被吃满。伪比对路线的优势是速度快、资源占用小Salmon比对一个样本往往只需几分钟缺点是它本质上是转录本水平的定量某些情况下对于新转录本或复杂基因结构的处理不如基因组比对路线灵活。我自己在实际项目中怎么选呢如果项目以差异表达分析为主要目标且样本量大我推荐直接用Salmon定量之后再在R里面用tximport整合。如果项目不仅要做表达定量还要看可变剪接、融合基因或者需要出BAM文件做IGV可视化那就老老实实走STAR featureCounts的路线。两种方案最后都能得到表达矩阵但要注意不同定量工具得到的counts数据不能混用这一点一定要在项目开始时确定好。我个人的默认方案是STAR featureCounts DESeq2组合。原因很朴素社区资料多、报错好查、审稿人认可度高。毕竟做科研项目稳定性比炫技重要得多选一套大多数人都验证过的流程能省掉大量踩坑时间。2. 从环境搭建到数据准备2.1 conda管理分析环境RNA-seq分析涉及的工具非常多STAR、fastp、samtools、featureCounts、R、DESeq2等等每个工具又有自己的依赖库。如果一股脑全装在系统里版本冲突是迟早的事。我强烈建议从一开始就用conda或者mamba建一个独立分析环境把所有工具装在一起互不干扰。安装环境的命令很简单# 建议用mamba比conda快很多 mamba create -n rnaseq -c bioconda -c conda-forge \ fastp star samtools subread multiqc r-base # 激活环境 conda activate rnaseq有一点需要提醒用conda安装R之后之后用R包时可能会遇到系统依赖缺失的问题最常见的包括libcurl、libxml2、openssl。如果你用的是Windows最好装WSL之后再跑流程macOS的M系列芯片在安装某些生物信息学软件时偶发兼容性问题遇到的话优先查conda-forge源。2.2 参考基因组与注释文件的选择很多人会忽略这一步但恰恰是这一步决定了后续所有分析的质量。参考基因组版本不统一或者注释文件和参考基因组版本不匹配是RNA-seq分析中最常见的低级错误而且一旦到了比对之后才被发现代价非常大。我的建议是人类数据优先使用GENCODE版本的GTF注释文件因为它的基因注释质量比Ensembl或UCSC更适合转录组而且版本号与Ensembl对齐。下载参考基因组时文件命名里带上版本号比如Homo_sapiens.GRCh38.dna.primary_assembly.fa和Homo_sapiens.GRCh38.108.chr.gtf然后在一个genome/目录里统一管理。一个GTF文件只用在一个版本上。GRCh37的GTF配上GRCh38的基因组比对率虽然不会崩但后续定量出来的基因会非常奇怪。在拿到FASTA和GTF之后先做一个索引完整性检查。用grep -c ^看FASTA里有多少条染色体用grep -v ^#统计GTF里有多少个转录本和基因做到心中有数再进入比对环节。2.3 建立STAR索引STAR比对的第一步是给参考基因组建索引这一步只需要做一次后续所有同版本样本都能复用。STAR建索引时特别吃内存人类基因组建议分配30GB以上内存否则容易中途被杀掉。具体命令如下STAR --runMode genomeGenerate \ --genomeDir /path/to/star_index_hg38 \ --genomeFastaFiles /path/to/Homo_sapiens.GRCh38.dna.primary_assembly.fa \ --sjdbGTFfile /path/to/Homo_sapiens.GRCh38.108.chr.gtf \ --sjdbOverhang 149 \ --runThreadN 16这里--sjdbOverhang的取值通常设为读长减1。比如测序读长是PE150那就设成149。这个参数影响剪接位点在索引中的表现形式虽然不是非设不可但设置正确能让比对更准确地处理跨越剪接点的reads。如果你的读长不是150bp比如有些平台是100bp记得改成99。3. 原始数据质控别急着上手干活先看清楚数据长什么样3.1 FASTQ格式基础与质量值拿到测序下机的原始数据一般是*.fastq.gz文件。FASTQ文件每四行为一个单位第一行是序列标识符第二行是碱基序列第三行是加号分隔符第四行是对应每个碱基的质量值字符串。质量值是用ASCII码表示Phred分数具体公式是Q -10 * log10(P)其中P代表该碱基被测序错误的概率。Q20表示错误率1%Q30表示错误率0.1%。现在的主流平台下机数据Q30通常在85%以上如果低于这个水平就要留意是不是测序过程有问题。拿到FASTQ文件后千万不要直接开始比对。先跑一遍fastp或者MultiQC看整体质量这一步花不了几分钟但能帮你提前发现大量问题。我看过一个真实案例某样本的reads里接头残留率高达30%几乎可以确定是建库环节出了问题这种情况如果不做质控直接比对最终差异分析结果里会出现大量假阳性。3.2 fastp过滤参数怎么设fastp是目前最常用的质控工具它把质控报告、接头去除、低质量碱基过滤、polyG尾裁剪整合在一起一条命令就能完成。我常用的命令是这样fastp \ -i sample_R1.fastq.gz \ -I sample_R2.fastq.gz \ -o sample_clean_R1.fastq.gz \ -O sample_clean_R2.fastq.gz \ --detect_adapter_for_pe \ --cut_front \ --cut_tail \ --cut_front_window_size 1 \ --cut_front_mean_quality 20 \ --cut_tail_window_size 1 \ --cut_tail_mean_quality 20 \ --length_required 36 \ --thread 16 \ --html sample_fastp.html \ --json sample_fastp.json其中--cut_front和--cut_tail分别是去除5’端和3’端低质量碱基窗口大小为1意味着逐个碱基判断质量如果质量低于Q20就剪掉。--length_required 36表示如果read被剪成不足36bp就直接丢弃。--detect_adapter_for_pe用于自动检测并剪切双端测序的接头序列。有人会纠结参数要不要开得太猛。我的建议是对于常规RNA-seq质量过滤的标准不必过高尤其是不要为了追求高clean率把过短的reads都留着。短reads在STAR比对里很大概率比对不上或比对到错误位置属于“留着添乱”的类型。过滤之后用grep -c ^SRR或者fastp自带的报告查看剩下的reads数量心里就有底了。3.3 质控报告怎么看质控报告的解读有一个关键顺序先看clean reads占总reads的比例通常在90%以上是正常的低于85%就要从建库或测序端找原因。再看接头残留率严格说应在1%以下。然后是GC含量分布。mRNA-seq在建库时经过oligo-dT富集GC分布会有一定偏向但如果出现完全偏离正常范围的双峰分布或尖峰说明可能存在污染或者PCR扩增偏差。最后看duplication水平。RNA-seq本身的duplication会比较明显尤其对于高表达基因但如果整体duplication超过70%且不是由于测序深度过低导致很可能建库时扩增循环数太多需要留意下游定量是否存在偏倚。做了多个样本的项目建议把fastp的结果都收集起来最后统一跑一次MultiQC汇总这样所有样本的质量指标就能在一份报告里横向对比一眼就能发现哪个样本异常。4. 比对与定量从FASTQ到表达矩阵4.1 STAR比对的核心参数STAR比对RNA-seq数据是目前的黄金标准。它最出色的能力是处理跨越剪接位点的reads通过它可以找到reads的剪接位点并正确比对到基因组上。其核心原理是“种子搜索扩展”先用高质量连续匹配的种子序列快速定位候选位置再在候选区域做延伸打分最终输出最优比对结果。我常用的STAR比对命令如下STAR --genomeDir /path/to/star_index_hg38 \ --readFilesIn sample_clean_R1.fastq.gz sample_clean_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outSAMunmapped Within \ --outFileNamePrefix sample_ \ --outFilterMismatchNmax 10 \ --outFilterMultimapNmax 10 \ --outFilterScoreMinOverLread 0.3 \ --outFilterMatchNminOverLread 0.3 \ --runThreadN 16--outSAMtype BAM SortedByCoordinate直接输出坐标排序后的BAM文件省去下游samtools sort的步骤。--outFilterMultimapNmax 10表示允许输出最多10个比对位置的多比对reads对于RNA-seq来说这些通常来自同源基因家族建议保留一份在后续分析中按默认处理方式排除。--outFilterMismatchNmax 10是允许的最大错配数如果测序质量很好把这个值调低到6-8能显著减少错误比对。比对耗时和资源同样需要关注。STAR单样本人类转录组16线程、PE150数据通常30分钟能跑完2亿条reads内存消耗一般在25-30GB左右。如果你的服务器内存只有16G建议换用Salmon或Hisat2。跑STAR之前先ulimit -n 65535或者用--limitBAMsortRAM适当提高内存上限防止BAM排序阶段崩掉。4.2 featureCounts定量的细节STAR比对得到的BAM文件下一步就是用featureCounts统计每个基因的reads数。featureCounts是subread软件包里的组件速度快内存占用少是目前RNA-seq定量的首选。这里有一个非常重要的参数细节如果你在比对时使用了--quantMode GeneCountsSTAR本身就能输出一个简单的基因counts文件但那个计数的规则比较粗糙对双端数据和多外显子基因的处理不如featureCounts细致。所以我还是建议用featureCounts单独跑一遍。featureCounts -a /path/to/Homo_sapiens.GRCh38.108.chr.gtf \ -o sample_counts.txt \ -T 16 \ -p \ --countReadPairs \ -s 0 \ sample_Aligned.sortedByCoord.out.bam-p表示输入的是双端数据--countReadPairs表示以read pair为单位计数即一对reads比对成功算1个计数如果不加这个参数featureCounts会按单端方式对每条read单独计数双端数据就会翻倍后续差异分析会出大问题。-s 0表示非链特异性建库。如果你的文库是链特异性需要根据实际建库方式设置-s 1或-s 2这个信息务必在项目开始时向建库方确认。featureCounts输出的文件会附带gene length、start、end等注释信息实际用于差异分析的就是最后一列counts数据。你可以用自己的脚本提取基因名和counts两列也可以用R的read.delim直接读取。4.3 counts矩阵生成以后所有样本跑完featureCounts后就会得到一张原始counts矩阵。行是基因列是样本单元格是reads数。这张矩阵是后续所有统计分析的输入但要注意这时的数据还“不干净”因为不同样本的测序深度不一样高表达基因和长基因天然会被测到更多reads所以直接拿原始counts做比较是错的。也不要直接用TPM或FPKM做差异分析。TPM和FPKM做了长度归一化它的设计初衷是用于基因表达水平的相对比较而DESeq2的统计模型需要的是整数型的原始counts。这一点我见过太多人弄错了——拿FPKM跑DESeq2跑完自己都不知道结果为什么那么离谱。所以正确的输入数据是featureCounts输出的原始counts矩阵。归一化、统计学检验都交给DESeq2内部处理你只需要保证数据格式正确。5. 差异表达分析理解DESeq2背后的原理5.1 为什么DESeq2要“原始整数counts”DESeq2是目前差异表达分析中使用最广的工具之一。它的核心思路是假设一个基因在一个样本里的counts服从负二项分布然后用所有基因的数据来估计每个样本的文库大小size factor和基因自身的离散度最后通过广义线性模型判断基因在不同分组之间是否存在显著差异。为什么要用负二项分布而不是泊松分布因为RNA-seq的counts数据存在“过离散”现象即方差大于均值。泊松分布假设均值等于方差对真实数据来说太理想化了负二项分布多了一个离散度参数能更好地描述生物重复之间的波动。而输入必须是原始整数counts原因也很直接DESeq2内部要自己估计size factor这个步骤用到的是所有基因的原始read数。如果你拿FPKM进去相当于告诉它文库大小都一样它内部会基于数值的大小做归一化结果自然是错的。5.2 完整启动一个DESeq2分析假设你有一个表达矩阵是count_matrix.txt列名是样本名行名是基因名。还有一个实验设计表colData.txt里面有样本名和分组信息。下面是启动分析的基础代码# 读取表达矩阵和样本信息 count_matrix - read.delim(count_matrix.txt, row.names 1, check.names FALSE) col_data - read.delim(colData.txt, row.names 1) # 保证表达矩阵的列与样本顺序一致 count_matrix - count_matrix[, rownames(col_data)] # 构建DESeq2对象 library(DESeq2) dds - DESeqDataSetFromMatrix( countData count_matrix, colData col_data, design ~ condition ) # 设定对照组水平 dds$condition - relevel(dds$condition, ref control) # 运行差异分析 dds - DESeq(dds) # 提取结果 res - results(dds, contrast c(condition, treatment, control)) res - as.data.frame(res)design ~ condition意思是差异只与分组有关。如果你的实验有多个变量比如包含性别或批次设计公式可以写成~ batch condition这样DESeq2会把这些因素作为协变量纳入模型减少它们对差异表达结果的干扰。这是很多人没注意到的地方——有批次效应不放进design里等于把批次效应当成了随机噪声结果可能被隐藏在批次差异里的假信号带偏。5.3 从结果表格里挑差异基因DESeq2的结果表里有几列核心信息baseMean该基因在所有样本中归一化后的平均表达量。log2FoldChange处理组相对对照组的表达倍数变化log2尺度。pvalue和padj显著性检验的p值和多重假设检验校正后的p值。挑选差异表达基因的常规标准是padj 0.05且|log2FoldChange| 1也就是表达量翻倍或减半。这个标准没有绝对的对错如果你的结果里差异基因太少可以放宽到padj 0.1如果差异基因太多也可以把阈值收紧到log2FoldChange 2。但我要提醒一句只看显著性是不够的。很多基因虽然padj显著但baseMean非常低比如在几乎所有样本里都只有1-2条reads。这类基因的差异往往不可靠最好在筛选时加上baseMean 10的条件。这就相当于生物信息学里的“看得见的表达量”可以滤掉大量低表达噪声。5.4 差异倍数收缩LFC shrinkageDESeq2里还有一个重要概念是“倍数变化收缩”。对于低表达基因或重复数较少的基因它的log2FoldChange波动会很大可能因为一个样本的微小变化就跑到5以上。为了解决这个问题DESeq2提供了一种经验贝叶斯收缩方法把这些不稳定的倍数变化向0收缩。许多文章里会推荐用lfcShrink函数生成收缩后的结果再用于后续的火山图或基因功能注释res_shrink - lfcShrink(dds, contrast c(condition, treatment, control), type apeglm)这种收缩处理对排名影响非常大。尤其是在样本量少、生物学重复只有2-3个的实验里收缩后的log2FoldChange排序比原始结果要稳健得多。做功能富集之前如果用的是自己挑选的差异基因列表我建议用收缩后的结果来排序。6. 功能富集分析与背后的统计原理6.1 从差异基因列表到GO/KEGG富集拿到差异表达基因列表后大家最常做的就是GOGene Ontology富集分析和KEGG通路富集分析。这个分析本质上做的是“超几何检验”给定一组差异基因看在某个功能类别比如“免疫应答”或“p53信号通路”里差异基因出现的比例是不是显著高于随机预期。超几何检验可以理解为“摸球问题”假设一共有N个基因其中M个基因属于某个通路你选了K个差异基因其中恰好有x个落在这个通路里那么计算这个事件发生的概率是多少。如果概率极小就说明这个通路和差异基因显著相关。在R里最常用的工具是clusterProfilerlibrary(clusterProfiler) library(org.Hs.eg.db) # 把差异基因的Symbol转成ENTREZID deg_symbol - read.table(deg_symbol.txt, header FALSE)$V1 deg_entrez - bitr(deg_symbol, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db) # GO富集 go_results - enrichGO( gene deg_entrez$ENTREZID, OrgDb org.Hs.eg.db, keyType ENTREZID, ont BP, pAdjustMethod BH, qvalueCutoff 0.05 ) # KEGG富集 kegg_results - enrichKEGG( gene deg_entrez$ENTREZID, organism hsa, pvalueCutoff 0.05, qvalueCutoff 0.05 )6.2 基因ID转换的丢失问题做富集分析遇到的一个高频问题是ID转换“丢失”——用bitr转换基因Symbol到ENTREZID时总会有一部分基因查不到对应条目然后被丢弃。这种情况太常见了尤其是那些非标准Symbol或者新型基因ID。我遇到过最夸张的一次3000个差异基因里只有1800个成功转换丢失率高达40%。排查下来发现原因是上游GTF版本比较老而org.Hs.eg.db数据库是新的导致一部分历史别名查不到。解决方法很简单转换完成后做个反向检查把丢失的Symbol用AnnotationDbi::mapIds再查找一次别名如果还找不到就接受现实毕竟富集分析基于的是可注释基因子集只要丢失率低于20%影响可以接受。6.3 ORA和GSEA怎么选常规的GO/KEGG富集属于ORAOver-Representation Analysis过表达分析它只关心“差异基因列表”里有哪类功能富集不看表达量的高低趋势差异。这种方式简单直观但缺点很明显你只用了显著差异的基因大量表达量有变化但不显著的基因信息被丢弃了。GSEAGene Set Enrichment Analysis的思路不一样。它把全体基因按某种指标排序比如fold change然后看预先定义的基因集是显著富集在排序的顶部还是底部从而推断哪些通路被整体激活或抑制。GSEA不需要你事先挑差异基因所以更适合处理那些差异不剧烈但功能协调变化的生物学过程。如果你做的是临床样本或者发育时间序列这类连续型数据我强烈建议把GSEA作为ORA的补充。用clusterProfiler里的gseKEGG或gseGO时输入需要的是一个按fold change排序的基因名向量注意保留gene symbol的顺序不要把排序后的数据再打乱。7. 实战中的高频问题与排查经验7.1 比对率过低很多人跑完STAR发现比对率只有60%甚至更低第一反应是怀疑比对参数不对。我遇到的情况里这种问题往往是样品污染或参考基因组版本不一致导致的。首先用fastp看质控报告里的GC含量和duplication如果GC异常、duplication异常高很可能是建库时混入了其他物种的RNA。其次检查参考基因组版本某些基因的mRNA序列和方法在不同版本之间差异很大比对率自然会下降。如果排除这些因素比对率仍然低可以检查一下reads长度是否和STAR的--sjdbOverhang参数匹配。建索引时设了149但后来发现实际读长是100bp数值差得不多也有一定影响重新建索引即可。7.2 批次效应的典型表现与处理批次效应可以理解为“实验操作顺序不同带来的系统性差异”。它不一定来自测序更常见的是同一个项目的样本分了两批提取RNA、两批建库结果在PCA图上第一批样本和第二批样本整齐地分成两个簇。处理批次效应的首选方案不是用工具“去批次”而是在实验设计阶段就做好平衡把对照和处理样本混匀在同一批建库和测序里。如果已经晚了就在DESeq2的design公式里加入批次变量例如design ~ batch condition。极端情况下如果批次效应过于严重再考虑使用ComBat-seq之类的工具但这类工具可能过度矫正应作为最后手段。我个人的经验是先画PCA图。如果同一处理的样本没有聚在一起而是按批次聚类那你后续解释任何差异结果之前都要先把批次问题解决。7.3 服务器内存不够怎么办STAR比对非常吃内存我见过不少人在小型工作站上跑STAR时直接OOM。解决办法有三个方向如果人类数据STAR必须用就是换内存更大的机器用Salmon替代STAR做定量它几乎不占什么内存速度还快通过--genomeLoad和共享内存优化一次只加载一个样本到内存适合同一台机器跑多个样本的场景。7.4 实用小工具让流程更顺畅跑全流程时不要完全靠肉眼跟踪进度建议把fastp、STAR、featureCounts这些步骤串成snakemake规则或nextflow流程。哪怕一开始只写最简单的shell循环也要把每一步的输出目录和日志文件固定下来。有了统一的命名规范比如results/20240101_version1/sample_fastp.log以后排查问题会节省大量时间。7.5 常见问题速查表现象可能原因推荐措施比对率70%样本污染/参考基因组版本不匹配检查GC含量、duplication核对FASTA与GTF版本duplication70%建库扩增过多确认后续定量结果是否异常必要时重做文库差异基因过多(上万)分组差异明显或存在批次效应检查PCA图考虑加入批次变量差异基因过少(几个)筛选阈值太严/样本间差异太小放宽padj或log2FC阈值检查重复数PCA中样本不按分组聚类批次效应或分组定义错误检查colData分组信息考虑批次处理GTF与基因版本不匹配转换ID时大量丢失换对应版本的GTF重新定量R包装不上系统依赖缺失用conda安装R或单独装libcurl、libxml2等我自己在实际操作里还有一个习惯每跑完一个样本就在日志里记录下该样本的关键指标——raw reads、clean reads、比对率、gene assigned比例。当整个项目的样本都跑完后把这些指标汇总成一张表一方面方便检查异常样本另一方面在写Methods或回复审稿人时可以直接引用这些质控数据。RNA-seq的分析流程虽然标准化程度高但每个项目都会遇到自己独特的麻烦。上面这些内容是我多次跑完整套流程之后沉淀下来的经验也是我每次带新人时都会重点强调的知识点。数据不会说谎但前提是你要用对工具、设对参数、并且每一步都理解为什么要这么做。如果你也是刚接触高通量测序数据分析不妨从这条路一步步走下来整个分析逻辑会非常清晰。
阅读完成 · 觉得有帮助？
