企业官网建设流程全解析

在 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，不够），+80=180（超过 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）基本能把一个组装结果的"骨架"描述清楚。

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 # 粗看是否存在长串 N（scaffold 的典型特征） awk '/^>/{next}{n=gsub(/N/,"N"); if(n>50) c++} END{print "含长N串的行数:", c+0}' 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：从某些系统拷过来的文件会带\r，length($0)会多算 1 个字符。一条序列多算 1bp 影响不大，但如果你是几百万条，总长度就偏了几 Mb，N50 也可能正好卡在临界点上错位。用file assembly.fasta或grep -c $'\r'检查。
• 存在空行：中间夹空行不会改变长度（空行 length 为 0），但会让某些解析器误判序列结束，用 awk 累加则不受影响。
• header 里带空格和注释：像>contig_1 length=1234 cov=56.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 条序列的 header（seqkit 已安装的情况下） seqkit sort -l -r assembly.fasta | grep '^>' | head -5

3. 纯 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_ALL=C sort -nr，加LC_ALL=C是为了避开 locale 带来的排序开销，在几百万行的数据上能快不少。-n按数值排序，-r降序。这里有个细节：如果长度都是整数，sort -nr完全够用；不要用默认的字典序排序，否则 "9" 会排在 "100" 后面。

定位这一步用 awk 一次遍历完成，把所有 Nx 值和 Lx 值一起算出来：

LC_ALL=C 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.sh，chmod +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) reader="pigz -dc" ;; *) reader="cat" ;; esac $reader "$fa" \ | awk '/^>/ { if (n > 0) print n; n = 0; next } { n += length($0) } END { if (n > 0) print n }' \ | LC_ALL=C 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 '{s+=length($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 三个工具结果对不上时的排查顺序

真遇到数字打架，别急着怀疑工具，按这个顺序查基本能定位：

按这个表走一遍，九成以上的"结果不一致"都能解释清楚。剩下那一成，多半是某一方用的不是标准的 Nx 定义（比如把阈值序列本身的长度和累计长度搞混了），这时候让对方报一下算法细节就能对上。

5. 批量处理：把单样本脚本改成流水线

5.1 循环处理目录下所有样本

一个项目动辄几十上百个样本，逐个手跑不现实。用一层 for 循环套住脚本，把结果追加到同一个文件里：

out="nxx_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并行起来。注意并行时不能直接往同一个文件追加，输出会交错。正确做法是每个样本写一个临时文件，最后统一 cat：

ls ./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 'NR>1{print $5"\t"$1}' nxx_summary.tsv | sort -n | head -10 # 按总长度排序，看有没有明显偏小的 awk 'NR>1{print $3"\t"$1}' nxx_summary.tsv | sort -n | head -10

总长度明显偏小的样本，通常是测序量不足或者质控过严导致大部分数据被丢；N50 特别低的样本，可能是组装参数过于保守，或者样本本身的物种复杂度极高、重复区多。这两种情况处理方式不同，不能一概而论。

还可以算一个衍生比值，看一眼长度分布形态：

awk 'NR>1{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_ALL=C sort -nr -S 2G --parallel=4 lengths.txt

-S 2G是给排序分配的内存上限，--parallel是多线程。这两个参数在条数上千万时能明显提速。

另一个容易忽略的是 IO。用pigz -dc代替gzip -dc解压，多核并行解压能快好几倍，尤其在 gz 压缩率很高的情况下。如果同一个文件要跑多次分析，先解压到本地磁盘再跑，比每次解压划算。

还有一点：不要在 NFS 挂载的网络盘上直接跑统计。我踩过一次，20GB 的文件读了两小时还没读完，换到本地盘只要两分钟。网络盘的顺序读带宽和延迟都不适合这种流式大文件扫描。

6.3 常见问题速查表

提示：任何对外报告的 N50 数字，建议在表格脚注里写清三件事——层级（contig/scaffold）、过滤阈值（如 ≥500bp）、总长度。这三项写清楚，别人才能复现你的数。

最后分享一个我自己一直在用的小习惯：每个项目的组装目录里放一个stats.md，把每次统计的命令、参数、结果、日期都记一行。看起来是额外工作，但当三个月后有人问你"当时那个样本的 N50 是多少、用的什么阈值"，你能十秒内翻出来，而不是重新跑一遍再猜当时怎么设的参数。做生信这行，可复现性比聪明更重要，一个随手记的习惯能省掉大量返工。另外，如果你手上的组装结果还带覆盖度信息，可以在长度统计的基础上一并把每条序列的平均深度也汇总出来，配合 N50 一起看，能更快判断这份组装是"整体偏碎"还是"只有低丰度物种碎"，这个思路往下延伸还能直接接到分箱前的质控环节。