企业官网建设流程全解析

宏基因组组装跑完之后，最常被问到的一句话就是：“这批样本的 N50 是多少？”组装的 FASTA 文件动辄几十 GB，用现成工具跑一遍要等，换了工具数字还对不上，第三方又突然要一个 N90，这些情况在实际项目里几乎每周都会遇到。这篇内容就聊一件事：在 Linux 命令行里，把宏组装数据的 N50、N90 连同 L50、auN 这些衍生指标一次算准、算快，并且把计算口径固定下来，避免不同人算出不同结果。涉及的代码可以直接抄走用，awk 和 Python 两个版本都给，读完就能在自己的数据上跑出结果。

1. 宏组装报告里的 N50/N90：这两个数字为什么总要自己复核

1.1 N50 在累积长度曲线上的几何意义

N50 的定义其实不复杂，但很多人背下来的是“长度加权中位数”这句口诀，真到要解释给别人听的时候卡壳。更准确的说法是：把组装结果里所有 contig 的长度从大到小排成一列，依次累加，当累加的碱基数刚好达到全部碱基总数一半的时候，当前这条 contig 的长度就是 N50。N90 同理，把 50% 换成 90%。

拿一组小数据手算一遍更直观。假设有 5 条 contig，长度分别是 1000、800、500、300、200，总长 2800 bp。降序累积依次是 1000、1800、2300、2600、2800。总长的一半是 1400，累积到第二条 800 的序列时超过 1400，所以 N50 等于 800。总长的 90% 是 2520，累积到第四条 300 的序列时才超过 2520，所以 N90 等于 300。

这个定义背后其实是一条“累积长度曲线”：横轴是按长度降序排列的序列编号，纵轴是累积碱基数。N50 就是在纵轴走到一半时，横轴对应那条序列的高度。所以 N50 天然偏向长 contig——只要有几条超长序列撑住曲线，N50 就会被拉得很高，哪怕整个组装的碎片化程度很严重。这也是为什么宏组装报告里 N50 一定要和 L50、总长、contig 条数一起看，单独一个 N50 几乎说明不了问题。

1.2 现成工具为什么还不够用

有人会问，现成的统计工具不是一条命令就出 N50 吗，何必自己写。这话对一半。现成工具的问题集中在三个地方。

第一是口径不透明。不同工具对“是否过滤短 contig”“N90 用首次达到还是插值”“是否把 N 碱基计入总长”的处理并不一致，默认参数下算出来的数字可能差几个百分点。在内部看趋势无所谓，一旦要写进交付报告或者和合作方对标，就必须知道自己报的到底是哪一套口径。

第二是输出格式不受控。项目里经常需要把几十个样本的统计结果拼成一张 TSV 表，直接喂给下游的绘图脚本或数据库。现成工具的默认输出是给人看的文本，要转成规整的表格还得再写一层解析，反而更麻烦。

第三是批量场景下的额外开销。每个样本都启动一次完整的工具进程，几百个样本跑下来，进程启动和依赖加载的时间可能比真正统计的时间还长。而统计长度这件事，本质上只是扫一遍 FASTA 文件、累加每行的字符数，用 awk 或一个几十行的 Python 脚本就够了，没有额外的外部依赖，扔到任何一台 Linux 机器上都能跑。

1.3 宏组装和单基因组组装在统计口径上的差别

宏基因组组装和单基因组组装在统计上有一个很实际的差别：宏组装的 contig 数量多出两三个数量级。一个细菌基因组组装可能只有几百条 contig，而一份土壤或肠道样本的宏组装动辄几十万到上百万条，其中大量是几百 bp 的短碎片。

这个差别带来两个直接后果。一是短 contig 对总长的贡献不可忽略，是否设置最小长度过滤会明显改变分母，进而改变 N50 和 N90。二是内存和 IO 的压力变大，把全部长度读进内存再做排序，在百万级序列上还能接受，在千万级就要开始算计内存了。所以宏组装的统计脚本在设计上要额外考虑“过滤阈值可配置”和“内存占用可控”这两点，这也是后面代码里会专门处理的地方。

2. 从 FASTA 到长度数组：Linux 下最容易被忽略的解析细节

2.1 只靠 awk 打通 gz 文件到长度列表

统计 N50 的第一步是把 FASTA 变成一列长度数字。这一步看起来最没技术含量，实际上出问题最多。一个稳妥的写法是：用解压流把文件送进 awk，awk 按 header 切分记录、累加序列行长度，输出每条的碱基数，最后交给 sort 排序。

zcat sample.assembly.fa.gz \ | awk '/^>/{ if (n) printf "%d\n", len n++; len = 0; next } { sub(/\r$/, "") gsub(/[ \t]/, "") len += length($0) } END{ if (n) printf "%d\n", len }' \ > sample.len

这段 awk 的逻辑值得逐句拆开看。遇到以>开头的行，说明上一条记录结束了，如果之前已经出现过记录（n大于 0），就把上一条累积的长度打出来，然后把长度计数器清零；遇到非 header 行，先去掉行尾的\r，再去掉可能存在的空格或制表符，把剩余字符数累加到len上。最后的END块负责补打最后一条记录，因为文件结尾不会有下一个 header 来触发打印。

这个实现有个隐含前提，就是序列可能跨多行。很多组装软件默认会把每条 contig 的序列折成 60 或 80 个字符一行，如果脚本假设“一条序列只占一行”，长度就会被算成行数，结果全错。用累积的方式处理就完全不受折行影响，这也是为什么要用len += length($0)而不是直接取某一行的长度。

2.2 多行序列、CRLF、空行与注释行

上面那段代码里特意处理了\r和空格，原因很实际。FASTA 文件从 Windows 环境转过来、或者经过某些文本处理工具之后，行尾会带上回车符，length($0)会把\r也数进去，每条序列多算一个碱基。单条看不出来，几十万条累计下来就是几十万 bp 的偏差，总长和 N50 都会轻微偏移。这种偏差不会让程序报错，只会让数字悄悄不对，最难查。

空行的问题类似。有些文件的记录之间存在空行，空行的长度是 0，正常情况下不影响累加结果，但如果文件开头有空行，第一个 header 之前就会有一段“无主”的序列内容，被算进第一条记录里。稳妥的做法是在解析时判断当前是否已经进入某条记录，没进入就跳过。Python 版本里用in_record这个标志位就是为了这个。

还有一种容易踩的情况是 header 行本身带描述文字，比如>contig_1 length=12345 cov=8.7。统计长度时完全不需要解析这些字段，只要认准行首是>就够了。但有人会图省事，直接从 header 里抓length=后面的数字当作实际长度，这在文件完整时结果一样，可一旦文件被截断、或者序列被某个工具重新折行改动过，header 里的数字就和实际序列对不上了。我的习惯是以实际序列字符数为准，header 里的长度只用来做交叉验证。

2.3 过滤阈值必须在统计前定死

宏组装数据里短 contig 的数量占比极高，不设阈值直接统计，N50 会被大量短序列拉低，而且不同批次之间因为短序列数量波动，N50 的可比性很差。所以业界常见做法是先设定一个最小长度，比如 500 bp 或 1000 bp，只统计达到阈值的序列。

这里的关键不是选哪个阈值，而是同一批数据、同一个项目里必须用同一个阈值，并且在报告里写清楚。过滤在哪个环节做也有讲究：可以在 FASTA 层面先用工具筛一遍再统计，也可以把阈值作为参数传给统计脚本，在累加长度的阶段直接跳过。后者的好处是不产生中间文件，几十 GB 的数据不用多占一份磁盘。

用 awk 实现过滤很简单，把打印那一行改成条件判断即可：if (n && len >= MIN_LEN) printf "%d\n", len，其中MIN_LEN通过awk -v MIN_LEN=500传入。Python 版本里就是一个if l >= args.min_len: yield l。别看改动只有一行，它能省掉一次全量数据写盘，在磁盘吃紧的机器上很实用。

3. 把累积和算明白：N50/N90 的判定规则与两套口径

3.1 降序累积的判定过程

拿到长度列表之后，排序和累加是核心。排序必须按数值降序，不能按字符串排序，这是最经典的坑之一。假设长度列表里有 1000、800、200，直接sort（默认字典序）得到的是 1000、200、800，因为字符1小于2小于8。用错排序方式，N50 会在某些数据上完全跑偏，而且因为小数据测试时可能碰巧没暴露，很容易带到生产环境里。

正确的做法是sort -nr，其中-n表示按数值比较，-r表示降序。排序完成后，用 awk 再扫一遍做累积：

awk -v p50=0.5 -v p90=0.9 ' { a[NR] = $1; total += $1 } END{ t50 = total * p50; t90 = total * p90; acc = 0 for (i = 1; i <= NR; i++) { acc += a[i] if (!h50 && acc >= t50) { n50 = a[i]; l50 = i; h50 = 1 } if (!h90 && acc >= t90) { n90 = a[i]; l90 = i; h90 = 1 } if (h50 && h90) break } printf "%d\t%d\t%d\t%d\t%d\t%d\n", NR, total, a[1], n50, n90, l50 }' sample.len

这里用数组把长度存下来，是因为累积到 N50 和 N90 需要的目标值依赖总长，而总长要先扫完一遍才知道。百万级序列的数组在 awk 里大概占几十 MB，可以接受；千万级就要考虑换成两遍扫描文件的方式，避免数组膨胀，具体在后面的性能章节展开。

3.2 “首次达到”与“线性插值”的差异有多大

这是最容易引起争议的一点。同样是 N50，不同工具可能给出两个不同的数。

“首次达到”口径就是我们上面手算的过程：累积碱基数第一次大于等于目标值时，取当前那条 contig 的长度。这是绝大多数工具和文献默认采用的定义，简单、可复现。

“线性插值”口径则假设在单条 contig 内部，碱基是随长度线性分布的。当累积到第 i 条时超过了目标，不直接取 L_i，而是从 L_i 里减去超出的那部分，得到L_i - (acc_i - target)。听起来更“精确”，实际上把 contig 当成了一个可以任意切分的连续体，而 contig 是离散对象，切一半没有生物学意义。

差异有多大？在真实宏组装数据上，因为最长的那几条 contig 占总长的比例通常不高，两种口径的差距一般在零点几个百分点，报告里几乎看不出区别。但在极端情况下差距可以很大：如果整个组装只有一条 1000 bp 的序列，总长 1000，目标 500，首次达到口径给出 N50 = 1000，插值口径给出 500，差了一倍。所以这件事的重要性不在于数值本身，而在于必须在报告里写清楚用的哪一套，否则跨项目、跨工具的数字没法比较。

我个人的做法是默认用“首次达到”，同时在脚本里保留一个--interp开关，遇到合作方明确要求插值口径时直接切换，输出文件里额外加一列口径标记，避免几个月后自己也记不清当时报的是哪个数。

3.3 L50、auN、NGx 这些衍生指标顺手一起算

既然已经把长度数组排好序了，多算几个指标几乎不增加成本，但信息量提升明显。

L50是累积到 50% 时用掉的 contig 条数，也就是上面代码里的l50。它和 N50 是一对：N50 高但 L50 也高，说明靠少数几条超长序列撑起来的，组装其实还是碎的。

auN的定义是sum(L_i^2) / sum(L_i)，可以理解为长度加权的平均长度。它比 N50 更平滑，不会因为某一条 contig 的长度跳变而剧烈波动，适合做批次之间的趋势对比。计算时只需要在遍历中累加平方和，注意数值范围，长度是几十万时平方是几百亿，用 64 位整数或双精度都没有问题。

NGx系列是用一个给定基因组大小（而不是实际总长）做分母来算的 Nx。宏组装一般不知道真实基因组大小，所以用不上，但在做单菌 bins 的评估时会遇到，知道这个概念有助于看懂别人报告里的数字。

一个实用的输出表长这样：

4. 三套实现：awk 流式、Python 脚本、现成工具对照

4.1 awk 两遍扫描的极简版本

如果只是临时看一个数，不想写脚本，awk 加上外部排序就是最短路径。但要注意 awk 数组在超大输入下的内存问题，可以改成两遍扫描同一个长度文件：第一遍只求总长，第二遍边读边累加，读到目标值就退出。长度文件本身很小（每条一行数字，百万条约几 MB），读两遍的成本可以忽略。

# 第一遍：总长 total=$(awk '{s+=$1} END{print s}' sample.sorted.len) # 第二遍：分别求 N50 和 N90 awk -v total="$total" ' { acc += $1 } !h50 && acc >= total*0.5 { n50 = $1; l50 = NR; h50 = 1 } !h90 && acc >= total*0.9 { n90 = $1; l90 = NR; h90 = 1 } h50 && h90 { exit } END{ printf "N50=%d\tL50=%d\tN90=%d\tL90=%d\n", n50, l50, n90, l90 }' sample.sorted.len

这个版本的思路是“能早停就早停”，N90 一般在一半左右的序列位置就能达到，后面的数据完全不用读。在千万级序列的场景下，这个早停能省下相当一部分时间。不过要注意：exit触发后END块仍然会执行，所以END里的打印是安全的，这点 awk 的行为和某些脚本语言不同，容易记混。

4.2 一个能批量出 TSV 的 Python 脚本

日常项目里我更倾向用 Python，因为要批量处理样本、要输出规整表格、要加各种口径开关，用脚本维护起来比一堆 shell 拼接的命令清爽得多。下面这个版本支持 gz 直读、最小长度过滤、插值口径开关和批量输入。

#!/usr/bin/env python3 import argparse import gzip import os import sys def iter_lengths(path, min_len=0): opener = gzip.open if path.endswith(".gz") else open cur, in_rec = 0, False with opener(path, "rt", encoding="ascii", errors="replace") as fh: for line in fh: if line.startswith(">"): if in_rec and cur >= min_len: yield cur cur, in_rec = 0, True else: cur += len(line.strip()) if in_rec and cur >= min_len: yield cur def nx(lengths, pct, interp=False): total = sum(lengths) if total == 0 or not lengths: return 0, 0 target = total * pct / 100.0 acc = 0 for i, L in enumerate(lengths, 1): acc += L if acc >= target: if interp: return int(round(L - (acc - target))), i return L, i return lengths[-1], len(lengths) def stats_for(path, min_len, interp): lens = list(iter_lengths(path, min_len)) if not lens: return None lens.sort(reverse=True) total = sum(lens) n50, l50 = nx(lens, 50, interp) n90, l90 = nx(lens, 90, interp) aun = sum(x * x for x in lens) / total name = os.path.basename(path) for suf in (".gz", ".fa", ".fasta", ".fna", ".fa.gz", ".fasta.gz"): if name.endswith(suf): name = name[: -len(suf)] break return [name, len(lens), total, lens[0], n50, l50, n90, l90, round(aun, 1), min_len] def main(): ap = argparse.ArgumentParser(description="assemble N50/N90 stats") ap.add_argument("fasta", nargs="+") ap.add_argument("-m", "--min-len", type=int, default=0) ap.add_argument("--interp", action="store_true") ap.add_argument("-o", "--out", default="-") args = ap.parse_args() header = ["sample", "contigs", "total_bp", "longest", "N50", "L50", "N90", "L90", "auN", "min_len"] out = sys.stdout if args.out == "-" else open(args.out, "w") out.write("\t".join(header) + "\n") for p in args.fasta: row = stats_for(p, args.min_len, args.interp) if row is None: sys.stderr.write("skip empty: %s\n" % p) continue out.write("\t".join(str(x) for x in row) + "\n") out.close() if __name__ == "__main__": main()

用法很直接：

python3 asm_stats.py -m 1000 asm/*.fa.gz -o all_samples.n50.tsv

几个实现细节值得说明。用line.strip()而不是rstrip("\n")，是为了同时去掉换行和回车；errors="replace"是为了防止极少数非 ASCII 字节直接抛异常中断整个批处理。文件名去掉后缀那段写得有点啰嗦，实际用的时候可以直接用os.path.basename(p).split(".")[0]，前提是样本名里不含点。

内存方面，list(iter_lengths(...))会把所有长度读进内存。Python 的整数对象不算小，每条大约 28 字节加上列表指针 8 字节，一百万条约 36 MB，一千万条约 360 MB。如果机器内存紧张，可以换成array.array("q")，每条只占 8 字节，代价是排序前需要转换一次。

4.3 现成工具的命令与它的默认口径

不想写脚本时，通用序列工具也能直接给出 N50，作为交叉验证很方便：

seqkit stats -a -T sample.assembly.fa.gz

-a输出全部统计项，-T输出制表符分隔格式，里面包含 N50。但它默认不过滤短序列，也不直接给 N90，需要的话可以先用seqkit seq -m 1000过滤，再用seqkit fx2tab -nl导出长度列，交给前面的 awk 处理。

我把三条路线放在一起对比，方便按场景选：

三种方式的结果在口径一致时应该完全相同。如果不同，基本可以确定是过滤阈值或排序方式不一致，这就进入下一节的排查流程。

5. 结果对不上时的完整排查链路

5.1 从“谁的排序不对”开始

两个人算同一个文件的 N50，结果差了百分之几，第一件事是核对分母，也就是总碱基数。如果总长不一致，问题出在解析或过滤环节；如果总长一致但 N50 不同，问题几乎一定在排序或取点规则上。

排查顺序我一般这么走：

1. 打印长度列表的前 10 行和后 10 行，确认是数值降序，不是字典序。看到1000后面跟着200再跟着800，就是排序错了。
2. 打印总长，和上游组装报告里的 total length 对一遍。差得离谱说明有序列没被解析到，比如文件里混入了非 FASTA 格式的段落。
3. 检查是否有空行被当成一条长度 0 的记录。空行混进长度文件后，NR会偏大，L50 跟着错位。
4. 确认过滤阈值是否一致。一边用 0，一边用 1000，N50 差几个百分点很正常。
5. 确认是否把 N 碱基计入长度。绝大多数工具是计入的，但自定义脚本如果只统计 ACGT，总长会偏小。

这个顺序的逻辑是“从最容易验证的开始”，前两步基本能定位八成问题，后面几步是收尾。

5.2 边界输入：空文件、单条序列、全部短于阈值

脚本在真实数据上跑得顺，往往在这些边界上翻车。

空文件或者全是空行的文件，长度列表为空，累加时总长为 0，total * 0.5还是 0，任何一条都会满足条件，如果代码没做判空，可能输出一条长度为 0 的记录，看起来像正常结果。Python 版本里我专门加了if row is None跳过并打到 stderr，就是为了在批处理时能一眼看到哪些样本有问题。

单条序列的情况前面提过，两种口径差异巨大。另外还有一些组装失败的小样本，整个文件只有寥寥几条序列，这时候 N50 的数值意义不大，报告里最好加一句样本量说明，避免被下游误读。

全部短于阈值的情况更隐蔽。比如阈值设成 1000，而这批样本的组装质量很差，最长 contig 只有 800，过滤后列表为空。此时正确行为是明确报“无满足阈值的序列”，而不是输出 0 让下游以为组装成功。这类静默失败在批量跑几百个样本时特别危险。

5.3 一份可以贴在报告里的口径说明

最后分享一个我觉得最有价值的实践：把口径写成一句固定的话，附在每份统计结果的旁边。

我常用的模板是：统计基于长度不小于 1000 bp 的 contig，N50/N90 采用首次达到口径（累积碱基数首次大于等于过滤后总长的 50%/90% 时对应 contig 的长度），碱基计数包含 N，不区分 contig 与 scaffold。

这句话只有几十个字，但它解决了一个长期困扰：半年后回看数据、或者交接给别人的时候，不用再去翻当时的脚本猜参数。配合脚本输出的min_len列和--interp开关状态，整套统计就是自解释的。

6. 几十 GB 输入下的性能取舍

6.1 解压与 IO 是真正的瓶颈

宏组装的 FASTA 压缩后动辄二三十 GB，解压后可能上百 GB。整个统计流程里，耗时的大头是解压和读取，awk 或 Python 的解析速度反而是次要的。实测下来，单线程zcat的输出速度通常在每秒一两百 MB，而 awk 处理纯文本能跑到每秒几百 MB，所以链路是被解压卡住的。

如果机器核数够，可以用并行解压工具替换zcat，把解压速度提上来，让后面的解析成为新瓶颈。这一步的收益在机械硬盘上尤其明显，因为压缩流是顺序读，并行解压能把 CPU 和磁盘都喂饱。

另一个常被忽略的点是磁盘缓存。同一个文件如果刚被上游组装流程写过，很可能还在页缓存里，第二次读会快很多。批量跑样本时，把统计步骤安排在组装刚结束时执行，比隔一天再跑要省不少时间。

6.2 内存占用从 O(N) 降到可控

长度数组的内存随序列条数线性增长。百万级无所谓，千万级就要注意了。三个可行的降内存手段：

一是外部排序加两遍扫描。用sort -nr把长度落到磁盘，再用 awk 两遍读取，内存占用恒定在几 MB。代价是多一次磁盘往返，但对 SSD 来说可以接受。

二是直方图分箱。不保存每条长度，而是维护一个“长度区间到条数”的映射，比如按 100 bp 分箱。统计时从大到小遍历箱体，用箱体代表值乘以条数做累积。内存是常数级，代价是 N50 最多有 100 bp 的舍入误差。对报告用的 N50 来说这个误差通常可接受，但如果要精确到个位就不能用。

三是用定长整数数组。Python 的array("q")或者 numpy 的int64，每条 8 字节，一千万条 80 MB，排序性能也比纯 Python 列表好得多。这是我认为性价比最高的折中方案。

6.3 批量跑样本时的目录组织

几百个样本一起跑，最容易出问题的不是算法，是文件管理。我的习惯是把统计脚本、长度中间文件和最终结果分开放：

work/ ├── asm/ # 原始组装结果 ├── stats/asm_stats.py # 脚本 ├── stats/lengths/ # 可选：长度中间文件 └── stats/result/ # 最终 TSV

批量执行的循环里加上失败隔离，一个样本出错不影响后面的：

mkdir -p stats/result for f in asm/*.fa.gz; do s=$(basename "$f" .fa.gz) python3 stats/asm_stats.py -m 1000 "$f" -o "stats/result/${s}.tsv" \ || echo "FAILED: $s" >> stats/result/failed.log done

跑完之后，把每个样本单独的 TSV 拼成一张总表，再删掉表头重复行即可。这样组织的好处是任何一个样本的结果都可以单独复查，出问题不用重跑全量。我在实际使用中发现，把每个样本的中间长度文件保留下来（-m过滤后通常只有几十 MB），在后期需要临时换口径重算时特别省事，直接对着长度文件重跑累积逻辑，几秒钟就出结果，不用再解压一遍几十 GB 的原始文件。

最后再提一句容易被忽略的：统计结果出来后，花十秒钟扫一眼longest和contigs这两列。如果某个样本的longest明显偏离同批次其他样本，或者contigs数量突然翻倍，通常意味着那个样本的组装参数或者输入数据有问题，比 N50 本身更能说明批次质量。这个习惯帮我提前发现过好几次上游数据被截断的情况。