简介：一份面向生物信息学入门者的RNA-seq转录组数据分析流程项目代码包，以两个组织、六个样本的实例贯穿数据质控、序列比对、定量表达、差异表达分析与功能富集等环节，适合需要从零上手转录组分析的学生和科研人员。包内共3个文件，以HTML说明页、inscode项目配置及gitignore辅助文件为主，压缩包仅8KB，轻量简洁，便于快速查看分析框架并二次修改。资源发布以来已有518人浏览学习。代码详细展示了FastQC、Trimmomatic、HISAT2、featureCounts、DESeq2等工具的调用思路，同时给出GO与KEGG富集分析以及火山图绘制的实现方案，并特别注明数据已经特殊修改，仅供学习参考，不必担心版权与隐私问题。拿到代码后可直接对照项目结构梳理流程要点，并借助示例中的参数设置与输出解读，理解从原始数据到生物学结论的完整链路；相比翻阅零散文档，这份紧凑的代码包能有效降低入门RNA-seq数据分析的学习门槛。

1. RNA-seq转录组数据分析流程：拿到FASTQ之后，你真正要面对的是什么

如果你手里刚拿到一批RNA-seq转录组测序数据，头一件事往往不是打开软件，而是对着几个G的FASTQ文件发愣：文件名后缀乱七八糟、质量报告一大堆、不知道从哪一步开始才算是“正式分析”。这个标题下的完整流程，本质上只回答一个问题——如何把一批下机读段变成一张差异基因表，并让这个结果经得起同行质疑。它适合两类人：一类是刚接触转录组、只会用在线工具点按钮的初学者，另一类是被上游流程折腾过、想确认自己每个参数是否选对的从业者。整条链路并不依赖某个黑匣子平台，而是由一组命令行工具串联：质控、剪接感知比对、基因计数、差异表达、通路解释，每一步都有明确的输入输出和可复现的日志。

2. 流程拓扑：为什么是“比对—定量—差异”三板斧

2.1 剪接感知比对为什么是RNA-seq的第一道分水岭

转录组数据与DNA重测序在比对策略上的根本差异，在于真核生物的成熟mRNA已经把内含子剪掉了。一段跨越外显子连接处的读段，如果直接丢给普通比对工具，会被当成“跨过大段缺失的异常比对”，轻则被标记为剪接位点冲突，重则直接被丢弃。所以RNA-seq的标准做法，是使用剪接感知比对器，它在比对时允许读段被拆成两个或多个片段，分别落在两个外显子上，中间用一段较长的“内含子跨度”连接起来。

常见做法是优先考虑STAR、HISAT2这类比对器。STAR之所以在长读长、高表达动态范围的转录组项目里几乎是默认选择，是因为它的核心思想不是逐碱基延伸，而是先建立参考基因组的后缀数组索引，把读段的种子序列快速定位到基因组位置，再向外扩展。这一策略让它在处理全长转录本、跨外显子junction时速度非常快，而且对双端读段的一致性要求更严格。代价是内存占用偏高，构建一次人类基因组索引通常需要30GB左右的内存，这一点在选型时必须有心理准备。

另一类方案是salmon或kallisto这类伪比对工具，它们不把读段逐碱基比对到基因组，而是直接把k-mer映射到转录本上，本质上是“跳过比对、直接定量”。它们的速度比STAR高一个数量级，但代价是拿不到基因组层面的信息，例如剪接位点、新转录本、变异信息全都没有。如果你的分析终点只是差异表达，伪比对是省时间的合理选择；如果你的流程下游还要做剪接分析、融合基因或者变异检测，老老实实走STAR更稳妥。

我一般会这样选型：有参考基因组、样本量在几十个以内、服务器内存足够——直接上STAR；样本量到了上百个、只关心基因表达量变化、机器配置一般——考虑salmon配合tximport。这个决定不是“哪个工具更先进”，而是“你的下游问题到底需要什么粒度的信息”。伪比对给不出基因组坐标，这是它的天然边界。

2.2 从BAM到基因表达矩阵：featureCounts的定量模型与参数选择

比对完成之后，每一个读段在BAM文件里都带上了基因组坐标。接下来要做的是把这些坐标折算回基因：一个读段落在某个基因的外显子区域上，就给它记上一笔。featureCounts是目前最主流的计数工具之一，它不依赖比对时是否设置了--quantMode GeneCounts，而是直接读取BAM文件和GTF注释，把重叠关系统计成“每个基因有多少个读段”。它的统计模型是meta-feature级别的：GTF里属于同一个基因的所有外显子会被合并成一个区间集合，读段只要与其中任何一段有重叠，就计入这个基因。

这里有一个新手最容易翻车的点：双端测序的一条读段对（read pair）本质上来自同一个cDNA片段，两条read分别测了同一段分子的两端。如果两条read都落在这个基因上，直接按“两条read”计数，相当于把一个片段数了两次。featureCounts对双端数据的默认行为是只统计“read对”中能与特征成功配对的那个端，加上-p参数后才会把两端配对视为一个计数单元。还有一个容易被忽略的参数--countReadPairs，它决定了在-p模式下是“用所有配对成功的片段”还是“仅用至少一端重叠的片段”进行统计。

参数选择上，链特异性是另一个必须明确的决定。常见建库策略有三种：未链特异性、正向链、反向链。大多数市售链特异性试剂盒是反向链策略，也就是第二条链的序列才是mRNA的反向互补。featureCounts里用-s 0表示不区分链，-s 1表示正向链，-s 2表示反向链。如果建库是反向链而你用了-s 0，结果通常“也能出数”，但计数噪声会明显增加；如果用了-s 1，大量基因的计数会趋近于零。我的建议是：拿到数据先问清楚建库方式，不要靠猜，不然下游差异分析里会出现一整批看起来“不表达”的假阴性基因。

实际跑一次featureCounts的命令通常是这样的：

featureCounts -a genes.gtf \ -o counts.txt \ -T 8 \ -p --countReadPairs \ -s 2 \ *.bam

-a指定GTF注释文件，-o指定输出文件前缀，-T是并行线程数，-p启用paired-end模式，--countReadPairs让每个片段只算一次，-s 2对应反向链建库。输出会生成一个counts.txt和同名的.summary文件，后者会按“已分配”“未分配”“多重比对”等类别统计每个BAM的读段去向，这是你判断数据是否干净的第一张体检表。如果“未分配”的比例超过30%，先别急着跑差异分析，回头检查GTF版本和比对参数，这个习惯能省掉后面一整轮的返工。

3. 上游质控与比对实操：从原始数据到BAM

3.1 参考基因组与GTF的准备流程：版本对齐是第一原则

整个流程里最不性感但回报率最高的步骤，是构建一套版本互相匹配的参考基因组和注释文件。很多项目做到一半发现“比对率正常但差异基因全是假象”，回头一查多半是基因组序列来自一个版本，GTF注释来自另一个版本，染色体命名规则甚至都不一致。常见做法是从同一个发布渠道同时下载基因组FASTA和对应的GTF注释，并且确认二者属于同一构建号。如果你用的参考物种有多个命名版本，例如UCSC风格和Ensembl风格，需要统一成一种，否则STAR构建索引时就会把注释里的染色体名全部落空。

拿到文件之后，养成给文件重命名的习惯：genome.fa、genes.gtf、star_index/三个目录，不要保留原始下载的乱码文件名。这样做的目的不是整洁强迫症，而是所有下游工具只认你给的路径，一旦某个文件的软链接断了，STAR会报一个含糊的“cannot open file”错误，浪费大量排查时间。另外，GTF文件建议先做一次简单的格式检查，确认第九列包含gene_id属性，因为STAR和featureCounts都依赖这个字段做聚合。

STAR索引的构建参数里，--sjdbOverhang是最常被问到的：它表示构建剪接位点数据库时，在注释的剪接位点两侧各取多长的序列。经验值是“读段长度减一”，比如150bp读长就填149。这个值的本质是模拟一条从剪接位点一侧跨到另一侧的最长读段，让STAR建立一个基于注释的剪接位点索引。取小了，长读段跨接点时无法完全匹配；取大了，索引构建的内存和时间会无谓增加。构建命令大致如下：

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

这个命令会生成SA、SAindex、genomeParameters.txt等一组文件。如果后续GTF注释更新了，需要重建索引，而不是把新GTF硬塞给旧索引。一个隐藏坑是--genomeDir指向的目录必须为空，否则STAR会拒绝覆盖已有索引。

3.2 用STAR跑通标准双端比对：核心参数与两遍比对的必要性

比对这一步是最容易出“看起来成功但实际无效”结果的地方。STAR一次典型运行会产生Log.final.out、SJ.out.tab、BAM文件多个输出。双端数据的标准命令如下：

STAR --genomeDir ./star_index \ --readFilesIn sample_R1.fastq.gz sample_R2.fastq.gz \ --readFilesCommand zcat \ --runThreadN 16 \ --outSAMtype BAM SortedByCoordinate \ --outFilterMultimapNMax 10 \ --outFilterMismatchNMax 2 \ --outSAMattributes All \ --quantMode GeneCounts \ --twopassMode Basic \ --outFileNamePrefix sample_

逐项说明：--readFilesCommand zcat告诉STAR输入文件是gzip压缩格式；--outSAMtype BAM SortedByCoordinate直接输出按坐标排序的BAM，省去之后samtools sort这一步；--outFilterMultimapNMax 10允许一个读段最多比对到10个位置，超过就视为不可唯一比对，这个值设置得越大，唯一比对比例越低，但能保留更多重复基因家族的信号，一般5～10是折中区间；--outFilterMismatchNMax 2限制读段与参考序列的错配数，转录组数据不是DNA重测序，错配数限制不建议设为0，因为RNA编辑和测序错误都是真实存在的；--quantMode GeneCounts会让STAR顺带输出一个ReadsPerGene.out.tab文件，这是STAR自己的基因计数，脚本层面可以直接用来做快速评估，不过正式定量我仍然建议单独跑featureCounts，因为后者的注释处理逻辑更透明。

--twopassMode Basic是全流程里最值得花时间理解的一个参数。它做的事情是：第一遍先用注释提供的剪接位点做常规比对，把比对结果中新发现的剪接位点收集起来；第二遍把这些新位点合并进索引重新比对一次。这样做的价值在于，注释文件不可能覆盖真实样本中所有剪接事件，某些组织特异的剪接位点在注释里根本不存在，如果不用两遍模式，这些读段会被直接判为比对失败。代价是运行时间几乎翻倍。对于差异表达分析，两遍模式不是必选项，但对于新转录本发现、剪接定量这类下游任务，它是必需品。

比对结束后不要急着删中间文件。Log.final.out每一行都有明确含义，重点看四个数字：输入的读段总数、唯一比对比例、多位点比对比例、未比对比例。其中“多位点比对比例”如果超过20%，说明样本很可能有大量重复基因家族表达，或者参考基因组里存在未masked的重复区域。

3.3 比对后的质量筛查：看哪些指标才不算玄学

比对完成之后，很多人只看一眼“mapping rate 90%”就长舒一口气，这其实远远不够。我一般会做三层检查。第一层看STAR的Log.final.out里的唯一比对比例，通常期望在80%以上，低于这个数就要怀疑样本污染或参考基因组不匹配；第二层用samtools flagstat看BAM的比对flag分布，重点确认paired reads比例和properly paired比例，如果properly paired低于70%，说明插入片段分布异常或两条read比对到不同染色体；第三层把未比对和多位点比对的读段比例加起来，与前面featureCounts的summary做交叉验证，确认“未分配”不是发生在比对之后。

对链特异性数据还有一个检查习惯：随机挑几个已知基因，比如看一个只在正向链表达的基因和一个只在反向链表达的基因，用samtools view抽几行BAM记录，人工确认它们的strand flag与建库方式是否一致。这一步虽然手动，但能在一小时内暴露建库信息被记反的问题，比跑完整个下游再回头排查要划算得多。RNA-seq的质控本质不是“跑一遍软件”，而是“让每一步输出都能自圆其说”，这个原则贯穿整个流程。

4. 下游差异表达分析：DESeq2最小可复现流程

4.1 从featureCounts输出到DESeq2对象：矩阵整理与样本表设计

差异表达分析从这里开始才真正进入统计学范畴。在拿到featureCounts输出的counts.txt之后，我不会直接把它喂给R，而是先做一次轻度整理：把第一列基因ID设为行名，去掉长度、区间的辅助列，只保留每个样本的count列，列名改成可读的样本名。这个动作既是为了DESeq2的输入格式，也是为了让下游可视化时图例不至于变成一串乱码。

DESeq2要求两样东西：count矩阵和样本信息表。count矩阵里的数值必须是非负整数，这是转录本计数的最原始形态；任何已经做过FPKM、TPM、log2转换的数据都不能直接喂进去，这是新手最常见的报错来源。样本信息表至少要有两列：样本名和分组条件，分组条件在分析里叫condition，它决定了设计公式。设计公式~condition表示最简单的两组比较；如果你的数据有批次信息，可以用~batch + condition把批次效应纳入模型，不要默认忽略它，否则批次和生物学差异混在一起，后面谁也说清不差异基因到底是批次的还是处理的。

实际构建DESeq2对象的代码如下：

library(DESeq2) # 读入featureCounts输出，行名是基因ID countData <- read.table("counts.txt", header=TRUE, row.names = 1, check.names = FALSE) # 只保留count列，去掉chr/start/end/strand countData <- countData[, 6:ncol(countData)] # 手动准备样本表，不要直接在脚本里硬编码分组 colData <- data.frame( row.names = colnames(countData), condition = c("control", "control", "treated", "treated") ) # 确保condition是因子，且对照组是第一水平 colData$condition <- factor(colData$condition, levels = c("control", "treated")) dds <- DESeqDataSetFromMatrix(countData = countData, colData = colData, design = ~ condition)

这段代码里有一个容易忽略的细节：check.names = FALSE。如果不加，R会把样本名里的横杠转换成点，后续所有输出里的列名都与原始文件名对不上，排查起来相当痛苦。condition因子的水平顺序直接决定了差异比较的方向，在R里，因子默认把第一个水平当作基线，所以要让对照组排在最前面，否则最后结果的fold change方向是反的。等你会用DESeq2之后会发现，这个“方向反了”的错误在真实项目里频率高得吓人，而且往往要等到做通路富集时才察觉。

4.2 核心分析与结果解读：DESeq、results与关键参数的取舍

构建好dds对象之后，DESeq2的核心分析只需要两行R代码，但背后的统计模型远比这两行复杂。它会先估计每个基因的size factor来消除测序深度差异，再估计基因的离散度参数，最后用负二项分布拟合每个基因的表达量，并做组间比较。这个过程的输出是一个results对象，本质上是一个包含六列的表：baseMean、log2FoldChange、lfcSE、stat、pvalue、padj。

dds <- DESeq(dds) res <- results(dds, alpha = 0.05, contrast = c("condition", "treated", "control")) res <- lfcShrink(dds, contrast = c("condition", "treated", "control"), res = res, type = "apeglm")

这里用了lfcShrink对log2FoldChange做收缩处理。为什么要收缩？低表达基因的count数大都很小，log2FoldChange的方差会非常大，一个在处理组只有1个count、对照组为0的基因，不加收缩时fold change可能显示为无穷大或几十倍，这显然不是真实生物学效应。apeglm是目前推荐的收缩方法，它把低表达基因的效应量朝0收缩，保留高表达基因的稳定估计。这里注意一个次序问题：先跑DESeq()，再跑results()，最后用lfcShrink()导出收缩后的值，顺序乱了会报错或得到完全不同的数值。

关于差异基因的筛选阈值，最保守且约定俗成的组合是padj < 0.05且|log2FoldChange| > 1。但不要把这当成铁律。如果你的后续实验验证成本高，可以收紧到padj < 0.01；如果样本量很小、生物学差异不强，适当放宽log2FC阈值到0.58（也就是1.5倍）也能接受。关键是一定要在筛选前先看padj的分布，如果padj普遍偏大，问题往往不在阈值，而在样本间离散度太高或组内存在一个离群样本。

4.3 从结果表到判断依据：MA图、PCA与差异基因表

拿到结果后，先不要急着筛选差异基因。第一个动作是画MA图，横轴是baseMean（平均表达量），纵轴是log2FoldChange。这个图的价值在于一眼看出收缩效果是否正常：高表达区域应该有稳定的离散分布，低表达区域的点应该被压缩向0。如果MA图看起来像一把竖起来的刷子，说明你的lfcShrink被跳过了，或者某个样本的比对出了问题。绘制MA图只需要plotMA(res)一行命令。

第二个动作是PCA。DESeq2提供了vst变换功能，把count数据变换到近似正态分布，再计算主成分。PCA的主要目的是看样本是否按分组聚类。如果对照组和处理组明显分开，那说明处理效应强于样本间噪声；反过来，如果同组样本没有聚在一起，差异分析结果的可信度就要打个问号。这一步骤不要用原始count直接跑，因为测序深度差异会主导前两个主成分，掩盖真实的生物学差异。

vsd <- vst(dds, blind = TRUE) plotPCA(vsd, intgroup = "condition")

blind = TRUE表示在变换时不使用分组信息，这样的PCA可以更客观地反映样本间的真实距离。PCA图上如果发现单个样本远离同组，不要直接删样本，先检查它的总count数是否明显低于其他样本，以及比对率是否正常。一个低质量样本可能拖垮整个差异分析，但一个“看起来离群”的样本也可能是真实的生物学异质性，删样本要记录依据，不能凭感觉。

5. 避坑记录：转录组流程里最常踩的五个坑

5.1 STAR构建索引时内存溢出

现象：运行STAR --runMode genomeGenerate时，报错信息中出现“FATAL ERROR: not enough memory for BAM sorting”或系统直接kill进程。这在人类基因组和大型基因组上是高频问题。

原因：STAR索引构建确实需要大量内存，但更隐蔽的原因是系统可用内存被其他任务占用，或者--runThreadN开得太大导致排序阶段内存需求叠加。

解决：构建索引前用free -g确认可用内存；建议基因组索引构建内存预留不少于基因组大小的25倍；如果确认系统内存充足仍然报错，查看是否有ulimit -v对虚拟内存做了限制。另一个可行方案是把--runThreadN从16降到8，内存占用会显著下降，速度损失可以接受。

5.2 双端数据在featureCounts里计数翻倍

现象：featureCounts输出的总分配数几乎是预计读段数的两倍，差异分析时所有基因的表达量都偏高，并且DESeq2的size factor估计异常。

原因：没有加-p参数。featureCounts默认把双端数据的每一端各当成独立读段计数，一条read pair被统计了两次。这是未加-p时的典型表现，在双端测序项目里几乎是必踩的坑。

解决：加上-p --countReadPairs。注意--countReadPairs只在-p模式下生效。验证方法很简单：用samtools view -c统计BAM里总共有多少个比对记录，和featureCounts summary里的总读段数对比，正常情况下两者应该接近，翻倍则说明参数漏了。

5.3 GTF版本与参考基因组不匹配，比对率正常但计数异常

现象：STAR比对率正常，Log.final.out里unique mapping有85%以上，但featureCounts的分配率只有40%，大量读段被计入“未分配”。

原因：GTF注释和参考基因组来自不同构建版本，最常见的是染色体命名不一致，例如一个用chr1，一个用1。读段能比对到染色体上，但GTF里的注释坐标落不到这些染色体条目上，导致重叠检测失败。

解决：确认基因组和注释来源一致。用grep检查GTF的第一列染色体名，用samtools view看BAM的染色体名，两边直接对比。如果出现不一致，用sed批量替换GTF里的染色体前缀，或者重新下载与基因组同版本的注释，后者明显更稳妥。

5.4 把FPKM矩阵当count矩阵喂给DESeq2

现象：DESeq2运行时报错“strange input”或者“estimateSizeFactors failed”，提示存在负数或非整数值，也可能直接提示所有基因的size factor估计失败。

原因：DESeq2要求输入非负整数的原始count。如果Input序列已经做过FPKM/TPM标准化，包含小数甚至零值，负二项模型的离散度估计完全没有意义。这个错误在新手场景中出现频率极高，因为很多人习惯先在上游输出表达量，再换工具做分析。

解决：回到featureCounts输出，保留下游标准化之前的原始count。DESeq2内置的vst和rlog函数会自动处理标准化，不需要你提前转换。如果实在没法回到原始数据，只能换用接受表达量矩阵的工具，比如limma的voom流程，它至少能处理log2-CPM输入。

5.5 链特异性设错导致一大半基因“不表达”

现象：差异分析结果里出现数量惊人的零表达基因，或者某些已知管家基因的count极低，featureCounts的summary显示“Chimeric fragments”和“Unassigned Ambiguity”偏高。

原因：链特异性参数-s设置错误。反向链建库被当成正向链或当成非链特异性处理时，反义链的读段全部落空，表达量统计的系统性偏差直接摧毁后续分析。

解决：先查建库试剂盒说明确认链方向，再用小样本试验。一个更聪明的验证方法是：随机挑一个已知的单链表达基因，比如某些lncRNA，在BAM里人工检查其正负链信号；或者用工具计算已知链特异基因的表达比例。修改-s参数后重新跑featureCounts，不要用旧结果继续下游。

6. 进阶验证：从差异基因表到通路层面的可解释信号

拿到差异基因表之后，如果只读基因名，这一轮分析的价值会立刻减半。差异基因不是终点，它只是把“哪些基因变了”这个问题交到了你手上，“这些变化共同指向了哪些生物学过程”才是转录组分析真正要交付的结论。这一章介绍的验证路径，是用GSVA把基因表达矩阵转换成通路活性矩阵，从通路层面确认差异信号是否有结构性的生物学解释。

GSVA的思路是：对每个样本、每个基因集，先用核密度估计得到基因表达水平的排序信息，再计算该样本在这个基因集上的富集分数。它不依赖你已经设置了阈值，而是使用全部基因的表达信息，因此比“对差异基因做超几何检验”保留更多弱效应信息。操作上，需要先准备一个经vst或logCPM变换的表达矩阵，再用基因集数据库构建gsva参数对象：

library(GSVA) library(GSEABase) expr <- assay(vsd) # vst变换后的表达矩阵 gs <- getGmt("h.all.v7.5.1.symbols.gmt") gsva_res <- gsva(expr, gs, method = "ssgsea", kcdf = "Gaussian", verbose = FALSE)

method = "ssgsea"是相对稳健的选择，它对单样本计算富集分数，本身就适合样本量较小、组内异质性较高的转录组数据。得到通路活性矩阵后，再用普通的t检验或秩和检验比较两组样本在每个通路上的活性差异。这个结果的解释力比单个差异基因强得多：如果一个通路在处理组的活性分数整体升高，并且这种升高在PCA上也能看到方向一致性，那这个结果就具备生物学可解释性，而不是一堆基因名的随机堆叠。

最后想分享一个习惯：我在每次转录组分析结束时，都会把三条核心证据放在同一个HTML报告里——差异基因火山图、关键通路的GSVA热图、目标基因的表达趋势图。三者能互相印证，才说明流程从数据到结论是自洽的。如果GSVA给出的方向与差异基因列表明显冲突，我不会去修改阈值“凑答案”，而是回头检查样本分组和比对环节。RNA-seq分析没有某个灵光一现的魔法步骤，它的可靠性来自每一步都有日志、每个结论都有交叉验证。希望帮到你。

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