手写 FASTA 解析器踩了 9 个坑:从 0 条记录到 3 条的完整调试实录
本文是《手写 FASTA 解析器》系列第 1 篇
我是生信方向研一新手,Python 零基础起步,研究方向是小鼠 VNTR(可变数目串联重复)变异分析。
为了搞懂序列文件是怎么被读进程序的,我从零手写了一个 FASTA 解析器 ——改了 9 版才对。
这个系列记录这 9 版每一次错在哪、现象是什么、病根在哪。所有性能数字都是实测值,不是估算。
本篇覆盖前 5 版,错误集中在语法层、数据结构层和环境层—— 这三类的共同点是
报错信息还愿意告诉你一些东西。到了中篇,错误就开始沉默了。
如果你正在学 Python 做生信,这些坑你大概率会一个不漏地踩一遍。
一、任务
写一个函数read_fasta(path)，读 FASTA 文件，返回{序列ID: 序列字符串}。
听起来十行代码的事。但真实的 FASTA 从来不长教科书那样。测试文件我故意埋了 5 个坑：
>VNTR_01 chr11:3100000-3100120 unit=CAG CAGCAGCAGCAGCAGCAGCAGCAGCAGCAG CAGCAGCAGCAGCAGCAGCAGCAGCAGCAG CAGCAGCAG >VNTR_02 chr2:9800500-9800560 unit=AATG aatgaatgaatgaatgaatgaatgaatg AATGAATGAATG >VNTR_03 chr19:1200000-1200045 unit=CCG CCGCCGCCGNNNNNCCGCCGCCG CCGCCGCCG| # | 坑 | 为什么会有 |
|---|---|---|
| 1 | 序列被折成多行（line-wrapped） | 习惯 60 或 80 列一折，不能假设"一行标题一行序列" |
| 2 | 记录之间有空行 | 各种工具产出的格式不统一 |
| 3 | 大小写混杂 | soft-mask（软屏蔽）：重复区/低复杂度区被标成小写 |
| 4 | 序列里有N | 测序或组装没定下来的位置 |
| 5 | 最后一行没有换行符 | 手写解析漏掉最后一条记录的头号原因 |
正确答案：3 条记录，长度 69 / 40 / 32 bp。
解析思路：缓冲累积
因为序列被折行，你不能读一行存一行。必须：
- 遇到
>→ 说明上一条记录结束了，先结账（把攒的碎片拼起来存档），再开新（取新 ID、清空篮子） - 遇到序列行 → 扔进篮子
- 文件读完 →还得再结一次账（最后一条后面没有下一个
>来触发保存）
代码骨架：
defread_fasta(path):records={}# 账本:ID -> 序列seq_id=None# 当前在读谁chunks=[]# 篮子:当前这条的碎片withopen(path,encoding="utf-8")asf:forlineinf:line=line.strip()ifnotline:continueifline.startswith(">"):# ← 第 1 处:结账 + 开新passelse:# ← 第 2 处:扔进篮子pass# ← 第 3 处:收尾结账returnrecords三个待填处，我填错了 9 次。
二、九次迭代实录
第 1 版：把 C/Java 语法带进来了
ifline.startswith(">"):{records={chunks};line++;chunks=0;}全是别的语言的习惯。Python 的对应写法：
| 写的 | Python |
|---|---|
if x:{ ... }花括号划代码块 | 用缩进（4 个空格），没有花括号 |
chunks=0;行末分号 | 不写分号 |
line++ | n = n + 1或n += 1 |
open(D:/data.fa) | 路径是字符串，必须加引号 |
records={chunks} | records[键] = 值 |
核心观念：缩进即语法（significant whitespace）。行首空格数不是排版，是"这句归谁管"。
第 2 版：整本账本被覆盖成一张纸条
records="".join(chunks)# 我写的跑出来：
共读到 9 条记录 Traceback (most recent call last): File "read_fasta.py", line 65, in <module> AttributeError: 'str' object has no attribute 'items'读报错的标准动作：从最后一行往上读。
- 最后一行 = 病名：
'str' object has no attribute 'items'→ 对一个字符串调用了.items()，而.items()是字典专属的 - 倒数第二行 = 案发地点：第 65 行，
for rid, seq in recs.items(): - 然后自己推理：
recs本该是字典，现在是字符串 → 函数返回的不是字典 → 回去找谁把records变成了字符串
但真正的线索更早就出现了：报错前那句共读到 9 条记录。9 不是 3，而 9 恰好是文件最后一行CCGCCGCCG的长度 —— 说明records早就变成了那串字符，len()数的是字符数。
记住这条：报错行常常不是病根，报错前那句可疑输出才是。
病根：
records="".join(chunks)# 错:把整本账本扔了,换成一张纸条records[seq_id]="".join(chunks)# 对:在账本里写上一页records = 新东西——整个变量改换门庭，原来的字典被丢弃records[键] = 新东西——往字典里塞一条，字典依然是字典
这也是新手对方括号的三种语义没分清：
| 写法 | 名字 | 含义 |
|---|---|---|
line[1:]、seq[3:9] | 切片 slice | 取一段，有冒号 |
parts[0] | 索引 index | 取一个，没冒号，编号从 0 开始 |
records["VNTR_01"] = seq | 字典赋值 | 往字典里存 |
口诀：有冒号 = 切一段，没冒号 = 挑一个。
第 3 版：拼错了对象
records[seq_id]="".join(line)# 错:line 是当前这一行records[seq_id]="".join(chunks)# 对:chunks 是攒了一路的篮子"".join()的作用是"用空字符串当胶水，把列表里每块粘起来"。喂给它一个单行字符串，它会把每个字符重新粘一遍 —— 结果还是那一行。
第 4 版：找不到文件（其实文件在）
PS D:\Program Files\Microsoft VS Code> python d:/work/read_fasta.py Traceback (most recent call last): File "d:\work\read_fasta.py", line 26, in read_fasta with open(path, encoding="utf-8") as f: FileNotFoundError: [Errno 2] No such file or directory: 'test_vntr.fa'关键线索不在报错里，在提示符里：终端站在D:\Program Files\Microsoft VS Code。
黑话：当前工作目录（cwd, current working directory）。大白话：你敲命令时终端"站"在哪个文件夹。写相对路径（"test_vntr.fa"这种不带盘符的），Python 就从 cwd 开始找，跟脚本自己放在哪儿一点关系都没有。
importos os.chdir("D:/somewhere_else")print(os.path.abspath("test_vntr.fa"))# D:/somewhere_else/test_vntr.faprint(os.path.exists("test_vntr.fa"))# False -> FileNotFoundError三种修法：
# 修法 1(最标准):先 cd 过去再跑。路径带中文或空格时一定加引号cd"D:\work\fasta_practice"python read_fasta.py# 修法 2(应急):传绝对路径recs=read_fasta(r"D:\work\fasta_practice\test_vntr.fa")# 修法 3(专业做法,上集群靠它):让脚本自己找到自己importos here=os.path.dirname(os.path.abspath(__file__))# __file__ = 脚本自己的路径fa=os.path.join(here,"test_vntr.fa")recs=read_fasta(fa)修法 3 的价值：不管从哪个目录运行都能跑通。以后把脚本和数据一起丢到 HPC 上、由 SLURM 在某个莫名其妙的目录里启动它，全靠这一招。
顺带一个坑：我在文件顶部写了path = r"D:/.../test_vntr.fa"，以为这样就设好了路径 —— 它是死代码。
黑话：作用域（scope）与遮蔽（shadowing）。大白话：def read_fasta(path):括号里的path是函数自己的本地盒子，装的是调用时传进来的值。它跟外面那个同名全局变量毫无关系，只是撞名了。
第 5 版：钥匙锁在门里，读出 0 条
ifseq_idisnotNone:# 这道门的含义:"之前已经有过记录吗"records[seq_id]="".join(chunks)seq_id=line[1:].split()[0]# 错:取 ID 被关在门里结果：0 条，或者 1 条键为None的鬼记录。
逐行追踪就明白了：
| 读到 | seq_id状态 | 门开吗 | 后果 |
|---|---|---|---|
>VNTR_01 | None（初始值） | 关 | 取 ID 被跳过 →seq_id还是None |
| 69 bp 序列 | None | — | 进篮子 |
>VNTR_02 | 还是None | 还是关 | 又跳过；chunks=[]在门外照样执行 →69 bp 蒸发 |
| … | 永远None | 永远关 | 每条序列都被清空丢掉 |
开门的钥匙被锁在门里面了—— 门要求"必须已经有 ID 才能进"，而"拿 ID"这个动作在门内。
判断法则（这条最值钱，可反复用）：
问自己一句：读到第一个标题行时，这个动作需不需要做？
- 结账→ 不需要（第一条前面没有"上一条"）→放门里
- 取新 ID→ 需要（第一条也得有名字）→放门外
- 清空篮子→ 需要 →放门外
正确层级：
ifline.startswith(">"):ifseq_idisnotNone:# 门:之前有过记录吗records[seq_id]="".join(chunks)# 门内:结账seq_id=line[1:].split()[0]# 门外:每次都必须取新 IDchunks=[]# 门外:每次都必须清空顺带说line[1:].split()[0]这三道工序：
| 写法 | 干什么 | 结果 |
|---|---|---|
line | 原始标题行 | >VNTR_01 chr11:3100000-3100120 unit=CAG |
line[1:] | 切片：从第 2 个字符起，把>甩掉 | VNTR_01 chr11:3100000-3100120 unit=CAG |
.split() | 按空白切块，得到列表 | ['VNTR_01', 'chr11:...', 'unit=CAG'] |
[0] | 取第 0 块 | VNTR_01 |
为什么必须切掉>：它是 FASTA 格式的标记符号，不属于 ID。BED 文件、samtools faidx、pysam 里的序列名都不带>，带着它永远匹配不上。
为什么split()不填参数：默认按任意空白切，且连续空白视为一个 —— FASTA 和 SAM 的标题行有时空格有时 Tab，不填参数正好两种都吃。
上篇小结
前 5 版的错误可以归成四层：
| 层 | 错什么 | 报错帮不帮你 |
|---|---|---|
| 语法层 | 花括号、分号、line++、路径没加引号 | 帮:直接SyntaxError |
| 数据结构层 | records =覆盖字典、join(line)拼错对象 | 半帮:报错在下游,得自己往上推 |
| 环境层 | 工作目录不对、参数被全局变量遮蔽 | 半帮:说"文件不存在",但文件其实在 |
| 控制流层 | 取 ID 被锁在if门里 | 不帮:完全不报错 |
注意最后一行 —— 从第 5 版开始,错误不再报错了,只是安静地给你一个错结果。
中篇会讲最阴的四版:序列整体错位但assert照样通过、空文件冒出鬼记录、
三个测试全绿却慢 125 倍、else被抢走导致append变成死代码。
另外附一份报错速查表和两个能反复用的 debug 套路。
