简介这是一份面向生物信息学初学者的RNA-seq转录组数据分析完整流程项目代码覆盖从原始测序数据到功能富集的核心环节。内容以两个不同组织的六个样本为例串联FastQC质控、Trimmomatic剪接、HISAT2比对、featureCounts定量、DESeq2差异表达分析并包含GO与KEGG富集分析及差异基因火山图绘制方法。压缩包共3个文件以可运行的inscode项目文件为主附带了项目说明html及gitignore配置文件整体约8KB便于直接查看代码结构和启动演示。已有518人学习下载适合希望系统掌握转录组上游分析流程并动手实操的入门学习者。资源中的数据经过特殊修改可在学习过程中放心使用帮助理解每一步工具参数设置、结果解读与常见分析思路为后续独立开展类似项目打下基础。1. RNA-seq转录组分析为什么完整流程值得跑通做转录组的人应该都遇到过这种场景拿着测序公司返回的FASTQ文件脑子里很清楚要的是差异表达基因列表但真到动手时才发现从原始数据到火山图之间隔着质控、过滤、比对、定量、差异分析、富集整整六道工序。我去年帮某课题组复现一套RNA-seq分析时发现他们三个人分别跑了三套流程比对软件一个用STAR一个用Hisat2定量一个用featureCounts一个用salmon最后差异基因列表对上的不到一半。问题不是工具不行而是流程里每一步的选型和参数会层层传导任何一处不严谨都会让结果变成一纸空谈。这份项目代码的价值就在于把一整套可靠的RNA-seq转录组分析流程固化成了可复现的脚本新手拿过去能按顺序跑通熟手也能在关键步骤上看到有据可依的参数配置。2. 全流程选型与路径从FASTQ到差异基因列表的六个环节2.1 流程的六个环节与每步的产出物RNA-seq分析无论数据来自人、小鼠还是其他物种主干路径大致是固定的。搞清楚每个环节输入什么、输出什么、验证什么比急着敲命令更重要。我按实际执行顺序把这条链路拆成六个环节环节输入输出需要验证的点质量评估原始FASTQ质量报告HTML碱基质量、GC含量、接头残留过滤清洗原始FASTQ干净FASTQ过滤后数据量是否够用比对干净FASTQ 参考基因组索引BAM文件比对率、重复率是否正常定量BAM GTF注释基因表达计数矩阵基因数是否在预期范围差异分析计数矩阵 样本分组差异基因列表表达变化方向是否和实验一致下游解读差异基因列表富集结果、可视化图通路结果能否对应生物学背景我在实际复现这份项目代码时它把前四个环节固化成了Shell脚本第五和第六环节用R脚本承接。这个拆分思路我认为是对的上游是重计算任务适合写在Shell里丢给服务器跑下游是交互式探索放在R里调整阈值和绘图更方便。值得注意的一个细节是每个环节之间通过中间文件衔接——过滤后的FASTQ、比对完的BAM、定量的count矩阵这样任何一个环节出问题都可以只重跑该环节不需要从头再来。2.2 软件选型为什么是STAR加featureCounts加DESeq2这套流程选了STAR做比对、featureCounts做定量、DESeq2做差异分析这个组合反映了目前最主流且相对稳定的技术路线。STAR的核心优势是剪接比对能力强它把基因组切分成suffix array索引搜索reads时能同时处理跨外显子的junction read对真核生物转录组来说这是刚需。featureCounts的优势体现在两点一是速度极快二是它对多映射read的处理策略是默认丢弃这符合大多数差异分析对保守性的要求。DESeq2选它不是因为功能最花哨而是它基于负二项分布建模对生物学重复数较少的实验设计比如三组重复有比较稳健的表现而且它对低表达基因的收缩估计能有效减少假阳性。这套组合里需要预先确认的是STAR的索引构建需要参考基因组FASTA和GTF注释文件featureCounts只接受BAM格式输入DESeq2的输入需要是原始整数计数而非归一化后的值。如果手头数据来自链特异性建库三条工具都要注意链方向参数的一致性STAR里是--outSAMstrandFieldfeatureCounts里是-sDESeq2不直接感知链信息但定量结果已经受前两者影响。这套流程把这几个参数的默认选择写在配置里跑之前花十分钟确认一遍比跑完再返工划算得多。3. 上游实操质量控制、STAR比对与featureCounts定量的参数细节3.1 环境准备与数据质量控制先用conda把环境固定住做生信分析的人最怕的就是软件版本漂移。两个月前能跑的脚本两个月后因为某个依赖库升级就报错这是最常见的翻车方式。这套流程代码里用conda环境锁定了所有依赖的版本我复现时按它提供的environment.yml直接建环境conda env create -f environment.yml conda activate rnaseqenvironment.yml里把fastp、STAR、featureCounts、samtools、multiqc这些工具的精确版本号固定住了这样做的好处是整个环境可复现换一台机器也能得到一致结果。我一般还会在激活环境后运行一次conda list | grep -E star|featurecounts|fastp确认关键工具确实加载的是指定版本。质控环节我重点看fastp参数。下面这个命令是这套流程里的核心过滤逻辑fastp \ -i sample_1.fastq.gz \ -I sample_2.fastq.gz \ -o clean_1.fastq.gz \ -O clean_2.fastq.gz \ --detect_adapter_for_pe \ --trim_poly_g \ --cut_front \ --cut_front_window_size 1 \ --cut_front_mean_quality 3 \ --cut_tail \ --cut_tail_window_size 1 \ --cut_tail_mean_quality 3 \ --cut_right \ --cut_right_window_size 4 \ --cut_right_mean_quality 15 \ --qualified_quality_phred 15 \ --length_required 36 \ --compression 6 \ --thread 8 \ --json sample.json \ --html sample.html逻辑说明--detect_adapter_for_pe让fastp自动识别并切除双端测序的接头序列这比手动指定接头序列更省心因为不同建库试剂盒的接头可能存在差异。--trim_poly_g针对NextSeq平台常见的polyG尾进行切除如果测序平台是HiSeq或Novaseq环境不同这条建议保留。三段式质量裁剪分别处理reads的头部、尾部和滑动窗口之所以分成三部分是因为reads不同位置的碱基质量分布差异很大。参数说明qualified_quality_phred 15表示碱基质量低于Q15会被当作低质量处理这是比较常见的阈值选择如果数据整体质量很好可以放宽到Q20但你会在reads保留率上付出代价。length_required 36过滤掉裁剪后长度不足36bp的reads这个阈值跟后续STAR比对的最小匹配长度有关设得太高会损失真实现有转录本的信息。跑完后我通常会看一眼HTML报告里的reads通过率和平均碱基质量通过率低于70%的话建议回到原始数据去确认建库质量是否有问题。3.2 STAR比对索引构建与两遍比对的核心参数STAR是这套流程里最吃内存也最需要细调参数的一步。第一步是构建参考基因组索引代码如下STAR \ --runMode genomeGenerate \ --genomeDir ./star_index \ --genomeFastaFiles genome.fa \ --sjdbGTFfile genes.gtf \ --sjdbOverhang 149 \ --runThreadN 16 \ --genomeSAindexNbases 14 \ --limitGenomeGenerateRAM 31000000000sjdbOverhang的值应该设置为读长减1这套流程的测试数据是150bp双端测序所以配置里写149。这个值影响剪接位点数据库的构建如果读长变了必须跟着改。genomeSAindexNbases 14是基因组suffix array索引的规模参数基因组越大这个值越大人基因组默认14够用基因组特别大时可以适当调大到15但我不建议随便改因为改得不好会显著增大内存占用。limitGenomeGenerateRAM设到31GB是为了确保索引构建过程不被内存限制中断。比对阶段用的是两遍比对策略第一遍生成剪接位点信息第二遍利用这些信息做更精准的比对STAR \ --runMode alignReads \ --genomeDir ./star_index \ --readFilesIn clean_1.fastq.gz clean_2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outSAMunmapped Within \ --outSAMattributes Standard \ --outFilterMismatchNmax 2 \ --outFilterMultimapNmax 10 \ --outFilterScoreMinOverLread 0.33 \ --outFilterMatchNminOverLread 0.33 \ --alignSJoverhangMin 8 \ --alignSJDBoverhangMin 1 \ --runThreadN 16 \ --outFileNamePrefix sample_逻辑说明outSAMtype BAM SortedByCoordinate直接输出按坐标排序的BAM文件省掉后续再单独跑一次samtools sort。readFilesCommand zcat告诉STAR输入文件是gzip压缩格式需要先解压读取如果不加这个参数STAR会把压缩文件当普通文本处理导致比对率惨不忍睹。outFilterMismatchNmax 2限制每个read最多只能有2个错配这个阈值对一致序列数据来说是合理的如果样本是肿瘤或含有较高变异率的细胞系可以放宽到3到4但需要警惕错配容忍度过高带来的多定位问题。参数说明outFilterMultimapNmax 10允许一个read最多比对到基因组的10个位置超过10个位置的read会被标记为多映射而丢弃这个阈值实际上决定了最终比对率里包含还是排除重复序列区域的信号。alignSJoverhangMin 8要求剪接位点两侧至少各有8bp的比对锚定长度这个值直接影响junction read的识别质量。两遍比对的好处是第一遍把新发现的junction位点合并进索引第二遍时这些位点已经可以作为已知剪接位点用于比对对非模式物种或注释不完整的基因组来说差距尤其明显。3.3 featureCounts定量从BAM到计数矩阵featureCounts这一步相对简单但它决定了最终输入DESeq2的计数矩阵长什么样。命令里有个值得注意的细节是-s参数它表示链特异性模式featureCounts \ -a genes.gtf \ -o counts.txt \ -t exon \ -g gene_id \ -p \ -s 2 \ -T 8 \ -B \ -C \ sample_1.bam sample_2.bam sample_3.bam逻辑说明-t exon -g gene_id指定了计数单元——exon是GTF注释文件第三列的特征类型gene_id是第九列的属性名意思是把所有比对到同一个基因所有外显子上的read累加起来作为该基因的表达量。-p表示输入的是paired-end数据featureCounts会把成对的read当作一个片段来计数避免同一片段的正反两条read被重复计算。-B要求片段的两端都能成功比对到参考基因组上-C要求一个片段的两端比对位置一致不允许一端比对到一条染色体、另一端比对到另一条染色体这两个选项都用于提高计数质量。参数说明-s 2适用于stranded reverse文库比如dUTP法建库这是目前主流链特异性建库方式如果数据是非链特异性文库这里应该改成s 0设错会导致一半的read计数丢失表现是最终计数矩阵里大量基因表达量接近0。-T 8是线程数这个参数和STAR的线程数保持一致即可。跑完后我建议看一眼counts.txt.summary文件它会列出每个样本成功分配的read数量如果所有样本的分配率都非常低先回头检查BAM的染色体命名是否与GTF一致这是最常见的静默错误来源。4. 差异表达分析实操DESeq2的输入、设计公式与结果解读4.1 DESeq2的输入格式与样本设计拿到featureCounts输出的计数矩阵后接下来桌面切换到R环境。DESeq2对输入数据有一个硬性要求必须是整数计数不能是FPKM、TPM这些已经归一化过的值。这是因为DESeq2内部会对原始计数做自己的归一化处理输入归一化值等于进行二次归一化结果会失真。读取数据并构建DESeq2对象的代码如下library(DESeq2) count_data - read.table(counts.txt, headerTRUE, sep\t, row.names1, skip1) count_matrix - count_data[, 6:ncol(count_data)] colnames(count_matrix) - c(ctrl_1, ctrl_2, ctrl_3, treated_1, treated_2, treated_3) sample_info - data.frame( row.names colnames(count_matrix), condition factor(c(ctrl, ctrl, ctrl, treated, treated, treated), levels c(ctrl, treated)) ) dds - DESeqDataSetFromMatrix( countData count_matrix, colData sample_info, design ~ condition )逻辑说明skip1是因为featureCounts输出文件的第一行是命令参数信息真正的内容从第二行开始。count_data[, 6:ncol(count_data)]的作用是把前五列注释信息Geneid、Chr、Start、End、Strand去掉只留下各样本计数列。设计公式~ condition告诉DESeq2我们只关心处理组与对照组之间的差异如果实验设计包含批次信息可以在design里加上批次变量比如写~ batch condition。参数说明factor(c(...), levels c(ctrl, treated))这一行的顺序决定了后续差异分析中的比较方向。DESeq2默认把factor的第一个水平作为参考组这里ctrl排在第一个所以最终结果里log2FoldChange的正负号表示处理组相对于对照组的变化方向。这个顺序非常关键写反了会导致所有基因的上调和下调方向互换而且不会产生任何报错属于最危险的错误之一。4.2 差异基因筛选标准与结果表的列解读构建好DESeq2对象后一行命令完成差异分析dds - DESeq(dds) res - results(dds, alpha 0.05, lfcThreshold 0) res_df - as.data.frame(res) res_df$gene - rownames(res_df) sig_genes - subset(res_df, padj 0.05 abs(log2FoldChange) 1)逻辑说明DESeq(dds)内部会依次完成大小因子估计、离散度估计和负二项分布模型拟合这一步不需要人工干预但如果样本数量太少会在此处报错至少需要两个生物学重复。results()里设置alpha 0.05表示显著性阈值lfcThreshold 0表示不设置表达变化倍数的最小门槛。筛选差异基因时我习惯在padj 0.05之外加一个abs(log2FoldChange) 1即表达量变化超过两倍的标准这个组合是当前文献中最常见的默认标准。结果表中每一列的含义需要解释清楚baseMean是该基因在所有样本中的平均归一化计数log2FoldChange是变化倍数以2为底的对数正值表示处理组上调负值表示处理组下调。lfcSE是log2FoldChange的标准误stat是Wald检验的统计量pvalue是未校正的P值padj是经过Benjamini-Hochberg方法校正后的P值。实际筛选时永远以padj为准而不是pvalue因为当基因数量在两万左右时单纯按pvalue 0.05筛选会出来上千个假阳性基因。排序和导出差异基因res_ordered - sig_genes[order(sig_genes$padj), ] write.csv(res_ordered, file DE_genes_padj0.05_log2FC1.csv, quote FALSE)我通常会导出一份中英文注释信息附在后面方便直接用Excel查看。还有个习惯是顺手跑一下plotMA(res, mainMA plot)看整体分布如果图形呈对称喇叭状说明数据质量可靠如果出现大面积偏移或形状异常说明定量环节可能有问题。5. 避坑排查五个让分析重跑的常见问题与处理办法5.1 参考基因组版本与GTF不匹配的沉默错误现象STAR比对率正常featureCounts也正常跑完但最终计数矩阵里大量基因的count值为0或异常低。 原因这是我见到的频率最高的翻车现场。参考基因组FASTA和GTF注释文件来自不同版本或不同来源比如基因组是ENSEMBLrelease 110的GTF却是GENCODE的两者对同一基因的ID命名方式不同featureCounts按gene_id聚合时自然匹配不上。 解决下载参照文件时把基因组和GTF的版本号记录下来尽量从同一版本同一来源下载。用下面的命令快速核对染色体命名是否一致grep ^ genome.fa | head -5 cut -f1 genes.gtf | sort -u | head -5如果一边是chr1另一边是1说明命名不兼容。需要一个脚本把GTF里的染色体名换成基因组里的写法或者反过来。从那以后我拿到新参考文件第一件事就是跑这两行命令花不了一分钟但能省掉一整天返工。5.2 双端测序比对率低的真实原因现象STAR比对率只有60%出头看着明显偏低但又不完全算失败。 原因这种情况最常见的原因是两条reads的序列方向信息在比对时没有被正确处理或者reads里残留了未切除干净的接头序列。另一个容易忽视的原因是--outFilterMismatchNmax设置得太严格对变异率较高的样本来说错配容忍度不够导致大量reads因为两三个错配就被判定为未比对。 解决先检查clean数据里接头残留是否已经清除fastp的HTML报告里有明确指标。再把outFilterMismatchNmax从2放宽到4如果比对率从60%提到85%以上说明原数据变异率比预期高放宽阈值是合理的。最后用samtools stats检查BAM文件里paired properly的比例这个数字应该在90%以上才算正常。5.3 STAR内存不足导致的任务卡死现象--runMode genomeGenerate或比对阶段运行到一半进程被OOM kill直接杀掉。 原因STAR对内存的需求超出预期尤其是构建索引时--genomeSAindexNbases设置过大或者--limitGenomeGenerateRAM低于实际需要。人基因组索引构建通常需要30GB左右内存而比对阶段每个线程也要消耗几GB。 解决在提交任务时直接申请足够的内存不要用服务器默认配置。如果基因组很大把genomeSAindexNbases降到12到13也可以显著降低内存占用代价是比对速度稍慢。还可以用samtools view -H检查一下BAM文件头部里的STAR版本信息确认比对时用的索引和当前基因组目录一致。5.4 featureCounts计数丢失多映射与配对选项检查现象featureCounts跑完日志显示某个样本成功分配的片段数只有二三十万明显少于BAM文件里的总片段数。 原因这种情况通常是一个或多个选项没设对。-s链方向设反会让一半reads因链信息不符被丢弃-C选项会把那些一端比对到线粒体另一端比对到核基因组的嵌合片段全部排除这些片段在很多样本里占比不小。 解决先把-s 2改成-s 0跑一次如果计数大幅回升说明是链方向问题。如果改了-s还是没改善再检查BAM里比对到线粒体上的比例用samtools idxstats sample.bam | grep MT看数值线粒体reads比例超过5%时需要单独处理常见做法是在比对前过滤掉线粒体序列或者专门加一条线粒体染色体参数。5.5 DESeq2设计公式与因子水平的方向陷阱现象差异分析结果里log2FoldChange方向整体和实验预期相反。 原因factor的levels顺序决定参考组而参考组决定了变化方向。比如实验是处理组vs对照组对照组必须是参考组但很多人写代码时习惯把处理组写在前面导致每个基因的上下调方向全部反转。这个错误不报错、不警告如果不拿已知marker基因验证方向很容易一路错到论文里。 解决运行差异分析之前用下面两行代码主动确认参考组colData(dds)$condition - relevel(colData(dds)$condition, ref ctrl) results(dds)$log2FoldChange在写文章之前选几个你明确知道该上调的基因比如经典marker基因看一眼log2FoldChange是不是正的。从那次之后我形成了个条件反射——每一组新的分析跑完DESeq2先看参考组定义对不对再往下走绝不动摇。6. 进阶用法把流程封装成脚本并验证批次效应6.1 流程封装与参数配置分离跑通一遍只是开始真正能提升效率的是把整套流程封装成一份可用配置文件驱动的脚本。我复现这套流程代码后做的最重要一件事是把所有变量抽出来放进config.sh脚本主体只负责按顺序执行。配置文件大致长这样GENOME_DIR/path/to/star_index GTF_FILE/path/to/genes.gtf FASTQ_DIR/path/to/fastq_files OUTPUT_DIR/path/to/output THREADS16 READ_LENGTH150 STRANDNESS2 REFERENCE_GROUPctrl TREATMENT_GROUPtreated然后在主流程脚本里通过source config.sh方式导入这些变量后续所有命令都用$THREADS、$STRANDNESS这样的变量代替硬编码。这样做的价值在于换一批数据时只需要改一个配置文件而不是钻进脚本里逐行替换路径。对批次效应明显的项目我还会在config.sh里额外加两个字段BATCH变量和CONDITION变量DESeq2的design就可以写~ batch condition把批次信息纳入统计模型。6.2 用PCA验证样本整体关系拿到计数矩阵后不要急着跑差异分析我现在的习惯是先做一次PCA快速看看全局。这段R代码是这套流程里验证批次效应最直接的方法library(DESeq2) vsd - vst(dds, blind TRUE) pca_data - plotPCA(vsd, intgroup c(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_bw()逻辑说明vst()对计数矩阵做方差稳定变换目的是让高表达基因和低表达基因的方差变得可比这一步对PCA和样本间相关性计算很重要。blind TRUE表示变换时不考虑分组信息避免在探索性分析时引入实验设计的先验影响。plotPCA返回PC1和PC2坐标把它们画出来就能直观看出各组的聚类情况。判断维度有三条第一同一组的样本应该聚在一起如果对照组里有一个样本明显游离在外说明它可能是批次异常样本需要回到原始数据检查第二处理组和对照组之间沿着PC1方向分开说明处理效应非常强差异分析结果大概率干净可靠第三如果PCA图上样本完全按批次聚类而不是按处理组聚类说明批次效应占据了主导地位这种情况下DESeq2的design里必须加入batch项否则结果会被批次效应污染。这是我从一次失败项目中吸取的教训——当时有一组RNA-seq数据由于建库时间不同PCA图里两个批次分别聚成两团起初没在意后来差异基因列表里全是批次差异而非处理差异。从那以后我每次拿到数据都强制把PCA走一遍确定样本关系正常才继续往下做希望帮到你。本文还有配套的精品资源点击获取