Linux下RNA-seq STAR比对与cufflinks拼接避坑指南

跑完质控之后的下一步，很多人会习惯性甩一句“直接上 STAR 比对吧”，然后被 30G 内存、几小时的索引构建、一堆看不懂的参数直接劝退。Linux系统下RNA-seq分析里，STAR比对和cufflinks拼接这两步，是整条流程中最吃硬件、也最容易在细节上翻车的环节：索引建错版本、线程开太多把机器拖死、BAM 没排序就往下走、拼接出来一堆假阳性转录本，最后回过头查问题时已经过去两天。这篇内容讲的就是从 FASTQ 到有参转录本拼接这一段路的完整落地方法，包含索引构建的内存测算、比对参数的逐项取舍、日志指标怎么看、以及 cufflinks 目前在实践中的真实定位和替代方案。适合已经能在 Linux 下敲命令、跑过质控、但还没把比对拼接流程真正跑通跑稳的人，也适合那些流程能跑但结果总感觉不对劲、想把每个参数的理由搞清楚的人。文中所有命令和参数都可以直接抄到自己服务器上改一改就用，踩过的坑我会一并写出来。

1. 把读段放回基因组：splice-aware 比对解决的到底是什么问题

1.1 RNA-seq 读段天生跨越外显子边界

DNA 测序的读段来自连续的基因组片段，比对时只要在参考序列上找到一段高度相似的区间就行。RNA-seq 不是这样。成熟 mRNA 经过剪接，外显子在转录本上是首尾相连的，但映射回基因组时中间隔着一整段内含子。一条 150bp 的双端读段，完全可能前半段落在外显子 A 的末尾，后半段落在外显子 B 的开头，而在基因组坐标上这两个位置相隔几千甚至十几万碱基。

这就意味着普通的、不允许大段缺口的比对算法在这里是失效的。你把这样的读段丢给一个只做错配和短 indel 的比对器，结果要么是比对不上被丢掉，要么是强行劈成两半各比对一段但拼不回一个转录本。真正的做法是让比对器知道“这里有剪接”这件事，允许在读段中间开一个巨大的 gap，并且这个 gap 的边界要尽量贴合已知的剪接位点基序。这类算法统称 splice-aware alignment，STAR 就是其中用得最广的一个。

理解这一点很重要，因为它决定了后面所有参数的默认值为什么长那样。比如--alignIntronMin默认是 20，因为小于 20bp 的间隔基本可以认为是 indel 或者比对噪声，不可能是真正内含子；--alignIntronMax默认给了很大的值，因为有些基因的内含子确实能长到几十万碱基。这些默认值不是拍脑袋来的，是对生物学事实的编码。

1.2 STAR 的两步走策略：种子搜索加剪接拼接

STAR 的全称是 Spliced Transcripts Alignment to a Reference，它的核心思路可以拆成两层。第一层是把整个参考基因组建成一个压缩的后缀数组结构（suffix array），再用最大可映射前缀（Maximal Mappable Prefix，MMP）去搜索。你可以把它想象成一本按字母顺序排好序的字典，任何一段序列都能在极短时间内找到它在这本字典里最长的那段前缀位置。第二层是当一条读段的前缀匹配到某个位置后，剩余的部分接不上了，STAR 会把剩下的序列当作一条新读段重新搜索，然后把两段结果拼起来，如果拼接点符合剪接信号，就认为这是一条跨剪接位点的比对。

这个策略的好处是快，代价是索引本身非常大。后缀数组加上各种辅助结构，人类基因组建完索引通常要占 27 到 30 GB 磁盘，运行时还要把这套结构读进内存。所以 STAR 对内存的胃口，本质上是算法换来的，不是实现得粗糙。

还有一点值得单独提：STAR 在建索引时可以把注释文件里的剪接位点一起编进索引（--sjdbGTFfile），这样比对阶段遇到已知剪接位点时优先级更高、更准。这也是为什么建索引时必须准备好 GTF，而不是只给一个 FASTA 就完事。

1.3 为什么很多人最终选了 STAR 而不是其他方案

转录组比对工具不止一个，HISAT2、TopHat 都是常见选择。TopHat 现在基本不用了，因为它底层依赖 Bowtie，速度慢，而且对多外显子转录本的处理方式比较绕。HISAT2 内存占用小很多，人类基因组大概 8 GB 左右就能跑，对硬件紧张的机器很友好。

那为什么还是大量流程用 STAR？主要是两个原因。一是速度快，在足够内存和多核的条件下，STAR 处理一个 3000 万读对的样本通常在一小时以内，比 HISAT2 快不少。二是它对剪接位点的检测更激进，能发现一些注释里没有的新剪接事件，这对做可变剪接分析的人很重要。代价就是内存。如果你的服务器只有 32 GB 内存，要跑人类样本的 STAR，基本上只能单样本串行，而且还得把 BAM 排序的内存控制住，后面会详细讲。

选择逻辑其实很朴素：机器内存够（≥64 GB）、样本量大、追求速度，用 STAR；机器内存紧、样本不多、能接受慢一点，HISAT2 更稳妥。两者输出的 BAM 在下游处理上完全通用，所以这个选择不影响后面的拼接步骤。

2. 建索引之前的账：内存、注释版本和参数三件事

2.1 参考基因组和注释文件必须是同一个版本

这是最容易被跳过、也最容易导致后面全盘返工的一步。参考基因组 FASTA 和 GTF 注释必须来自同一个版本发布。你可能觉得“染色体名字对得上就行”，实际上不同版本之间染色体坐标会有细微调整，外显子边界会变，有些基因名会合并或拆分。如果 FASTA 用的是某个版本的组装，GTF 用的是另一个版本，STAR 在建索引阶段可能不报错，比对也能跑完，但剪接位点会错位，拼接出来的转录本坐标全是偏的，你拿去和别人的结果对比就是对不上。

更隐蔽的一种情况是染色体命名风格不一致。有的 FASTA 里染色体叫chr1，有的叫1。这两种命名混用，STAR 会在日志里警告说注释里有大量序列找不到对应染色体，然后默默忽略掉整份注释。你如果只看最后“比对完成”几个字，根本发现不了。

我的习惯是在建索引前先做两个检查：用grep ">" genome.fa | head看一眼染色体命名，再用cut -f1 annotation.gtf | sort -u | head看一眼注释的命名，两边对一遍。花三十秒，能省两天。

grep ">" reference.fa | head -5 awk '$3=="exon"{print $1}' annotation.gtf | sort -u | head -5

2.2 genomeGenerate 的内存需求和参数匹配

建索引这一步的内存消耗，经验值是基因组大小的 8 到 10 倍左右。人类基因组约 3.1 Gb，需要 30 GB 上下；小鼠类似；植物里小麦这种十几 Gb 的基因组，就需要上百 GB 内存，普通服务器直接放弃，得换分染色体策略或者用内存更省的工具。

一个常被忽略的参数是--genomeSAindexNbases。默认值 14 是给人类这个量级的基因组准备的。如果你做的是真菌、细菌或者某个小型生物的基因组，基因组只有几十 Mb，用 14 会导致索引文件异常大、比对速度反而下降。官方给的调整公式是取min(14, log2(GenomeLength)/2 - 1)的整数部分。比如基因组 30 Mb，log2(30e6) 约等于 24.8，除以 2 减 1 约等于 11.4，那就设成 11。

另一个是--sjdbOverhang，默认 100。这个参数的含义是“剪接位点两侧用于构建注释序列的序列长度”，理论上应该设成 read length 减 1。你如果做的是 150bp 双端，就设 149。如果建索引时用了 100，而实际读段长 150，STAR 会照常跑，但剪接位点附近的比对灵敏度会略降一点。这个差异在小样本上不明显，在大样本差异表达分析里可能会体现在低表达基因的检出率上。既然知道了，就顺手改对。

STAR --runMode genomeGenerate \ --runThreadN 16 \ --genomeDir ./star_index \ --genomeFastaFiles reference.fa \ --sjdbGTFfile annotation.gtf \ --sjdbOverhang 149 \ --genomeSAindexNbases 14

建索引的耗时也和线程强相关。16 线程跑人类基因组，大概 40 分钟到 1 小时；8 线程可能要 1.5 小时以上。这一步是一次性的，建完就能反复用，所以值得多给点线程。

2.3 索引目录里的产物和复用策略

--genomeDir指定的目录里跑完会出现一批文件，常见的有Genome、SA、SAindex、chrName.txt、chrLength.txt、genomeParameters.txt。这几个文件是一个整体，缺一不可，拷贝索引时必须整个目录一起拷，单独拷SA过去是用不了的。

一个实用技巧：索引建好后把整个目录打个压缩包存着，同时在目录里放一个 README 记录这次用的 FASTA 和 GTF 的版本号、md5、建索引的完整命令。半年后你回头跑新样本，看到一堆叫star_index的目录，没有这份记录根本分不清哪个对应哪个版本。我踩过这个坑，最后把三份索引全删了重建。

还有一点，同一次项目里所有样本必须用同一份索引。中途觉得“注释更新了，重新建一份”，然后新旧样本混着用，下游拼接阶段会直接报错或者给出无法解释的结果。

提示：建索引前先确认可用内存，用free -g看一眼。如果可用内存不到基因组大小的 10 倍，先把其他任务停掉，或者换到内存更大的节点上做。索引建到一半被 OOM killer 杀掉，日志里通常只有一行 Killed，很容易误判成程序 bug。

3. 单样本比对实操：每个参数为什么这么设

3.1 输入处理：压缩文件直接喂给 STAR

大部分测序数据交付时都是.fastq.gz，STAR 支持直接读压缩文件，靠--readFilesCommand zcat这一行。这样省掉了先解压再比对的磁盘占用和 IO 时间，尤其是样本多的时候差别明显。如果是.bz2就换成bzcat，如果是不压缩的直接去掉这个参数。

双端数据的输入写法是把两个文件依次放在--readFilesIn后面，第一个是 R1，第二个是 R2。顺序搞反不会报错，STAR 会照常跑完，但比对率会掉得很惨，因为它默认 R1 和 R2 来自同一片段的两端、方向相反。如果你的文库是特殊建库方式（比如某些单链建库），还需要额外加--outSAMstrandField intronMotif或者调整--readFilesIn的顺序。普通双端文库不用管这些。

输出前缀--outFileNamePrefix建议带上样本名和斜杠，比如./aln/SAMPLE_A_，这样所有输出文件都会以这个前缀开头并落到指定目录。不加的话默认输出到当前目录，文件名是Aligned.sortedByCoord.out.bam这种通用名，样本一多必然互相覆盖。

3.2 关键参数逐项拆开看

真正影响结果的参数就那么十来个，其余大部分用默认值即可。下面这张表是我在项目里固定使用的配置，每一项后面写了理由。

这几个比例类参数（0.66、0.04）是 STAR 作者给的推荐值，在常规真核转录组数据上表现稳定。什么时候需要改？如果你的物种和参考基因组亲缘关系较远（比如用近缘物种的基因组做参考），错配比例得放宽到 0.06 甚至 0.08，否则比对率会低得没法看。反过来，如果做的是同一品系的回交群体，可以收紧到 0.03 以降低假阳性。

--alignSJoverhangMin这个参数值得单独说。它控制的是“一条读段跨越剪接位点时，两侧各需要多少碱基锚定才算数”。默认 8 是一个平衡点。如果你做的是可变剪接精细分析，想找更多新剪接事件，可以降到 5，但要做好心理准备：假阳性会明显增加，需要下游用更严格的过滤条件去筛。反过来做常规差异表达，保持 8 就好，甚至提到 10 更稳。

3.3 排序内存：STAR 最容易 OOM 的地方

--outSAMtype BAM SortedByCoordinate这一行很方便，但它背后是 STAR 在内存里做排序。默认情况下 STAR 会尽可能多占内存来做这件事，如果在共享服务器上，很容易把整台机器拖垮或者被系统杀掉。

解决办法是显式限制排序内存：--limitBAMsortRAM 8000000000，单位是字节，8 GB。这个值和基因组索引占用的内存是分开算的。人类基因组下，索引本身占 27 GB 左右，排序再给 8 GB，单样本总共 35 GB 上下。服务器 64 GB 内存的话，理论上能并行两个样本，但实际不建议，因为系统缓存和其他进程还要占，并行跑两个很容易互相挤压导致两个都变慢。稳妥做法是一次跑一个，用多线程把单个样本的速度拉满。

如果内存实在紧张，有个替代方案：改成--outSAMtype BAM Unsorted，让 STAR 只输出未排序的 BAM，然后单独用 samtools 排序，samtools sort -@ 8 -m 2G。samtools 的排序是外部排序，会把中间结果写到临时文件，内存占用可控，代价是多一次磁盘读写。数据量大的时候这一步可能多花十几分钟，但换来的是不会 OOM。

STAR --runThreadN 12 \ --genomeDir ./star_index \ --readFilesIn ./clean/SAMPLE_A_R1.fq.gz ./clean/SAMPLE_A_R2.fq.gz \ --readFilesCommand zcat \ --outFileNamePrefix ./aln/SAMPLE_A_ \ --outSAMtype BAM SortedByCoordinate \ --outSAMattributes NH HI AS NM MD \ --outFilterMismatchNoverReadLmax 0.04 \ --outFilterMultimapNmax 10 \ --alignSJoverhangMin 8 \ --alignSJDBoverhangMin 1 \ --alignIntronMax 1000000 \ --limitBAMsortRAM 8000000000

4. 比对日志怎么读，以及比对率低时从哪查起

4.1 Log.final.out 里真正值得看的几行

STAR 跑完会在输出前缀后面生成Log.final.out，里面有二三十个指标，大部分可以一眼扫过。真正需要盯的是这几项：

这里要注意“一般参考范围”的措辞。不同物种、不同组织、不同建库方式差异很大。同一个实验室同一批样本，比对率应该彼此接近；如果其中一个样本突然比别的低 20 个百分点，那才是真正需要查的信号。绝对数值本身参考意义有限。

还有一个经常被忽略的指标是Number of input reads。这个数应该和你质控后保留的读段数对得上。如果对不上，说明输入文件有问题，可能是 R1 和 R2 行数不一致，或者压缩包损坏导致 STAR 提前停止读取。

4.2 比对率异常时的逐层排查顺序

我一般按这个顺序查，从最省时间的开始：

第一层，看 R1 和 R2 是不是搞反了。这个错误极其常见，把两个文件对调再跑一个样本，比对率如果从 40% 跳到 85%，问题就找到了。检查方法不难，看 R1 文件的前几条序列，正常 R1 的起始碱基分布比较集中（取决于建库方式），但这需要经验。最快的验证还是直接对调跑一遍。

第二层，看接头和低质量碱基残留。用 fastqc 看一眼质控报告里的 adapter content 和 per base sequence quality。如果接头占比明显，说明质控那一步的修剪没做干净，回去重跑 trimming。

第三层，看 rRNA 残留。总 RNA 建库如果没有做 rRNA 去除或者 polyA 富集，rRNA 能占到六七成。这种情况下比对率可能还行，但有效数据量少得可怜。判断方法是从 BAM 里统计比对到 rRNA 区域的读段比例，或者直接用比对率乘以“非 rRNA 比例”估算有效数据。

第四层，看参考基因组版本。如果是新物种或者非模式生物，参考组装质量可能不高，N 碱基多、组装碎片化，这种情况下比对率天然偏低。可以考虑换一个组装版本，或者用转录组本身做参考。

第五层，看文库类型。链特异性建库如果参数没对应上，或者是单细胞数据、降解严重的 FFPE 样本，比对率都会受到影响，需要针对性调整参数而不是一味怀疑数据。

4.3 BAM 的排序、索引和快速核查

STAR 输出排序好的 BAM 之后，还差一个索引文件才能被下游工具快速访问。这一步用 samtools：

samtools index ./aln/SAMPLE_A_Aligned.sortedByCoord.out.bam

同时跑一下samtools flagstat做交叉验证，它给出的是从 BAM 层面统计的比对情况，可以和 STAR 日志互相印证。两者数字如果差得很多，说明中间某个环节出了问题。

samtools flagstat ./aln/SAMPLE_A_Aligned.sortedByCoord.out.bam samtools idxstats ./aln/SAMPLE_A_Aligned.sortedByCoord.out.bam | head -20

idxstats特别有用，它能按染色体列出比对的读段数。人类样本里线粒体染色体上的读段比例如果特别高（比如超过 20%），说明样本可能有降解或者线粒体污染，下游分析时要考虑是否剔除。这个信息在 STAR 日志里是看不到的。

5. 拼接这一步的现实选择：cufflinks 还是 StringTie

5.1 cufflinks 在今天的定位要说清楚

cufflinks 曾经是转录本拼接的标准工具，配套的 cuffcompare、cuffmerge、cuffdiff 组成了一整套分析链路。它的工作原理是基于比对上的读段构建一个“重叠图”，然后在这个图上寻找最小路径来推断转录本结构。在有参考注释（-G参数）的情况下，它会倾向于沿用已知的转录本结构，只对读段支持充分的新结构做补充。

问题是这套工具在 2016 年之后基本停止更新了。停止维护带来的直接后果是：遇到新版本的 BAM 格式或者新的注释文件格式，可能出现各种兼容性问题；在多线程和大数据量下的表现也不如后来的工具。更关键的是，cufflinks 对低表达转录本不敏感，无参拼接（不给-G）时的假阳性率相当高，经常给你输出几万个“新转录本”，其中大部分是噪音。

但 cufflinks 并没有完全消失，还有一些老项目和发表过的流程在用它，所以知道怎么跑、怎么理解它的输出，仍然是有价值的。如果你是新起一个项目，我建议直接用 StringTie，它出身于同一批人，思路类似但实现更现代。如果是复现老流程或者接手别人的代码，那就照着原来的跑，但心里要清楚它的局限。

5.2 用 cufflinks 跑一遍有参拼接的完整流程

先给最基础的命令，还是接在 STAR 输出的 BAM 后面：

cufflinks -p 8 \ -G annotation.gtf \ -o ./cufflinks_out/SAMPLE_A \ -b reference.fa \ ./aln/SAMPLE_A_Aligned.sortedByCoord.out.bam

-b参数给参考 FASTA，是为了做偏差校正（bias correction），它会统计读段在转录本不同位置的分布，纠正一些建库和测序引入的位置偏好。这个参数在样本量小的时候效果有限，但加了没有坏处。

-p是线程数。这里要提醒一句，cufflinks 的多线程效率和 STAR 完全不是一个量级，开 8 线程和开 16 线程的差距可能只有 20%，主要瓶颈在算法本身。所以别指望靠加线程把速度拉起来。

多个样本跑完后需要合并。cufflinks 的做法是先为每个样本生成一个 GTF，然后用cuffcompare和cuffmerge合成一个统一的转录本集合：

ls ./cufflinks_out/*/transcripts.gtf > assemblies.txt cuffmerge -g annotation.gtf -s reference.fa -p 8 assemblies.txt

cuffmerge会过滤掉一些明显不可信的转录本，比如那些没有对应 reads 支持的、长度过短的、或者只在单一样本出现的。但过滤规则比较黑盒，输出结果需要自己再检查一遍。

内存方面，cufflinks 跑人类基因组大概要 6 到 8 GB，比 STAR 客气不少，但样本多的时候并行跑也会把内存吃满。建议同时最多跑三到四个样本。

5.3 “有参拼接”和“从头组装”是两个概念

这里必须澄清一个常见混淆。STAR 加 cufflinks 这条路，常常被叫做“无参拼接”或者“de novo assembly”，这个说法是不准确的。STAR 需要参考基因组才能比对，cufflinks 是在比对结果的基础上推断转录本结构，所以整条路是基因组引导的（genome-guided），属于有参分析的范畴。

真正的从头组装（de novo assembly）是不依赖任何参考基因组，直接从读段之间的重叠关系重建转录本，常用工具是 Trinity 这类。这条路用在完全没有参考基因组的物种上，计算量比有参分析大一到两个数量级，内存需求也高得多，而且组装出来的序列没有染色体坐标，下游注释全靠自己。

区分清楚这一点很关键，因为它决定了你能不能直接用别人的流程。如果某篇文献说“我们做了 de novo assembly”，然后方法里写的是 STAR + cufflinks，那是写错了。如果你的物种有靠谱的参考基因组，就走有参这条路，成本和可靠性都好得多。

5.4 换成 StringTie 之后的流程差异

StringTie 的命令结构几乎和 cufflinks 一一对应，迁移成本很低：

stringtie -p 8 \ -G annotation.gtf \ -o ./stringtie_out/SAMPLE_A.gtf \ -l SAMPLE_A \ -b ./ballgown_out/SAMPLE_A \ ./aln/SAMPLE_A_Aligned.sortedByCoord.out.bam

-l给样本前缀，会写进输出的 GTF 里，方便后面区分来源。-b是可选参数，输出 ballgown 格式的中间文件，如果后面要做可变剪接分析可以用来做可视化。

多样本合并用--merge：

ls ./stringtie_out/*.gtf > merge_list.txt stringtie --merge -G annotation.gtf -p 8 \ -o ./stringtie_out/merged.gtf \ merge_list.txt

两者在实际项目中的差异，我整理成下面这张表，仅代表我自己的使用体感，不是绝对结论：

换成 StringTie 之后最大的体会是速度。同一个样本，cufflinks 要跑四十分钟的，StringTie 十来分钟就出结果，而且结果的一致性看起来更好。当然这只是体感，真要比较灵敏度需要设计对照实验。

6. 批量跑和踩过的那些坑

6.1 用 shell 循环把样本串起来

样本一多，逐个敲命令肯定不现实。最简单的办法是写一个样本列表文件，然后 for 循环读进去。关键在于从文件名里正确提取样本名，这决定了后面输出的文件名会不会互相覆盖。

#!/bin/bash set -euo pipefail INDEX=./star_index CLEAN=./clean OUT=./aln THREAD=12 SORT_RAM=8000000000 mkdir -p "$OUT" for R1 in "$CLEAN"/*_R1.fq.gz; do SAMPLE=$(basename "$R1" _R1.fq.gz) R2="$CLEAN/${SAMPLE}_R2.fq.gz" if [[ ! -f "$R2" ]]; then echo "缺少 R2: $SAMPLE" >&2 continue fi STAR --runThreadN $THREAD \ --genomeDir $INDEX \ --readFilesIn "$R1" "$R2" \ --readFilesCommand zcat \ --outFileNamePrefix "$OUT/${SAMPLE}_" \ --outSAMtype BAM SortedByCoordinate \ --limitBAMsortRAM $SORT_RAM samtools index "$OUT/${SAMPLE}_Aligned.sortedByCoord.out.bam" done

几个细节值得说。set -euo pipefail让脚本在任何一个命令失败时立刻停下，而不是默默往下跑导致后面的样本全错。basename "$R1" _R1.fq.gz用后缀剥离的方式取样本名，比用cut或sed稳妥，不会因为样本名里带下划线而出错。循环里先检查 R2 是否存在，因为交付的数据里偶尔会缺几个 R2 文件，不检查的话 STAR 会收到一个空路径然后报一堆看不懂的错。

6.2 内存溢出时的现象和判断方法

被系统 OOM killer 杀掉的进程，在 shell 里的表现通常是Killed或者Exit code 137。137 这个数字是 128 加 9，9 是 SIGKILL 的信号编号。看到这个基本可以断定是被系统强行终止的，不是程序自己崩的。

判断到底是索引内存不够还是排序内存不够，可以看 STAR 日志里的最后一行。如果日志里已经开始输出比对进度然后突然停住，多半是排序阶段的内存问题；如果一开始就停住，那可能是索引加载阶段的问题。相应地，前者调小--limitBAMsortRAM或者换BAM Unsorted，后者就得换机器了。

还有一种情况是没被杀掉，但系统开始疯狂使用 swap，整台机器慢得像卡死。这时候free -g会看到 swap 使用量飙升，vmstat 1能看到 si 和 so 两个列持续非零。这种状态比直接杀掉更糟，因为会拖慢所有并行任务。解决办法是监控起来，跑大任务前先看一眼内存余量，留出至少 20% 的缓冲。

6.3 输出目录的组织方式

批量跑完之后最容易乱的是输出目录。每个样本 STAR 会生成十来个文件，cufflinks 或 StringTie 又有各自的目录结构，几十个样本堆在一起，找东西全靠ls加grep。

我的组织方式是三级结构：原始数据一个目录、比对结果一个目录、拼接结果一个目录。比对结果里所有样本平铺，靠前缀区分；拼接结果里每个样本一个子目录。另外单独建一个logs目录，把每个样本的 STAR 日志和拼接日志收集一份进去，方便后续统一 grep 检查比对率。

再进一步，可以用一个简单的汇总脚本把所有样本的Log.final.out里关键指标抽出来做成一张表，一次性看完全部样本的比对情况：

for f in ./aln/*_Log.final.out; do SAMPLE=$(basename "$f" _Log.final.out) UNIQ=$(grep "Uniquely mapped reads %" "$f" | awk -F'|' '{print $2}' | tr -d ' ') MULTI=$(grep "% of reads mapped to multiple loci" "$f" | awk -F'|' '{print $2}' | tr -d ' ') echo -e "${SAMPLE}\t${UNIQ}\t${MULTI}" done

把这个输出重定向到文件里，用表格软件打开，一眼就能看出哪个样本是离群的。这比逐个打开日志文件效率高得多。

6.4 时间预期和资源规划的参考值

最后给一组我这边实际跑出来的经验值，供规划资源时参考。人类样本，3000 万读对，150bp 双端：

这套数字在小鼠、斑马鱼这类基因组小一些的物种上会成比例下降，在植物大基因组上会显著上升。真正规划的时候，把单样本的峰值内存乘以你想并行的样本数，再留 20% 余量，就是服务器的内存需求。硬盘方面，一个样本的 BAM 大概 2 到 4 GB，加上索引和其他中间文件，一百个样本的项目准备 1 TB 空间是比较稳妥的。

我在实际项目里最深刻的体会是，比对和拼接这两步的成败，八成取决于前期准备而不是命令本身。参考基因组的版本核对、索引参数和读段长度的匹配、内存的提前测算，这三件事在开工前花半小时搞清楚，能避免后面两天的返工。剩下两成才是参数微调，而且大部分情况下默认值就够用，不要为了调参而调参。真遇到比对率异常，按 R1/R2 顺序、接头残留、rRNA 污染、参考版本这个顺序往下查，基本都能定位到原因。