1. 项目全景复盘一条贯穿RNA-seq分析的完整链路高通量测序数据分析做到后面你会发现真正决定一个项目生死的往往不是某个高深的算法而是你有没有在动手前把整条逻辑链条想清楚。RNA-seq项目尤其如此它不像基因组测序那样“测完就等于大半完成”RNA-seq的每一个分析环节都在回答一个具体的生物学问题而这些问题之间又环环相扣。1.1 破题RNA-seq到底在测什么、算什么RNA-seq本质上是测量细胞内所有RNA分子的“快照”。实验人员把细胞里的总RNA提取出来用oligo(dT)富集带poly(A)尾巴的mRNA或者用rRNA去除试剂盒去除核糖体RNA然后进行片段化、反转录成cDNA、加上测序接头最后上机测序。你拿到手里的fastq文件就是这些cDNA片段两端或单端几十到几百个碱基的序列。高通量测序数据分析的宏观目标就是把不计其数的短读段reads重新“拼”回它们原本属于的基因数一数每个基因有多少条reads落到它上面。基因覆盖的reads数量就构成了表达量。所以RNA-seq数据分析和WGS全基因组测序最大的不同在于WGS重点在于发现变异位点而RNA-seq的核心矛盾在于“比对到剪接后的转录本”和“对转录本进行准确量化”。你做的每一步从质控、比对到定量都在围绕“准确量化”这四个字展开。我当时接手第一个正式项目时测序公司给了我大概30个样本的RNA-seq数据双端150bp每个样本6-8G的原始fastq。面对这些数据新手最容易犯的错误就是直接拿一个工具链从头跑到尾跑到最后发现结果一塌糊涂却不知道问题出在哪一步。正确的做法是建立清晰的阶段意识。1.2 分析链路拆解从FASTQ到差异基因的路线图RNA-seq数据分析按惯例可以切成三段上游原始数据清理、中游比对与定量、下游差异表达与功能解读。这三级你可以理解成上游把机器的原始产物fastq变成干净、可信的序列。中游把干净的序列精确映射到基因组/转录组并统计每个基因的计数counts。下游基于counts矩阵做统计建模筛选差异表达基因并映射到生物学通路pathway上。有人会问为什么一定要区分上下游其实这是工程化思维边界清晰才能阶段化验收。比如上游质控不能达标下游差异分析做得再花哨都是沙上建塔。每个阶段都有明确的输入输出文件格式我在动手前一定会先确认这些fastq文件命名是否规范R1/R2是否对应。参考基因组和注释文件的版本是否匹配。链特异性建库方向strandedness是什么这直接影响比对后的定量参数。如果你使用的是常见的Illumina平台一般会有两种情况链特异性文库dUTP法和非链特异性文库。这一步信息不明确后面featureCounts和HTSeq-count会算错表达量而且错得非常“自然”不仔细看很难发现。另外我建议在项目启动前把工具环境搞定。Conda是目前最省心的方案新手可以直接用Miniconda创建一个独立的分析环境避免不同工具依赖冲突的坑。我自己常用的环境构建命令是conda create -n rnaseq python3.9 conda activate rnaseq conda install -c bioconda -c conda-forge fastqc multiqc trim-galore hisat2 subread samtools这个命令会把后续用到的上游分析工具一次性装齐。我在刚开始做生信时经常因为环境问题耗费大量时间后来养成了“先配环境再做分析”的习惯确实省了很多无意义的折腾。2. 从海量读长中抽取信号的预处理实操拿到测序公司交付的fastq文件后我的习惯是先别急着做比对。为什么因为测序仪偶尔会给出一些质量不达标的读段或者接头没有完全去除干净。这些“脏数据”如果不清理会在比对时产生大量错配导致下游定量偏低甚至引入假阳性。2.1 用FastQC和MultiQC建立质控基线结果怎么看FastQC是最常用的测序数据质量评估软件相信做这一行的人都听过。运行方式非常简单fastqc -t 8 *.fastq.gz -o ./qc_raw注意-t参数指定线程数如果你的服务器有32个CPU核心跑8~16线程都不会有太大压力。跑完后每个样本会生成一个html报告你可以用浏览器打开逐项检查。但当你手里有二三十个样本时每个都点开看效率太低所以我会再跑一次MultiQC把目录下所有结果汇总成一个报告multiqc ./qc_raw -o ./qc_raw_multiqcFastQC报告里需要重点关注的项目大概四个Per base sequence quality看每个碱基位置的测序质量分数。如果红线中位数低于Q20说明这一段测序质量差下游要重点考虑低质量碱基的切除。Per sequence GC content看所有读段的GC含量分布。正常应该近似正态分布如果出现一个异常峰可能意味着文库有污染或者PCR偏好性严重。Adapter Content看接头序列的残留程度。如果这个模块大面积超标就说明建库时接头去除不彻底必须修剪。Overrepresented sequences如果某条序列占比异常高可能是接头二聚体也可能是某些高表达基因的产物。我碰到过一种常见情况样本的Per base sequence quality在末端位置第140~150bp明显下滑Q30比例从90%跌到70%。这种时候我会倾向于做一次末端修剪因为下游的比对算法在低质量区域容易产生错误匹配从而干扰定量。2.2 用trim_galore去除接头与过滤低质量碱基参数怎么定trim_galore是包装了Cutadapt的常用软件适合处理Illumina双端数据。最常见的运行命令如下trim_galore --phred33 -q 25 --length 36 --stringency 3 --paired -o ./clean Sample_R1.fastq.gz Sample_R2.fastq.gz这里每个参数都在回答一个实际分析问题--phred33声明质量分数体系是Phred33。这是目前Illumina的标准格式但如果你的数据比较老可能是Phred64格式需要手动调整。-q 25以滑窗方式切除质量值低于25的碱基。为什么是25Phred质量值25对应99.7%的碱基正确率这个标准平衡了准确性和保留率是我们实验室验证过的经验值。--length 36如果一条读段经过修剪后长度小于36bp就直接丢弃。36是最小匹配阈值太短了比对软件无法可靠地定位到转录本上。--stringency 3允许接头匹配的错配数达到3个。实际测试中--stringency 1过严格容易漏掉接头--stringency 3在去除效率与误切之间是比较理想的平衡点。修剪之后我会把原始和修剪后的fastq各自跑一遍FastQC对比一下前后质量变化。例如某批牙髓干细胞的样本经过trim后Q30比例从84%提高到94%这让我对接下来的比对结果很有底气。这里要特别提醒一句修剪不是越狠越好。质量太差我们可以切但如果你把读段切到30bp以下反而会让比对软件更难处理。过度修剪还会造成序列信息量不足特别在转录本异构体的区分上非常不利。我一般会把保留的读段长度阈值设在36~50bp之间具体看测序读长和插入片段分布少于这个阈值的数据就当作废了。3. 以精准比对撬动表达定量的上游关键实操预处理做完下面进入让无数新手头疼的比对环节。RNA-seq比对与DNA-seq比对有本质差异基因组是连续的而转录组的成熟mRNA由于剪接事件的存在一个读段可能会跨过外显子连接处。如果直接用普通比对软件比如BWA去比对这些跨越剪接位点的读段就会被当成比对不上的“孤儿”从而大量丢失。3.1 参考基因组版本选择与Hisat2参数细节RNA-seq比对工具的主流是STAR和Hisat2。STAR的优点是速度快、准确度高但吃内存非常夸张一个人基因组大概需要30GB左右的内存如果你的服务器内存不够大就会跑得很痛苦。Hisat2则基于FM索引算法对内存更友好同时专门针对剪接位点进行了优化是我在小服务器上做常规转录组分析时的首选工具。比对前需要准备好参考基因组索引。这一步最需要注意的是基因组版本和注释版本的匹配。我一直推荐使用Ensembl或GENCODE发布的参考基因组和GTF注释文件比如当前最新的人基因组版本是GRCh38或叫hg38。千万别从不同来源拼凑比如基因组用UCSC的hg38注释用Ensembl的GRCh38虽然在大部分位点能对上但在基因命名空间上会有差异导致特征计数的时候出现基因ID匹配不上的情况。Hisat2建索引hisat2-build -p 20 genome.fa genome_indexHisat2比对命令hisat2 -p 16 --dta -x ./genome_index/genome \ -1 clean_R1.fq.gz -2 clean_R2.fq.gz \ -S output.sam 2 alignment_log.txt离最终比对结果可用还要再过一道关键的手续——按坐标排序并转成BAM。这一步直接用samtools完成samtools sort - 16 -O BAM -o sample_sorted.bam output.sam samtools index sample_sorted.bam这里为什么强调加--dta参数是因为它为下游Hisat2官方推荐的转录本组装工具StringTie做了优化设计。不过如果你用的是featureCounts来做定量有没有--dta影响不大。真正影响后续featureCounts的是排序方式samtools sort按基因组坐标排序是标准操作很多新手不小心只做了-b转BAM没排序后面featureCounts就会报错。比对完成后马上要检查比对率。我一般用cat alignment_log.txt看summary整体比对率在90%以上正常放心往下走。整体比对率在80%90%还能接受注意检查是否有过多多比对multi-mapping。整体比对率低于80%强烈警惕可能是样本污染、参考基因组选择错误或者建库问题。此时必须往回查FASTQ质控报告而不要带着这个不靠谱的BAM继续向下做。我曾经收到一份样本比对率只有62%。当时一度怀疑是参考基因组版本不匹配后来查了一圈才发现客户送测的样本物种其实是小鼠而他们提交信息时写的是人。这件事给我上了很深刻的一课拿到数据先确认物种、确认参考基因组永远是第一步。3.2 链特异性与BAM质控的隐性工程文库方向strandedness是RNA-seq分析中最容易踩的暗坑。RNA-seq建库分为链特异性stranded和非链特异性unstranded。链特异性文库能保留转录本原来的方向信息因此在定量时能够提供更精确的正反义链区分这对检测反义RNA转录本时尤其重要。具体到操作层面使用featureCounts做定量时必须指定-s参数-s 0非链特异性文库-s 1链特异性文库且read1与基因链方向相同-s 2链特异性文库且read1与基因链方向相反。如果你不确定文库类型可以用RSeQC的infer_experiment.py脚本对BAM文件做一次快速判断。这个方法能帮你在第一次flow中就确定方向而不用等到定量结果出来后才发现基因表达量整体掉了一半那时候挽回的代价就大多了。比对环节的另一个重要操作是reads distribution检查。你可以使用RSeQC的read_distribution.py统计reads在外显子、内含子、基因间区的分布比例。一个理想的RNA-seq文库reads落在外显子的比例应该在60%以上。如果内含子比例高得离谱说明样本有基因组DNA污染或者annotation版本不对。有时候实验操作或RNA提取过程中混入基因组DNA会导致这种情况单靠生信手段很难完全消除唯一建议是重新提取RNA再建库。4. 定量工具比的不是想象力而是看得懂的“数值”定量是RNA-seq从“序列”走向“数量”的关键一环。很多新手在这时候会纠结到底用RPKM、FPKM还是TPM其实这几个指标都是在回答同一个问题基因表达量有多大但它们的归一化原理各有不同。4.1 使用featureCounts构建基因水平的counts矩阵在差异表达这个步骤之前我强烈建议所有原始计数raw counts都用featureCounts这类软件来生成。即使你是从Salmon这类伪比对工具拿到定量结果也建议保留一个基因水平的counts版本以备后用。featureCounts属于subread软件包安装非常简单conda install -c bioconda subread基础运行命令featureCounts -a genome.gtf \ -o counts.txt \ -T 8 \ -p --countReadPairs \ -s 2 \ -t exon \ -g gene_id \ sample1_sorted.bam sample2_sorted.bam sample3_sorted.bam各参数含义-aGTF注释文件的位置。-T线程数按服务器核心数调整。-p --countReadPairs处理双端测序时按read pair计算避免把同一个片段计两次。-s 2链特异性方向需要提前用RSeQC判断。-t exon只统计落在外显子上的reads。-g gene_id将外显子计数聚合到基因水平。GTF里可以按gene_id聚合也可以按transcript_id聚合。默认经验是差异表达用gene_id可变剪接分析才需要transcript_id。完成后counts.txt文件里有几列关键信息Geneid、Chr、Start、End、Strand、Length以及每个样本的原始计数。别嫌它格式蠢这些全都有用。还有很常见的一个问题是大家把raw counts聊到和RPKM/TPM混为一谈。我只想提醒一件事DESeq2的输入是raw counts不是RPKM或者TPM。RNA-seq的差异表达检验本质是对counts进行负二项分布建模降低文库大小和基因长度对测序深度的影响需要靠软件内在的normalization方法。如果你手动提供RPKM只会让统计模型失效这是一条走不通的死路。4.2 Salmon与伪比对的场景化选择除了传统“比对基因组再计数”的路线现在很多团队也喜欢用Salmon/Kallisto这类伪比对工具。这类工具不直接把reads比对到基因组上而是比对到一个包含所有已知转录本序列的“转录本索引”上然后用期望最大化算法EM估计转录本丰度。好处是速度极快、内存占用低而且不需要GTF文件来“数”基因——你直接得到的是转录本水平transcript level的定量结果。我一般在客户要求做转录本异构体分析或者对样本量大、时间紧的筛选项目时会优先考虑Salmon。骨架流程# 构建转录本索引 salmon index -t gentrome.fa -i salmon_index -p 8 --gencode # 定量 salmon quant -i salmon_index \ -l A \ -1 clean_R1.fq.gz -2 clean_R2.fq.gz \ -p 8 \ --validateMappings \ -o sample_quant/其中-l A让Salmon自动判断文库方向。这一步相当于把前面人工判断链特异性的环节省掉了适合自动化率高的流程。但在这里必须强调一个上下游衔接问题Salmon输出的是转录本水平的TPM和counts。要做基因水平的差异分析需要先用R包tximport把转录本计数汇总到基因水平。如果这一步直接偷懒用sum的方式聚合counts那些多对一的关系处理起来就不稳妥。我一般这样处理library(tximport) library(readr) files - file.path(sample_quant, quant.sf) names(files) - c(sample1, sample2, sample3) txi - tximport(files, typesalmon, tx2genetx2gene, countsFromAbundanceno)countsFromAbundanceno表示保留原始counts供DESeq2做后续统计这样既保留了伪比对在区分异构体上的灵敏度又不丢失下游分析的严谨性。5. 揭示生物学意义的差异表达与富集套路定量完成之后你就得到了一个行为准则清晰的“基因×样本”矩阵。接下来最重要的任务是从一堆数字中找出统计学上有显著差异的基因并把它们放到生物学功能层面去解释。5.1 用DESeq2做差异表达分析的正确姿势在edgeR、limma和DESeq2三巨头里如果只挑一个作为主力工具我强推DESeq2。它内置了估算基因间离散度dispersion和独立过滤independent filtering的方法在小样本量实验中表现得相对稳健而且用户基数大遇到问题好搜解决方案。启动DESeq2的第一步是读取原始counts矩阵和样本设计信息。这里特别提醒设计矩阵的分组列必须是因子而且对照组在第一水平。否则你会得到完全相反的符号log2FoldChange这在判断“上调/下调”时会致命。下面这段代码是我在几乎所有RNA-seq项目里都能稳定运行的基础模板library(DESeq2) # 1. 读取counts矩阵 countData - read.csv(counts_clean.csv, row.names 1, check.names FALSE) countData - round(countData) # 确保是整数 # 2. 设置样本信息对照组在前 colData - data.frame( condition factor(c(ctrl, ctrl, treated, treated), levels c(ctrl, treated)) ) # 3. 构建DESeqDataSet dds - DESeqDataSetFromMatrix( countData countData, colData colData, design ~ condition ) # 4. 标准化与差异检验 dds - DESeq(dds) # 5. 提取结果 res - results(dds, alpha 0.05) res - lfcShrink(dds, coef condition_treated_vs_ctrl, type apeglm)为什么在第5步用lfcShrink做一次收缩因为对于低表达基因原始log2FoldChange非常不稳定可能只是多了一条reads就显示10倍差异。使用apeglm方法对倍数变化进行收缩可以让低表达基因的差异倍数变得更保守可靠从而减少后续验证的假阳性率。接下来筛选显著差异基因我的常规阈值为sig_genes - subset(res, padj 0.05 abs(log2FoldChange) 1)注意padj是校正后的P值这比pvalue要严格得多。因为差异检验动辄检测两万多个基因如果不做多重假设检验校正假阳性数量会非常可观。顺手应该做一个PCA图检查样本聚类情况。DESeq2提供了plotPCA函数或者你也可以先vst(dds)获取方差稳定变换后的表达矩阵再做PCA。如果同组样本完全分不开那就说明批次效应严重或分组逻辑有问题这时候强行找差异基因容易得到一堆没有意义的假“噪声”。批次效应是另一个绕不开的问题。如果一批样本在测序仪上分了两条lane或者建库时间不同做样本聚类时往往会看到批次聚在一起。我可以接受轻微的同组内分散但如果对照和处理组之间因为批次效应直接重叠我建议引入批次变量到design公式中design ~ batch condition额外提一嘴很多新手喜欢直接在excel里算p值再自己挑差异基因。这种操作请务必停手——不经收缩、不估算离散度、不考虑模型的差异分析是无根的浮萍发出去会被审稿人批斗得体无完肤。5.2 从GO到KEGG富集分析的正确解读方式差异基因拿到手之后不能只看着几百个基因发呆。接下来要回答的问题是这些基因在哪些生物学过程BP、细胞组分CC、分子功能MF或者信号通路中显著富集我推荐使用clusterProfiler这个R包是目前公认API设计比较顺手的富集分析工具。经典代码library(clusterProfiler) library(org.Hs.eg.db) # 差异基因列表最好是Entrez ID gene_list - rownames(sig_genes) # GO富集 ego - enrichGO(gene gene_list, universe rownames(countData), OrgDb org.Hs.eg.db, keyType ENSEMBL, ont BP, pAdjustMethod BH, pvalueCutoff 0.05, qvalueCutoff 0.05) # KEGG富集 ekegg - enrichKEGG(gene gene_list, organism hsa, pvalueCutoff 0.05)这里有两个非常关键的细节。第一是必须设置universe背景基因集。有些教程会省略掉这个参数但这样富集分析的结果会被所有被检测过的基因干扰导致富集结果虚高。背景基因一般等于counts矩阵里所有检出基因也就是所有在你这批样本中表达量大于零的基因。第二是ID类型要跟注释匹配。如果你的差异基因是ENSEMBL格式而org.Hs.eg.db默认接受ENTREZID那你在传进去之前最好先bitr()转换一下否则会得到大量NA。KEGG富集本身依赖联网接口。如果在分析环境中无法正常访问KEGG API装不上或者超时就可以尝试用gson包本地化KEGG数据或者退而求其次选用ReactomePA和Msigdb的GSEA分析。GSEA与GO/KEGG最大的不同在于GSEA不做简单的“有/无”二值化筛选而是保留整个表达谱的排序信息去综合分析一个基因集合的协同变化趋势因此在挖掘生物学功能时往往更能还原真实调控模式。写富集结果时我建议把通路名称、富集基因数、p.adjust值这些关键指标放在表格里再配一个气泡图或条形图。千万别只看一个TOP30气泡图就下结论一定要看具体富集到了哪些基因以及这些基因在你的差异表达列表里的方向和倍数。6. 多年实战凝结的避坑指南与高效技巧这部分价值密度最高。我把过去几年在RNA-seq分析中踩过的坑和修复方法全部整理出来希望能帮你省下一年摸索时间。6.1 新手最容易忽略的五个典型毛病第一忽略链特异性导致定量混乱。你可以在比对结束后马上用RSeQC的infer_experiment.py判断文库方向这一步两分钟就能完成但能救回一个项目。我曾经见过一个伙伴因为没判断链方向直接用-s 0跑featureCounts结果样本内基因表达量普遍掉到一半以下差异分析做出来一片狼藉。第二使用旧版本的参考基因组和GTF注释。全新的基因注释版本往往包含更多新验证的转录本信息基因命名规范也更完整。如果是两三年都没更新的项目数据建议重新下载最新版本并重建索引不然基因ID可能对不上后续数据库更新。第三把标准化表达量TPM/RPKM输入DESeq2。这是致命的概念错误。DESeq2需要的是raw counts因为它要基于read计数和文库大小去做离散度估计。TPM已经经过基因长度归一化强行套进负二项分布模型会破坏模型假设结果毫无统计意义。第四忽视比对率异常就一路硬跑。无论哪个环节的比对率低于预期都应该停下来排查原因。样本是否搞混了参考基因组是否为对应物种建库是否有污染这些问题早发现早解决千万不要等到差异分析结束才回头查。第五单样本组做差异分析。这是很多刚起步的实验室容易踩的大坑。每组至少3个生物学重复是底线不然DESeq2算不出离散度检验结果毫无可靠性。经费实在紧张数量也不能少于2个否则只能当探针筛选项目发内部报告不能作为严谨结论。我把常见问题整理成速查表方便随时查阅现象可能原因排查手段比对率低于80%物种混淆、参考基因组版本错误、建库污染检查fastq物种来源确认参考基因组RNA-seq QC报告GC含量异常双峰基因组DNA污染或PCR扩增偏差检查文库片段分布必要时重新建库所有样本PCA完全分不开严重批次效应或分组逻辑错误检查建库批次、上机lane引入批次变量差异分析结果几乎全是0counts矩阵中有大量NA或0未做低表达基因过滤用rowSums(counts 10) 3过滤低表达基因某基因在样本A高表达但在B完全为0比对或定量错误、链方向参数错误检查该基因的BAM覆盖情况反向验证链方向6.2 流程工程化让分析跑得更快更稳的细节RNA-seq分析不是只会敲几行命令就行它本质上也是一个工程问题。在一个数十样本的项目中手动一句一句跑命令会让你浪费大量时间在等待和复制粘贴上。我强烈建议用Snakemake或者Nextflow把流程管理起来即使不用完全容器化也可以把每一步的规则rule定义清楚保证等一下批样本进来时自动重跑。举个简单的Snakemake逻辑rule all: input: results/counts.txt rule fastqc: input: data/{sample}_R1.fastq.gz output: qc/{sample}_R1_fastqc.html shell: fastqc -o qc {input} rule trim: input: r1data/{sample}_R1.fastq.gz, r2data/{sample}_R2.fastq.gz output: r1clean/{sample}_R1.fastq.gz, r2clean/{sample}_R2.fastq.gz shell: trim_galore --paired --phred33 -q 25 -o clean {input.r1} {input.r2}如果你现在还是单步手动操作十几个样本把这些规则文件建立起来并不难却能让你从重复劳动中解放出来。另一个能直接提升速度的小技巧是给samtools sort和featureCounts都设置恰当的线程数同时留意IO带宽瓶颈。如果服务器磁盘是机械硬盘跑上百G的数据时会明显拖慢速度有条件的话尽量把数据放在SSD或高速临时目录中。关于资源估算我一般按经验值准备内存Hisat2比对一个样本约需10~15GB内存STAR则需要30~50GB。Salmon和featureCounts相对轻量8~16GB就够用了。如果服务器内存就32GB强行跑STAR请准备好CPU被卡顿和等待的姿势。写在最后如果非要说做RNA-seq数据分析到底什么最要紧我会毫不犹豫地回答耐心和验证精神。我在一次次踩坑后逐渐养成了在每一步都留一个检查点的习惯——质控看一眼Q30比对看一眼比对率定量看一眼链方向和counts分布差异分析看一眼PCA和收缩效果。这些检查点加起来不会超过半小时却能把一个可能跑偏数日的项目及时拉回正轨。最后分享一个每次项目必做的小动作分析完成后把关键中间文件clean fastq、BAM、counts矩阵备份到独立目录并写一份MD5清单。这个做法让我在反复调参时始终有稳定可靠的版本可供回滚也让同事之间协作时不会互相覆盖文件。数据分析本身没有绝对标准但只要你每一步都踩实结果自然经得起推敲。
