前面七课，我们处理的数据基本都直接写在 Python 代码里：
seq = ”ATGCGT”reads = [”ATGC”, ”GGAA”, ”CCGT”]
然后再把它们交给 function：
gc_content(seq)generate_kmers(seq, 3)find_motif(seq, motif)
但真正做生物信息学时，数据当然不会这样存在。
你面对的是：
genome.fareads.fastqvariants.vcfannotation.gtfsample.bed
所以从 Lesson 08 开始，我们要把前面学到的 Python 世界和磁盘上的真实数据连接起来。
这一课最重要的流程是：
文件路径 ↓open ↓file object ↓读取 text ↓解析 text ↓Python object ↓分析 ↓写回文件
也就是说，今天真正学习的不是几个 open() 方法，而是：
怎样让数据从磁盘进入 Python，再经过处理后回到磁盘。
这会直接为 Lesson 09 的 FASTA / FASTQ parser 做准备。
1. open() 与 with：Python 怎样访问文件
1.1 从一个最简单的文本文件开始
假设当前目录中有：sequence.txt
内容是：ATGCGTACGT
Python 可以：
with open(”sequence.txt”, ”r”) as f: sequence = f.read()
现在：
可能得到：
这里第一次出现几个新东西：
”sequence.txt” → 文件路径”r” → read modeopen(...) → 打开文件f → file objectf.read() → 读取文件内容
其中真正需要建立的概念是：
open() 并不是“把文件变成字符串”。
而是首先创建一个file object，通过这个对象和磁盘上的文件交互。
可以想成：
sequence.txt ↑ │ ffile object
然后：
才把其中的数据读取出来。
1.2 为什么总是推荐 with open(...) as f
你当然可以写：
f = open(”sequence.txt”, ”r”)sequence = f.read()f.close()
但是更推荐：
with open(”sequence.txt”, ”r”) as f: sequence = f.read()
原因是：
离开 with 代码块后，Python 会自动关闭文件。
也就是说：
负责：打开文件→使用文件→自动关闭。
这叫context manager。
现在不需要深入研究它的实现机制。
对于你目前来说，直接形成习惯：
就很好。
尤其科研代码中，一个脚本可能会处理成千上万个文件，忘记关闭文件并不是一个值得保留的习惯。
1.3 "r"、"w"、"a"：你准备怎样使用文件？
最常见的三个 mode：
r → readw → writea → append
with open(”sequence.txt”, ”r”) as f: ...
实际上 "r" 是默认值，所以：
with open(”sequence.txt”) as f: ...
也可以。
但初期写出来有助于理解。
with open(”output.txt”, ”w”) as f: ...
非常重要：
"w" 如果发现文件已经存在，会覆盖原文件。
例如原来output.txt里面有 1000 行结果。
执行：
with open(”output.txt”, ”w”) as f: f.write(”hello”)
原来的内容就没有了。
科研分析里一定要对这一点有意识。
with open(”log.txt”, ”a”) as f: ...
表示：保留旧内容 + 把新内容加到最后。
不过数据分析结果是否适合用 "a"，取决于任务本身，不应该为了“不覆盖”就全部使用 append。
2. 读取文件：整个读取，还是逐行处理？
这是 genomics 中非常重要的区别。
假设sequences.txt内容：
有几种读取方式。
2.1 read()：整个文件一次读进来
with open(”sequences.txt”, ”r”) as f: text = f.read()
得到类似：
'ATGCGT\nGGCAAT\nTTAACC\n'
注意：\n 代表 newline。
也就是说磁盘中看起来是：
ATGCGT GGCAAT TTAACC
Python string 实际可以理解成：
ATGCGT\nGGCAAT\nTTAACC\n
这就是为什么：
可以把它拆开。
如果文件很小：几KB，几MB，整个读入内存通常没有问题。
但你的研究以后很可能面对：
5 GB FASTA 30 GB FASTQ 100 GB sequencing data
这时候f.read()的思路就值得警惕。 因为它意味着：
把整个文件内容一次加载到内存。
2.2 for line in f：逐行处理
对大文件，更重要的是：
with open(”sequence.txt”, ”r”) as f: for line in f: print(line)
这里的思维完全不同：
文件 ↓读取第 1 行 → 处理 ↓读取第 2 行 → 处理 ↓读取第 3 行 → 处理 ↓...
而不是：
这其实是你以后处理 FASTA / FASTQ 时最重要的编程习惯之一。
对于 genomics：
能逐条处理，就不要默认把整个大文件读进内存。
2.3 readline() 和 readlines()
你也会在别人代码里看到：
它读取一行。
例如：
with open(”sequences.txt”) as f: first_line = f.readline()
而：
会把所有行读取成一个 list：
[ ”ATGCGT\n”, ”GGCAAT\n”, ”TTAACC\n”]
所以可以简单理解：
read()→ 整个文件 → 一个 stringreadline()→ 一次一行 → 一个 stringreadlines()→ 所有行 → 一个 listfor line in f→ 一行一行迭代
对于你后面的大型 genomics 文件，我最希望你习惯的是：
而不是：
因为后者仍然会一次把所有行放进内存。
3. 文本解析：strip() 和 split() 才是真正的主角
打开文件只是第一步。
真正的数据处理往往发生在：
raw text↓parse↓structured data
这个过程中。
3.1 为什么经常看到 line.strip()
假设文件：
执行：
with open(”sequences.txt”) as f: for line in f: print(repr(line))
每一行后面有：\n。因此，经常写：
with open(”sequences.txt”) as f: for line in f: line = line.strip() print(repr(line))
前面 Lesson 02 已经认识过 strip()。
现在终于看到它为什么在真实文件处理中如此常见。
一个非常典型的模式就是：
with open(”sequences.txt”) as f: for line in f: line = line.strip() if line == ””: continue print(line)
意思是：
逐行读取↓去掉换行等首尾空白↓跳过空行↓处理真正的数据
3.2 split()：把一行文本拆成字段
真实文件通常不是：
这么简单。
比如一个 tab-separated 文件：
seq1 ATGCGT seq2 GGCAAT seq3 TTAACC
这里实际上是：
seq1\tATGCGT seq2\tGGCAAT seq3\tTTAACC
\t 表示 tab。
逐行读取：
with open(”sequences.tsv”) as f: for line in f: line = line.strip() seq_id, seq = line.split(”\t”) print(seq_id, seq)
对于：
执行：
得到：
然后利用 Lesson 03 学过的 unpacking：
seq_id, seq = line.split(”\t”)
于是：
seq_id → ”seq1”seq → ”ATGCGT”
到这里你应该会发现：
前面那些看起来零散的 Python 基础，现在开始真正组合起来了。
string+split+list+unpacking+for+function+file IO
3.3 文件中的数字首先也是 string
这是 R 用户尤其值得注意的一点。
假设文件：
seq1 ATGCGT 6seq2 GGCAAT 6
读取：
seq_id, seq, length = line.strip().split(”\t”)
此时length不是6：这个integer，而是："6"。也就是type(length)得到str。
如果后面需要数学计算：
length = int(length)
这正好回到 Lesson 01：
文件本身只保存文本 ↓ Python 读取成 str ↓ 根据数据含义做类型转换
以后读：
VCF BED GFF
时你会不断遇到这个过程。
4. 写文件：分析结果怎样保存下来？
读取：f.read() 写入：f.write()
例如：
with open(”output.txt”, ”w”) as f: f.write(”ATGCGT”)
文件就会包含：
4.1 write()不会自动加换行
假设：
with open(”output.txt”, ”w”) as f: f.write(”seq1”) f.write(”seq2”)
结果是：
如果需要换行：
with open(”output.txt”, ”w”) as f: f.write(”seq1\n”) f.write(”seq2\n”)
得到：
这一点非常简单，但也是文件输出最常见的小错误之一。
4.2 输出 tab-separated table
例如：
seq_id = ”seq1”seq = ”ATGCGT”length = len(seq)
可以写：
with open(”output.tsv”, ”w”) as f: f.write(f”{seq_id}\t{seq}\t{length}\n”)
结果：
这里第一次真正把前面学过的：
连接起来。
4.3print(..., file=f) 也可以写文件
Python 还有一种非常方便的写法：
with open(”output.txt”, ”w”) as f: print(”seq1”, file=f) print(”seq2”, file=f)
得到：
因为 print() 默认自动换行。
甚至：
with open(”output.tsv”, ”w”) as f: print(seq_id, seq, length, sep=”\t”, file=f)
也可以生成 tab-separated 行。
所以：
和：
都很常见。
初期你可以这样理解：
write→ 对具体写入的字符串控制更直接print(..., file=f)→ 输出人类可读的文本非常方便
5. pathlib：不要靠字符串拼文件路径
到目前为止我们一直写：
但真实项目可能是：
project/├── data/│ ├── raw/│ │ └── reads.fastq│ └── genome.fa├── results/└── scripts/
过去常见的写法可能是：
filename = ”reads.fastq”path = ”data/raw” + filename
现代 Python 更推荐：
然后：
path = Path(”data”) / ”raw” / ”reads.fastq”print(path)
这也是这节课非常值得养成的习惯。
5.1 Path 表示一个路径
from pathlib import Pathpath = Path(”data/genome.fa”)type(path)
现在 path 不是普通 string。
它表示：
文件系统中的一个路径。
可以：
非常适合科研脚本。
5.2 / 在 Path 中表示拼接路径
这是第一次看到/不是除法。
例如：
data_dir = Path(”data”)path = data_dir / ”genome.fa”print(path)
如果：
可以：
path = Path(”results”) / f”{sample}.txt”path
以后批量处理几十个 sample 时非常自然。
5.3 判断文件是否存在
返回True或False。
判断它是否是文件:
判断是否为目录。
from pathlib import Pathpath = Path(”data/genome.fa”)if not path.exists(): raise FileNotFoundError(f”File not found: {path}”)
5.4 当前工作目录：一个非常重要的概念
运行：
可以查看：
current working directory
也就是：
Python 当前把哪里当作“当前位置”。
假设：
是：
那么：
实际指向：
/home/user/project/data/genome.fa
所以 relative path：
本质上依赖：
current working directory
这也是为什么初学 Python 时经常会遇到：
明明“文件就在那个文件夹里”。
问题往往不是文件不存在，而是：
Python 当前工作的目录不是你以为的目录。
遇到这种情况，第一件事应该是：
from pathlib import Pathprint(Path.cwd())
确认自己究竟在哪里。
这个习惯以后非常有用。
6. 把前七课全部串起来：写一个真正的 sequence 文件处理程序
现在假设：
内容：
seq1 ATGCGTseq2 GGCCAATTseq3 ATATAT
我们的目标是：
- 得到 sequence ID 和 sequence
先复用 Lesson 07：
def gc_content(seq): if len(seq) == 0: raise ValueError(”Sequence is empty”) gc = seq.count(”G”) + seq.count(”C”)
然后：
from pathlib import Pathdef analyze_sequences(input_path, output_path): input_path = Path(input_path) output_path = Path(output_path) with open(input_path, ”r”) as infile, open(output_path, ”w”) as outfile: outfile.write(”seq_id\tlength\tgc_content\n”) for line in infile: line = line.strip() if line == ””: continue seq_id, seq = line.split(”\t”) length = len(seq) gc = gc_content(seq) outfile.write( f”{seq_id}\t{length}\t{gc}\n” )
调用：
analyze_sequences( ”sequences.tsv”, ”sequence_summary.tsv”)
得到：
seq_id length gc_contentseq1 6 0.5seq2 8 0.5seq3 6 0.0
这段代码值得认真看，因为它已经不是“Python 语法练习”了。
它拥有一个真正科研脚本的雏形：
input file ↓read ↓parse ↓function ↓analysis ↓output file
而且几乎所有东西你都已经学过。
String
line.strip()line.split(”\t”)
Tuple unpacking
Condition
Loop
Function
Exception
f-string
f”{seq_id}\t{length}\t{gc}\n”
File IO
Path
这就是为什么前几课一直强调：
不需要孤立地背方法。
真正重要的是知道这些工具什么时候组合起来。
科研文件不一定永远完美。
比如：
seq1 ATGCGTseq2 GGCCAATTbad_lineseq3 ATATAT
如果直接：
seq_id, seq = line.split(”\t”)
处理bad_line时就会报错。
我们现在已经有能力主动给出更清楚的信息：
def read_sequences(path): sequences = {} with open(path) as f: for line_number, line in enumerate(f, start=1): line = line.strip() if line == ””: continue fields = line.split(”\t”) if len(fields) != 2: raise ValueError( f”Invalid format at line {line_number}” ) seq_id, seq = fields sequences[seq_id] = seq return sequences
这里有一个很值得注意的新组合：
for line_number, line in enumerate(f, start=1):
Lesson 05 的 enumerate() 现在不再是在 list 上做练习，而是在真正的文件上发挥作用。
如果第 37 行坏掉：
ValueError: Invalid format at line 37
显然比：
ValueError: not enough values to unpack
更适合科研程序。
因为你马上知道应该去检查输入文件的哪一行。
你可以暂时这样建立迁移关系：
| |
|---|
readLines() | f.read() |
writeLines() | f.write() |
strsplit() | .split() |
trimws() | .strip() |
file.exists() | Path.exists() |
getwd() | Path.cwd() |
file.path() | Path(...) / ... |
但有一个区别很重要。
R 用户很容易直接想到：
read.delim()read.table()readr::read_tsv()
然后整个文件变成 data frame。
Python 后面当然也有：
但我们现在故意先不使用 Pandas。
因为对于：
理解：
本身就是很重要的能力。
尤其你后面需要写 genome assembly 和 sequence-processing code，这种思维比“会调用 read_csv()”重要得多。
假设reads.fastq有40GB。不要默认：
text = f.read()lines = f.readlines()
更应该想到：
with open(path) as f: for line in f: ...
也就是data streaming的思维。
下一课我们就会进一步把它升级为：
def read_fasta(path): ... yield ...
也就是：
file iteration+parser+generator
从而可以做到：
for seq_id, seq in read_fasta(”genome.fa”): ...
即使文件非常大，也不必一次全部加载到内存。
不过generator 是 Lesson 09的内容，这一课先把：
pathopenfile objectlinestripsplitwrite
真正掌握好。
Exercises
Exercise A｜预测 read() 的结果
假设文件：seq.txt
内容：
ATGCGT GGCAAT
代码：
with open(”seq.txt”) as f: x = f.read()print(repr(x))
假设文件最后也有一个正常换行符。
问：
x 是什么？type(x) 是什么？ len(x) 是多少？
x是：ATGCGT\nGGCAAT\n。type(x)是str，len(x)是14。
Exercise B｜逐行读取并过滤空行
假设：read.txt。
内容：
ATGCGT GGCAAT ATNNGC TTAACC
要求读取所有非空行，并且过滤掉含 "N" 的 sequence。
最后得到：
[ ”ATGCGT”, ”GGCAAT”, ”TTAACC”]
def read_reads(infile): sequences = [] with open(infile, ”r”) as f: for line in f: read = line.strip() if read == ””: continue if ”N” in read: continue sequences.append(read) return sequencesread_reads(”reads.txt”)
reads = []with open(”reads.txt”) as f: for line in f: read = line.strip() if read == ””: continue if ”N” in read: continue reads.append(read)
reads = []with open(”reads.txt”) as f: for line in f: read = line.strip() if read != ”” and ”N” not in read: reads.append(read)
两种都可以。
第一种对复杂数据清洗通常更容易读。
Exercise C｜用 Path 构建 sample 路径
给定：sample = "sample01"
希望建立：project/data/sample01/reads.fastq
要求使用 pathlib
然后分别得到：
文件名 suffix parent directory
sample = ”sample01”from pathlib import Pathfile = ( Path(”project”) / ”data” / sample / ”reads.fastq”)print(file.name)print(file.suffix)print(file.parent)
Exercise D｜解析一个 TSV 文件
假设：sequences.tsv
内容：
seq1 ATGCGTseq2 GGGAAAseq3 ATATGC
写：def read_sequences(path):
要求返回：
{ ”seq1”: ”ATGCGT”, ”seq2”: ”GGGAAA”, ”seq3”: ”ATATGC”}
def read_sequences(path): sequences = {} with open(path, ”r”) as f: for line in f: line = line.strip() if line == ””: continue seq_id, seq = line.split(”\t”) sequences[seq_id] = seq return sequencesread_sequences(”sequences.tsv”)
Exercise E｜给 parser 加上错误检查
继续使用：read_sequences()
现在要求：
- 每一行必须恰好有两个 tab-separated fields
def read_sequences(path): sequences = {} valid_bases = set(”ATCGN”) with open(path, ”r”) as f: for line_number, line in enumerate(f, start=1): line = line.strip() if line == ””: continue fields = line.split(”\t”) if len(fields) != 2: raise ValueError(f”Invalid format at {line_number}”) seq_id, seq = fields if not set(seq) <= valid_bases: raise ValueError(f”Invalid sequence at {line_number}”) sequences[seq_id] = seq return sequencesread_sequences(”sequences.tsv”)
Exercise F｜综合题：生成 sequence summary
这是本课最值得完整写一遍的题。
输入：sequences.tsv
seq1 ATGCGTseq2 GGGAAAseq3 ATATGCseq4 CCCGGG
写：def summarize_sequences(input_path, output_path):
生成summary.tsv
要求内容：
seq_id length gc_contentseq1 6 0.5seq2 6 0.5seq3 6 0.3333333333333333seq4 6 1.0
要求复用：gc_content()
而不是在主函数中重新写一遍 GC 计算逻辑。
from pathlib import Pathdef gc_content(seq): if len(seq) == 0: raise ValueError(”Sequence is empty”) gc = seq.count(”G”) + seq.count(”C”) return gc / len(seq)def summarize_sequences(input_path, output_path): input_path = Path(input_path) output_path = Path(output_path) if not input_path.exists(): raise FileNotFoundError(f”Input file not found: {input_path}”) with open(input_path,”r”) as infile, open(output_path,”w”) as outfile: outfile.write(”seq_id\tlength\tgc_content\n”) for line_number, line in enumerate(infile, start=1): line = line.strip() if line == ””: continue fields = line.split(”\t”) if len(fields) != 2: raise ValueError(f”Invalid format at {line_number}”) seq_id, seq = fields length = len(seq) gc = gc_content(seq) outfile.write(f”{seq_id}\t{length}\t{gc}\n”)summarize_sequences(”sequences.tsv”, ”summary.tsv”)
10个月宝宝每天需要喝多少奶粉？
