简介转录组测序RNA-seq是研究基因表达调控的重要手段而规范的生信分析流程是获得可靠结果的前提。NGS数据从原始FASTQ文件到差异表达基因需要经历质量控制、序列比对、转录本定量与统计检验等关键步骤每一步的工具选型与参数设置都会影响结论。本文围绕一套完整的RNA-seq转录组分析流程代码系统讲解trimmomatic质控、STAR比对、featureCounts定量以及DESeq2差异分析的实现细节并结合常见错误场景说明如何避免双端数据误统计、参考基因组命名不一致等问题。该流程具备高可复现性与批量处理能力适合生信新手快速上手也能帮助从业者建立模块化的分析思路。1. RNA-seq转录组分析流程代码别再手动拼shell命令了RNA-seq转录组分析听上去就是把测序数据比对到参考基因组上数一数reads但真正跑过一遍的人都知道从fastq到差异基因表格中间每一步都可能翻车。我最近在带A同学做模拟项目X时发现他手动敲命令把双端reads按单端计数结果整个差异表全乱了。这份RNA-seq转录组数据分析流程项目代码最大的价值是把质控、比对、定量、差异分析四段串成一条完整可复现的流程你只需要准备样本表和参考基因组路径就能从原始数据跑到DESeq2的差异表。适合刚入门的生信新手也适合要批量处理样本、不想每次改脚本的从业者。2. 流程全貌为什么把trimmomatic、STAR、featureCounts串成一条snakemake链拿到一批转录组原始数据之后最核心的问题是选择哪套工具链。跑通流程不难难的是每个环节的选型要匹配数据特征和机器配置。这份项目代码把整套流程拆成四个环节质控清洗、比对、定量、差异分析。下面把这套选型逻辑讲透你后面自己在别的项目里换工具时也知道从哪下手。2.1 质控与清洗trimmomatic比cutadapt稳但参数要按测序平台调拿到fastq后的第一件事不是直接比对而是先跑FastQC。FastQC会生成一个html报告重点看Per base sequence quality和Adapter Content两页。如果你发现reads的3端质量分掉到20以下或者出现了大段接头序列那么不清理直接比对会让后续比对器花大量时间去匹配含接头的序列最终uniquely mapped比例低到你不敢相信。质控清洗工具的选择上常见做法是cutadapt配合fastp但这里推荐trimmomatic。原因在于trimmomatic把接头去除和质量滑动窗口合并到一个命令行里双端reads只需要写一条命令cutadapt去接头很干净但滑动窗口修剪需要单独配参数在流程脚本里写起来容易漏项尤其是R1/R2两个文件都要分别指定稍不注意就变成只修了R1。固定在一套流程里trimmomatic的写法更稳定。trimmomatic PE \ ${SAMPLE}_R1.fastq.gz ${SAMPLE}_R2.fastq.gz \ ${SAMPLE}_trim_R1.fastq.gz ${SAMPLE}_unpaired_R1.fastq.gz \ ${SAMPLE}_trim_R2.fastq.gz ${SAMPLE}_unpaired_R2.fastq.gz \ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \ SLIDINGWINDOW:4:20 \ MINLEN:36这行命令里PE表示双端模式后面四个输出文件按顺序对应R1保留、R1未配pair、R2保留、R2未配pair。ILLUMINACLIP的三个数字是seed错配数、palindrome剪切阈值、simple剪切阈值TruSeq3-PE.fa要放在你trimmomatic安装目录下的adapters文件夹里。SLIDINGWINDOW:4:20表示从read头开始每4个碱基算一个窗口窗口平均质量低于20就从这里截断。MINLEN:36是保留最短长度比36还短的read直接丢掉。血泪经验是这些参数不是写死就完事。第一次跑完修剪后要把输出的trim_R1文件再丢进FastQC看一遍如果输出文件大小掉了一半说明SLIDINGWINDOW裁得太狠或者是文库本身质量差。这种情况下我会把窗口质量阈值降到15同时把MINLEN降到30多保留一些reads用于下游分析。2.2 比对环节STAR吃内存但junction识别和稳定性更值得比对器是流程里第一个分岔路。Hisat2基于图比对速度更快、内存需求低适合只有16G内存的机器但它对剪接位点的预测算法在某些新数据上不如STAR稳定尤其是非模式生物。STAR使用后缀数组搜索比对准确率高对跨外显子junction的支持很成熟所以这份项目代码默认选了STAR。你可能担心内存实际上STAR的内存消耗大头在基因组索引建好索引后每个样本比对时峰值内存会降下来。STAR --runMode genomeGenerate \ --genomeDir /data/reference/star_index \ --genomeFastaFiles /data/reference/genome.fa \ --sjdbGTFfile /data/reference/genes.gtf \ --sjdbOverhang 149 \ --runThreadN 16建索引时--sjdbOverhang的值一般是读长减1比如PE150测序就是149。这个参数不是随便填的填0就代表不用剪接位点信息短外显子识别会打折。GTF文件要和参考基因组来自同一版本GRCh38配GRCh38的注释混用不同版本大概率会在下游特征计数时发现基因覆盖率极低。比对阶段的命令要顺便产出排序后的BAMSTAR \ --genomeDir /data/reference/star_index \ --readFilesIn ${SAMPLE}_trim_R1.fastq.gz ${SAMPLE}_trim_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outFileNamePrefix ${SAMPLE}_ \ --runThreadN 8--outSAMtype BAM SortedByCoordinate直接输出coordinate排序BAMfeatureCounts最喜欢这种格式。--quantMode GeneCounts是STAR的附加功能会顺手生成一个ReadsPerGene.out.tab文件这个文件里就是每个基因的raw count可以拿来做快速检查但正式定量仍然走featureCounts保证流程输出统一。如果你手头的机器内存只有16GSTAR建索引时很容易进程被操作系统杀掉。这种场景下我一般会改用Hisat2或STAR的shared memory模式但这份流程代码默认你有一个像样的计算节点。小基因组例外细菌、真菌的基因组索引很小16G内存完全够用。2.3 定量环节featureCounts不搞花哨-p -B -C组合最关键定量是这个流程里最容易被低估的一步。从比对后的BAM直接统计read重叠是经典路线代表工具是featureCounts和HTSeq另一条路线是用kallisto或Salmon做伪比对快、省内存但依赖转录组序列可复现性要差一些。这份流程选了featureCounts原因是它的输入是BAM和GTF输出能明确对应到每个基因号所有统计逻辑都在文档里出了问题容易排查。featureCounts \ -a /data/reference/genes.gtf \ -o ${SAMPLE}_counts.txt \ -p -B -C \ -T 8 \ -t exon -g gene_id \ ${SAMPLE}_Aligned.sortedByCoord.out.bam这里最关键的是-p、-B、-C三个参数。--p表示数据是双端让featureCounts把一条read的两端合并计算一个count-B要求两对read都比对到基因组上才统计-C用于排除包含嵌合片段的reads。很多人在这一步漏掉-p结果双端reads被当成两个单端read统计count数整体翻倍后面所有差异分析都会失真。这是整个流程里出现频率最高的翻车点。参数-t exon指定统计外显子区域-g gene_id指定把外显子计数汇总到基因号。如果你的GTF里用的是transcript_id改成转录本水平也行但差异分析一般先汇总到基因。是否统计multi-mapping reads是另一个值得想清楚的问题。featureCounts默认只统计uniquely mapped加上-M才会把多重比对纳入如果不小心处理反而会引入定位模糊的噪声所以默认不加M。3. 动手复现样本表、脚本、R代码一次跑通这一章我把流程拆成三个可操作步骤先准备好样本表和参考文件再跑主流程脚本最后进入差异分析。每一步输出什么、该检查什么都会写清楚。3.1 样本表与参考文件先定命名规范别让条件因子中间出乱子在跑流程之前先把样本表整理好。这份流程代码用CSV格式第一行是列名后面每行是一个样本sample,condition,fastq1,fastq2 ctrl_1,ctrl,raw/ctrl_1_R1.fastq.gz,raw/ctrl_1_R2.fastq.gz ctrl_2,ctrl,raw/ctrl_2_R1.fastq.gz,raw/ctrl_2_R2.fastq.gz treat_1,treat,raw/treat_1_R1.fastq.gz,raw/treat_1_R2.fastq.gz treat_2,treat,raw/treat_2_R1.fastq.gz,raw/treat_2_R2.fastq.gzsample列必须唯一condition列不要写大小写混合比如Ctrl和treat这种到最后R里做factor时很容易因为变量名匹配不上而报错。如果有批次信息建议在CSV里再加一列batch后面DESeq2设计公式会用到。参考文件方面准备好参考基因组的fasta文件和对应的GTF注释文件另外STAR索引要先建好路径固定下来后面脚本里不会改来改去。3.2 主流程脚本每一步输出都是检查点这里用一份bash脚本作为流程核心最好嵌入到snakemake里管理但如果你刚接触流程调度直接用bash跑也没有问题。脚本以样本名为参数输出目录按样本分隔。#!/usr/bin/env bash set -euo pipefail SAMPLE$1 THREADS8 GENOME_DIR/data/reference/star_index GTF/data/reference/genes.gtf echo [trim] ${SAMPLE} trimmomatic PE \ ${SAMPLE}_R1.fastq.gz ${SAMPLE}_R2.fastq.gz \ ${SAMPLE}_trim_R1.fastq.gz ${SAMPLE}_unpaired_R1.fastq.gz \ ${SAMPLE}_trim_R2.fastq.gz ${SAMPLE}_unpaired_R2.fastq.gz \ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \ SLIDINGWINDOW:4:20 \ MINLEN:36 echo [star] ${SAMPLE} STAR \ --genomeDir ${GENOME_DIR} \ --readFilesIn ${SAMPLE}_trim_R1.fastq.gz ${SAMPLE}_trim_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outFileNamePrefix ${SAMPLE}_ \ --runThreadN ${THREADS} echo [count] ${SAMPLE} featureCounts \ -a ${GTF} -o ${SAMPLE}_counts.txt \ -p -B -C -T ${THREADS} \ -t exon -g gene_id \ ${SAMPLE}_Aligned.sortedByCoord.out.bam echo [done] ${SAMPLE}脚本里set -euo pipefail是必须的意思是任何一个命令失败或者管道里发生问题就立刻退出避免后续步骤在错误状态下继续跑。如果你不习惯至少要保留set -e。STAR的--readFilesCommand zcat用来解压.gz文件如果fastq没有压缩这一行要去掉。注意STAR的输出文件名由--outFileNamePrefix决定这里设置成${SAMPLE}_最终BAM名字就是${SAMPLE}_Aligned.sortedByCoord.out.bamfeatureCounts必须对准这个文件。跑完每个样本先看三样东西修剪后的fastq大小是否异常STAR的Log.final.out里uniquely mapped比例以及featureCounts输出文件里reads assignment summary。如果uniquely mapped低于70%第一步去检查Sample ID是不是搞混了或者参考基因组是不是匹配错了物种。3.3 DESeq2差异分析跑之前先看这三行主流程跑完后把所有样本的counts合并成一个矩阵然后进入R环境。这一步很多新手会直接把featureCounts的输出bind在一起但featureCounts输出的列名是${SAMPLE}_Aligned.sortedByCoord.out.bam不去掉后缀R会把它当成一长列列名后面脚本会卡住。library(DESeq2) # 1. 读入合并后的counts矩阵 counts - read.table(all_counts.txt, header TRUE, row.names 1, check.names FALSE) counts - counts[, grep(_count$, colnames(counts))] colnames(counts) - sub(_count$, , colnames(counts)) # 2. 读入样本信息 samples - read.csv(samples.csv, stringsAsFactors FALSE, row.names 1) samples$condition - factor(samples$condition, levels c(ctrl, treat)) # 3. 构造DESeq2对象 dds - DESeqDataSetFromMatrix(countData round(counts), colData samples, design ~condition) dds - DESeq(dds) res - results(dds, contrast c(condition, treat, ctrl)) write.csv(res, diff_treat_vs_ctrl.csv)第一个注释里的round()很关键DESeq2不接受带小数的count数据如果你在合并时用某种归一化方法先处理过会报错。第二步samples的行名要和counts列名顺序一致DESeq2会根据样本名自动匹配所以这一步错不了。design~condition表示只有处理条件一个因子效果等同于两组比较如果样本表里有batch列改用design~batch condition。contrast三个参数分别是变量名、比较组、参考组顺序写反差异方向就反过来。跑完DESeq后要留意res里很多低表达基因的padj是NA这是DESeq2默认过滤掉低counts基因的结果不是程序崩溃。如果你发现padj列全NA去查design和colData的匹配情况通常是某个组只有一两个样本离散度估计不出来。4. 避坑指南RNA-seq流程跑挂的五个常见现场平时在群里看到最多的问题都是基础配置或者参数漏项不是算法本身出了问题。下面这五条是真实项目里反复出现过的按环节分开写方便你对照排查。4.1 上游质控和比对环节的翻车临时目录与染色体命名第一条FastQC命令没有报错但没有输出html报告。现象是终端显示出fastqc进程完成但结果目录里只有fastqc_data.txt没有fastqc.html。原因是服务器/tmp目录空间紧张或权限不足FastQC在生成html时没有拿到临时文件的读写权限。解决方法是给FastQC显式指定临时目录fastqc -d /data/tmp -o /data/fastqc ${SAMPLE}_R1.fastq.gz并先用mkdir -p /data/tmp建目录。并行跑多个样本时建议每个样本一个独立输出目录避免临时文件互相覆盖。第二条STAR比对后uniquely mapped比例是0%Log文件里也没有报错。现象是Log.final.out显示Uniquely mapped reads数位0或接近0所有reads被标在unmapped: other。最常见的原因是参考基因组fasta的染色体名是1、2、3而GTF文件里是chr1、chr2、chr3STAR比对到的read染色体位置在GTF对应的转录本里找不到结果就是比对率异常低。解决方法是比对前先查两者命名samtools faidx genome.fa查看fasta染色体cut -f1,2 genes.gtf | head查看GTF染色体。如果确实不一致用脚本统一参考基因组命名后重新建STAR索引不要直接在GTF里手动加chr那个坑更大。4.2 下游定量和差异分析环节的翻车参数漏项与因子顺序第三条featureCounts的计数结果是预期的两倍。现象是某个已知基因的表达量明显翻倍所有count都偏大。原因是双端数据没有加-p参数featureCounts默认把mate1和mate2当成两个独立read计数。解决方法是给featureCounts加-p -B -C另外确保输入BAM是coordinate排序。如果你拿到的是STAR默认输出的name-sorted BAM先samtools sort再跑featureCounts否则会看到大量read pairs are not adjacent的警告计数结果同样不可信。第四条DESeq2输出表里pvalue和padj全是NA。现象是差异表能生成但所有p值都是空log2FoldChange也是NA。原因可能有两个colData里列名和design公式不一致比如design写了~batch condition但colData里没有batch列或者某个条件因子忘了设成factor数字被当成连续变量。还有一个高发原因是一个组只有一个样本DESeq2没法估计离散度全部样本被当作离群值排出。解决方法是先用resultsNames(dds)看看实际可用的系数名再决定用results(dds, namecondition_treat_vs_ctrl)还是contrast写法。如果组内只有单样本换成edgeR或limma-voom并调整min.group.size。第五条流程跑一半磁盘爆满无法续跑。现象是STAR或featureCounts报No space left on device但前面的中间文件都还在。原因是STAR输出的BAM虽然排序但并未自动删除中间SAM文件加上trimmomatic的unpaired文件、临时fastq堆积很容易把家目录塞满。解决方法是脚本里在BAM生成后删除${SAMPLE}_Aligned.sortedByCoord.out.sam并在每跑完一个样本后用df -h看剩余空间。另外bash脚本可以加断点续跑逻辑if [ -f ${SAMPLE}_Aligned.sortedByCoord.out.bam ]; then echo skip star; else STAR ...; fi。以后重跑同一批样本时已经完成的步骤会自动跳过省时间也省磁盘。5. 进阶技巧用multiqc和PCA给结果上两道保险主流程跑完拿到差异表先别急着去注释那些p0.05的基因先做两个低成本的检查。第一道保险是整合质控报告第二道保险是看样本聚类。这两步能过滤掉一批因为样本混淆、批次效应造成的假阳性。第一道保险用multiqc。它能把所有FastQC报告、STAR的Log.final.out、featureCounts的summary汇总成一个html页面整体看一遍uniquely mapped比例、3端质量分布、基因计数分配到各区域的比例。如果所有样本都在80%左右只有某一个样本只有40%那基本可以断定这个样本出了问题建议在后续分析里直接剔除而不是保留下来让差异分析强行补偿。multiqc results/ -o results/multiqc/ samtools quickcheck -v results/star/*.bam echo all bam oksamtools quickcheck是验证BAM完整性的最轻量方式如果BAM文件在比对中途被截断它会返回错误。这个命令在批量分析里很好用配合multiqc的bam统计可以提前发现问题。不要等到DESeq2报错再回头查文件。第二道保险是用PCA验证样本关系。从counts矩阵做VST变换后直接画PCA看处理组和对照组是否各自成群。vsd - vst(dds, blind FALSE) plotPCA(vsd, intgroup condition)这里blindFALSE的意思是让DESeq2在估计离散度时考虑分组信息PCA的结果更能体现组间差异。如果画出来的PCA里对照组和处理组完全交错而批次信息又是按另一列分明说明设计公式里要加入批次变量。另一种常见情况是同一处理组的几个样本分散得很远这时先回fastqc报告看是不是测序深度差异太大再考虑是否要用更严格的过滤处理。这两道保险花不了十分钟却能避免你在后续GO分析里埋头扒一堆假阳性。从那以后我每次拿到新一批RNA-seq数据都会强制先跑一遍multiqc和PCA确认样本调用和分组有效才敢看diff.csv希望帮到你。本文还有配套的精品资源点击获取