1. 项目拆解hisat2 在 RNA-seq 流程里到底做了什么1.1 为什么在这篇里单独聊“比对”这个环节RNA-seq 流程走到序列比对这一步多数人纠结的不是“他会不会比对”而是“怎么在样本很多的时候把 hisat2 跑得又快又稳”。我见过两种极端一种是一行命令一个样本写了二十行几乎一模一样的代码样本少还能忍样本一多就完全失控另一种是直接上 for 循环但对循环里的变量、日志、断点续跑完全没有概念跑挂了只能在 nohup.out 里痛苦地翻日志。先说清楚比对在整体流程里的位置。一条标准的 RNA-seq 分析主线大概是fastq 原始数据 → 质控与去接头 → 比对到参考基因组 → 转录本组装或定量 → 差异表达分析。比对是整个流程的枢纽后面所有结论都建立在比对结果上。如果这一步出的 bam 文件本身有问题后续 stringtie 组装的转录本、ballgown 或 deseq2 算出来的差异基因全都不可信。hisat2 属于 splice-aware aligner也就是剪接感知比对工具。它和普通 DNA 比对工具最大的区别在于转录组数据里大量 reads 会跨越外显子连接区exon-exon junction如果比对工具不处理剪接位点这些 reads 要么被错误拆分要么直接被丢弃最终导致基因定量严重偏低。hisat2 通过索引里预存的剪接位点信息和局部比对策略能把这些跨剪接位点的 reads 正确放回基因组上这也是它在转录组流程里成为主流选择的核心原因。1.2 hisat2 和 STAR、bowtie2 怎么选选工具之前先看看几个常见选项的差异。很多人第一次跑转录组都会纠结到底用 STAR 还是 hisat2也有同学误把 bowtie2 拿来比对 RNA-seq 数据那基本是给自己挖坑。工具是否支持剪接比对内存占用速度适用场景hisat2支持较低人类基因组索引约 20 GB 上下快普通转录组差异表达分析工作站即可跑STAR支持高人类基因组索引动辄 30 GB 以上内存中高精度转录本定量服务器资源充足时推荐bowtie2不支持剪接很低快DNA 测序、小 RNA、ChIP-seq不适合常规 RNA-seq从实践经验来看做普通的转录组差异表达分析hisat2 的准确率完全够用而且它对服务器的要求低得多。STAR 在某些特殊场景下更稳比如需要检测融合基因或者对转录本异构体组装精度要求极高但代价是内存和运行时间都明显上涨。如果手里的服务器只有 64 GB 内存跑十几个样本的 STAR 比对任务很容易因为并行度过高直接 OOM同样的数据丢给 hisat2开 16 个线程也基本没压力。bowtie2 就更加不建议了它本身不支持剪接比对跑出来的比对率会惨不忍睹而且跨外显子的 reads 全部碎掉后面定量结果没有任何参考价值。我在实际项目中见过有人用它跑 RNA-seq比对率只有 30% 左右排查了半天才发现是工具选错了。所以这篇博文以 hisat2 为主线后面提到的所有脚本、参数、踩坑经验全部围绕这个工具展开。2. 动手之前索引构建与关键参数2.1 从注释文件提取剪接位点再构建索引很多新手第一次用 hisat2走的流程是下载基因组 fasta直接 hisat2-build 一把梭然后拿默认索引开始比对。这个流程跑通没问题但比对质量其实打了折扣因为你没有把基因注释里的剪接位点信息放进去。hisat2 的官方推荐做法是先利用 GTF/GFF 注释文件提取剪接位点splice sites和外显子exons信息再把这两个文件传给 hisat2-build最后构建出的索引里就包含了已知剪接位点。这样比对时跨越 junction 的 reads 能被更准确地放回基因组尤其是那些只跨越一个剪接位点的短 reads。# 1. 从 GTF 提取剪接位点 hisat2_extract_splice_sites.py genome.gtf genome.ss # 2. 从 GTF 提取外显子信息 hisat2_extract_exons.py genome.gtf genome.exon # 3. 把两个信息文件传给 hisat2-build hisat2-build -p 8 --ss genome.ss --exon genome.exon genome.fa genome_tran这里有几个细节值得注意。第一基因组 fasta 和 GTF 的染色体命名必须一致比如都是 chr1、chr2或者都是 1、2不能一个带 chr 前缀一个不带否则提取出来的剪接位点对不上基因组比对时反而可能出现错误。第二GTF 不是越全越好建议先过滤掉那些不标准的基因类型只保留我们关心的类型。有些物种的 GTF 里包含大量预测基因和假基因提取出的 splice site 文件会非常臃肿建索引的时间和磁盘占用都会明显增加。第三hisat2-build 的 -p 线程数不必开太高索引构建阶段并不是线性加速的4 到 8 个线程就够用了。索引构建完成后观察一下生成的 .ht2 文件大小心里有个底。人类基因组 grch38 的索引大约 20 GB 左右视注释版本和剪接位点数量略有浮动。如果是某物种或者特殊参考基因组索引会小很多跑起来也快得多。2.2 两个影响下游结果的参数--dta 与 --rna-strandnesshisat2 参数看着多实际每个样本跑下来常用的没几个。但有两个参数直接决定下游流程能不能顺利走通我必须单独拿出来讲。第一个是 --dta。这个参数的本意是让 hisat2 输出的比对结果更适合下游转录本组装工具典型的就是 stringtie。加了 --dta 之后hisat2 会倾向于输出更简洁的比对信息减少不必要的 soft-clip 和长 CIGAR 拼接避免因为单个 reads 的复杂比对方式干扰后续组装。如果你计划用 stringtie 做转录本定量这个参数一定要加。不加也能跑但 stringtie 处理时可能会遇到很多奇怪的比对记录组装出来的转录本质量也会受影响。还有一个同名参数 --dta-cufflinks那是给老一代 cufflinks 流程准备的。现在大部分人已经不用 cufflinks 了所以不要两个都加。这两个参数在 hisat2 里不能同时使用同时指定会直接报错。第二个是 --rna-strandness。这个参数设置的是文库链特异性方向常见值有三个unstranded、FR、RF。很多人搞不清 FR 和 RF 的区别简单理解RNA-seq 文库制备时如果使用了链特异性建库方法比如 dUTP 法reads 的方向就带有链信息。双端测序时如果 read1 比对到基因组正链、read2 比对到负链这种叫 FR如果反过来叫 RF。具体用哪个值取决于你用的建库试剂盒和实验设计我一般在分析前先查清楚这批数据是哪种文库而不是凭感觉猜。如果这个参数设置错误最直接的后果是比对本身不会报错但基因覆盖度会明显偏向错误链下游转录本定量时大量基因的表达量会被低估或高估。在项目里遇到“所有基因都只有正链覆盖”“比对率正常但 stringtie 结果乱七八糟”这类情况十有八九是链方向搞错了。可以用 RSeQC 的 infer_experiment.py 去检查 bam 文件的链方向它会给出一个链特异性判断用来核对你的设置是否符合预期。infer_experiment.py -r genes.bed -i sample.sorted.bam输出结果里会显示 reads 在正链和负链上的比例如果结果和预期不一致回头改 --rna-strandness 再跑一遍比对。3. 不利用循环的写法能跑但真的不好维护3.1 逐样本处理的常规操作先看看最直接的写法不利用循环每个样本单独写一段比对命令。这也是很多刚接触 RNA-seq 的同学第一反应会采用的方案。假设你手头有 4 个双端样本分别是 sample1 到 sample4数据放在 ~/data 下参考基因组索引在 ~/ref/genome_tran输出目录是 ~/aln。逐样本写出来的命令大概是这个样子# 样本1 hisat2 -p 16 --dta \ -x ~/ref/genome_tran \ -1 ~/data/sample1_1.fastq.gz \ -2 ~/data/sample1_2.fastq.gz \ -S ~/aln/sample1.sam \ --summary-file ~/aln/sample1_summary.txt # 样本2 hisat2 -p 16 --dta \ -x ~/ref/genome_tran \ -1 ~/data/sample2_1.fastq.gz \ -2 ~/data/sample2_2.fastq.gz \ -S ~/aln/sample2.sam \ --summary-file ~/aln/sample2_summary.txt # 样本3、样本4 依此类推如果是单端数据把 -1 和 -2 换成 -U 指定一个 fastq 文件即可。命令本身的逻辑很简单指定索引、输入文件、输出文件加一个 summary-file 把比对统计信息单独保存下来。这种写法直观最不容易让新手产生“脚本是不是有问题”的焦虑因为每一条命令都是完整的、可独立运行的。样本只有两三个的时候这么写完全没问题处理起来也很清晰。3.2 真实项目里最容易翻车的三个问题但一旦样本量上升到十几个、几十个逐样本写法的弊端就会集中爆发。我总结了一下自己在真实项目里遇到过的三个主要问题。第一个问题是维护成本失控。假设你已经写了 20 个样本的比对命令跑了一半发现参考基因组注释版本需要更新索引要重新构建比对命令里的 -x 路径全都要改。这时候你得在脚本里翻出全部 20 段一段一段手动替换稍有不注意就会漏掉一个样本。如果用的是循环脚本你只需要改一个索引路径变量剩下的全自动生效。第二个问题是错误定位困难。逐样本写法没有任何统一的日志管理每个样本只输出一个 sam 文件和 summary 文件。如果某个样本比对中途报错比如内存不足、输入文件损坏、磁盘空间满了你只能在终端里往上翻日志或者在目录里挨个找哪个样本缺输出文件。真实项目里你通常是晚上挂起任务第二天早上来看如果不做统一日志记录排查哪个样本失败了都是一件非常痛苦的事情。第三个问题是新增样本时容易复制粘贴出错。这是最隐蔽的坑。你复制了样本5的命令改成样本6但输出文件忘了改写成 sample5_1.sam或者样本6的数据文件名是 sample6_R1.fastq.gz你还在用 sample6_1.fastq.gz。这些错误在跑命令的时候不会立刻暴露等到比对完成你发现 sample6 的输出其实是空文件或者样本6因为找不到输入文件报了错才会回过头来检查命名。而循环脚本通过统一的命名规则拼接路径这类低级错误发生概率会低很多。所以我现在的观点很明确样本数小于等于 4 个逐样本写法还勉强可以接受样本数超过 6 个强烈建议直接上循环。不要觉得循环难几分钟就能写好但接下来能帮你省下大量排错时间。4. 利用循环批量跑一个脚本管完所有样本4.1 基础循环模板比对、排序、删 SAM 一气呵成用循环处理批量样本核心思路就是把样本名放进一个数组里for 循环遍历所有样本名然后通过样本名去拼接输入文件路径、输出文件路径和日志文件路径。这样不管你有多少个样本脚本主体只需要写一遍。下面这个脚本是实际项目中我经常用的基础模板包括了 hisat2 比对、samtools 排序转 bam、删除中间的 sam 文件以及保存比对统计信息#!/bin/bash # 配置区域 data~/data # fastq 存放目录 idx~/ref/genome_tran # hisat2 索引前缀 aln~/aln # 比对输出目录 log~/aln/logs # 日志目录 threads16 # 每个样本的线程数 mkdir -p $aln $log # 样本名列表按实际项目修改 samples(sample1 sample2 sample3 sample4 sample5 sample6) # 循环处理每个样本 for s in ${samples[]}; do echo processing $s at $(date) # 比对 hisat2 -p $threads --dta \ -x $idx \ -1 $data/${s}_1.fastq.gz \ -2 $data/${s}_2.fastq.gz \ -S $aln/${s}.sam \ --summary-file $aln/${s}_summary.txt \ 2 $log/${s}_hisat2.log # 排序并转 bam samtools sort - 8 -m 2G \ -o $aln/${s}.sorted.bam \ $aln/${s}.sam # 建索引 samtools index $aln/${s}.sorted.bam # 删除中间的 sam 文件释放磁盘空间 rm -f $aln/${s}.sam echo $s done at $(date) done脚本里最关键的是 -1 和 -2 文件路径用${s}_1.fastq.gz和${s}_2.fastq.gz拼接。这意味着你的原始数据文件名必须遵循统一的命名规则。如果样本名是 sample1双端文件就应该是 sample1_1.fastq.gz 和 sample1_2.fastq.gz如果实际文件名是 sample1_R1.fastq.gz那你就把拼接规则改成${s}_R1.fastq.gz总之保证格式一一对应。排序这一步我单独说明一下。hisat2 默认输出的是 SAM 格式文本文件体积大概是 BAM 的好几倍。十几个样本的 sam 文件轻轻松松占掉几百 GB磁盘不够就直接卡死。所以在比对完成后立刻转成 coordinate-sorted 的 BAM 并删除 SAM是基本操作。stringtie 等下游工具需要的正是 coordinate-sorted BAM所以不用担心排序格式问题。4.2 日志、断点续跑和资源分配这些工程化细节基础的 for 循环能跑通但真实项目不能只满足于“能跑”。如果这个脚本要挂在一台服务器上跑一夜你回来却不知道它执行到哪一步了那就非常被动。我在实践中会额外加上几个工程化处理。日志管理是最基本的一条。上面的模板里已经写了2 $log/${s}_hisat2.log把 hisat2 运行时输出到 stderr 的进度信息单独保存到每个样本自己的日志文件。为什么要单独存因为 if 你只用一个 nohup.out 收集所有输出几十个样本的日志混在一起等出问题时找一条报错信息能把人折磨疯。每个样本一个日志文件哪个样本出事就直接看哪个文件定位快得多。断点续跑是我强烈建议加上的功能。循环脚本一旦中途失败比如第 4 个样本因为内存不足挂了等修复问题后从头再跑前面已经跑完的样本又会被重跑一遍白白浪费时间。解决办法很简单在循环体开始处检查这个样本是否已经处理完成。for s in ${samples[]}; do if [ -f $aln/${s}.sorted.bam ]; then echo skip $s, already done continue fi # 比对和排序... done判断条件可以选“sorted bam 文件是否已经存在”因为一个样本只要走到排序完成这一步就说明前期比对没问题。更严谨的做法是额外生成一个 done 标记文件全部流程成功结束后才 touch 一下但实际使用中检查 bam 文件是否存在已经够用了。资源分配是另一个容易被忽视的问题。很多人拿到 64 核的服务器就想着把 30 个样本全部同时丢进去每个样本用 4 个线程看起来并行度拉满。但 hisat2 比对和 samtools sort 都会占用大量内存30 个样本同时跑内存很可能直接爆掉OOM 之后连失败的日志都来不及写。我个人的经验是先用串行的 for 循环把整个流程跑通确保每个样本都能正确完成如果时间来不及再考虑样本级并行。并行时一般控制 3 到 4 个样本同时跑每个样本开 8 个线程就能稳定地吃满一台 32 核、128 GB 内存的服务器。如果你想更灵活地控制并发数可以考虑 GNU parallel。但新手我不建议一上来就折腾这个先把 for 循环版本跑稳理解每一步在做什么再考虑提速的问题。5. 跑完后的质量检查与常见报错5.1 高频问题速查表比对的常见问题其实翻来覆去就那么几类。我把它们在真实项目中出现的频率、表现和解决办法整理成一张速查表方便你排查问题的时候直接对照。现象可能原因解决办法总体比对率只有 30% 或更低物种选错了参考基因组确认测序物种重新下载对应基因组构建索引比对率 80% 左右但明显低于预期接头污染或低质量 reads 太多用 fastp 或 trimmomatic 先做质控再比对70000 运行到一半内存不足直接退出并行样本过多或线程过多减少同时运行的样本数降低 -p 参数samtools sort 报错或卡住磁盘空间不足sam 文件太大及时删除中间 sam或改用 -m 参数控制排序内存所有基因都只有正链或负链覆盖--rna-strandness 设置反了用 infer_experiment.py 检查后重新比对样本列表里重复执行同名任务数组里有重复样本名循环前用 sort -u 对样本列表去重比对结果 bam 文件 reads 数比预期少很多fastq 双端配对信息损坏或不匹配用 seqkit stats 检查 -1 和 -2 的 reads 数是否一致这里面有两个问题特别容易“带偏节奏”。一个是物种和参考基因组不匹配。真实项目里偶尔会混入来源不明的样本或者样本标签写错了你拿着人的参考基因组去比对外源物种的数据比对率当然低得离谱。所以每次拿到新样本我的习惯是先随机抽一个 fastq用 seqkit 统计一下 reads 长度和数量再跑一个 fastqc 看看 GC 含量分布心里有数了再进入批量比对阶段。另一个是接头污染。hisat2 本身不会帮你切除接头它只是一个比对工具。如果原始 fastq 里大量 reads 末端带着接头序列这些 reads 要么比对不到基因组上要么比对到错误的位置最终影响的是整体比对率和定量的准确性。所以标准的 RNA-seq 流程里质控这一步不能省fastp 跑完再进入 hisat2 是很多项目里的标准操作。5.2 从比对率看数据质量别只顾埋头跑命令hisat2 每个样本运行结束后会在 summary 文件里输出一个 Overall alignment rate。这个数字是性价比最高的质量指标。正常的人类或小鼠转录组数据双端 150 bp经过合理质控后比对率通常在 85% 到 95% 之间。如果低于 80%就应该停下来查一查是参考基因组版本不对、样本污染、还是数据质量太差。批量跑完后我的建议是把所有样本的 summary 文件集中起来用一条命令快速扫一眼比对率分布grep overall alignment rate ~/aln/*_summary.txt如果大部分样本都在 90% 左右突然有一两个样本只有 60%那大概率是这几个样本的数据有问题单独拿出来重新走一遍质控和比对而不是对全部样本推倒重来。看活一个样本有没有跑好还可以配合 samtools flagstat 对 BAM 文件做一次快速统计。这个命令会告诉你总 read 数、比对上的 read 数、双端均比对上的比例等。比对率正常的样本flagstat 输出的 paired mapped 比例一般也都在 85% 以上。如果这个比例偏低说明虽然 read 比对上了但双端配对关系异常同样值得关注。样本量多的情况下我建议把所有 flagstat 结果汇总到一起看趋势。比如用循环把所有样本的 flagstat 输出汇总到一个文件里for s in ${samples[]}; do samtools flagstat $aln/${s}.sorted.bam $aln/${s}.flagstat done再配上一个简单的观察逻辑如果所有样本的统计指标分布都差不多说明这批数据整体没有明显问题如果某个样本的指标和群体差异很大那这个样本需要特别留意。不要觉得这些检查多余真实项目的经验是那些最后被判定为“异常样本”的数据往往就在这里露出了端倪。6. 后续流程怎么衔接6.1 给下游 stringtie 的输入准备hisat2 比对完成后你手里是一堆 coordinate-sorted BAM 文件和它们的索引。接下来最常见的走向有两种一是用 stringtie 做转录本组装和定量二是用 featureCounts 直接统计基因表达量。如果你计划走 stringtie 路线前面强调的 --dta 参数就是为这一步铺路的。stringtie 读取 BAM 文件后会依据参考注释进行转录本组装并估算每个转录本的表达量。它要求输入的 BAM 是 coordinate-sorted 格式这也是我在脚本里直接用 samtools sort 排序的原因之一。stringtie 的基本调用方式大约是stringtie sample.sorted.bam \ -p 8 \ -G annotation.gtf \ -o sample.gtf \ -B \ -e其中的 -G 指定参考注释-o 指定输出转录本 GTF-B 会额外输出每个样本的覆盖度文件供后续 ballgown 做差异分析使用-e 表示只估计参考注释中已知转录本的表达量。不过 stringtie 的参数细节不是今天这篇的重点等我后续专门写组装和定量流程时再展开。这里要强调的是不要在比对这一关偷懒。BAM 文件名最好和样本名保持一致sorted.bam 和对应的 .bai 索引文件应该放在同一个目录。很多下游分析工具在读取多个样本时会自动按样本名匹配 BAM 文件如果你的命名混乱后面跑差异分析时光是整理样本映射表就要费不少功夫。6.2 想提速时并行策略要稳妥如果你这批样本特别多比如超过 30 个串行跑确实会让时间线拉得很长。这时候再考虑并行也不迟但要讲究策略。我最推荐的稳妥并行方式是这样的先确认服务器的 CPU 核心数和内存总量然后估算每个 hisat2 样本需要的线程数和内存。以一台 64 核 256 GB 内存的服务器为例每个样本开 8 线程、预留 16 GB 内存一次并行跑 4 个样本比较稳妥既能充分利用资源又不会内存爆掉。如果你的服务器内存小一些就相应减少并行数量。如果你用的是 for 循环串行版本想过渡到并行可以用 xargs 指定并发数比如# 生成样本名列表 ls ~/data/*_1.fastq.gz | sed s/_1.fastq.gz// | sort -u \ | xargs -P 4 -I {} sh -c bash run_one_sample.sh $1 _ {}这里的关键是把每个样本的处理过程封装进一个独立的脚本再由 xargs 以 4 个并发作业的方式调度。这样比你自己在后台手动开多个 for 循环要规范得多。不过这种封装方式适合你已经把单样本流程跑通的情况千万不要第一次跑比对就这么整排错时会非常痛苦。我个人在实际操作中的体会是RNA-seq 比对这个环节脚本写得简单清晰比炫技更重要。循环也好不循环也罢核心目标是让每个样本的比对结果可追踪、可复现、可检查。先把一个样本的完整流程走通再把所有样本纳入批量管理最后才谈并行提速——这是最稳的路线。这篇主要把 hisat2 比对阶段讲透了后续的 stringtie 组装、转录本定量、差异表达分析我会继续更新在这个系列的下一篇里。如果你现在正要处理一批真实的 RNA-seq 数据建议先照着上面第 4 节的循环脚本把流程跑通再根据自己的数据特点调整细节。