前四课我们一直在建立 Python 的“对象和数据结构”：
str → 一条序列list → 一批对象tuple → 一组固定数据dict → key → valueset → unique elements
到现在为止，我们大多是在手动操作一个对象。
这一课加入两个能力：
if → 根据条件决定做什么for → 对一批对象重复做什么
1. if：根据条件决定是否执行代码
1.1 if / elif / else
例如根据 GC content 做简单判断：
gc_content = 0.62if gc_content > 0.6: print(”High GC”)
基本结构是：
这里最需要注意的是两个东西：
Python 不使用 R 的 {}：
if (gc_content > 0.6) { print(”High GC”)}
而是使用：
if gc_content > 0.6: print(”High GC”)
缩进本身就是代码结构的一部分。
gc_content = 0.52if gc_content > 0.6: print(”High GC”)elif gc_content >= 0.4: print(”Medium GC”)else: print(”Low GC”)
1.2 在条件中组合多个判断
seq = ”ATGCGTAA”if seq.startswith(”ATG”) and seq.endswith(”TAA”): print(”Stars with ATG and ends with TAA”)
如果只需要满足其中一个条件：
stop_codon = ”TAG”if stop_codon == ”TAA” or stop_codon == ”TAG” or stop_codon == ”TGA”: print(”Stop codon”)
后面学到 set 之后，其实有一个更清晰的写法：
stop_codons = {”TAA”, ”TAG”, ”TGA”}if stop_codon in stop_codons: print(”Stop codon”)
1.3 in 与 not in 在条件判断中的使用
检查 motif：
seq = ”ATGCGTACGT”if ”CGT” in seq: print(”Motif found”)
检查 DNA 是否有 ambiguous base：
if ”N” in seq: print(”Sequence contains N”)
if ”N” not in seq: print(”No ambiguous base”)
2. for：逐个处理 sequence 中的元素
for 是这一课最重要的部分。
最简单的例子：
seq = ”ATGC”for base in seq: print(base)
2.1 遍历 string 和 list
String:
seq = ”ATGC”for base in seq: print(base)
List:
genes = [”TP53”, ”EGFR”, ”KRAS”]for gene in genes: print(gene)
再比如一批 reads：
reads = [ ”ATGCGT”, ”GGCTAA”, ”ATCCGT”]for read in reads: print(read)
现在我们第一次可以真正表达：
对每一条 read 做某件事。
例如计算每条 read 的长度：
for read in reads: print(len(read))
for read in reads: gc_count = read.count(”C”) + read.count(”G”) gc_content = gc_count / len(read) print(f”{read}: {gc_content:.2%}”)
到这里，你已经开始写真正的批量 sequence-processing code 了。
2.2 用循环累积结果
例如，我们不用 .count()，自己统计 G 的数量：
seq = ”ATGCGGT”g_count = 0for base in seq: if base == ”G”: g_count = g_count + 1print(g_count)
你以后会大量看到这种模式：
count = 0for something in data: if condition: count += 1
其中：
就是：
的简写。
类似：
等价于：
2.3 用 dictionary 做碱基计数
现在可以把 Lesson 04 的 dictionary 加进来。
先建立：
counts = { ”A”: 0, ”C”: 0, ”G”: 0, ”T”: 0}seq = ”ATGCGGTAA”for base in seq: counts[base] += 1print(counts)
3. 如何遍历 Dictionary
Dictionary 和 list 不太一样，因为它同时有key和value, 所以有几种不同的遍历方式。
3.1 直接遍历 dictionary：得到 key
sequences = { ”seq1”: ”ATGCGT”, ”seq2”: ”GGGAAA”, ”seq3”: ”ATCCGT”}for seq_id in sequences: print(seq_id) print(sequences[seq_id])
3.2 .items()：同时得到 key 和 value
items() 用来在循环中同时取得 dictionary 的 key 和 value。
这是以后最常用的 dictionary 遍历方式之一：
for seq_id, seq in sequences.items(): print(seq_id, seq)
你应该能看到上一课 tuple unpacking 的影子：
seq_id, seq = (”seq1”, ”ATGCGT”)
因为：
本质上提供的是一组：
for seq_id, seq in sequences.items(): gc_count = seq.count(”G”) + seq.count(”C”) gc_content = gc_count / len(seq) print(f”{seq_id}: {gc_content:.2%}”)
这已经非常接近以后处理 FASTA的代码形式
3.3 遍历 set
bases = {”A”, ”C”, ”G”, ”T”}for base in bases: print(base)
可以逐个取得四种 base。
但是不要写出这种逻辑：
“set 的第一个一定是 A，第二个一定是 C。”
因为 set 的目的本来就不是保存 position。
所以：
list → position 很重要set → membership 更重要
4. range()、enumerate() 与 zip()：需要位置时怎么循环？
普通：
只能直接得到base
但 sequence algorithms 很多时候还需要：
这时会用到这一组工具
4.1 range()：生成整数序列
range() 用来产生一系列整数，最常用于按 index 循环。
for i in range(5): print(i)
range() 可以指定起点和终点：
for i in range(2, 6): print(i)
也可以加入 step：
for i in range(0,10,2): print(i)
所以基本结构是：
类似于：
seq = ”ATGC”for i in range(len(seq)): print(i, seq[i])
4.2 enumerate()：同时得到 index 和元素
刚才我们写：
for i in range(len(seq)): print(i, seq[i])
Python 更常写成：
for i, base in enumerate(seq): print(i, base)
for position, base in enumerate(seq, start = 1): print(position, base)
这里有一个非常重要的生信应用。
Python index：
Rosalind 或某些 biological position：
于是：
有时可以让输出位置更符合题目要求。
但要记住：
改的是循环产生的编号，不是 Python string 自己的 indexing 规则。
seq[0] 仍然永远是第一个字符。
4.3 zip()：同时遍历多个对象
genes = [”TP53”, ”EGFR”, ”KRAS”]chroms = [”chr17”, ”chr7”, ”chr12”]for gene, chrom in zip(genes, chroms): print(gene, chrom)
5. for 在序列算法中的真正用途
现在进入最重要的应用部分。
5.1 找所有 motif，而不是只找第一处
Lesson 02 我们学习了：
它只能直接告诉我们第一处。
现在我们可以自己检查所有位置。
seq = ”GATATATGCATATACTT”motif = ”ATAT”for i in range(len(seq) - len(motif) + 1): if seq[i:i+len(motif)] == motif: print(i)
为什么是：
len(seq) - len(motif) + 1
假设：
sequence length = 10motif length = 3
合法的起始 index 是：
总共：
个位置。
因此：
恰好生成：
如果继续到：
剩下的 sequence 已经不够一个完整的 3-mer。
这条公式以后会不断回来：
长度为 n 的 sequence长度为 k 的 window共有：n - k + 1个完整 window
如果题目要求 1-based position：
for i in range(len(seq) - len(motif) + 1): if seq[i:i+len(motif)] == motif: print(i + 1)
5.2 生成所有 k-mer
现在 genome assembly 的一个核心操作已经可以正式写出来。
seq = ”ATGCGT”k = 3for i in range(len(seq) - k +1): kmer = seq[i:i+k] print(kmer)
如果不只是打印，而是希望以后继续使用，可以先建立空 list：
kmers = []for i in range(len(seq) - k + 1): kmer = seq[i:i+k] kmers.append(kmer)print(kmers)
所以非常经典的模式是：
results = []for item in data: result = ... results.append(result)
6. while 与循环控制
for 的含义通常是：
对这一批东西逐个处理。
while 的思维不同。
count = 0while count < 5: print(count) count += 1
为什么 while 要特别小心？
如果我们忘记count += 1
那么count永远是0，条件一直为true，于是形成infinite loop，也就是死循环。
所以写 while 时必须问自己：
什么东西最终会让这个条件变成 False？
这是一条非常好的检查原则。
break 与 continue
这两个属于“改变循环正常执行流程”，所以放在一起理解。
break：立刻结束整个循环
例如找到第一个 N 就停止：
seq = ”ATGCNTA”for base in seq: if base == ”N”: print(”N found”) break
一旦碰到N
后面的字符不再继续检查。
适合：
找到第一个满足条件的对象后就已经完成任务。
continue：跳过当前这一轮
seq = ”ATGCNTA”for base in seq: if base==”N”: continue print(base)
所以：
break→ 整个循环停止continue→ 只跳过当前这一轮
7. 一个综合例子：DNA sequence QC
给定：
sequences = { ”seq1”: ”ATGCGT”, ”seq2”: ”ATGCNT”, ”seq3”: ”GGCCAA”, ”seq4”: ””}
我们想检查：
可以写：
for seq_id, seq in sequences.items(): if len(seq) == 0: print(f”{seq_id}: empty seuence”) elif ”N” in seq: print(f”{seq_id}: contains N”) else: print(f”{seq_id}: valid”)
8. 本课练习
这次重点不是写很多题，而是把 if + for + 数据结构真正连起来。
Exercise A｜检查 DNA sequence
给定：
seq = "ATGCNTAG"
逐个检查每个 base。
如果发现：
A / C / G / T
以外的字符，打印：
Invalid base: N
提示：
valid_bases = {"A", "C", "G", "T"}
思考如何组合：
for if not in
seq = ”ATGCNTAG”valid_bases = {”A”, ”T”, ”G”, ”C”}for base in seq: if base not in valid_bases: print(”Invalid base:”, base) print(”Vaild bases = ”, valid_bases)
Exercise B｜自己实现 nucleotide counting
不要使用：
seq.count()
给定：
seq = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCT"
先建立：
counts = { "A": 0, "C": 0, "G": 0, "T": 0 }
然后只用：
for dictionary
完成计数。
最后输出：
A C G T
的数量。
这题非常重要。
seq = ”AGCTTTTCATTCTGACTGCAACGGGCAATATGTCT”counts = { ”A” : 0, ”T” : 0, ”C” : 0, ”G” : 0}for base in seq: if base in counts: counts[base] += 1print(counts)
Exercise C｜所有 motif positions
给定：
seq = "GATATATGCATATACTT" motif = "ATAT"
找出 所有 motif 出现的位置。
要求输出：
1-based positions
seq = ”GATATATGCATATACTT”motif = ”ATAT”for i in range(len(seq) - len(motif) + 1): if seq[i:i + len(motif)] == motif: print(i+1)
Exercise D｜生成 k-mer
给定：
seq = "ATGCGTAC" k = 4
生成所有 4-mer，并保存进：
kmers = []
seq = ”ATGCGTAC”k = 4kmers = []for i in range(len(seq) - k + 1): kmer = seq[i:i+k] print(kmer) kmers.append(kmer)
Exercise E｜小型 genomics QC
给定：
reads = [ ”ATGCGT”, ”NNNNNN”, ”ATGC”, ””, ”GGCCAATT”]
依次处理每条 read：
· 如果是空字符串，打印 "empty"
· 如果含 N，打印 "contains N"
· 否则，如果长度小于 5，打印 "too short"
· 其他情况打印 "pass"
reads = [ ”ATGCGT”, ”NNNNNN”, ”ATGC”, ””, ”GGCCAATT”]for read in reads: if read == ””: print(”empty”) elif ”N” in read: print(”contains N”) elif len(read) < 5: print(”too short”) else: print(”pass”)
