title: Lesson 09｜FASTA / FASTQ Parser 与 Generator：逐条处理大型序列文件jupyter: python3
Lesson 08 我们已经完成了一个非常关键的转换：
磁盘上的文件↓open↓逐行读取↓strip / split↓Python object↓分析
我们处理的文件还是比较简单的 TSV：
seq1 ATGCGTseq2 GGCAATseq3 TTAACC
每一行就是一个完整 record：
所以：
seq_id, seq = line.split(”\t”)
就能完成解析。
但真实的生物信息学文件通常不会这么简单。
例如 FASTA：
>seq1ATGCGTACGT>seq2GGCCAATT>seq3ATATGC
看起来似乎也不复杂。
可是很快就会遇到：
>seq1ATGCGTACGTGGTTAACCCCGG>seq2GGCCAATTATGC
这时候：
一个 sequence record 可以跨很多行。
于是你不能再用：
seq_id, seq = line.split(”\t”)
这种“每一行就是一个 record”的思路。
这一课我们真正要学习的是：
raw lines↓识别 record boundary↓积累属于同一个 record 的内容↓组装一个完整 biological record↓逐条交给后续分析
这就是：
parser
而进一步，如果文件有几十 GB，我们又不希望：
解析整个文件↓全部保存到 memory↓最后才开始分析
更希望：
解析 seq1↓交给程序解析 seq2↓交给程序解析 seq3↓交给程序...
这就会引出这一课最重要的新工具：
以及：
generator
今天我们会完成三个真正的 bioinformatics parser：
FASTA → dictionary parserFASTA → generator parserFASTQ → generator parser
1. FASTA 的真正结构：不要把它理解成“两行一组”
先看最简单的 FASTA：
直觉上可能会觉得：
第 1 行 → header第 2 行 → sequence第 3 行 → header第 4 行 → sequence
于是想写：
header = ...sequence = ...
甚至认为：
这是一个非常危险的理解。
因为合法 FASTA 完全可能是：
>seq1ATGCGTACGTAACCGG>seq2GGGAAATT
这里：
真正的 sequence 是：
也就是说 FASTA 的逻辑不是：
headersequenceheadersequence
而是：
headersequence linesequence linesequence line...下一个 header
所以真正的 record boundary 是：
遇到下一个以 > 开头的 header。
1.1 FASTA record 的心智模型
例如：
>seq1 some descriptionATGCGTAACCGGTTAA
可以理解成：
record├── header│ └── seq1 some description└── sequence ├── ATGCGT ├── AACCGG └── TTAA
Parser 要做的是：
组合成：
所以我们需要：
当前 sequence ID+当前 sequence 的所有行
两个“状态”。
2. 第一个 FASTA parser：先返回 dictionary
先暂时不考虑超大文件。
我们希望：
内容：
>seq1ATGCGTAACCGG>seq2GGGAAA>seq3ATATGCGC
最后得到：
{ ”seq1”: ”ATGCGTAACCGG”, ”seq2”: ”GGGAAA”, ”seq3”: ”ATATGCGC”}
2.1 先观察每一行
with open(”sequences.fa”) as f: for line in f: line = line.strip() print(repr(line))
我们可以利用：
判断当前是不是 header。
例如：
if line.startswith(”>”): ...
否则：
就是 sequence line。
3. Parser 为什么需要“记住当前状态”
假设现在程序读到：
我们需要记住：
接下来看到：
不能立即结束 record，因为后面可能还有：
所以还需要一个容器：
然后：
sequence_parts.append(”ATGCGT”)sequence_parts.append(”AACCGG”)
得到：
等到遇到：
我们终于知道：
seq1 已经结束了。
于是：
得到：
然后保存：
sequences[current_id] = ””.join(sequence_parts)
再开始新的：
current_id = ”seq2”sequence_parts = []
整个状态变化可以画成：
>seq1↓current_id = ”seq1”sequence_parts = []ATGCGT↓[”ATGCGT”]AACCGG↓[”ATGCGT”, ”AACCGG”]>seq2↓先保存 seq1↓current_id = ”seq2”sequence_parts = []
这就是 parser 的核心。
4. 完整的 FASTA dictionary parser
def read_fasta(path): sequences = {} currend_id = None sequences_parts = [] with open(path) as f: for line in f: line = line.strip() if line == ””: continue if line.startswith(”>”): if current_id is not None: sequences[current_id] = ””.join(sequences_parts) currend_id = line[1:] sequences_parts = [] else: sequences_parts.append(line) if current_id is not None: sequences[current_id] = ””.join(sequences_parts) return sequences
这一段需要认真理解。
它第一次真正涉及：
parser state
4.1 为什么是 line[1:]
如果：
那么：
是：
而：
就是：
所以：
把 > 去掉。
如果 header 是：
>seq1 chromosome_1 gene_x
那么：
得到：
现在我们暂时把整个 header 保存下来。
后面需要时，可以进一步拆：
seq_id = line[1:].split()[0]
得到：
5. 一个非常容易漏掉的问题：最后一个 record
看看核心部分：
if line.startswith(”>”): if current_id is not None: sequences[current_id] = ””.join(sequence_parts) current_id = line[1:] sequence_parts = []
你会发现：
当前 record 是在“遇到下一个 header”时才被保存的。
比如：
>seq1AAAA>seq2CCCC>seq3GGGG
当看到：
保存：
看到：
保存：
但是文件结束之后：
后面没有新的 header。
所以如果我们不额外处理：
if current_id is not None: sequences[current_id] = ””.join(sequence_parts)
最后一个 sequence 就会消失。
这是 parser 中非常典型的：
end-of-file edge case
以后写 parser 时应该经常问自己：
6. 为什么使用 list + "".join()，而不是不断 +=
我们现在用了：
sequence_parts = []sequence_parts.append(line)sequence = ””.join(sequence_parts)
你可能会问：
为什么不直接：
sequence = ””sequence += line
对于很小的 sequence：
当然也能运行。
但 string 是：
immutable
也就是说：
本质上可能需要不断产生新的 string。
因此对于大量 sequence fragments，Python 中更典型的模式是：
parts = []parts.append(...)parts.append(...)parts.append(...)result = ””.join(parts)
类似：
这个模式以后会反复出现。
7. 给 FASTA parser 加验证
真实 parser 不应该默认输入永远正确。
例如错误文件：
sequence 出现在任何 header 之前。
或者：
header 是空的。
我们可以加强 parser：
def read_fasta(path): sequences = {} current_id = None sequence_parts = [] with open(path) as f: for line_number, line in enumerate(f, start=1): line = line.strip() if line == ””: continue if line.startswith(”>”): if current_id is not None: sequences[current_id] = ””.join(sequence_parts) current_id = line[1:].strip() if current_id == ””: raise ValueError( f”Empty FASTA header at line {line_number}” ) sequence_parts = [] else: if current_id is None: raise ValueError( f”Sequence before header at line {line_number}” ) sequence_parts.append(line) if current_id is not None: sequences[current_id] = ””.join(sequence_parts) return sequences
这里 Lesson 07 的：
开始真正参与 parser 设计。
8. Dictionary parser 有什么问题？
我们的：
最终返回：
{ ”seq1”: ”...”, ”seq2”: ”...”, ”seq3”: ”...”}
对于一个：
可能完全没问题。
可是如果是：
甚至更大的 metagenome assembly：
那么这个设计意味着：
读取 seq1↓放进 dictionary读取 seq2↓也放进 dictionary读取 seq3↓继续放...整个 FASTA↓全部留在 memory
但假设我们的目的只是：
其实没有必要把所有 sequence 同时留在内存。
我们真正希望：
读到 contig1↓算长度↓完成读到 contig2↓算长度↓完成
也就是：
一次产生一个 record。
这就是 generator 的工作。
9. yield：第一次认识 generator
先不看 FASTA。
从一个非常小的例子开始。
普通 function：
def get_numbers(): return [1, 2, 3]
调用：
此时：
已经是：
三个数字全部已经产生。
现在写：
def generate_numbers(): yield 1 yield 2 yield 3
调用：
numbers = generate_numbers()
你可能以为：
会是：
但实际上：
是：
真正使用：
for number in numbers: print(number)
才得到：
10. yield 和 return 的关键区别
普通：
意味着：
而：
意味着：
产生一个 value↓暂停 function↓保留当前位置和状态↓下次继续
例如：
def demo(): print(”A”) yield 1 print(”B”) yield 2 print(”C”)
执行：
此时甚至还没有打印：
然后：
输出：
并得到：
程序暂停在：
第二次：
从那里继续：
然后产生：
第三次：
继续：
然后 generator 结束。
11. Generator 的真正心智模型
可以把普通 function 想成：
function↓完成所有工作↓return 最终结果
而 generator：
function↓做到一部分↓yield 一个结果↓暂停 ↑ │调用者再次请求 │继续 ←─────↓yield 下一个结果
所以 generator 特别适合：
比如：
FASTA recordsFASTQ readsVCF variantsSAM alignments日志文件大型 TSV
12. 把 FASTA parser 改成 generator
现在我们不再返回：
{ ”seq1”: ”...”, ”seq2”: ”...”, ...}
而希望：
for seq_id, seq in read_fasta(”genome.fa”): print(seq_id, len(seq))
每解析完一个 sequence，就：
完整代码：
def read_fasta(path): current_id = None sequence_parts = [] with open(path) as f: for line_number, line in enumerate(f, start=1): line = line.strip() if line == ””: continue if line.startswith(”>”): if current_id is not None: sequence = ””.join(sequence_parts) yield current_id, sequence current_id = line[1:].strip() if current_id == ””: raise ValueError(f”Empty FASTA header at line {line_number}”) sequence_parts = [] else: if current_id is None: raise ValueError(f”Sequence before header at line {line_number}”) sequence_parts.append(line) if current_id is not None: sequence = ””.join(sequence_parts) yield current_id, sequence
这就是我们这一课最重要的一段代码。
13. Generator FASTA parser 究竟发生了什么？
假设：
>seq1AAAA>seq2CCCC>seq3GGGG
运行：
for seq_id, seq in read_fasta(”test.fa”): print(seq_id, seq)
程序大致经历：
打开文件读 >seq1↓记录 seq1读 AAAA↓保存到 sequence_parts读 >seq2↓seq1 已经完整↓yield (”seq1”, ”AAAA”)
此时 parser：
暂停。
调用者拿到：
执行：
然后下一轮 for 再向 generator 请求 record。
Parser 从暂停位置继续：
开始 seq2↓读 CCCC↓遇到 >seq3↓yield (”seq2”, ”CCCC”)
如此继续。
所以内存中通常只需要保留：
而不是：
14. Streaming：为什么 generator 特别适合 genomics
假设一个 FASTA：
包含：
如果：
sequences = read_all_fasta(...)
返回 dictionary：
可能根本不可行。
而 generator：
for seq_id, seq in read_fasta(path): ...
则近似：
磁盘 ↓当前 contig ↓分析 ↓丢掉 ↓下一个 contig
这就是：
streaming
要注意 generator 并不神奇地让 sequence 本身消失。
如果某一条 chromosome 本身：
你仍然需要保存这一条 sequence 的内容。
但是你不再需要同时保存：
这已经是巨大的差别。
15. 一个真正实用的 FASTA streaming 分析
假设我们只想统计：
sequence IDlengthGC content
完全不需要把整个 FASTA 保存成 dictionary。
可以：
def gc_content(seq): if len(seq) == 0: raise ValueError(”Sequence is empty”) gc = seq.count(”G”) + seq.count(”C”) return gc / len(seq)with open(”summary.tsv”, ”w”) as outfile: outfile.write(”seq_id\tlength\tgc_content\n”) for seq_id, seq in read_fasta(”genome.fa”): length = len(seq) gc = gc_content(seq) outfile.write(f”{seq_id}\t{length}\t{gc}\n”)
这里的数据流非常漂亮：
genome.fa ↓read_fasta() ↓一个 sequence ↓len / gc_content ↓写 summary.tsv ↓丢掉 sequence ↓下一个
整个分析不需要：
这种变量。
这就是大型 sequencing workflow 非常重要的一种设计方式。
16. Generator 是 iterator，所以可以直接 for
Lesson 05 我们已经用过：
以前的：
通常是：
现在 generator 也可以：
forrecordinread_fasta(path): ...
因为 generator 是一种：
iterator
你可以暂时这样建立关系：
iterable└── 可以被 for 遍历iterator└── 可以逐个产生下一个值generator└── 一种很方便的 iterator
现阶段不需要深入 Python iterator protocol。
只需要记住：
让 function 变成：
调用它后得到：
然后：
会不断向它请求下一个结果。
17. Generator 只能走一次
这是一个很重要的特点。
g = read_fasta(”genome.fa”)
第一次：
for seq_id, seq in g: print(seq_id)
会遍历所有 sequence。
然后再：
for seq_id, seq in g: print(seq_id)
通常什么都不会输出。
因为 generator 已经：
也就是已经走到结尾。
如果想重新读取：
g = read_fasta(”genome.fa”)
重新创建一个 generator。
或者直接：
for seq_id, seq in read_fasta(”genome.fa”): ...
以后需要第二遍时再次：
for seq_id, seq in read_fasta(”genome.fa”): ...
18. 不要随手 list(generator)，否则可能失去 streaming 优势
你当然可以：
records = list( read_fasta(”genome.fa”))
得到：
[ (”seq1”, ”...”), (”seq2”, ”...”), ...]
但是这意味着：
generator↓所有 records↓list↓全部进入 memory
所以如果你写 generator 本来是为了处理：
然后马上：
那就基本把 streaming 的优势又取消了。
当然对于测试小文件：
records = list(read_fasta(”tiny.fa”))
非常方便。
关键不是：
而是要知道：
你是在主动选择把所有结果 materialize 到 memory。
19. FASTQ：record boundary 又不一样了
FASTA 的 record：
FASTQ 的一个标准 record 通常是四行：
分别是：
line 1 → headerline 2 → sequenceline 3 → +line 4 → quality
连续 records：
@read1ATGCGT+IIIIII@read2GGCAAT+HHHHHH
所以 FASTQ parser 的逻辑和 FASTA 不完全相同。
FASTA：
FASTQ：
因此它反而更适合用：
一次读取四行。
20. 一个最基础的 FASTQ parser
def read_fastq(path): with open(path) as f: while True: header = f.readline() if header == ””: break sequence = f.readline() plus = f.readline() quality = f.readline() header = header.strip() sequence = sequence.strip() plus = plus.strip() quality = quality.strip() yield header, sequence, quality
这里第一次在文件 parser 中真正使用：
流程是：
读 header↓如果 EOF → break↓读 sequence↓读 +↓读 quality↓yield 一个 read↓下一轮
21. 为什么用 header == "" 判断 EOF，而不是 strip() == ""
这一点值得注意。
在真正到达文件结尾时返回：
所以：
表示：
没有更多数据。
但如果文件中只是一个空行：
它并不等于：
这两个概念不同：
”\n”→ 文件里真的有一行，只是内容为空””→ EOF，没有任何字符可读
这也是为什么文件处理时：
具有特殊意义。
22. FASTQ parser 必须做验证
上面的 parser 太信任输入。
假设文件坏了：
quality 行没了。
我们的 parser 可能仍然产生一个奇怪的 record。
所以应该检查：
def read_fastq(path): with open(path) as f: record_number = 0 while True: header = f.readline() if header == ””: break sequence = f.readline() plus = f.readline() quality = f.readline() recored_number += 1 if(sequence == ”” or plus == ”” or quality == ””): raise ValueError(f”Incomplete FASTQ record {record_number}”) header = header.strip() sequence = sequence.strip() plus = plus.strip() quality = quality.strip() if not header.startswith(”@”): raise ValueError(F”Invalid FASTQ header in recored {record_number}”) if not plus.startswith(”+”): raise ValueError(F”Invalid FASTQ separator in recored {record_number}”) if len(sequence) != len(quality): raise ValueError(f”Sequence/quality length mismatch in recored {record_number}”) yield header, sequence, quality
现在 parser 已经相当像真正的数据工具了。
23. 为什么 FASTQ 的 sequence 和 quality 长度必须相等？
例如：
sequence：
一共：
quality：
也是：
因为每一个 nucleotide 都对应一个 quality score。
所以：
len(sequence) == len(quality)
应该成立。
如果：
对应：
只有 5 个 quality character，就说明 record 有问题。
因此：
if len(sequence) != len(quality): raise ValueError(...)
是非常自然的 validation。
24. 为什么这里 yield header[1:], sequence, quality
原始 header：
和 FASTA 类似，我们往往不需要：
所以：
得到：
于是：
for read_id, seq, qual in read_fastq(”reads.fastq”): ...
非常自然。
25. FASTA 和 FASTQ parser 的差别
现在可以比较：
| | |
|---|
| > | @ |
| | 常见格式中 record 内 sequence 为一行 |
| | |
| | |
| 当前 header + sequence parts | |
| (seq_id, seq) | (read_id, seq, qual) |
所以：
Parser 不是背一个固定模板。
真正重要的是：
先理解文件格式↓找出 record boundary↓决定程序需要保存什么 state↓逐条产生 structured record
26. Parser 和 analysis 应该分开
一个非常值得现在建立的软件设计习惯是：
不要写：
def read_fasta_and_calculate_gc_and_filter_and_write(...): ...
更推荐：
def read_fasta(path): ...
只负责：
然后：
只负责：
最后主程序组合：
for seq_id, seq in read_fasta(path): gc = gc_content(seq)
为什么？
因为同一个：
以后可以被不同任务复用。
比如：
for seq_id, seq in read_fasta(path): print(seq_id, len(seq))
或者：
for seq_id, seq in read_fasta(path): print(seq_id, gc_content(seq))
或者：
for seq_id, seq in read_fasta(path): kmers = generate_kmers(seq, 31)
或者：
for seq_id, seq in read_fasta(path): if len(seq) >= 1000: ...
Parser 不需要知道：
“用户之后想对 sequence 做什么。”
它只负责：
这叫：
separation of concerns
也是以后写科研软件非常重要的设计原则。
27. 一个完整 pipeline：过滤短 contig 并写 FASTA
假设：
我们想保留：
的 contig。
现在已经非常简单：
def filter_fasta(input_path, output_path, min_length): with open(output_path, ”w”) as outfile: for seq_id, seq in read_fasta(input_path): if len(seq) < min_length: continue outfile.write(f”>{seq_id}\n”) outfile.write(f”{seq}\n”)
调用：
filter_fasta( ”assembly.fa”, ”assembly.filtered.fa”, 1000)
数据流：
assembly.fa↓read_fasta↓一个 contig↓length filter↓写入 filtered.fa↓下一个 contig
注意这里甚至不需要：
因为每一个通过过滤的 record 可以立即写入文件。
这就是：
stream input+stream processing+stream output
28. FASTA sequence 输出时为什么有时要 wrap
我们刚才输出：
outfile.write(f”{seq}\n”)
于是 sequence 全部在一行。
例如：
>chr1ATGCGTACGTACGTACGTACGTACGT...
这在很多情况下是合法的。
但是很多 FASTA 文件会把 sequence 按：
等长度换行。
例如：
>seq1ATGCGTACGTAACCGGTTAAGGCC...
现在我们已经有能力自己完成。
例如每 80 个字符：
for i in range(0, len(seq), 80): outfile.write( seq[i:i + 80] + ”\n” )
这里重新使用了：
以前学到的内容又回来了。
29. 写一个 write_fasta_record()
为了不让写 FASTA 的代码到处重复，可以单独封装：
def write_fasta_record(outfile, seq_id, seq, width=80): outfile.write(f”>{seq_id}\n”) for i in range(0, len(seq), width): outfile.write(seq[i:i + width] + ”\n”)
然后过滤程序变成：
def filter_fasta(input_path, output_path, min_length): with open(output_path, ”w”) as outfile: for seq_id, seq in read_fasta(input_path): if len(seq) < min_length: continue write_fasta_record(outfile, seq_id, seq)
结构明显更清楚：
read_fasta→ 负责读filter_fasta→ 负责决定留不留write_fasta_record→ 负责 FASTA 输出格式
30. 不要默认 header 只有一个单词
真实 FASTA header 可能是：
>NC_000001.11 Homo sapiens chromosome 1
如果：
得到：
NC_000001.11 Homo sapiens chromosome 1
如果只想要 identifier：
current_id = line[1:].split()[0]
得到：
这两个设计都可能是正确的。
关键是：
Parser 应该明确自己保留的是完整 header，还是第一个 identifier field。
不要无意识地丢掉 metadata。
例如更清楚的设计可能是：
header = line[1:].strip()seq_id = header.split()[0]
然后决定：
还是：
31. 不要在 parser 里过度“清洗”数据
假设 FASTA sequence：
你当然可以：
sequence = sequence.upper()
但这属于一个设计决定。
Parser 是否应该：
还是：
没有永远唯一答案。
对于科研程序，更重要的是：
行为应该明确，而不是偷偷修改输入。
例如：
def read_fasta(path, uppercase=False): ...
然后：
if uppercase: sequence = sequence.upper()
会比无条件修改更清楚。
现在不必立刻把 parser 设计得这么复杂，但需要开始有这种意识。
32. 一个更完整的 FASTA generator
到目前为止，我们可以写成：
from pathlib import Pathdef read_fasta(path): path = Path(path) if not path.exists(): raise FileNotFoundError(f”FASTA file no found: {path}”) current_id = None sequence_parts = [] with open(path, ”r”) as f: for line_number, line in enumerate(f, start=1): line = line.strip() if line == ””: continue if line.startswith(”>”): if current_id is not None: sequence = ””.join(sequence_parts) yield current_id, sequence current_id = line[1:].strip() if current_id == ””: raise ValueError(f”Empty FASTA header at line {line_number}”) sequence_parts = [] else: if current_id is None: raise ValueError(f”sequence before FASTA header at line {line_number}”) sequence_parts.append(line) if current_id is not None: sequence = ””.join(sequence_parts) yield current_id, sequence
调用
for seq_id, seq in read_fasta(”genome.fa”): print(seq_id, len(seq))
这已经是一个真正有实用价值的 parser。
33. 一个更完整的 FASTQ generator
from pathlib import Pathdef read_fastq(path): path = Path(path) if not path.exists(): raise FileNotFoundError(f”FASTQ file no found: {path}”) with open(path, ”r”) as f: record_number = 0 while True: header = f.readline() if header == ””: break sequence = f.readline() plus = f.readline() quality = f.readline() record_number += 1 if(sequence == ”” or plus == ”” or quality == ””): raise ValueError(f”Incomplete FASTQ record {record_number}”) header = header.strip() sequence = sequence.strip() plus = plus.strip() quality = quality.strip() if not header.startswith(”@”): raise ValueError(f”Invalid FASTQ header in record {record_number}”) if not plus.startswith(”+”): raise ValueError(f”Invalid FASTQ separator in record {record_number}”) read_id = header[1:] yield read_id, sequence, quality
for read_id, seq, qual in read_fastq( ”reads.fastq” ): print( read_id, len(seq) )
34. FASTQ quality 现在先当作 string
例如：
目前：
只是：
我们暂时不解释：
ASCIIPhred+33quality score conversion
因为那会把这一课的重点从：
parser + generator + streaming
带到 sequencing quality encoding。
现在只需要知道：
每一个 quality character↔一个 nucleotide
所以：
len(sequence) == len(quality)
是一个重要的格式检查。
35. yield 后面的值也可以是 tuple
我们写：
其实等价于产生：
所以：
forrecordinread_fasta(path): print(record)
得到：
('seq1', 'ATGC...')('seq2', 'GGAA...')
然后：
for seq_id, sequence in read_fasta(path):
只是再次使用：
tuple unpacking
也就是说：
yield+tuple+for+unpacking
又组合起来了。
36. next()：手动看看 generator 怎样工作
假设：
g = read_fasta(”tiny.fa”)
可以：
得到：
再：
得到：
如果没有更多结果，再：
会出现：
通常你不需要自己处理它。
因为：
已经自动帮你处理。
也就是说：
for record in generator: ...
本质上不断做类似：
直到：
37. Generator 让 pipeline 非常自然
有了：
之后，可以直接：
for seq_id, seq in read_fasta(”assembly.fa”): if len(seq) < 1000: continue gc = gc_content(gc) print(seq_id, len(seq), gc)
甚至可以把不同处理阶段继续拆开。
例如：
def long_sequences(records, min_length): for seq_id, seq in records: if len(seq) >= min_length: yield seq_id, sqe
于是：
records = read_fasta(”assembly.fa”)records = long_sequences( records, min_length=1000)for seq_id, seq in records: print(seq_id)
数据像水一样流过：
FASTA↓read_fasta↓long_sequences↓analysis
这种设计以后在：
genomicsdata engineeringcommand-line pipelines
里都会不断遇到。
38. Generator 并不意味着“代码一定更快”
需要纠正一个容易产生的误解：
它最大的优势通常是：
memory efficiency+lazy evaluation+streaming
例如：
需要把大量 integer 全部放进 memory。
而：
本身不会生成一个装满一亿个 integer 的 list。
Generator 也类似：
这叫：
lazy evaluation
它可能改善整体 workflow，但重点不是“神奇加速”，而是：
不提前计算和储存所有结果。
39. R → Python：generator 可以怎样理解？
R 里很多常见 workflow 更偏：
读取完整 object↓vector / dataframe↓操作整个 object
Python generator 的思路更接近：
所以不要强行寻找一个：
R generator = Python generator
的一对一对应。
真正值得建立的新思维是：
不是所有数据都必须先完整存在内存里，然后程序才能开始分析。
这对 genomics 尤其重要。
40. 本课综合程序：FASTA summary
现在把 Lesson 01–09 串起来。
目标：
input:assembly.faoutput:assembly_summary.tsv
输出：
seq_id length gc_contentcontig1 18342 0.4832contig2 9271 0.5011...
代码：
from pathlib import Pathdef gc_content(seq): if len(seq) == 0: raise ValueError(”Sequence is Empty”) gc = seq.count(”C”) + seq.count(”G”) return gc / len(seq)def read_fasta(path): path = Path(path) if not path.exists(): raise FileNotFoundError(f”FASTA file not found: {path}”) current_id = None sequence_parts = [] with open(path) as f: for line_number, line in enumerate(f, start=1): line = line.strip() if line == ””: continue if line.startswith(”>”): if current_id is not None: yield current_id, ””.join(sequence_parts) current_id = line[1:].strip() if current_id == ””: raise ValueError(f”Empty FASTA header at line {line_number}”) sequence_parts = [] else: if current_id is None: raise ValueError(f”Sequence before header at line {line_number}”) sequence_parts.append(line) if current_id is not None: yield(current_id, ””.join(sequence_parts))def summarize_fasta(input_path, output_path): with open(output_path, ”w”) as outfile: outfile.write(”seq_id\tlength\tgc_content\n”) for seq_id, seq in read_fasta(input_path): outfile.write( f”{seq_id}\t” f”{len(seq)}\t” f”{gc_content(seq)}\n” )
summarize_fasta( ”genome.fa”, ”genome_summary.tsv” )
这段程序已经非常接近一个小型真实 bioinformatics tool。
这里串起了 Lesson 01–09
String
line.strip()line.startswith(”>”)””.join(...)
List
sequence_parts = []sequence_parts.append(line)
Tuple
Unpacking
for seq_id, seq in read_fasta(...):
Condition
Loop
和：
continue / break
Function
gc_content()read_fasta()summarize_fasta()
Exception
raise ValueError(...)raise FileNotFoundError(...)
File IO
Path
Generator
这已经不是孤立语法了。
它们开始形成：
Exercises
Exercise A｜预测 FASTA parser 的输出
文件：
内容：
执行：
for seq_id, seq in read_fasta(”test.fa”): print(seq_id, seq, len(seq))
问输出是什么？
第一个 record：
需要拼接成：
长度：
第二个：
长度：
因此：
seq1 ATGCGTAA 8seq2 CCGG 4
这里最重要的是：
FASTA 的换行通常不是 sequence 本身的一部分。
Exercise B｜找出 parser bug
下面 parser 有问题：
def read_fasta(path): current_id = None parts = [] with open(path) as f: for line in f: line = line.strip() if line.startswith(”>”): if current_id is not None: yield current_id, ””.join(parts) current_id = line[1:] parts = [] else: parts.append(line)
对于：
它只会产生：
为什么？
应该怎样修改？
因为 record 是在：
时才：
所以：
在遇到：
时会被产生。
但是：
后面没有第三个 header。
文件直接结束了。
因此需要在 loop 之后：
if current_id is not None: yield current_id, ””.join(parts)
完整修复：
def read_fasta(path): current_id = None parts = [] with open(path) as f: for line in f: line = line.strip() if line.startswith(”>”): if current_id is not None: yield current_id, ””.join(parts) current_id = line[1:] parts = [] else: parts.append(line) if current_id is not None: yield current_id, ””.join(parts)
这是 parser 中非常典型的：
Exercise C｜只打印长度超过 1000 bp 的 contig
假设已经有：
要求：
读取 assembly.fa只输出 length >= 1000 的 sequence ID 和 length
不要把所有 sequence 放入 list。
for seq_id, seq in read_fasta(”assembly.fa”): if len(seq) > 1000: print(seq_id, len(seq))
Exercise D｜计算 FASTA 中 sequence 的数量
要求写：
def count_sequences(path): ...
利用：
不要再重新解析一次 FASTA 格式。
答案
def count_sequences(path): count = 0 for seq_id, seq in read_fasta(path): count += 1 return count
调用：
n = count_sequences(”genome.fa”)
也可以：
def count_sequences(path): return sum( 1 for _ in read_fasta(path) )
不过第二种涉及 generator expression。
现在能看懂即可，不要求优先这样写。
第一种更适合当前阶段。
Exercise E｜找到最长 contig
要求：
def longest_sequence(path): ...
返回：
例如：
>seq1AAAA>seq2CCCCCCCC>seq3GGG
返回：
不能把整个文件先变成 dictionary。
答案
def longest_sequence(path): longest_id = None longest_seq = ”” for seq_id, seq in read_fasta(path): if len(seq) > len(longest_seq): longest_id = seq_id longest_seq = seq if longest_id is None: raise ValueError( ”FASTA file contains no sequences” ) return longest_id, longest_seq
这里虽然整个 FASTA 可以很大，但程序只保留：
当前 sequence+目前发现的最长 sequence
而不是所有 sequence。
这就是 streaming algorithm 的一个典型例子。
Exercise F｜写一个 FASTA filter
要求：
def filter_fasta( input_path, output_path, min_length): ...
只保留：
的 sequence。
输出要求 sequence 每行最多：
def filter_fasta(input_path, output_path, min_length, width = 60): with open(output_path, ”w”) as outfile: for seq_id, seq in read_fasta(input_path): if len(seq) >= min_length: outfile.write(f”>{seq_id}\n”) for i in range(0, len(seq), width): outfile.write(seq[i:i + width] + ”\n”)
def write_fasta_record(outfile, seq_id, seq, width=60): outfile.write(f”>{seq_id}\n”) for i in range(0, len(seq), width): outfile.write(seq[i:i + width] + ”\n”)def filter_fasta(input_path, output_path, min_length): with open(output_path, ”w”) as outfile: for seq_id, seq in read_fasta(input_path): if len(seq) < min_length: continue write_fasta_record(outfile,seq_id,seq)
Exercise G｜判断下面代码是否仍然是 streaming
代码 1：
for seq_id, seq in read_fasta( ”genome.fa”): print(seq_id)
代码 2：
records = list( read_fasta(”genome.fa”))for seq_id, seq in records: print(seq_id)
哪一个保持了 generator 的主要 memory 优势？
代码 1：
for seq_id, seq in read_fasta(...):
保持 streaming。
大致：
代码 2：
records = list(read_fasta(...))
会首先把 generator 中的所有结果收集成：
因此：
所以对于大型 FASTA：
通常会失去 generator 最重要的 memory advantage。
但对于小测试文件：
仍然非常方便。
Exercise H｜FASTQ parser
假设：
内容：
@read1ATGCGT+IIIIII@read2GGCA+HHHH
要求：
for read_id, seq, qual in read_fastq(...): ...
依次产生：
(”read1”, ”ATGCGT”, ”IIIIII”)(”read2”, ”GGCA”, ”HHHH”)
并检查：
header 以 @ 开头第三行以 + 开头sequence 与 quality 长度一致
def read_fastq(path): path = Path(path) with open(path) as f: record_number = 0 while True: header = f.readline() if header == ””: break sequence = f.readline() plus = f.readline() quality = f.readline() record_number += 1 if (sequence == ”” or plus == ”” or quality == ””): raise ValueError(f”Incomplete FASTQ record {record_number}”) header = header.strip() sequence = sequence.strip() plus = plus.strip() quality = quality.strip() if not header.startswith(”@”): raise ValueError(f”Invalid header in FASTQ record {record_number}”) if not plus.startswith(”+”): raise ValueError(f”Invalid separator in FASTQ record {record_number}”) if len(sequence) != len(quality): raise ValueError(f”Sequence/quality length mismatch in FASTQ record {record_number}”) yield (header[1:], sequence, quality)
for read_id, seq, qual in read_fastq(”reads.fastq”): print(read_id, seq, qual)
Exercise I｜综合题：统计 FASTQ reads
写：
def summarize_fastq( input_path, output_path): ...
输入：
@read1ATGCGT+IIIIII@read2GGGG+HHHH@read3ATATAT+JJJJJJ
输出：
read_id length gc_contentread1 6 0.5read2 4 1.0read3 6 0.0
要求：
复用 read_fastq()复用 gc_content()逐条读取不要把全部 reads 放入 list
def summarize_fasta(input_path, output_path): with open(output_path, ”w”) as outfile: outfile.write(f”read_id\tlength\tgc_content\n”) for read_id, seq, qual in read_fastq(input_path): length = len(seq) gc = gc_content(seq) outfile.write(f”{read_id}\t{length}\t{gc}\n”)
Exercise J｜本课最重要的思考题
下面两种程序设计：
sequences = read_all_fasta(”assembly.fa”)for seq_id, seq in sequences.items(): ...
for seq_id, seq in read_fasta( ”assembly.fa”): ...
什么时候设计 1 可能更合适？
什么时候设计 2 更合适？
如果：
例如：
或者同一批 sequences 要被重复使用很多次，那么 dictionary：
可能非常方便。
因为数据已经全部存在 memory。
如果：
例如：
那么 generator：
for seq_id, seq in read_fasta(...):
通常更适合。
所以并不是：
generator 永远比 dictionary 好
真正的判断标准是：
这其实已经开始进入：
程序设计与算法选择
而不仅仅是 Python syntax。
# Lesson 09 最重要的几个心智模型## ① Parser 的核心不是 `split()`Lesson 08 的 TSV：```text一行 = 一个 record
所以：
就够了。
但 FASTA 告诉我们：
真正的 parser 必须：
识别 record boundary+保存 state+组装完整 record
