手写 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 unitCAG CAGCAGCAGCAGCAGCAGCAGCAGCAGCAG CAGCAGCAGCAGCAGCAGCAGCAGCAGCAG CAGCAGCAG VNTR_02 chr2:9800500-9800560 unitAATG aatgaatgaatgaatgaatgaatgaatg AATGAATGAATG VNTR_03 chr19:1200000-1200045 unitCCG 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_idNone# 当前在读谁chunks[]# 篮子:当前这条的碎片withopen(path,encodingutf-8)asf:forlineinf:lineline.strip()ifnotline:continueifline.startswith():# ← 第 1 处:结账 开新passelse:# ← 第 2 处:扔进篮子pass# ← 第 3 处:收尾结账returnrecords三个待填处我填错了 9 次。二、九次迭代实录第 1 版把 C/Java 语法带进来了ifline.startswith():{records{chunks};line;chunks0;}全是别的语言的习惯。Python 的对应写法写的Pythonif x:{ ... }花括号划代码块用缩进4 个空格没有花括号chunks0;行末分号不写分号linen n 1或n 1open(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, encodingutf-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 过去再跑。路径带中文或空格时一定加引号cdD:\work\fasta_practicepython read_fasta.py# 修法 2(应急):传绝对路径recsread_fasta(rD:\work\fasta_practice\test_vntr.fa)# 修法 3(专业做法,上集群靠它):让脚本自己找到自己importos hereos.path.dirname(os.path.abspath(__file__))# __file__ 脚本自己的路径faos.path.join(here,test_vntr.fa)recsread_fasta(fa)修法 3 的价值不管从哪个目录运行都能跑通。以后把脚本和数据一起丢到 HPC 上、由 SLURM 在某个莫名其妙的目录里启动它全靠这一招。顺带一个坑我在文件顶部写了path rD:/.../test_vntr.fa以为这样就设好了路径 —— 它是死代码。黑话作用域scope与遮蔽shadowing。大白话def read_fasta(path):括号里的path是函数自己的本地盒子装的是调用时传进来的值。它跟外面那个同名全局变量毫无关系只是撞名了。第 5 版钥匙锁在门里读出 0 条ifseq_idisnotNone:# 这道门的含义:之前已经有过记录吗records[seq_id].join(chunks)seq_idline[1:].split()[0]# 错:取 ID 被关在门里结果0 条或者 1 条键为None的鬼记录。逐行追踪就明白了读到seq_id状态门开吗后果VNTR_01None初始值关取 ID 被跳过 →seq_id还是None69 bp 序列None—进篮子VNTR_02还是None还是关又跳过chunks[]在门外照样执行 →69 bp 蒸发…永远None永远关每条序列都被清空丢掉开门的钥匙被锁在门里面了—— 门要求必须已经有 ID 才能进而拿 ID这个动作在门内。判断法则这条最值钱可反复用问自己一句读到第一个标题行时这个动作需不需要做结账→ 不需要第一条前面没有上一条→放门里取新 ID→ 需要第一条也得有名字→放门外清空篮子→ 需要 →放门外正确层级ifline.startswith():ifseq_idisnotNone:# 门:之前有过记录吗records[seq_id].join(chunks)# 门内:结账seq_idline[1:].split()[0]# 门外:每次都必须取新 IDchunks[]# 门外:每次都必须清空顺带说line[1:].split()[0]这三道工序写法干什么结果line原始标题行VNTR_01 chr11:3100000-3100120 unitCAGline[1:]切片从第 2 个字符起把甩掉VNTR_01 chr11:3100000-3100120 unitCAG.split()按空白切块得到列表[VNTR_01, chr11:..., unitCAG][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 套路。