做转录组测序（RNA-seq）分析的人，几乎每天都要和“序列比对”这一步打交道。拿到下机数据之后，质控、过滤、比对、定量一套流程走下来，最耗时也最容易出问题的就是比对环节。这个环节把几千万条reads放回参考基因组上，决定了后续基因表达定量和转录本组装的质量。而hisat2这几个工具，很多做生信的同学都不陌生，但我发现不少教程只是给了命令，没讲清楚为什么这么写、不同写法之间差异有多大。
我这次就把自己的完整操作经验拆开来写。标题里提到的“不利用循环”和“利用循环”，其实是批处理两种思路：一种是逐条样本“人肉机械式”敲命令，另一种是用Shell循环批量搞定十几个甚至几十个样本。很多人问过我这俩到底差在哪，更想让我把代码贴出来直接能用。所以我打算分两部分说，同时把比对环节里涉及的核心概念、参数选型、常见报错一起整理了。
文章适合正在入门RNA-seq分析的同学，或者已经跑通流程但只会照抄命令、想弄明白原理的从业者。我会用自己实际跑过的一套双端测序数据来举例，配合每一步操作的说明，尽量让你看完之后不仅会跑，还能改，还能排错。
1. RNA-seq比对这一步到底在做什么，为什么非用hisat2
先花点篇幅把“RNA-seq比对”的定位讲清楚。整个RNA-seq标准分析流程，大概是：原始下机数据（FASTQ）→ 质量评估和过滤（fastqc/fastp/trimmomatic）→ 序列比对（hisat2/STAR/bowtie2）→ BAM文件处理（samtools）→ 基因表达定量（featureCounts/HTSeq）→ 差异表达分析（DESeq2/edgeR）→ 功能富集。比对是这整条链路里承上启下的关键节点，后面定量基因表达量、找差异基因、拼接新转录本都要回归到BAM文件里比对上的reads上。
为什么RNA-seq不能用普通的DNA比对工具？因为真核生物的基因组里有内含子，转录组测序得到的reads往往跨过外显子和外显子的剪接连接区域。普通的比对工具只允许连续比对，遇到跨越剪接位点的read就直接比对不上或者被截断。hisat2之所以在转录组领域被广泛使用，用的是一种叫做“层级图FM索引”的比对算法，它把基因组的索引分成全局FM索引和局部FM索引，能高效处理跨内含子的read。配合已知剪接位点信息，它在中间价位上能找到外显子跳跃、可变剪接等事件。
我对比过几款常用工具，简单列个选型感受：
| 工具 | 算法特点 | 适用场景 | 资源占用 | 我的使用建议 |
|---|---|---|---|---|
| hisat2 | 基于BWT/FM索引，支持剪接位点识别 | 常规转录组定量、转录本组装 | 较低（4线程跑4G数据约15分钟） | 默认首选，平衡速度和内存 |
| STAR | seed搜索，长read拼接 | 高精度、检测融合基因 | 高（需30G以上内存，建立索引磁盘占用大） | 样本量大或服务器内存充裕时可以用 |
| bowtie2 | 不支持跨剪接位点 | 只做miRNA、ncRNA等非剪接转录本 | 低 | 不适合RNA-seq |
| BWA | 不支持剪接，为DNA测序设计 | DNA重测序、变异检测 | 中 | RNA-seq里几乎不用 |
我个人的经验是：如果只是做常规的转录组差异表达分析，没特殊需求，hisat2在大多数情况下的性价比最高。它的索引文件比STAR小得多，运行速度快，加上官方提供人类、小鼠等常见物种的预构建索引，下载解压就能用，节省掉自己费内存构建索引的时间。这种“开箱即用”的感觉在生信工具里真的不多。
2. 两个关键前置动作：参考基因组索引和输入数据的准备
跑hisat2比对之前，有两件事闲着也是闲着，早点做完能省很多事。一个是准备好参考基因组并构建索引，另一个是把FASTQ里的reads“收拾整齐”，别拿着带一堆接头和低质量碱基的数据直接去喂比对工具。
2.1 参考基因组和注释文件从哪里来
这一步我见过很多人卡住，本质上是“不知道去哪下载对的东西”。常见的数据源是GENCODE、Ensembl或UCSC。以人类为例，推荐去GENCODE下载全套FASTA和GTF注释，版本号要对得上。比如GRCh38就对应gencode的v44左右版本（具体版本每年更新，以官网为准）。有人会图省事，从不同平台拼凑了一份基因组和一份注释，比对跑完才发现基因注释里有大量ID对不上，这会给后面的featureCounts定量挖坑。
需要下载的两类核心文件：一级是参考基因组FASTA（一条条染色体的序列），另一级是GTF注释文件（记录每个基因、转录本、外显子坐标）。构建hisat2索引只需要FASTA，但如果想在比对时利用已知剪接位点，就要从GTF里提取剪接位点信息。有人问不用已知剪接位点会怎样，当然可以，hisat2照样能通过GTF索引去发现跨内含子的read，只是准确率和速度稍逊，所以标准流程里我通常都会把剪接位点文件准备出来。
2.2 用hisat2-build构建索引
hisat2构建索引的命令其实很直接：
# 假设FASTA文件为 genome.fa hisat2-build -p 8 genome.fa genome_index这里genome_index是索引的前缀名，运行完会生成genome_index.1.ht2、genome_index.2.ht2等一批后缀文件，这些都是后续比对时-x后面要指定的前缀名。如果要对人类或小鼠做全基因组索引，内存要求不低，我实测32GB内存的机器构建人类全基因组索引依然吃力，这时我的建议是直接去hisat2官网下载预构建索引包，里面存储了基因组索引、转录组索引以及snp索引，省心很多。
如果基因组太大或想进一步压缩，还可以考虑用--snp附带常见SNP信息来构建，这样比对时不会因为个别SNP位点降低比对率。这一步属于优化项，不是必需。
2.3 比对前要把reads处理成什么状态
我见过不少人为了省事直接跳过质控就拿原始FASTQ跑比对，结果比对率惨不忍睹，回头查半天发现是接头污染加上前几个循环低质量碱基过多。这一步建议至少跑一下fastp，做去接头和低质量过滤。fastp一条命令搞定，而且自动输出JSON/HTML报告：
fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz -o sample_clean_R1.fastq.gz -O sample_clean_R2.fastq.gz -h sample_fastp.html -j sample_fastp.jsonfastp处理完之后，比对输入的数据就应该是质量均匀、几乎不含接头序列的干净reads。这里有个常见误区：fastp默认会trim到前90个碱基，如果你的测序读长是PE150，本来150bp的reads会被截到90bp，丢失部分有效信息。所以我习惯加一个参数--max_len1 150 --max_len2 150关掉自动截断，或者另外指定。总之，看报告、理解每个参数再动手，比盲目套命令可靠得多。
3. hisat2的常用参数你真的搞懂了吗
很多人拿到一条命令就直接复制粘贴，参数是干什么的完全没概念。等哪天条件变了要改，就手足无措。我在这儿把自己的用法和参数背后的逻辑说明白。
3.1 最常用的一组命令
假设我们有一份经过fastp处理的样本，双端reads分别叫sample_clean_R1.fastq.gz和sample_clean_R2.fastq.gz，同时我们下载了官方的hisat2索引，前缀是grch38/genome。那么基础命令是：
hisat2 -p 8 --dta -x grch38/genome -1 sample_clean_R1.fastq.gz -2 sample_clean_R2.fastq.gz -S sample.sam-p 8：使用8个线程并行计算。线程越多越快，但注意不是线性加速，而且内存随之增加。-x：指定索引前缀。它的后面不要加.1.ht2这类后缀，否则会找不到。-1/-2：双端reads。如果是单端数据用-U。-S：指定输出SAM文件的文件名。如果不加-S，结果会直接哗啦啦打印到屏幕上，你只能重定向到一个文件。--dta：为下游转录本组装设计的参数。加了之后，比对结果会包含更适合组装工具的比对记录。如果你后续直接用featureCounts定量，这个参数的影响有限，但也不会造成伤害，所以我一般统一加上。
3.2 链特异性参数的选择
做链特异性建库的样本，需要告诉hisat2测序链的方向。常见的取值是RF和FR。怎么判断用哪个？看你的文库构建试剂盒说明：目前Illumina的TruSeq Stranded Total RNA基本都是RF，早期的dUTP法也多数是RF；很多用SMARTer建库的自由链是FR（其实更准确应该是fr）。他这个参数会直接影响后续定量的正负链基因分配。
我曾经有一个样本，怎么跑怎么觉得基因表达量分布异常，后来发现就是链特异性设置反了，改完--rna-strandedness RF才恢复正常。怎么验证呢？可以用RSeQC的infer_experiment.py脚本自动判断文库链方向：
infer_experiment.py -r genes.bed -i sample.sorted.bam输出里的1++,1--,2+-,2-+比例如果接近RF，那链特异性参数就没问题。
3.3 已知剪接位点和新剪接位点的处理
hisat2不仅能利用已知的GTF注释来辅助剪接位点检测，还能把发现的新剪接位点输出到独立文件。我的习惯是：比对已知物种时，从GTF提取剪接位点文件，用--known-splicesite-infile传入。同时打开--novel-splicesite-outfile保留新发现的剪接位点，这样再跑第二轮时可以叠加这部分信息。
提取剪接位点文件用官方脚本：
hisat2_extract_splice_sites.py genome.gtf > genome.splice_sites然后在比对命令里加：
hisat2 -p 8 --dta --known-splicesite-infile genome.splice_sites \ --novel-splicesite-outfile sample.novel.splice_sites \ -x grch38/genome -1 sample_clean_R1.fastq.gz -2 sample_clean_R2.fastq.gz \ -S sample.sam这么做的好处是，比对时优先在已知剪接位点附近搜索，速度更快、跨剪接位点read的命中率更高。我在一个含大量可变剪接事件的组织样本上实测，加入已知剪接位点文件后比对率提升约2%到3%，别小看这个数字，差异表达分析后期它可能就是显著性边缘的稻草。
3.4 输出格式和压缩技巧
默认输出是纯文本SAM，非常大。一个20G的FASTQ测序样本，SAM文件轻松上60G。我通常不会保留SAM，直接用samtools边比对边压缩：
hisat2 -p 8 --dta -x grch38/genome -1 sample_clean_R1.fastq.gz -2 sample_clean_R2.fastq.gz 2> sample.hisat2.log | samtools view -bS - > sample.bam注意这里2> sample.hisat2.log是关键的。hisat2的比对日志（reads数、比对率）会输出到标准错误，如果不用2>接住，它会和BAM二进制数据混在一起（实际不会混到文件，但屏幕会刷屏，且你看不到统计信息）。所以我把日志单独保存，后续要报告比对率就直接看这个log文件，非常方便。
4. 不利用循环：一条一条命令硬扛到底
现在进入标题里的正式议题。先说“不利用循环”的方式。这种方式适合样本数量很少、临时调试、或者想逐个确认每个样本比对是否正常的情况。
4.1 单样本逐条执行的完整流程
还是拿一个叫sample_A的样本举例。比对前的所有数据都在clean_data/目录下，后续输出统一放在hisat2_out/。第一步，逐条把比对跑完：
cd hisat2_out hisat2 -p 8 --dta \ -x ~/ref/grch38/genome \ -1 ../clean_data/sample_A_clean_R1.fastq.gz \ -2 ../clean_data/sample_A_clean_R2.fastq.gz \ -S sample_A.sam \ 2> sample_A.hisat2.log跑完检查log：
cat sample_A.hisat2.log里面会写总的reads数量、整体比对率、每类比对情况（比如多重比对、唯一比对）。确认没什么离谱的异常，再继续第二步，把SAM转成排序后的BAM：
samtools view -bS sample_A.sam > sample_A.bam samtools sort -@ 8 -o sample_A.sorted.bam sample_A.bam samtools index sample_A.sorted.bam这里顺序是view转格式、sort排序、index建索引。排序后的BAM是后续featureCounts和IGV可视化都需要的标准形态。当然你也可以直接管道一步到位：
samtools view -bS sample_A.sam | samtools sort -@ 8 -o sample_A.sorted.bam -然后rm sample_A.sam sample_A.bam释放空间。每次跑到这里，都得多留一个心眼：确认samtools sort成功退出之后再删除中间文件，否则你可能连重新跑比对的机会都没有，数据只能从头再来。
4.2 这种方式的坑在哪里
如果只有3个样本，这么敲没问题。样本到10个以上就会开始痛苦。我的一个师弟曾经拿着一份50个样本的数据，老老实实复制粘贴替换文件名，中间不小心把某个样本的R1和R2搞反了，最后定量出来那个样本所有基因表达量全乱了，排查了一下午。这就是“不利用循环”最大的问题：人不是机器，重复性操作中极易出错，而且出错后极难察觉。
另一个让我彻底放弃这种方式的原因是时间碎片化。每跑完一个样本，你得盯着屏幕等结果，再手动敲下一条命令。一旦晚上人不在，这个流程就断掉了。手工逐个跑，还意味着没法顺手汇总所有样本的报告，前一个比对完了你要手动记录，后一个还要再记一遍。把时间花在无意义的敲回车键上，对干我们这行是完全不划算的。
所以我一般只有两个场景才会用逐条的方式：一是只处理一两个快速测试样本；二是某个样本报错反复调试参数时，我单拎出来单独跑。不然一律走循环。
5. 利用循环：批量处理几十个样本的优雅姿势
下面重点聊聊“利用循环”的方式，这也是我实际工作中默认选用的方法。在Linux环境下，bash for循环就够用，不需要python脚本，也不需要在服务器上装额外的工具。核心只有一个思想：把“样本名列表”提取出来，然后在循环体内用变量去拼接所有相关的文件名。
5.1 先理清文件命名结构
我接手的数据通常是一个目录里放着sample_A_clean_R1.fastq.gz、sample_A_clean_R2.fastq.gz、sample_B_clean_R1.fastq.gz……一直到sample_Z。注意关键点：样本名要作为文件名的前缀，且前缀之间不能有歧义。样本名最好规范一致，比如SampleA就全是SampleA，不同样本不要有的叫Sample001有的叫Sample1，不然循环正则匹配容易出现误匹配。
写循环的关键是把变量提取出来。比如sample_A_clean_R1.fastq.gz这个文件，我想提取出sample_A这个前缀。最简单是这样：
for fq1 in ../clean_data/*_clean_R1.fastq.gz; do sample=$(basename "$fq1" _clean_R1.fastq.gz) echo "$sample" donebasename的作用是去掉路径和末尾的_clean_R1.fastq.gz，剩下的就是样本名。这一招是循环批处理的核心，理解了它，循环的骨架就立住了。
5.2 完整的多样本循环比对脚本
把比对、格式转换、排序整合进一个循环里，这是我长年在用的模板：
#!/bin/bash set -euo pipefail INDEX=~/ref/grch38/genome CLEAN_DIR=../clean_data OUT_DIR=./hisat2_out THREADS=8 mkdir -p "$OUT_DIR" for fq1 in "$CLEAN_DIR"/*_clean_R1.fastq.gz; do sample=$(basename "$fq1" _clean_R1.fastq.gz) fq2="$CLEAN_DIR/${sample}_clean_R2.fastq.gz" echo "===== Processing $sample =====" hisat2 -p "$THREADS" --dta \ -x "$INDEX" \ -1 "$fq1" -2 "$fq2" \ 2> "$OUT_DIR/${sample}.hisat2.log" \ | samtools view -bS - \ | samtools sort -@ "$THREADS" -o "$OUT_DIR/${sample}.sorted.bam" samtools index "$OUT_DIR/${sample}.sorted.bam" done这个脚本里有几个地方我特别说一下：
set -euo pipefail是bash脚本的“安全气囊”。-e让任何命令出错立刻退出，-u提醒未定义变量（这个对初学者非常友好，比如文件名拼错了，变量为空，脚本直接报错而不是瞎跑），-o pipefail确保管道中任何一个环节出错（不只是最后一条命令）都会让整个管道返回失败。别小看这三行，没有它们，你可能比对失败了却还在继续跑后面的排序，最后拿一个空BAM去做定量，等到结果异常白白耗掉几天时间。
为什么要用管道直接把hisat2输出接到samtools，而不是先生成SAM文件？因为SAM太大了。每个样本几十GB的纯文本真的可能把磁盘挤爆，管道流式处理可以避免中间文件的产生。注意管道里我给hisat2只留了-1 "$fq1" -2 "$fq2"，没加-S，因为没写-S时比对结果会输出到标准输出，正好喂给管道。我把标准错误用2>单独写到日志文件，这样既能看到进度和比对率，又不污染管道数据流。
循环跑完后，可以通过MultiQC汇总所有样本的hisat2日志，一键生成质量报告：
multiqc ./hisat2_out --filename RNAseq_hisat2_reportMultiQC的输出长得既直观又省事，比对率、reads数一目了然，写报告或者组会汇报时直接贴图。
5.3 别让循环卡死在“无样本”情形
有时候CLEAN_DIR目录写了通配符但是没有任何匹配文件，for循环仍然会执行一次，变量被设置成类似../clean_data/*_clean_R1.fastq.gz的原文，然后脚本就莫名其妙地拿一个不存在的文件名去跑。虽然-e会让它停下来，但错误信息会让新手摸不着头脑。我习惯在前面加一个检查：
shopt -s nullglob fq1_files=("$CLEAN_DIR"/*_clean_R1.fastq.gz) if [ ${#fq1_files[@]} -eq 0 ]; then echo "No R1 fastq files found in $CLEAN_DIR" >&2 exit 1 fishopt -s nullglob可以让没有匹配时通配符展开成空，然后检查数组长度。这是个小细节，但实际跑批处理时真的能救命。
5.4 循环过程中断如何断点续跑
另一个高频痛点：跑第15个样本时集群掉线了，或者内存不够被kill了，前面的14个样本已经产出sorted.bam，重头跑一遍又浪费时间。我的习惯是循环里加上判断，如果目标文件已存在，直接跳过：
for fq1 in "$CLEAN_DIR"/*_clean_R1.fastq.gz; do sample=$(basename "$fq1" _clean_R1.fastq.gz) out_bam="$OUT_DIR/${sample}.sorted.bam" if [ -s "$out_bam" ]; then echo "$sample already done, skip" continue fi ... done-s判断文件存在且非空，稳妥一点还可以配合samtools quickcheck验证BAM完整性。这个方法让我在大规模批处理时可以随时中断、随时续跑，不用担心中途意外影响最终产出。
6. 比对结果怎么看：SAM字段、比对率和常见判断
跑完了不等于完事了。拿到BAM之后，先别急着冲进定量，花几分钟把结果质量盘一盘。这一步做扎实了，后面会少很多麻烦。
6.1 SAM格式的核心字段速览
SAM文件（以及排序后的BAM）每行是一条比对记录。关键字段从左到右依次是：QNAME（read名字）、FLAG（比对标志位）、RNAME（比对到的染色体/contig）、POS（起始位置）、MAPQ（比对质量值）、CIGAR（比对操作串）。FLAG是个二进制数字的组合值，比如0表示双端read都比对上了，4表示这条read没比对到，16表示负链比对，99通常是read1正链且在正确配对方向上。CIGAR串则记录了连续的match、mismatch、insertion、deletion。初学不必深究每个数字，但至少要能看懂41这类异常flag提示可能存在非冗余或嵌合比对。
我排查比对问题时第一件事就是less -S打开BAM看看，或者用samtools view加过滤条件提取特定read。比如想看看某个基因区域里都哪些reads覆盖：
samtools view sample.sorted.bam chr1:1000000-2000000 | head -20这会快速展示该区域的对齐情况，阅读FLAG和CIGAR基本能判断是正常外显子连接还是跨内含子比对。
6.2 用samtools flagstat统计比对质量
最常用的快速统计命令：
samtools flagstat sample.sorted.bam输出里会看到total reads、mapped reads、properly paired、secondary、supplementary等数字。其中mapped的比例（比对率）是最直观的整体质量指标。以人全转录组普通polyA文库为例，干净数据比对到GRCh38后总比对率（包括多重比对）在90%以上比较正常，严谨一点看uniquely mapped占比。比对率低于80%就要重点排查。
另外，properly paired的比例也很关键，它表示双端reads正确配对且距离合理。这个指标如果偏低，大概率是文库插入片段异常、样品降解或参考序列错误。我见过有样本比对率还行但properly paired只有六成，后来发现是样本建库时RNA降解导致插入片段变短，read有一端比对超出片段范围而低质量，这个在数据预处理阶段肉眼是看不出来的。
6.3 比对率低最常来自哪些原因
总结一下我实际排查过的案例：
参考基因组与样本物种不匹配。这个听着低级但真有人犯过，跑完漂亮地输出“0% overall”，然后才意识到样本是斑马鱼，参考索引用的却是人类。解决方式是引入MultiQC前确认物种信息，或者比对前用
kraken2做一个快速的物种分类注释。至少要在拿到log时就怀疑，而不是三个月后才发现。接头污染和高比例低质量reads。没用fastp或者参数用错了，比对率会低好几个点。检查fastp报告里passed filter reads的比例，如果低于70%，大概率是原始数据本身质量问题或者质控参数不合理。
rRNA污染。未做rRNA去除的文库会有大量reads比对到rDNA区域，看起来整体比对率不低，但实际有效信息很少。可以在比对后用featureCounts统计rRNA区域的reads占比，超过10%就该怀疑建库过程有问题。
数据库版本不匹配。GTF和基因组版本对不上，跨外显子的reads找不准剪接位点，mapped里会出现很大比例的“mapped but not properly paired”。这种情况换官方配套版本即可。
嵌合reads/融合转录本。部分肿瘤样本会发生大量基因融合，导致reads无法唯一比对到单一区域，这属于生物层面的正常现象，不用纠结。
7. 常见问题与排查技巧实录
最后集中整理一些我在实操里踩过并且花费时间不下两小时的报错和怪异行为。
7.1 hisat2-build构建索引报错或内存不足
内存不足几乎是上服务器时最容易碰到的情况。人类全基因组索引建议至少40GB可用内存，低于这个阈值就很容易被Out of Memory杀死。解决办法很直接：直接去官方下载预构建索引，不要自己build。如果研究的物种非模式生物，没有预构建索引，可以只用“转录组外显子序列”来构建，或分染色体分批构建然后合并，这个操作更复杂，但至少能跑完。
常见的另一个报错是“chromosome names end with \r”或“Encountered a non-ACGT character”。这通常是因为下载的FASTA有Windows换行符，或者包含U（RNA基因组序列）和N等模糊字符。办掉它的方式是：
tr -d '\r' < genome.fa > genome_unix.fa然后再跑hisat2-build。
7.2 索引版本不匹配引发“Error: Invalid index file”
这个报错常见于把GRCh37的数据和GRCh38的索引混用，或者用旧版本hisat2去加载新版本索引。虽然提示信息很长很吓人，但本质就是版本冲突。我的习惯是每次都把索引的下载链接和版本号记在preprocessing脚本备注里，确保所有样本都在同一个版本环境下比对。
7.3 循环里面变量名拼接出错
我见过有人把样本名带出来了，但是在循环里写着$fq1_dir/${sample}_clean_R1.fastq.gz却忘加引号，文件名一旦有空格就碎成多段。所以我的建议：
"$CLEAN_DIR/${sample}_clean_R1.fastq.gz"所有变量引用都加双引号，这是一个能避开无数隐形崩溃的好习惯。
7.4 管道中退出码被吞掉
如果你没有set -o pipefail，hisat2因为某种原因没成功（比如索引路径错误），但samtools view可能还是生成了一个几乎空的BAM，最终管道退出码为0。脚本不会报错，你也不会注意到。我自己的习惯是不仅加set -euo pipefail，还会在循环里加一行判断BAM大小：
if [ ! -s "$out_bam" ]; then echo "ERROR: $sample BAM is empty" >&2 exit 1 fi这样问题在第一时间暴露，而不是留到定量后才发现。
7.5 单端reads怎么处理
不是所有数据都是双端。如果实验是单端测序（比如某些miRNA或small RNA），hisat2命令只要把-1/-2换成-U：
hisat2 -p 8 --dta -x grch38/genome -U sample_clean.fastq.gz -S sample.sam其余处理流程相同。只是后续定量时，featureCounts也要记得不设-p参数，否则默认统计双端时会把片段计数逻辑搞错。
把比对这关打通之后，后面定量和分析就水到渠成了。之前提过这系列“未完待续”，是因为后面还有samtools处理细节、featureCounts定量、DESeq2差异分析以及结果可视化，都可以单独拉出来写好几篇。我自己跑完一个完整项目，最大的感受是：在这个环节上花的力气绝对不会白费。前期多花点时间理解参数，比后期对着一个谁也说不出哪来的异常表达矩阵挠头，要有用得多。
