3. FASTA 格式详解：从文件结构到命令行统计

FASTA 格式是我进入生信领域接触的第一个文件格式。说实话，当年我在实验室用 perl 写统计脚本的时候，遇到过一个坑：一个细菌基因组有 300 多条 scaffold，我直接把整个文件当作一条序列去算 GC 含量……结果出来的值怎么都不对。后来才搞清楚，FASTA 里每条序列以 > 开头，序列行可以跨多行，而且绝不能统计空白行里的东西。所以这一章，咱们就把 FASTA 的里里外外掰开揉碎讲清楚。

3.1 文件结构：> 号、标识符和序列行

一个标准的 FASTA 文件长这样：

>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGT
GTCGATAGCAGCTTCTGAACTGGTTACCTGCCGTGAGTAAATTAAAATTTTATTGAC
TTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGACAGATAAAAATTA

核心就两个部分：

• 标识符行（defline）：以 > 开头，后面紧跟序列 ID，空格后可加描述信息。整行一般不超过 80 个字符（FASTA 规范建议）。
• 序列行：由 A、T、C、G、N（模糊碱基）或小写字母组成。可以连续多行，每行长度通常为 60–80 个字符。序列结束后，要么是文件末尾，要么是下一个 >。

有些文件会在两条序列之间加空行，这其实不规范。我建议你用 awk 或 sed 预处理去掉空行。

3.2 质量检查：常见错误与避坑指南

我接手过不少公共数据库下载的细菌基因组 FASTA，里面常见的问题有：

我曾经帮同事排查一个注释项目，序列长度明明都是 2M，BLAST 结果却显示有 200 万条短片段——原因是 FASTA 里插入了大量空格和数字，导致 seqkit 把每行分开统计。后来我用下面这条命令清理，问题立刻解决：

sed -e 's/[^ACGTNacgtn>]//g' -e '/^$/d' input.fasta > clean.fasta

3.3 命令行快速统计序列长度和 GC 含量

统计这些信息，说白了就是为 BLAST 注释前的数据质控做准备。我一般用 awk 或 bioawk，一行命令搞定。

3.3.1 用 awk 手工统计

假设细菌基因组是单条染色体，或者你想要每条 scaffold 的详情：

awk '/^>/ {if (seqlen) print seqlen; seqlen=0; next} 
     {seqlen += length($0)} 
     END {print seqlen}' sample.fasta

如果你还想要 GC 含量：

awk '/^>/ {if (seqlen) print id, seqlen, gc; id=$0; seqlen=0; gc=0; next} 
     {seqlen += length($0); gc += gsub(/[GCgc]/,"",$0)} 
     END {print id, seqlen, gc}' sample.fasta

注意：gsub 会修改 $0，但这里我们用它的返回值（替换次数）。我建议在正式分析前先在小样本上测试一下。

3.3.2 使用 bioawk 两句话

# 统计每条序列长度和 GC 比例
bioawk -c fastx '{print name, length(seq), gc(seq)}' sample.fasta

bioawk 内置了 gc() 函数，返回的是 0-1 的小数值。如果想看 GC 百分比，乘 100 就行：

bioawk -c fastx '{printf "%s\t%d\t%.1f\n", name, length(seq), gc(seq)*100}' sample.fasta

3.3.3 统计整个文件的总长度和平均 GC

# 总碱基数
bioawk -c fastx '{tlength += length(seq)} END {print tlength}' sample.fasta

# 总 GC 数
bioawk -c fastx '{tgc += gc(seq)*length(seq)} END {print tgc}' sample.fasta

3.4 示例：手绘 FASTA 结构及质控流程

下面这张图把 FASTA 的构成和常见校验点串在一起了。你用它对照自己的文件，能快速定位问题。

你可以看到，从文件到统计结果之间，必须经过空行检测、ID 重复检查、非法字符过滤，最后才能用命令行安全地算出每一条序列的长度和 GC 含量。这一步如果跳过，后面的 BLAST 注释很有可能被错误的数据坑死。