RNA-seq转录组数据分析流程详解：从fastq到差异基因

简介：一个面向生物信息学初学者的RNA-seq转录组数据分析项目代码，以两个不同组织的6个样本为例，完整演示从原始测序数据到生物学结论的闭环流程，覆盖数据质控、序列比对、定量表达、差异分析及功能富集等核心环节。压缩包内共3个文件，包含一个.inscode代码文件、一个index.html说明页面及一个.gitignore配置文件，整体大小仅8KB，结构精简，便于快速上手。已有518人学习该资源。代码中整合了FastQC、Trimmomatic、HISAT2、featureCounts和DESeq2等常见工具的分析步骤，并给出GO与KEGG富集分析、差异表达基因火山图的绘制方法；同时附带样例数据（经特殊修改）与结果解读提示，读者可边看边练，将分析流程直接落地为可运行的项目。

1. RNA-seq转录组数据分析流程在解决什么问题：拿到fastq后真正要花心思的地方

做转录组分析一年多后，我最怕的不是数据量大，而是朋友发来一句话：RNA-seq转录组数据分析流程听上去很标准，但真正跑起来，很多人会卡在各式各样的参数选择和异常排查上。对方以为这是一条命令的事，但从fastq到差异基因，你要经过质控、比对、定量、差异表达统计、可视化五个环节，每一步的错误都会被下游结果放大。这篇笔记把我常用的项目代码和参数整理成一条可以直接照抄的路线，从原始测序数据一路跑到火山图和热图，并把我踩过的几个坑按“现象、原因、解决”写清楚。适合三类人：刚接触生信、手里有公共数据想自己分析的实验人员；需要把这套流程批量跑几十个样本的工程师；以及想弄明白比对率和DESeq2结果为什么总出问题的研究助理。核心原则是先立住数据质量，再做统计推断，只看图不看质控报告是拿不到可靠差异基因列表的。

2. 从fastq到比对结果：fastp质控与STAR两遍比对的完整命令

拿到测序数据后别急着上比对工具，先花10分钟做质控。这一步省下的时间会在后面几十倍地还回来。

2.1 为什么要先做质控：接头残留是最容易被忽略的比对率杀手

转录组文库上机前虽然有过一轮QC，但送到你手里的fastq里仍然可能藏着三类问题：接头序列没切干净、read末端的polyG尾巴（部分测序平台双色化学法特有的标记）、以及整体偏低的质量分数。这些问题不会让程序报错，只会悄悄压低唯一比对率，或者把多比对reads的比例推高。我在模拟数据里只添加一小段残留接头，STAR唯一比对率就从95%掉到了82%左右，这个幅度足以让你的差异基因列表大洗牌。

质控参数要针对转录组而不是基因组来设置。转录组reads主要落在外显子区域，polyG尾巴对带polyA尾的基因干扰尤其明显，如果不处理，后续featureCounts定量时会有大量reads数到基因间区或者基因上游区域，下游结果里出现一堆假阳性。常见做法是把fastp放在样本循环里逐对reads处理，跑完后先看HTML报告再决定是否调整参数。

2.2 fastp最小命令与参数解读

下面这段是我项目代码里最常用的fastp调用，直接复制到终端就能跑：

fastp -i data/raw/S1_R1.fastq.gz -I data/raw/S1_R2.fastq.gz \ -o data/clean/S1_R1.fastq.gz -O data/clean/S1_R2.fastq.gz \ --trim_poly_g --length_required 30 \ --qualified_quality_phred 20 --cut_front --cut_tail \ -h results/fastp/S1.html -j results/fastp/S1.json

参数说明：

• --trim_poly_g：切除read末尾的polyG。这个参数对照转录组数据建议默认开启，尤其是测序平台较新的数据。
• --length_required 30：过滤掉长度低于30bp的read。太短的read很难定位到基因组上，还容易比对到重复区域产生假计数。
• --qualified_quality_phred 20：平均质量Q20以下整条read丢掉。转录组还可以再严格一点设Q25，但要看整体数据质量，强行调高会删掉太多reads。
• --cut_front --cut_tail：在read两端做滑动窗口质量裁剪。对转录组有用，但别开得太狠，否则大量读长不均一的reads会堆到比对文件的末尾，影响后续splice junction识别。

跑完fastp后别急着跑下一步，先打开S1.html看两个位置：per base sequence content的尾部曲线，以及adapter污染统计。如果尾部曲线出现明显的A/T分离，说明polyG没切干净，需要回头检查--trim_poly_g；如果曲线正常但adapter比例高于5%，说明接头切除还有遗漏。

2.3 STAR索引构建与两遍比对

比对工具选择上，我一般直接用STAR而不用其他比对器。它的剪接位点支持好，两遍比对模式能重新学习新剪接位点，速度上对常规人类样本也够快。代价是内存占用大，人类基因组索引需要30GB左右内存。如果你的机器内存只有16GB，老老实实换其他方案，但接受剪接位点检测略粗糙的代价。

索引构建先做一次，后续所有同参考基因组样本共用：

STAR --runMode genomeGenerate \ --genomeDir ref/star_index \ --genomeFastaFiles ref/GRCh38.fa \ --sjdbGTFfile ref/genes.gtf \ --sjdbOverhang 149 \ --runThreadN 16

参数说明：

• genomeDir：索引输出目录。建议单独建一个ref目录，和结果目录分开，避免不小心把索引文件混进结果里。
• sjdbOverhang：一般设为读长减1，PE150就写149。这个值影响剪接位点数据库的构建，写大了增加索引体积，写小了降低剪接reads的比对灵敏度，只在建索引时设置一次。
• runThreadN：线程数，按机器核数设，不要超过物理核心数。

比对命令，默认开两遍模式：

STAR --genomeDir ref/star_index \ --readFilesIn data/clean/S1_R1.fastq.gz data/clean/S1_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outSAMattributes NH HI AS NM MD \ --quantMode GeneCounts \ --twopassMode Basic \ --outFileNamePrefix results/bam/S1_ \ --runThreadN 16

参数说明：

• --outSAMtype BAM SortedByCoordinate：直接输出按坐标排序的BAM，下游featureCounts需要这个顺序。
• --outSAMattributes NH HI AS NM MD：保留比对标签。featureCounts和samtools后续都会用到NH和HI，AS和NM在检查比对质量时有用，建议默认带上。
• --twopassMode Basic：两遍比对。先做一遍快速比对收集新的剪接位点，再带着这些新位点做第二遍正式比对。只有研究可变剪接事件时需要开两遍；如果只做基因水平的差异表达，一遍比对就够了，速度能快不少。我一般会在基因差异表达项目里关掉这个参数。
• --quantMode GeneCounts：顺便输出每个基因的raw counts，但我不建议直接拿它替代featureCounts，它在处理多overlap和rRNA过滤上的默认行为不如featureCounts直观。

2.4 比对结果检查：比对率与染色体reads分布

STAR跑完第一步不是急着定量，先看results/bam/S1_Log.final.out文件。重点关注两行：Uniquely mapped reads number和Uniquely mapped reads %。正常转录组样本唯一比对率一般在85%以上，低于75%时先别往下走，查接头、查污染、查参考基因组版本是否匹配。

还可以顺手检查reads的染色体分布，看有没有线粒体或rRNA污染：

samtools idxstats results/bam/S1_Aligned.sortedByCoord.out.bam | \ sort -k3 -nr | head -10

输出里chrM占比如果比同批次样本明显偏高，说明文库里有较多线粒体表达量。组织样本里chrM占5%到20%都算正常，但单核细胞或某些细胞系会更高，这个比例需要和你的样本类型横向对比。chrM reads异常高时，可以在比对完成后再过滤掉MT上的reads，避免干扰后续内参基因的评估。这一步做完，BAM文件才算合格，可以进入定量环节。

3. 定量一步：featureCounts批量生成表达矩阵

比对完成后的BAM里记录的是每个read的基因组位置，但转录组分析需要知道哪些read属于哪个基因的外显子，最后汇总成基因表达矩阵。这一步工具很多，我基本固定用featureCounts，原因是参数少、内存占用低、输出格式直接能被DESeq2吃下去。

3.1 搞清楚featureCounts到底在数什么

featureCounts的原理是把BAM里的比对片段和GTF注释里的外显子区间做交集，一个片段与某个基因的任意外显子有overlap，就记到该基因名下。它默认只统计唯一比对的reads，一个read同时比对到多个位置会被直接丢弃。这也是前面为什么强调保留NH和HI标签的原因。

这里有个很容易被忽略的细节：配对测序时，featureCounts默认数的是fragment而不是read。也就是说一对paired-end reads共用一个fragment，只计一次。如果你的下游差异分析工具用的是raw counts，必须理解这个计数单位，否则在文库大小归一化时会高估或低估样本间的真实比例。定量这一步不要把工具当黑匣子，至少要清楚它输出矩阵里每一行每一列是怎么来的。

3.2 单条命令跑完所有样本并合并矩阵

featureCounts可以在一次调用里接收多个BAM文件，不需要自己写循环逐个跑。示例命令：

featureCounts -a ref/genes.gtf \ -o results/counts/counts_all.txt \ -T 8 -p -s 0 \ results/bam/S1_Aligned.sortedByCoord.out.bam \ results/bam/S2_Aligned.sortedByCoord.out.bam \ results/bam/S3_Aligned.sortedByCoord.out.bam \ results/bam/S4_Aligned.sortedByCoord.out.bam

参数说明：

• -a：GTF注释文件路径，必须和STAR建索引时用的一致。
• -T 8：线程数。featureCounts对线程扩容效果一般，设8到16就够。
• -p：配对末端样本启用fragment计数，不要漏掉这个参数。
• -s 0：链特异性设置。0表示未链特异性文库，1表示reads来自正链，2表示来自负链。你的文库是stranded还是unstranded，直接决定这个值，不确定时先看文库构建说明，或者用默认0跑完看gene body coverage判断。
• 后面的多个BAM文件会自动合并到一个矩阵里，输出文件counts_all.txt中每个样本占一列。

跑完以后用head看一眼输出文件，注意第一行以#开头，这是命令行记录；第二行才是表头，包含Geneid、Chr、Start、End、Strand、Length六个属性列，从第7列开始才是真正的样本counts。这个格式固定，读取时按这个结构处理就不会乱。

3.3 原始计数矩阵的读取与TPM转换

在R里读取featureCounts输出时，要跳过第一行注释：

raw_counts <- read.table("results/counts/counts_all.txt", header = TRUE, row.names = 1, skip = 1, check.names = FALSE) gene_length <- raw_counts$Length counts <- raw_counts[, 6:ncol(raw_counts)] colnames(counts) <- sub("_Aligned\\.sortedByCoord\\.out\\.bam$", "", basename(colnames(counts)))

这里skip = 1是因为第一行是参数注释，直接跳过；check.names = FALSE是为了让表达矩阵的列名保留原始BAM文件名，方便用basename提取样本名。从第6列开始取是因为默认输出里第2到第6列是Chr、Start、End、Strand、Length，这些是基因属性不是表达量。

RNA-seq差异分析用raw counts可以直接喂给DESeq2，不需要做长度归一化。但如果你要做样本间表达量比较、画PCA或者热图，一般用TPM。从raw counts转TPM的常见做法是先用基因长度换算成RPK，再按样本归一化：

rpk <- counts / (gene_length / 1000) tpm <- t(t(rpk) / colSums(rpk) * 1e6)

这段代码核心就两行，但有个地方经常出问题：colSums(rpk)的分母是总RPK，不是raw count总和，写错一步TPM就会整体偏移。我一般只用TPM做样本间表达量的横向展示，差异分析一律回退到DESeq2的median-of-ratios归一化，因为TPM对低表达基因的波动太敏感，作为统计检验输入并不合适。

4. 差异表达分析：DESeq2设计式、完整代码与结果解读

表达矩阵生成后，核心分析才算开始。差异表达分析我固定用DESeq2，它比edgeR的入门门槛略高一点，但优点是对小样本量和非平衡设计的处理更稳定，而且设计公式的设计思路能迫使你提前把实验结构想清楚。

4.1 把表达矩阵变成DESeq2能吃的对象

DESeq2需要两个输入：counts矩阵和样本信息表格。有个前置条件：counts矩阵必须是整数，千万不能传TPM或FPKM进去，DESeq2内部需要根据整数计数来估计离散度。如果从外部拿到的是已经归一化过的表达矩阵，就别再套DESeq2了。

样本信息通常是一个数据框，行名是样本名，列至少包括分组条件：

sample_info <- data.frame( row.names = colnames(counts), condition = factor(c("control", "control", "treat", "treat")), batch = factor(c("batch1", "batch1", "batch2", "batch2")) )

行名必须和counts矩阵的列名完全一致，DESeq2按行名对齐，顺序无关。condition和batch都要转成factor，否则DESeq2会报错。

design公式是DESeq2的核心设计决策。常见写法是~ batch + condition，意思是在扣除批次效应的前提下比较分组差异。如果你的实验没有批次问题，直接写~ condition即可。design公式里加进条件外的东西，都会改变结果的log2FC含义，所以别乱加变量，只有确定会影响表达的协变量才放进去。

4.2 完整DESeq2运行代码

这是整个流程里最标准化的一段，我每次都是复制同一个模板：

library(DESeq2) dds <- DESeqDataSetFromMatrix( countData = counts, colData = sample_info, design = ~ batch + condition ) dds <- DESeq(dds) res <- results(dds, contrast = c("condition", "treat", "control"), alpha = 0.05) res_shrink <- lfcShrink(dds, contrast = c("condition", "treat", "control"), res = res, type = "apeglm")

lfcShrink是获取稳定log2FC的关键一步。DESeq2原始输出的log2FC对低表达基因波动很大，尤其是那些counts接近0的基因，shrink之后基因排名更稳定，也更适合用来画火山图和排序。type = "apeglm"是常用的一种收缩算法，需要先安装对应包；如果没有装，可以跳过这行直接用results的输出，但注意两套log2FC数值会有差异，你最终写进文章里的那套要固定用同一个来源。

results输出表格里每列的含义可以用下面这个方式理解：

4.3 结果过滤、火山图与热图

拿到结果后，筛选显著基因的常见阈值是padj < 0.05且|log2FC| >= 1。写成R代码：

sig <- res_shrink[!is.na(res_shrink$padj) & res_shrink$padj < 0.05 & abs(res_shrink$log2FoldChange) >= 1, ]

过滤后别忘了输出一张csv备份，基因ID、padj、log2FC、baseMean都要留全，方便后面随时回溯。

火山图和热图是这个流程最后能直接发给合作者的产出。火山图推荐用shrink后的log2FC，否则低表达基因会在图上拖出一条假差异的长尾。热图我一般用VST变换后的数据而不是raw counts：

vsd <- varianceStabilizingTransformation(dds) top_genes <- head(order(res_shrink$padj), 50) pheatmap::pheatmap(assay(vsd)[top_genes, ], scale = "row")

scale = "row"是按基因做标准化，让每个基因在样本间的相对高低体现出来，不然高表达基因会把低表达基因的颜色差异压没。热图样本聚类如果和你预期分组不一致，先别急着改图，回去看原始数据是不是存在批次问题。

还有一个常见困惑：DESeq2与edgeR、limma给出的差异基因列表经常不一样，这不是谁错了，而是它们对低表达基因的处理和对离散度的估计方式不同。我的习惯是以DESeq2为主要结果，再用limma做一次一致性交叉验证。如果两套列表overlap率很低，优先怀疑数据质量或样本分组，而不是急着解释生物学意义。

5. RNA-seq流程里的五个常见坑：现象、原因与排查路径

这一章把我踩得比较深的五次翻车记录写出来，每条按“现象、原因、解决”展开，跑流程遇到异常时可以直接对照排查。

1. 比对率只有60%多

现象：Log.final.out里Uniquely mapped只有60%到70%，Unmapped比例明显偏高。
原因：最常见的是参考基因组版本和文库来源不一致；其次是接头没切干净，短read无法比对；还有一种情况是样本被rRNA污染，大量reads比对到rRNA区域后被过滤掉。
解决：先看Unmapped reads的长度分布，集中分布在30到50bp的话基本是短read污染；用fastp重新检查adapter；再用一个rRNA参考序列做比对确认污染比例。版本问题只能重新准备同一来源的基因组和GTF，这是没有后悔药的一步，项目开始前必须确认参考版本，中途换版本意味着全部样本重跑。

2. PCA把样本分开，但和实验分组无关

现象：画出PCA后，前两个主成分把样本按送测批次分成两团，组间差异反而不明显。
原因：批次效应是转录组分析里最玄学的部分。细胞培养批次、RNA提取时间、建库试剂盒、测序lane，任何一个因素都会贡献系统性噪声。
解决：在设计阶段尽量让不同分组均匀分布在同一次建库和测序中；无法避免时，在DESeq2的design里加入批次因子，比如~ batch + condition。也可以用limma的removeBatchEffect先看PCA是否变干净，但要记住批次因子只在design里参与建模，results里的log2FC才是校正后的结果。

3. padj全部是0.999，没有任何差异基因

现象：DESeq2运行不报错，但results表格里padj几乎全等于0.999。
原因：每组生物学重复太少，最常见是每组只有2个样本，离散度估计不出来；或者组内有一个离群样本把方差撑得过大。
解决：实验设计上至少每组3个重复，这是老生常谈但现实里最容易做不到。分析层面先看sizeFactors是否出现异常样本，再检查样本聚类，把离群样本剔除后重新跑。剔除样本要在方法部分写明依据，否则审稿人会有疑问。

4. featureCounts输出的counts全为0

现象：矩阵能正常生成，但部分样本整行或整列全是0。
原因：BAM里的染色体命名和GTF里的染色体命名对不上。常见情况是GTF用chr1而BAM的@SQ里是1，或者反过来。
解决：用samtools view -H查看BAM的@SQ列表，再用head查看GTF第一列，确认两者命名风格一致。如果不一致，重新准备来源匹配的参考文件，不要手动改GTF，手动替换染色体前缀很容易引入坐标错位。

5. 热图完全看不出分组结构

现象：差异基因热图花成一团，样本没有按condition聚类。
原因：没有做方差稳定化直接用了raw count；或者使用了全部基因而不是高变基因；又或者差异本来就少，硬画也画不出结构。
解决：改用VST或rlog后的数据画热图，并且用padj排序前50的基因，或按行方差排序取top2000高变基因。如果这样还是看不出聚类，回第4章重新检查design公式和重复数，问题大概率不在图上而在数据里。

排查完这些坑，整个流程的稳定性会提高一个档次。建议每次跑完分析保存一份sessionInfo输出，工具版本对不上时，这个文件能帮你省下大量定位时间。

6. 沉淀成项目代码：最小Snakemake流程与参数归档习惯

流程能跑通只是第一步，真正值得投入的是把它整理成随时可以重跑、可以交接的项目代码。我现在的做法是用Snakemake把质控、比对、定量串起来。最小示例：

SAMPLES = ["S1", "S2", "S3"] rule all: input: "results/counts/counts_all.txt" rule fastp: input: "data/raw/{sample}_R1.fastq.gz", "data/raw/{sample}_R2.fastq.gz" output: "results/clean/{sample}_R1.fastq.gz", "results/clean/{sample}_R2.fastq.gz" shell: "fastp -i {input[0]} -I {input[1]} -o {output[0]} -O {output[1]} " "--trim_poly_g --length_required 30" rule star: input: "results/clean/{sample}_R1.fastq.gz", "results/clean/{sample}_R2.fastq.gz" output: "results/bam/{sample}_Aligned.sortedByCoord.out.bam" shell: "STAR --genomeDir {config[genome]} " "--readFilesIn {input[0]} {input[1]} --readFilesCommand zcat " "--outSAMtype BAM SortedByCoordinate " "--outFileNamePrefix results/bam/{wildcards.sample}_ " "--twopassMode Basic --runThreadN 16" rule featureCounts: input: expand("results/bam/{sample}_Aligned.sortedByCoord.out.bam", sample=SAMPLES) output: "results/counts/counts_all.txt" shell: "featureCounts -a {config[gtf]} -o {output} -T 8 -p {input}"

有个血泪经验：STAR的--outFileNamePrefix会决定最终BAM文件名，Snakemake里声明的output必须和实际产物完全一致，否则任务永远提示缺输出，实际上文件已经躺在目录里了。跑之前先在终端手动执行一条STAR命令，确认完整文件名，再写进规则。

除此之外，我习惯在每个项目根目录放一份version.txt，记录工具名、版本号和参考文件路径，所有命令里的参考文件都通过配置文件传给Snakemake，而不是硬编码在脚本里。很多不可复现的结果，源头都是版本漂移。现在的我每接一个新项目，先建好目录结构、固定参考版本、写好规则，再开始处理样本，数据跑完半年后还能原样重跑。希望帮到你。

本文还有配套的精品资源，点击获取