上一课我们已经有了三类 sequence:
str → ordered + immutablelist → ordered + mutabletuple → ordered + immutable
这一课开始学习两种思维完全不同的数据结构:
dict → key → valueset → unique elements
它们的重点已经不再是:
“第几个元素是什么?”
而开始变成:
“根据一个 key 找到对应的数据”
以及:
“某个元素是否存在?”
这两个问题在生物信息学中非常常见。
例如:
sequences = { ”seq1”: ”ATGCGT”, ”seq2”: ”GGCTAA”, ”seq3”: ”ATCCGT”}
我们真正想表达的不是三个序列排在第 1、2、3 位,而是:
seq1 → ATGCGTseq2 → GGCTAAseq3 → ATCCGT
这就是 dictionary 的核心。
1. Dictionary:用 key 找 value
1.1 创建 dictionary 与读取数据
例如保存 gene 和 chromosome:
gene_chrom = { ”TP53”: ”chr17”, ”EGFR”: ”chr7”, ”KRAS”: ”chr12”}print(gene_chrom)type(gene_chrom)
gene_chrom = { ”TP53”: ”chr17”, ”EGFR”: ”chr7”, ”KRAS”: ”chr12”}gene_chrom[”TP53”]
1.2 Dictionary 的 value 可以是什么?
例如保存序列:
sequences = { ”seq1”: ”ATGCGT”, ”seq2”: ”GGCTAA”, ”seq3”: ”ATCCGT”}
保存 read counts:
read_counts = { ”sample01”: 12500000, ”sample02”: 9800000, ”sample03”: 14300000}
也可以:
gene_info = { ”gene”: ”TP53”, ”chrom”: ”chr17”, ”strand”: ”-”, ”exons”: [1, 2, 3, 4, 5]}print(gene_info[”gene”])print(gene_info[”exons”])
1.3 Dictionary 和 R 的对应关系
你可以暂时把:
gene_info = { ”gene”: ”TP53”, ”chrom”: ”chr17”, ”strand”: ”-”}
类比成 R 的 named list:
gene_info <- list( gene = ”TP53”, chrom = ”chr17”, strand = ”-”)
但 Python dictionary 更应该建立成:
key-value mapping / hash table
这个概念。
尤其以后出现:
或者:
时,你不应该再把它理解成“像 R list 一样取东西”。
而应该直接想到:
2. 添加、修改与删除 dictionary 内容
Dictionary 和上一课的 list 一样,是mutable
所以我们可以直接修改原来的 dictionary。
2.1 添加和修改 key-value
例如:
gene_chrom = { ”TP53”: ”chr17”, ”EGFR”: ”chr7”}gene_chrom[”KRAS”] = ”chr12”print(gene_chrom)
gene_chrom = { ”TP53”: ”chr17”, ”EGFR”: ”chr7”}print(gene_chrom)gene_chrom[”TP53”] = ”chromosome 17”print(gene_chrom)
所以同一句语法:
可能表示两件事情:
这非常重要
2.2 删除:pop() 与 del
例如:
gene_chrom = { ”TP53”: ”chr17”, ”EGFR”: ”chr7”, ”KRAS”: ”chr12”}chrom = gene_chrom.pop(”KRAS”)print(chrom)print(gene_chrom)
例如:
删除之后就不存在 "EGFR" 这个 key 了。
初期不用纠结什么时候一定要 pop()、什么时候一定要 del。
先记住:
pop() → 删除,并取得被删除的 valuedel → 直接删除
3. 查找 dictionary:in、get() 与 keys/values/items
这一部分是 dictionary 真正高频的操作。
3.1 in:判断 key 是否存在
例如:
gene_chrom = { ”TP53”: ”chr17”, ”EGFR”: ”chr7”}”TP53” in gene_chrom
这里有一个特别容易忽略的点:
对 dictionary 使用 in,默认检查的是 key。
例如:
虽然 "chr17" 明明是 dictionary 里的 value。
因为 Python 实际问的是:
有没有一个 key 叫 "chr17"?
没有。
所以:
和:
回答的是完全不同的问题。
3.2 get():安全地取得 value
例如:
得到:
相比之下:
如果 "BRCA1" 不存在,会产生:
所以:
和:
虽然都可以取 value,但含义稍有不同:
| | |
|---|
d[key] | | KeyError |
d.get(key) | | None |
还可以嵌套
gene_chrom = { ”TP53” : { ”chr7”: 2245 }}gene_chrom.get(”TP53”).get(”chr7”)
gene_chrom.get(”BRCA1”, ”unknown”)
这个功能以后做 counting 时会非常有用。
例如你以后会看到类似:
counts[base] = counts.get(base, 0) + 1
现在不要求你完全理解这句代码。
等学完 for loop 后我们会自己写出来。
3.3 keys()、values() 和 items()
这三个方法放在一起学最合理,因为它们只是从 dictionary 的三个角度看数据。
gene_chrom = { ”TP53”: ”chr17”, ”EGFR”: ”chr7”}gene_chrom.keys()
注意这里每一对:
下一课学习 for 后,你会经常看到:
for gene, chrom in gene_chrom.items(): ...
4. Dictionary 在生物信息学中的典型用途
4.1 Sequence ID → sequence
FASTA:
以后我们自己解析时,很自然可以得到:
sequences = { ”seq1”: ”ATGCGT”, ”seq2”: ”GGCTAA”}sequences[”seq1”]
4.2 碱基 → count
counts = { ”A”: 10, ”C”: 8, ”G”: 12, ”T”: 9}counts[”G”]
4.3 Node → neighbors
再往后进入 graph 时:
graph = { ”ATG”: [”TGC”], ”TGC”: [”GCG”, ”GCA”], ”GCG”: [”CGT”]}graph[”TGC”]
其生物信息学含义可以理解成:
于是:
dictionarykey → nodevalue → neighbors
这就是为什么 dictionary 最后可以非常自然地表示 adjacency list。
原课程大纲里强调 dictionary 对 assembly 很重要,核心正是这种 graph[node] = neighbors 的结构。
5. Set:只关心“有哪些”,不关心重复
Set 解决的则是另外一个问题:
有哪些不同的元素?
bases = {”A”, ”C”, ”G”, ”T”}type(bases)
Set 最重要的三个特点可以先记成:
无重复没有 index特别适合 membership testing
5.1 创建 set 与自动去重
例如:
kmers = {”ATG”, ”TGC”, ”ATG”, ”GCG”}kmers
set 实际上回答:
这条序列中出现过哪些不同字符?
而不是:
每个字符出现了几次?
Set 不表示元素的固定位置,因此不能像 list 一样:
这种写法会报错。
6. Set 的添加、删除与集合运算
6.1 add():加入一个元素
add() 用来向 set 中加入一个元素。
base = {”A”, ”C”, ”G”}base.add('T')base
base = {”A”, ”C”, ”G”}base.add('T')base.add('T')base
6.2 remove() 与 discard():删除元素
这两个方法功能非常接近,所以应该放在一起比较。
bases = {”A”, ”T”, ”G”, ”C”}bases.remove(”T”)bases
bases = {”A”, ”T”, ”G”, ”C”}bases.remove(”T”)bases.remove(”T”)
如果重复执行会产生:KeyError
bases = {”A”, ”T”, ”G”, ”C”}bases.discard(”T”)bases.discard(”T”)bases
6.3 Set 的交集、并集和差集
假设两个样本检测到的 genes:
sample_a = {”TP53”, ”EGFR”, ”KRAS”}sample_b = {”TP53”, ”BRCA1”, ”KRAS”}
这部分如果你有 R 的集合操作经验应该会很熟悉:
| | |
|---|
| intersect(a, b) | a & b |
| union(a, b) | a | b |
| setdiff(a, b) | a - b |
7. Dictionary 与 Set 的共同核心:快速判断“有没有”
现在可以把本课两个对象真正串起来。
对于 list:
genes = [”TP53”, ”EGFR”, ”KRAS”]
我们经常关心:
第几个元素是什么?
对于 dictionary:
gene_chrom = { ”TP53”: ”chr17”, ”EGFR”: ”chr7”}
我们关心:
TP53 这个 key 有没有?如果有,它对应什么?
对于 set:
genes = {”TP53”, ”EGFR”, ”KRAS”}
我们只关心:
TP53 有没有?
这一点和上一课的 list 连起来。
gene_chrom = { ”TP53”: ”chr17”}gene_chrom[”EGFR”] = ”chr7”gene_chrom
8. 本课练习
这次练习仍然只保留真正重要的几组。
Exercise A|FASTA-like dictionary
建立:sequences = { ”seq1”: ”ATGCGT”, ”seq2”: ”GGGAAA”, ”seq3”: ”ATCCGT”}请完成:1. 取得 seq2 的 sequence。2. 判断 ”seq4” 是否存在。3. 加入: seq4 → TTGGCA4. 把 seq1 修改成: ATGCGTAA5. 使用 len() 判断现在有多少条 sequence。
sequences = { ”seq1”: ”ATGCGT”, ”seq2”: ”GGGAAA”, ”seq3”: ”ATCCGT”}print(sequences[”seq2”])print(”seq4” in sequences)sequences[”seq4”] = ”TTGGCA”print(sequences)sequences[”seq1”] = ”ATGCGTAA”print(sequences)len(sequences)
Exercise B|[] 还是 .get()?
不要先运行,判断下面两句的区别:
sequences["seq10"]
sequences.get("seq10")
假设 "seq10" 不存在:
两句分别会发生什么?
进一步思考:
如果“找不到 seq10”本身就意味着数据有问题,你更希望用哪一种?
这已经开始涉及科研脚本中的设计选择。
第一个产生KeyError,第二个产生None。
如果“找不到 seq10”本身就意味着数据有问题,更希望用第一种,报告错误信息。
相反,如果需要考虑seq10不存在的情况,就用第二种,或者用get指定返回0
sequences.get("seq10", 0)
Exercise C|Set 与 DNA QC
给定:
seq1 = "ATGCGTAC" seq2 = "ATGCNTAC"
分别运行:
set(seq1) set(seq2)
思考:
seq1 中有哪些不同字符? seq2 中有哪些不同字符?
然后判断:
"N" in set(seq1) "N" in set(seq2)
结果分别是什么?
实际上这里还可以直接:
"N" in seq
所以再思考一个更重要的问题:
如果只是判断有没有 N,为什么其实完全不需要先 set(seq)?
不要为了“用了高级数据结构”而使用它,这是写科研代码时很重要的习惯。
seq1 = ”ATGCGTAC”seq2 = ”ATGCNTAC”print(set(seq1))print(set(seq2))print(”N” in set(seq1) )print(”N” in set(seq2) )
判断是否存在N,可以直接在str层面检测
Exercise D|两个 k-mer 集合
给定:
sample1 = {"ATG", "TGC", "GCG", "CGT"} sample2 = {"ATG", "TGC", "GCA", "CAT"}
请得到:
两个样本共同的 k-mer
两个样本所有不同的 k-mer
sample1 特有的 k-mer
sample2 特有的 k-mer
sample1 = {”ATG”, ”TGC”, ”GCG”, ”CGT”}sample2 = {”ATG”, ”TGC”, ”GCA”, ”CAT”}print(sample1 & sample2)print(sample1 | sample2)print(sample1 - sample2)print(sample2 - sample1)
最后,把前三课的数据结构连起来:
DNA sequence ↓ stringmany sequences ↓ listsequence ID → sequence ↓ dictionaryunique k-mers ↓ setnode → neighbors ↓ dictionary