简介一套面向生物信息学与转录组初学者的RNA-seq分析项目代码覆盖从原始测序数据到差异结果解读的完整流程包括数据质控、序列比对、定量表达、差异表达分析以及GO/KEGG功能富集和火山图绘制等核心环节。资源以两个不同组织共六个样本为具体示例串联FastQC、HISAT2、DESeq2等常用工具每个步骤均有可运行的代码逻辑能够帮助零基础学习者理解参数设置、中间文件处理与分析思路快速建立转录组分析的整体框架。压缩包体积约8KB仅含3个文件以HTML说明、inscode工程及gitignore配置为主结构精简便于直接查看项目全貌目前已有518人学习浏览获得一定关注。所有演示数据均经特殊修改仅用于教学参考读者可借此理解差异基因筛选、富集分析与结果可视化的具体实现并能举一反三迁移到自有转录组数据中。1. RNA-seq转录组数据分析流程一句话说清这条流水线能帮你解决什么RNA-seq转录组数据分析流程本质上是一条从下机数据到差异基因列表的加工流水线原始FASTQ先做质控过滤掉接头和低质量read然后比对到参考基因组再在基因或转录本层面计数得到表达矩阵最后做差异分析和功能富集。很多第一次接触的人以为难点在差异分析但实际上让整个流程翻车的往往是上游某个小环节参考基因组版本选错、GTF注释文件对不上、定量时参数口径不一致。这条流程适合三类人手上有转录组数据但不知道标准处理链路是什么的初学者已经在用别人的流程但想搞清楚每一步在干什么的从业者以及想自建一套稳定、可复现流程用于多个项目的课题组。看完这篇文章你至少能把从FASTQ到表达矩阵的完整套路搭起来并且知道每一步出问题之后从哪里排查。2. 流程设计与工具选型先想清楚分析链路再开跑数据2.1 一条标准流程的六个环节与数据流转录组数据分析流程的骨架其实是稳定的不同项目之间差异只体现在细节参数上。最常见的做法是把流程拆成六个环节原始数据质控、数据过滤、比对、定量、差异分析、下游注释与可视化。这六个环节不是随意排列的每一步的输出格式是下一步的输入文件格式的衔接是否顺畅直接决定了你后续要花多少时间填坑。第一个环节拿到的是下机FASTQ这是最原始的序列数据里面混杂着接头序列、低质量碱基和长度过短的read。第二个环节做过滤清洗产出的还是FASTQ但质量分布会明显改善。第三个环节把clean后的read比对到参考基因组上产出SAM或BAM文件这里面记录的是每条read在基因组上的坐标和比对质量。第四个环节把BAM里的比对信息汇总到基因或转录本上得到counts矩阵行是基因名、列是样本名。第五个环节对counts做统计建模筛出差异基因。最后一个环节做可视化、GO/KEGG富集把差异基因列表转化成生物学结论。我见过不少项目在第二个环节就出了问题质控报告显示接头率很低、测序质量也不错于是直接跳过过滤步骤。这在某些情况下确实可行但一旦某个样本出现了污染或接头残留后面所有样本都要重新跑返工成本反而更高。常规建议是每个样本都过一遍质控过滤之后再看一眼报告确认数据形态正常再进入比对。2.2 工具选型比对与定量的不同路线怎么定工具选型是流程设计里最容易被低估的一步。很多人直接拿现成流程里的工具用不做选择但如果替换输入物种或者数据规模变大性能表现会有明显差距。比对环节目前被用得最多的是STAR和HISAT2两者走的是完全不同的索引策略。STAR采用拼接式比对占用内存较大但速度非常快长read和跨外显子read的比对能力强是目前RNA-seq的主流选择。HISAT2用层级图索引内存占用只有STAR的三分之一左右适合内存受限的服务器在剪接位点识别上表现也不错。如果做的是small RNA或扩增子测序则可以选择bowtie2这类局部短序列比对工具没必要上STAR这种重型方案。定量环节Alignment-based方式用featureCounts或HTSeq-count输入BAM文件、输出基因计数。这种方式的好处是中间产物好排查比对率低、多比对、ambiguous都能在BAM层面看到原始原因。转录本定量的工具salmon和kallisto走的是伪比对路线不需要完整的基因组索引文件直接用转录组索引做快速定量速度快几倍但后续做异构体分析时需要注意兼容性问题。我的习惯是如果项目只需要差异基因层面用STARfeatureCounts简单直接、稳定、易排查如果项目需要异构体层面或数据量很大考虑salmon。差异分析环节DESeq2是最稳妥的选择它对小样本量、低表达基因和离散度估计做了比较充分的处理edgeR计算速度快对经典的两组比较也很可靠limma-voom则适合微阵列经验背景、想用线性模型框架的人。做常规转录组差异分析我一般直接用DESeq2它的方差稳定化变换对后续PCA和热图也比较友好。比对工具索引策略内存占用适用场景定量搭配STAR拼接式索引高30G左右RNA-seq标准流程、跨外显子featureCountsHISAT2层级图索引低内存受限服务器、转录组featureCountsbowtie2FM索引低small RNA、扩增子featureCountssalmon/kallisto伪比对索引低大规模转录本定量自带定量输出2.3 参考基因组与注释文件版本对应关系最容易翻车参考基因组的选择不只是一个下载哪个文件的问题它决定了后面所有比对和定量的正确性。常见的人选来自UCSC、Ensembl和GENCODE三个来源同一个物种、同一个版本号不同来源的染色体命名习惯都有差异。UCSC习惯在染色体名前加chr比如chr1、chr2线粒体写chrMEnsembl默认不加chr线粒体写MT。这个差异会直接导致后续比对和注释对不上比对率看起来正常但featureCounts定量出来全是零。注释文件的格式也值得区分。GTF和GFF3都用来描述基因结构但字段定义不同featureCounts和STAR对GTF的兼容性最好优先使用GTF格式。还有一点容易被忽略GTF文件版本必须和基因组版本完全匹配。比如参考基因组下载的是GRCm39注释文件却用了GRCm38的GTF坐标系统差了几个碱基到几十个碱基跨外显子计数会显著错乱。我的操作建议是一个项目内固定使用某一个来源的基因组和配套注释最好在流程目录里写清楚来源、版本号和下载日期不要凭文件夹名可能是hg38来判断。版本对应关系对初学者来说几乎是个黑匣子但一旦建立起基因组版本注释版本染色体命名三件事必须一致的认知很多神秘问题就自动消失了。3. 从FASTQ到表达矩阵STARfeatureCounts最小复现流程3.1 目录结构与环境准备动手之前先把目录结构定下来这一步虽然简单但能避免后期非常多的混乱。我的经验是每类中间产物放一个目录命名规则全流程统一至少在项目内不出现最终版最终版2这种文件夹名。mkdir -p {raw,clean,bam,counts,results,scripts} cd your_project_dir conda create -n rna-seq python3.10 -y conda activate rna-seq conda install -c bioconda fastp star subread -y目录里raw放原始FASTQclean放过滤后的FASTQbam放比对结果counts放计数矩阵results放差异分析输出scripts放后续要写的R脚本。conda环境名用了rna-seq一次性把fastp、star、subread装好。subread是featureCounts所属的软件包装它就能直接用featureCounts命令。3.2 质控与过滤fastp的常用参数配置fastp是目前用得最多的快速质控工具单条命令同时完成质量报告、接头过滤、低质量碱基修剪和read长度过滤。对于双端测序数据常见的最小命令如下fastp -i raw/sample1_R1.fastq.gz \ -I raw/sample1_R2.fastq.gz \ -o clean/sample1_R1.fastq.gz \ -O clean/sample1_R2.fastq.gz \ --detect_adapter_for_pe \ --trim_front1 4 \ --qualified_quality_phred 20 \ --length_required 36 \ --thread 16 \ -j clean/sample1_fastp.json \ -h clean/sample1_fastp.html参数说明里值得注意几个细节--detect_adapter_for_pe让fastp自动识别双端数据的接头序列省去手动指定接头序列的步骤--trim_front1 4把R1的5端前4个碱基统一修剪掉这是为了去除RNA建库过程中可能引入的引物残留偏倚--qualified_quality_phred 20表示质量值低于Q20的碱基会被处理这个阈值在常规转录组项目里是够用的--length_required 36表示过滤后长度小于36bp的read直接丢弃保留太短的read只会增加比对的噪音。输出文件里JSON和HTML是fastp生成的报告JSON的机器可读性更好后续做多样本质控汇总时可以直接用脚本解析。这一步跑完后建议不要只看一眼过滤率就继续而是抽查一两个样本的clean FASTQ大小确认输出文件不为空、大小合理。3.3 比对STAR索引构建与两遍比对STAR比对前必须先给参考基因组构建索引。索引只需要建一次后续所有样本都用同一个索引目录。构建索引时一定要同时提供基因组序列和GTF注释STAR会利用注释文件里的剪接位点信息来提高比对准确率。mkdir -p star_index STAR --runMode genomeGenerate \ --genomeDir star_index/ \ --genomeFastaFiles reference_genome.fa \ --sjdbGTFfile reference_annotation.gtf \ --sjdbOverhang 149 \ --runThreadN 16--sjdbOverhang的取值有讲究推荐设置成read长度减1。如果你用的是150bp双端测序这里就写149如果是100bp写99。这个值影响STAR对可变剪接位点的识别填错不会报错但会轻微降低跨外显子read的比对准确率。构建索引是一个比较耗时的过程人类基因组需要约30GB内存和二三十分钟属正常现象。比对阶段我常用的是两遍比对模式这是一个很多人忽略但很实用的设置。第一遍先快速比对找到所有样本共同的新剪接位点第二遍使用这些新剪接位点构建新的索引再重新比对一次。对数据质量一般或物种注释不完善的场景这能显著提升跨剪接位点的比对率。STAR --genomeDir star_index/ \ --readFilesIn clean/sample1_R1.fastq.gz clean/sample1_R2.fastq.gz \ --readFilesCommand zcat \ --outSAMtype BAM SortedByCoordinate \ --outSAMunmapped Within \ --outSAMattributes NH HI AS NM MD \ --twopassMode Basic \ --runThreadN 16 \ --outFileNamePrefix bam/sample1_--readFilesCommand zcat告诉STAR输入文件是gzip压缩格式不加这个参数STAR会直接把压缩文件当普通文本读然后报错。--outSAMtype BAM SortedByCoordinate输出按坐标排序的BAM文件这个格式是featureCounts直接能用的省掉了samtools sort这一步。--outSAMattributes里的NH HI AS NM MD是后续做变异检测或可视化时比较常用的标签保留它们不会漏掉信息。--outSAMunmapped Within把未比对read写入BAM文件而不是单独输出这样BAM文件更完整samtools flagstat统计结果也更有参考价值。比对完成后先看一眼每个样本目录下的Log.final.out文件里面有summary mapping rate之类的结果。比对率低于70%时建议停下来排查不要带着低质量比对结果往下走。3.4 定量featureCounts生成表达矩阵featureCounts的输入是排序后的BAM文件输出是一个counts矩阵。路径上STAR已经直接生成了sorted BAM所以不需要额外排序可以直接喂给featureCounts。featureCounts -T 8 -p --countReadPairs \ -a reference_annotation.gtf \ -o counts/expression_counts.txt \ bam/sample1_Aligned.sortedByCoord.out.bam \ bam/sample2_Aligned.sortedByCoord.out.bam \ bam/sample3_Aligned.sortedByCoord.out.bam-P是双端测序必须开的选项配合--countReadPairs表示按fragment而不是按read来计数。一对paired-end read如果比对到同一个基因算做一个fragment这样可以避免同一段RNA被两端的read各计一次导致计数翻倍。不开-p、用了-p但没写--countReadPairs、或者-p和--countReadPairs混用是初学者最常见的手误。具体可以看输出文件末尾的运行日志里面的Multi-mapping、Unassigned等分类计数能帮你判断参数是否正确。输出文件counts/expression_counts.txt里第一列是Geneid后面每列是一个样本的counts。第一行有注释信息用R读取时需要skip1。文件里还会显示unassigned的统计信息这些统计信息建议保留一份副本它们是后续排查比对质量问题的重要依据。4. 常见问题与排查报错之外更隐蔽的定量结果异常4.1 比对率正常但基因计数全为零染色体命名不一致现象STAR的summary比对率在90%以上看起来一切正常但featureCounts输出矩阵里几乎所有样本的counts都是零或接近零。原因参考基因组的FASTA和GTF注释文件来自两个不同的数据源。最常见的是参考基因组用了UCSC版本染色体命名是chr1、chr2GTF注释文件却来自Ensembl命名是1、2反过来也一样。染色体名对不上STAR比对到的坐标无法被featureCounts映射到基因上于是表现为比对成功但定量还是零。解决检查GTF文件里的染色体名格式再看BAM文件里染色体名格式。处理方式有两种一是换成配套的基因组和注释二是用sed批量替换GTF里的染色体前缀。我倾向于优先换文件因为批量替换命名有风险部分染色体命名差异不止chr前缀还可能涉及线粒体等特殊序列。换完之后重新构建STAR索引再跑一遍。4.2 featureCounts的assigned比例偏低参数与注释格式的坑现象featureCounts运行日志显示Overall assignment rate只有50%或更低大量reads被归入Unassigned_Ambiguity或Unassigned_NoFeatures。原因一是双端参数没开--countReadPairs没加或者漏了-p导致每个fragment被算成两条独立read计入ambiguous类别二是注释文件选得太粗或太细比如用了只包含变体级别的注释文件就会产生大量NoFeatures三是BAM里混入了比对到线粒体内的reads这部分在常规基因定量时容易被误算为ambiguous。解决先在featureCounts日志里找到unassigned分类统计看是哪一类占大头。如果是Ambiguity检查双端参数是否配对如果是NoFeatures检查GTF是否完整、基因注释版本是否过旧。线粒体reads的处理方式取决于项目需求一般不需要特意去除但如果比例异常高考虑是否线粒体污染可以加一个--ignoreDup或过滤掉比对到chrM的reads后再试一次。提示featureCounts日志最后部分会列出每个样本的assigned/unassigned统计明细这是排查定量异常的第一现场务必保留。4.3 DESeq2的PCA图显示处理组和对照组混在一起设计公式与批次效应现象做完方差稳定化变换后画PCA图发现样本没有按处理组和对照组分开反而是按测序批次或者文库制备日期聚成两堆或者干脆混在一起分不开。原因DESeq2的design公式没有包含批次信息或者包含的批次信息有误。另一种常见情况是样本信息表(coldata)里的行名和counts矩阵的列名不完全一致导致DESeq2读入时把样本排列错乱了这种情况表面看不报错但结果完全不可信。解决先核对coldata的行名是否和counts列名严格一致这一步可以通过colnames(counts)和rownames(coldata)直接比对。然后检查design公式如果有批次信息写成~batchcondition没有批次信息但有多组比较写成~condition。DESeqDataSetFromMatrix里design参数一旦设置错了后面再改会比重新构建整个dds更麻烦所以建议在构建之前就确认清楚。library(DESeq2) counts - read.table(counts/expression_counts.txt, headerTRUE, row.names1, skip1) coldata - read.csv(sample_info.csv, row.names1) all(colnames(counts) rownames(coldata)) dds - DESeqDataSetFromMatrix(countData counts, colData coldata, design ~ batch condition)重点说下最后一行代码里的design参数如果你有一个批次变量列名为batch、一个分组变量列名为condition那么design就写成~batchcondition先放批次再放条件。就算你对批次不感兴趣只要它实际存在不放进design就相当于让批次噪声混入组间差异里结果会异常。all(colnames(counts) rownames(coldata))这一行代码非常值钱它能在你被数据顺序坑到之前就拦住你。5. 从表达矩阵到结论差异基因可视化与快速人工复核技巧5.1 用火山图和已知marker基因做结果合理性验证DESeq2跑完差异分析后第一件事不是急着看富集结果而是先做一次人工合理性检验。我习惯先画一张火山图看看差异基因数量、up/down比例是否正常。如果处理和对照之间几乎没有差异基因、或者差异基因数量大到完全不正常这说明上游数据或者分组设计有问题没必要继续往下走。library(EnhancedVolcano) res - results(dds, contrast c(condition, treat, control)) EnhancedVolcano(res, lab rownames(res), x log2FoldChange, y pvalue, pCutoff 0.05, FCcutoff 1)第二步是验证已知marker基因。比如你在做药物处理转录组而该药物有文献报道的靶基因那就检查这些靶基因是否出现在差异基因列表里方向上是否合理。这一步不需要什么复杂方法就是把基因名和log2FoldChange列出来逐个看。很多流程在这个环节翻车过某个核心marker基因根本没出现在差异列表里结果是注释文件版本错位导致的漏统计。5.2 富集分析里的基因排序方式比你想的重要富集分析最容易做错的一个点是背景基因的选择。做GO富集时背景应该用本次实验检测到的所有表达基因而不宜用物种全部基因组基因。检测到的基因数量取决于表达阈值和过滤标准一般用counts矩阵里做过非零过滤的基因集作为背景。用错了背景富集结果会明显偏高富集到的通路会被无关基因稀释。另一个点是基因排序方式。GSEA类富集分析要求输入基因列表是排序的排序统计量一般用DESeq2结果里的stat列而不是log2FoldChange。用log2FoldChange排序会导致方差大的基因排在前面用stat列会更接近统计学意义上处理组相对于对照组的差异强度。library(clusterProfiler) library(org.Hs.eg.db) geneList - res[order(res$stat, decreasing TRUE), stat] names(geneList) - rownames(res) gse - gseGO(geneList geneList, OrgDb org.Hs.eg.db, keyType ENSEMBL, ont BP, minGSSize 10, maxGSSize 500)这段代码里geneList的names必须是ENTREZID或ENSEMBL格式具体取决于OrgDb的选择。如果直接拿SYMBOL往里传通常会报错可以先通过bitr做ID转换转换的时候注意丢失的基因不要超过百分之五超过说明转换方式有问题。在这个环节我自己有过一次印象很深的踩坑富集结果看起来非常漂亮所有通路都显著后来发现是background用了全基因组而没有用检测到的基因集整个结果虚高得离谱。从那以后我把背景基因列表列成了每次富集前必查的一个步骤希望这篇内容也能帮你避免类似的返工希望帮到你。本文还有配套的精品资源点击获取
