在 Linux 上算宏基因组组装结果的 N50 和 N90看着像件小事一个 awk 加一个 sort 就完事了。可真到了要拿去汇报、写方法学、跟合作方对齐数字的时候你会发现两个人的结果对不上是常态——有人说 N50 是 12kb有人跑出来是 9.8kb还有人报的是 scaffold 的 N50你手上这份是 contig 的。数字本身不难算难的是算对、算清楚、算得能被人复现。这篇就把这件事从定义到脚本到踩坑完整走一遍怎么在纯命令行环境里从 FASTA 抠出长度、怎么排序累加定位 N50/N90/L50/L90、怎么用现成工具交叉验证、批量样本怎么汇总成一张表以及在宏基因组这种物种混杂、丰度悬殊的场景下这些数字该怎么读。刚接触组装的同学可以照着抄脚本做过几轮的老手可以重点看交叉验证和排查那两节。1. 先把 N50/N90 的定义固定下来避免后续对数字的理解不一致1.1 从一堆长短不一的碎片说起宏基因组组装的产出本质上是一大堆长度参差的序列片段。测序读长只有几百 bp组装器把它们拼接成更长的 contig短的可能只有两三百 bp长的能到几十万 bp。你没法用平均长度来评价这批东西的质量因为平均长度会被极端值带偏一条 500kb 的长 contig 扔进十万条短片段里平均值立刻失真。所以领域内约定俗成的做法是看分位数——把序列按长度从大到小排好队然后从上往下累加长度累计到总长度的某个百分比时当前那条序列的长度就是对应的 Nx 值。N50 就是累计到总长度 50% 时那条序列的长度。举个例子假设有 6 条 contig长度分别是 100、80、60、40、20、10 kb总长 310kb一半是 155kb。从大到小累加100累计 100不够80180超过 155停。此时停下时的那条长度是 80kb所以 N50 80kb。注意不是 180也不是 155而是触发阈值的那一条自身长度。这一点是最容易写错的地方很多手写脚本在这里出错。这套逻辑的价值在于它对长序列敏感对海量短序列不敏感。组装质量好长 contig 多累计到 50% 需要的条数少N50 就高组装碎全是几千 bp 的短片段得累加几千条才够一半N50 就低。所以 N50 能比较真实地反映你拼出来的东西到底有多完整。1.2 L50 和 N90 各自回答什么问题N50 单独看有时不够一般会配三个伴生指标一起看。L50 是累计到 50% 时用掉了多少条序列它和 N50 是一体两面N50 高、L50 小说明少数几条长序列就撑起了半边天组装连贯性好N50 和 L50 同时变大说明你有很多中等长度的序列在贡献总量属于碎但均匀。N90 是累计到 90% 时的序列长度它衡量的是尾部——也就是覆盖绝大多数碱基需要多长的序列。N90 和 N50 之间的差距很有信息量两者接近说明长度分布集中没有特别离谱的长尾两者差得远说明组装结果里是少数超长 海量短碎片的混合形态这在宏基因组里非常典型因为高丰度物种能拼出较长的 contig低丰度物种就只能留下零碎片段。L90 同理是覆盖 90% 所需的序列条数。还有一个常被忽略的指标是最大序列长度和总长度这四个数放一起总长、最大长度、N50、N90基本能把一个组装结果的骨架描述清楚。指标含义数值大/小代表什么宏基因组中的典型量级总长度所有序列长度之和与数据量、组装冗余度相关几百 Mb 到几 Gb最大长度最长那条序列的长度大说明存在高丰度、低复杂度的优势物种几十 kb 到几百 kbN50累计 50% 时的序列长度大说明组装连贯通常 5kb~50kbL50累计 50% 用掉的序列条数小说明长序列贡献大数百到数千N90累计 90% 时的序列长度与 N50 差距小说明分布均匀通常 1kb~5kbL90累计 90% 用掉的序列条数反映碎片化程度数千到数万1.3 宏基因组组装为什么不能参照单基因组的评价标准这是我在实际项目里反复要解释的一件事。做单基因组组装的时候N50 动不动就是几 Mb因为整个基因组是同一个体的同一套序列重复区处理好就能拼得很长。宏基因组完全不同一个样本里可能有几百上千个物种丰度从 0.01% 到 30% 不等低丰度物种的测序深度可能只有 2×~5×根本不足以跨越重复区组装器只能保守地断开。所以你在宏基因组结果上看到 N50 8kb 就慌是完全没必要的。判断宏基因组组装好坏重点不在 N50 的绝对值而在于总长度是否合理跟预期物种总基因组大小比对、N50 与 N90 的比例关系、以及后续分箱binning能不能拿到高质量 MAG。有些项目为了追求高 N50把短于 1000bp 的 contig 全丢掉N50 立刻涨上去但总长度掉了一大截物种多样性信息也一起丢了——这种美颜式统计在评审里是要被追问的。我的习惯是任何一份宏基因组组装结果的 N50 数字必须同时附上三个限定条件——是 contig 还是 scaffold、过滤阈值是多少、总长度是多少。缺一个这个数字就没法横向比较。2. 数据准备进脚本之前先做三件事2.1 先确认手里是 contig 文件还是 scaffold 文件组装器通常会输出两个层级的文件contig 是仅由读长直接拼接、没有 gap 的连续序列scaffold 是在 contig 基础上利用配对信息把方向、顺序确定下来中间用一串 N 连接。两者算出来的 N50 差距可能很大因为 scaffold 的长度包含了 N 占位而且会把多条 contig 串成一条。确认方法很直接看 FASTA 的 header 和序列里有没有连续的 N# 看前两条序列的 header grep ^ assembly.fasta | head -2 # 粗看是否存在长串 Nscaffold 的典型特征 awk /^/{next}{ngsub(/N/,N); if(n50) c} END{print 含长N串的行数:, c0} assembly.fasta如果命名的后缀是_contig、_len、_cov这类一般就是 contig 级如果是_scaffold或者序列中间出现几十上百个连续 N那就是 scaffold 级。报数字的时候一定要标注清楚contig N50和scaffold N50是两个不同的指标混用会让整个表格失去意义。另外还要注意一个陷阱有些流程会把多个样本的 contig 合并到一个文件里或者把不同 k-mer 参数的组装结果 cat 到一起。这种文件算出来的 N50 是混合体没有生物学意义。跑之前用grep -c ^数一下条数心里有个数。2.2 FASTA 格式的五个常见脏数据我遇到过的格式问题基本集中在五类任何一个都会让长度统计悄悄出错序列跨多行这是正常情况但脚本必须把同一 header 下的多行长度累加不能只看第一行。很多随手写的脚本用length($0)在 header 行之后只取一行结果每条序列都被截断。行尾是 CRLF从某些系统拷过来的文件会带\rlength($0)会多算 1 个字符。一条序列多算 1bp 影响不大但如果你是几百万条总长度就偏了几 MbN50 也可能正好卡在临界点上错位。用file assembly.fasta或grep -c $\r检查。存在空行中间夹空行不会改变长度空行 length 为 0但会让某些解析器误判序列结束用 awk 累加则不受影响。header 里带空格和注释像contig_1 length1234 cov56.7这种用某些工具解析时会截断。自己写脚本只按^判断不受影响。文件被截断拷贝中断、磁盘满了FASTA 末尾少一段。算出来的总长度会比预期小很多且没有报错。用tail -2看最后几行是不是以序列内容结尾再和上游日志里的总长度对一下。注意\r导致的多算 1bp 在小规模数据上察觉不到但在过滤阈值恰好设成 500bp 时可能让一批实际 499bp 的序列混进统计这类偏差极难事后发现。2.3 一行命令先摸清数据底细正式跑指标之前我习惯先花十秒做个体检确认文件没坏、规模符合预期# 序列条数、文件大小、前两条 header echo 条数: $(grep -c ^ assembly.fasta) ls -lh assembly.fasta grep ^ assembly.fasta | head -2 # 如果是 gz注意用 zcat / pigz -dc pigz -dc assembly.fasta.gz | grep -c ^条数这个数很关键。如果组装结果只有几百条 contig那多半是用了极严格的参数或者组装器直接把大部分数据丢掉了如果有几百万条说明组装很碎后面统计和排序的耗时要提前评估。同时把条数和总长度、N50 一起记录下来形成自己项目的基线下次同类样本一对比就知道有没有异常。还有个实用技巧把最长的几条 header 打出来看看如果最长的那条远大于其他序列往往对应某个高丰度物种这种情况 N50 会明显被拉高不能单纯当作组装好的证据。# 最长的 5 条序列的 headerseqkit 已安装的情况下 seqkit sort -l -r assembly.fasta | grep ^ | head -53. 纯 awk 方案零依赖算出全套指标3.1 算法拆成四步就清楚了不管用什么语言实现N50 计算永远是这四步第一从 FASTA 中提取每条序列的长度得到一个纯数字列表第二把这个列表按降序排列第三求总长度并算出 50%、90% 的阈值第四从头累加第一次超过阈值时记录当前长度和已用条数。拆成四步的好处是可以分段调试。第一步错了后面全错第二步排序方向反了用了升序结果会离谱到不行第三步阈值算成总长度的一半再除以条数这种就彻底跑偏了。我见过最典型的一个 bug 是排序用了sort -n而不是sort -nr跑出来的 N50 只有几百 bp当事人以为组装彻底失败了其实是脚本升序排列累加永远从最短的开始。选择纯 awk sort 这套组合的原因很实际任何一台 Linux 机器上都有不需要装 conda、不需要联网、不需要担心版本冲突流式处理内存占用只跟序列条数有关不跟文件大小有关而且逻辑透明出了错一眼能看出在哪一步。相比之下把整个 FASTA 读进 Python 字典的写法遇到 30GB 的文件就直接卡死了。3.2 第一段把长度从 FASTA 里抠出来提取长度这一步的关键是处理多行序列。判断逻辑是遇到开头的行说明上一条序列结束了把累加器里的值输出并清零其他行就把当前行长度加到累加器里文件结束时别忘了输出最后一条。awk /^/ { if (n 0) print n; n 0; next } # 新 header输出上一条长度 { n length($0) } # 普通序列行累加长度 END { if (n 0) print n } # 输出最后一条 assembly.fasta lengths.txt wc -l lengths.txt # 条数应与 grep -c ^ 一致这里有个自检点必须做lengths.txt的行数要和grep -c ^ assembly.fasta的结果完全相等。如果少了说明有序列被漏掉最常见的原因就是 header 行和序列行之间夹了空行导致逻辑错乱或者文件末尾没有换行符导致最后一条没被 END 块捕获。如果多了通常是文件里混入了非 FASTA 内容。顺手把总长度也算出来对一下awk {s$1} END{printf 总长度: %d bp (%.2f Mb)\n, s, s/1e6} lengths.txt这个数应该和上游组装日志里的总长度基本一致差值一般在过滤掉的短序列范围内。对不上就停下来查别急着算 N50。3.3 第二段降序排序与累加定位排序用LC_ALLC sort -nr加LC_ALLC是为了避开 locale 带来的排序开销在几百万行的数据上能快不少。-n按数值排序-r降序。这里有个细节如果长度都是整数sort -nr完全够用不要用默认的字典序排序否则 9 会排在 100 后面。定位这一步用 awk 一次遍历完成把所有 Nx 值和 Lx 值一起算出来LC_ALLC sort -nr lengths.txt | awk { len[NR] $1 total $1 if ($1 max) max $1 } END { N NR if (N 0) { print 空文件; exit } for (i 1; i N; i) { cum len[i] if (!d50 cum total * 0.5) { n50 len[i]; l50 i; d50 1 } if (!d90 cum total * 0.9) { n90 len[i]; l90 i; d90 1 } } printf 序列条数\t%d\n, N printf 总长度\t%d\n, total printf 最长序列\t%d\n, max printf 平均长度\t%.0f\n, total / N printf N50\t%d\n, n50 printf L50\t%d\n, l50 printf N90\t%d\n, n90 printf L90\t%d\n, l90 }这里特意用了d50、d90两个标志位而不是用if (!n50)判断。原因是万一某条序列长度恰好被算成 0比如空行被误判成序列!n50会一直成立n50 会被反复覆盖成 0最后输出一个荒谬的结果。用独立的布尔标志位就不会有这个问题这是我从一次数据里混了空序列的事故中总结出来的。还有一点total * 0.5用的是浮点运算awk 默认双精度几十 Gb 的总长度完全在精度范围内不用担心。但如果你在别的语言里用 32 位整型存总长度超过 2.1Gb 就会溢出变成负数——这也是个经典坑。3.4 完整脚本一次跑完全部指标把上面两段拼起来加上 gz 支持和参数检查就是一个能直接用的脚本。存成nxx.shchmod x之后可以直接对任意 FASTA 调用#!/bin/bash # nxx.sh —— 计算 FASTA 的 N50/N90/L50/L90 # 用法: bash nxx.sh assembly.fasta[.gz] set -euo pipefail fa$1 [[ -f $fa ]] || { echo 文件不存在: $fa 2; exit 1; } case $fa in *.gz) readerpigz -dc ;; *) readercat ;; esac $reader $fa \ | awk /^/ { if (n 0) print n; n 0; next } { n length($0) } END { if (n 0) print n } \ | LC_ALLC sort -nr \ | awk -v sample$(basename $fa) { len[NR] $1; total $1 if ($1 max) max $1 } END { N NR if (N 0) { printf %s\t0\t0\t0\t0\t0\t0\t0\t0\n, sample; exit } for (i 1; i N; i) { cum len[i] if (!d50 cum total * 0.5) { n50 len[i]; l50 i; d50 1 } if (!d90 cum total * 0.9) { n90 len[i]; l90 i; d90 1 } } printf %s\t%d\t%d\t%d\t%d\t%.0f\t%d\t%d\t%d\n, sample, N, total, max, n50, total/N, l50, n90, l90 }跑一下看看输出printf 样本\t条数\t总长\t最长\tN50\t平均\tL50\tN90\tL90\n bash nxx.sh assembly.fasta在 20GB 左右的宏基因组 FASTA 上这套流程在普通 SSD 机器上的耗时大概是awk 提取长度 1~3 分钟瓶颈在磁盘 IO排序几秒到十几秒累加定位不到 1 秒。内存占用取决于序列条数几十万条的量级也就几十 MB完全不用担心。如果只想快速要两个数可以压缩成一行grep -v ^ assembly.fasta | awk {slength($0)} END{print s} # 仅适用于单行序列但我要提醒这行命令只对每条序列占一行的 FASTA 成立一旦遇到多行序列就全错。我见过有人拿它统计组装结果跑出来的总长度比实际小了三分之一就是因为组装器输出的 FASTA 每行 60 字符换行。别图省事。4. 现成工具交叉验证结果不一致时如何定位4.1 seqkit stats 的 -a 参数seqkit 是这几年用得最顺手的 FASTA/Q 工具集stats子命令默认只给条数、总长、最短、平均、最长要拿到 N50 得加-a全量统计seqkit stats -a -T assembly.fasta-T是纯 tab 分隔输出方便直接贴进 Excel 或喂给下游脚本。-a之后会多出 Q1/Q2/Q3、sum_gap、N50、GC 等列。需要注意的是seqkit 只给 N50不给 N90想凑齐还是得自己算。另外它对 gz 文件直接支持不用手动解压这点在大文件上很省事。一个容易被忽略的差异seqkit 统计的是文件里实际存在的序列包含所有长度不做任何过滤。如果你要的是过滤掉 500bp 以下之后的 N50得先seqkit seq -m 500过滤一遍再统计不能拿原始文件的数去跟别人过滤过的数比。4.2 assembly-stats 与同类工具的批量模式assembly-stats这类工具的好处是能一次接收多个文件输出对比视图assembly-stats -t sample1.fa sample2.fa sample3.fa-t是 tabular 输出一行一个样本字段包括总长、条数、N50、N90 等关键字。不同版本的字段顺序可能略有差异第一次用的时候先不加-t跑一个文件看清它打印的字段名再决定怎么取列别上来就cut -f6。这类工具的共性问题是它们的默认过滤阈值各不相同有的默认不过滤有的默认丢掉 500bp 以下。所以拿两个工具的结果对比之前先确认两者的输入是不是同一个文件、阈值是不是同一档。4.3 QUAST 报告里该看哪一行需要出正式评估报告的时候QUAST 是常用选择quast.py -o quast_out -t 8 --min-contig 0 assembly.fasta跑完看quast_out/report.txt。这里有个必须提醒的细节QUAST 有--min-contig参数不同版本默认值不一样常见默认是 500而这个参数会影响统计口径。报告表格里同时有多列比如contigs ( 0 bp)、contigs ( 500 bp)、contigs ( 1000 bp)而下方的N50、L50行默认对应的是经过阈值过滤的结果。也就是说如果你直接读N50那一行却没有注意列头的阈值说明很容易把过滤后的 N50当成全量 N50报出去。我的做法是跑 QUAST 时显式写上--min-contig 0让自己的报告口径和手算脚本完全一致如果要报过滤后的值就单独再跑一次--min-contig 500两份报告分开放文件名里带阈值杜绝混淆。4.4 三个工具结果对不上时的排查顺序真遇到数字打架别急着怀疑工具按这个顺序查基本能定位排查项检查方法常见结论输入文件是否同一份md5sum对比一个是过滤后的一个是原始的是否含 gz 且被重复解压看文件大小与条数条数翻倍总长虚高过滤阈值是否一致查工具参数与文档一方丢了 500bp 以下contig / scaffold 混用看 header 与 N 串报的是不同层级行尾 CRLFfile或grep -c $\r每条多算 1bp条数多时偏差明显是否合并了多样本grep -c ^对比单样本预期混合文件指标无意义按这个表走一遍九成以上的结果不一致都能解释清楚。剩下那一成多半是某一方用的不是标准的 Nx 定义比如把阈值序列本身的长度和累计长度搞混了这时候让对方报一下算法细节就能对上。5. 批量处理把单样本脚本改成流水线5.1 循环处理目录下所有样本一个项目动辄几十上百个样本逐个手跑不现实。用一层 for 循环套住脚本把结果追加到同一个文件里outnxx_summary.tsv printf 样本\t条数\t总长\t最长\tN50\t平均\tL50\tN90\tL90\n $out for fa in ./assemblies/*.fasta; do bash nxx.sh $fa $out done column -t $out | less -S有几个工程上的细节值得说。第一循环里不要用cd用相对路径或绝对路径避免中途路径错乱。第二样本多的时候建议加日志把每个样本的开始和结束时间打出来不然跑到一半卡住你都不知道卡在哪个文件上for fa in ./assemblies/*.fasta; do echo [$(date %F %T)] 处理 $fa run.log bash nxx.sh $fa $out 2 run.log done第三如果机器核数够可以用xargs -P并行起来。注意并行时不能直接往同一个文件追加输出会交错。正确做法是每个样本写一个临时文件最后统一 catls ./assemblies/*.fasta \ | xargs -P 4 -I{} bash -c bash nxx.sh {} {}.nxx.tmp cat ./assemblies/*.nxx.tmp nxx_summary.tsv rm -f ./assemblies/*.nxx.tmp-P 4是并发数别开太大。这类任务的瓶颈在磁盘读取开 4~8 个并发基本就能把 IO 打满再往上加只会让磁盘排队总耗时反而变长。我实测过在机械盘上开 16 并发比开 4 并发慢了将近一倍。5.2 汇总成 TSV 之后做初步判断拿到汇总表第一件事不是急着下结论而是横向扫一遍有没有异常值。我的习惯是先用 awk 把几个关键列拉出来排序看极值# 按 N50 升序看哪些样本最碎 awk NR1{print $5\t$1} nxx_summary.tsv | sort -n | head -10 # 按总长度排序看有没有明显偏小的 awk NR1{print $3\t$1} nxx_summary.tsv | sort -n | head -10总长度明显偏小的样本通常是测序量不足或者质控过严导致大部分数据被丢N50 特别低的样本可能是组装参数过于保守或者样本本身的物种复杂度极高、重复区多。这两种情况处理方式不同不能一概而论。还可以算一个衍生比值看一眼长度分布形态awk NR1{printf %s\tN90/N50%.2f\n, $1, $8/$5} nxx_summary.tsv | sort -t -k2 -n | head比值接近 1 的说明长度分布集中比值很小的比如 0.1 以下说明结果是少数超长 海量短碎片的组合这种样本在做分箱时要格外注意低丰度物种的信息可能大量丢在短片段里。5.3 结果解读的几条经验线跑了这么多项目我对宏基因组组装的 N50/N90 有几条自己的经验判断线仅作参考不能当硬标准contig N50 在 5kb 以下、同时总长度远低于预期说明组装质量确实有问题优先查数据质量和组装参数而不是急着往下游走。contig N50 在 5~20kb、N90 在 1kb 以上属于大多数土壤、肠道样本的正常范围可以正常进入分箱流程。N50 高但总长度明显偏小警惕过度过滤。先确认是不是丢掉了大量短 contig这会影响低丰度物种的回收率。同批次样本之间 N50 差异超过 3 倍先查是不是建库或测序环节的批次效应再看生物学原因。6. 常见问题与排查技巧实录6.1 排序与累加环节的四个经典错误这部分是我在帮别人看脚本时见得最多的几乎每个都能让结果偏离正确值第一种排序方向反了。sort -n是升序累加从最短的开始N50 会算出个几百 bp 的荒谬值。解决办法很简单看到 N50 小得离谱第一反应就是检查有没有-r。第二种用累计长度而不是当前序列长度。有些人写if (cum total*0.5) print cum输出的是累计长度比正确答案大。要输出的是len[i]。第三种对每条序列只取第一行长度。这在多行 FASTA 上会严重低估且各序列被低估的程度还不一样取决于行长N50 完全失真。判断方法算出来的总长度和上游日志对不上基本就是这个原因。第四种阈值算错。有人用total / 2 / N当阈值或者把 90% 写成 0.9 * N条数的 90%这完全是另一回事。记住 Nx 永远是对总长度取比例不是对条数。还有一个隐蔽的用 Python 的statistics.median之类去估算虽然思路不同但有人真这么干。N50 不是中位数两者在长度分布偏态时差别很大别混。6.2 大文件处理内存、IO 与耗时宏基因组 FASTA 动辄几十 GB处理时要注意几个点。纯 awk sort 的方案是流式的内存占用只和条数有关不随文件变大而线性增长这是它的最大优势。但sort在数据量超过内存时会自动落盘做外部排序速度会掉下来。如果你的序列条数上千万可以给 sort 指定缓冲区LC_ALLC sort -nr -S 2G --parallel4 lengths.txt-S 2G是给排序分配的内存上限--parallel是多线程。这两个参数在条数上千万时能明显提速。另一个容易忽略的是 IO。用pigz -dc代替gzip -dc解压多核并行解压能快好几倍尤其在 gz 压缩率很高的情况下。如果同一个文件要跑多次分析先解压到本地磁盘再跑比每次解压划算。还有一点不要在 NFS 挂载的网络盘上直接跑统计。我踩过一次20GB 的文件读了两小时还没读完换到本地盘只要两分钟。网络盘的顺序读带宽和延迟都不适合这种流式大文件扫描。6.3 常见问题速查表现象可能原因处理办法N50 只有几百 bp排序用了升序改用sort -nr加LC_ALLC总长度比日志小很多多行序列只取了首行用累加式 awk 提取长度总长度比预期大一点行尾 CRLF 多算 1bpdos2unix或sed -i s/\r$//条数与 grep 统计不符空行、末尾无换行用tail -c 1检查文件结尾和别人的数差一截过滤阈值不同明确500bp还是0bp报的是 scaffold N50文件层级搞混看 header 与连续 N 串跑得极慢网络盘 IO 或未并行解压拷到本地盘、用pigzgz 文件统计出错未解压直接读用pigz -dc或工具自带支持提示任何对外报告的 N50 数字建议在表格脚注里写清三件事——层级contig/scaffold、过滤阈值如 ≥500bp、总长度。这三项写清楚别人才能复现你的数。最后分享一个我自己一直在用的小习惯每个项目的组装目录里放一个stats.md把每次统计的命令、参数、结果、日期都记一行。看起来是额外工作但当三个月后有人问你当时那个样本的 N50 是多少、用的什么阈值你能十秒内翻出来而不是重新跑一遍再猜当时怎么设的参数。做生信这行可复现性比聪明更重要一个随手记的习惯能省掉大量返工。另外如果你手上的组装结果还带覆盖度信息可以在长度统计的基础上一并把每条序列的平均深度也汇总出来配合 N50 一起看能更快判断这份组装是整体偏碎还是只有低丰度物种碎这个思路往下延伸还能直接接到分箱前的质控环节。