Python学习 | 02 String 与 DNA 序列

这一课的目标

学完以后，你应该能熟练理解并使用：

seq = ”ATGCGTACGT”seq[0]seq[-1]seq[2:6]seq[::-1]len(seq)seq.count(”G”)seq.find(”CG”)seq.replace(”T”, ”U”)seq.upper()seq.lower()”GC” in seq

并且能够自己完成：

nucleotide counting

transcription

reverse sequence

reverse complement

motif searching

GC content

1. String 到底是什么？

Python 中：

seq = "ATGCGTACGT"

seq 是一个：

str

检查：

type(seq)

得到：

str

你可以暂时把 string 理解成：

一串有顺序的字符。

例如：

A T G C G T A C G T

它们有明确顺序，因此 Python 可以通过“位置”访问其中每一个字符。

这就是：

indexing

2. Python 最重要的规则：从 0 开始

这是 R 用户必须尽快形成肌肉记忆的地方。

假设：

seq = ”ATGCGT”

字符位置是：

字符:      A   T   G   C   G   TPython:    0   1   2   3   4   5

所以：

seq[0]

seq[1]

3. R 与 Python indexing 对照

假设：

ATGCGT

R 的思维：

1 2 3 4 5 6

Python：

0 1 2 3 4 5

所以：

以后凡是看到：

seq[i]

请先在脑中提醒自己：

i 是 Python index，不是人类习惯的第 i 个。

4. 从末尾开始取：负数 indexing

-1表示最后一个字符

seq = ”ATGCGT”seq[-1]

seq[-2]

位置可以想象成

字符:       A   T   G   C   G   T正向:       0   1   2   3   4   5反向:      -6  -5  -4  -3  -2  -1

这个在序列处理中非常好用。

比如判断末尾是不是 stop codon：

seq[-3:]

5. slicing：一次取一段序列

Python slice：

[start:end]

包含 start，不包含 end

即左开右闭：

[start, end)

seq = ”ATGCGTACGT”seq[2:6]

6. 为什么 Python 要这么设计？

一开始 R 用户通常觉得“不包含末端很奇怪”。

但其实很有逻辑。

例如：

seq[2:6]

长度正好：

6 - 2

len(seq[2:6])

这让很多序列和坐标计算非常方便。

7. slicing 的生信意义非常大

以后你会频繁写：

seq[start:end]

比如提取 genomic region、read 的一部分、k-mer 等。

而且一个特别重要的联系是：

Python slicing 的 [start:end) 思维，与 BED 文件常用的 0-based half-open coordinate 非常接近。

所以学会 Python slice 后，以后理解 BED 坐标体系其实会自然很多。

8. slice 可以省略 start

seq = ”ATGCGTACGT”seq[:3]

这对于 codon 很直观：

start_codon = seq[:3]print(start_codon)

9. slice 可以省略 end

seq = ”ATGCGT”seq[3:]

10. 取最后三个字符

seq = ”ATGAAATGA”seq[-3:]

stop_codon = seq[-3:]print(stop_codon)

11. slice 还有第三个参数：step

完整形式：

seq[start:end:step]

例如：

seq[::2]

表示：

每隔一个字符取一个。

seq = ”ATGCGT”seq[::2]

12. 一个非常经典的写法：反转字符串

seq = ”ATGCGT”seq[::-1]

13. 字符串是 immutable 的

seq = ”ATGC”seq[0] = ”G”

python的str是immutable的，创建之后不能直接修改其中某一个位置。

14. 为什么不能修改？

字符串对象：

"ATGC"

创建以后本身不能被原地修改。

所以：

seq[0] = "G"

不允许。

如果你想得到：

GTGC

需要创建一个新的字符串。

seq = ”G” + seq[1:]print(seq)

15. 这和 Lesson 01 的对象概念连起来了

之前我们说：

seq = ”ATGC”

可以理解为：

seq ─────→ ”ATGC”

执行：

seq = ”GTGC”

之后：

seq ─────→ ”GTGC”

seq 改变了它绑定的对象。

不是原来的 "ATGC" 被修改了。

这个概念以后理解 list 的 mutable 行为会非常关键

16. 字符串拼接

seq1 = ”ATGC”seq2 = ”GGTA”seq = seq1 + seq2print(seq)

17. 字符串重复

”AT” * 3

”N” * 10

这在生成 padding 或 mock sequence 时偶尔很好用

18. len()：序列长度

seq = ”ATGCGT”len(seq)

19. .count()：统计字符出现次数

seq = ”ATGCGCGTAA”print(seq.count(”A”))print(seq.count(”GC”))

20. 第一个 Rosalind 风格程序：碱基计数

seq = ”ATGCGCGTAA”a = seq.count(”A”)c = seq.count(”C”)g = seq.count(”G”)t = seq.count(”T”)print(a, c, g, t)

注意：

print(a, c, g, t)

默认会用空格隔开多个对象。

21. .upper() 与 .lower()

seq = ”atGcGTaa”seq = seq.upper()print(seq)

seq = ”atGcGTaa”seq = seq.lower()print(seq)

22. .replace()：字符串替换

一个非常直接的生信应用：

DNA → RNA transcription

dna = ”ATGCGT”rna = dna.replace(”T”,”U”)print(rna)

23. 注意 .replace() 不会修改原字符串

dna = ”ATGCGT”dna.replace(”T”, ”U”)print(dna)

24. in：判断 motif 是否存在

seq = ”ATGCGTACGT””CGT” in seq

”AAA” in seq

25. not in

判断 sequence 是否没有 ambiguous base N。

seq = ”ATGCGT””N” not in seq

实际数据 QC 时非常实用。

26. .find()：找到 motif 的位

seq = ”ATGCGTACGT”seq.find(”CGT”)

注意返回的是：

Python 的 0-based index。

即第一个 "CGT" 从 index 3 开始。

27. 找不到时会返回 -1

seq.find(”AAA”)

position = seq.find(”AAA”)if position == -1:    print(”Motif not found”)

28. .index() 和 .find() 很像

seq.index(”CGT”)

区别在于：

如果找不到：

seq.find("AAA")返回-1，

而

seq.index("AAA")会报错

seq.index(”AAA”)

所以对于探索性代码使用find()

29. .startswith() 与 .endswith()

seq = ”ATGCGTAA”seq.startswith(”ATG”)

seq.endswith(”TAA”)

30. .strip()：以后读文件时超级重要

seq = ”ATGCGT\n”

len(seq)

seq = seq.strip()print(seq)len(seq)

以后读 FASTA / FASTQ 文件时，.strip() 会频繁出现。

31. 特殊字符：\n 与 \t

print(”ATGC\nGGTA”)

print(”gene\tlength”)

32. GC content

seq = ”ATGCGCGTAA”gc_count = seq.count(”G”) + seq.count(”C”)gc_content = gc_count / len(seq)print(f”GC content: {gc_content:.2%}”)

33. 注意空序列问题

seq = ””len(seq)

如果计算：

gc_count / len(seq)

会发生：

ZeroDivisionError

以后写科研级函数时就需要处理这种情况。

不过等我们学到 function + if 后再正式解决。

34. Reverse sequence ≠ Reverse complement

reverse sequence

seq = ”ATGC”seq[::-1]

reverse complement

GCAT

陷阱：

35. 怎么计算 complement？

seq = ”ATGC”seq[::-1].replace(”A”, ”T”).replace(”T”, ”A”).replace(”G”, ”C”).replace(”C”, ”G”)

错在哪里？

36. 一个非常值得理解的字符串陷阱

seq = ”AT”seq[::-1].replace(”A”, ”T”).replace(”T”, ”A”)

37. 方法一：用 translate()

table = str.maketrans(”ATCG”,”TAGC”)seq = ”ATGC”complement = seq.translate(table)print(complement)print(complement[::-1])

38. 为什么这里先给你 translate()？

虽然这个方法稍微超前一点，但它非常适合：

字符到字符的一一映射。

所以它比连续 .replace() 更安全、更清晰。

以后我们还会用 dictionary 写一次 reverse complement，届时你会更理解其内部逻辑。

39. Reverse complement 推荐写法

seq = ”AAAACCCGGT”table = str.maketrans(”ACGT”, ”TGCA”)revcomp = seq.translate(table)[::-1]print(revcomp)

40. 处理小写序列

seq = seq.upper()table = str.maketrans(”ACGT”, ”TGCA”)revcomp = seq.translate(table)[::-1]

41. 一个很实用的 sequence cleaning

实际序列可能有空格或换行：

seq = ” ATGCGT\n”seq = seq.strip().upper()

意思就是：

去掉头尾空白和换行

转成大写

以后从文本文件读取序列时，你会大量看到这种 chaining。

42. Method chaining

seq.strip().upper()叫做 method chaining。

因为：

seq.strip()

返回一个 string。

string 又可以继续：

.upper()

所以：

seq.strip().upper()

非常自然。

例如：

dna = ”atgc\n”rna = dna.strip().upper().replace(”T”,”U”)print(rna)

43. 一个 R → Python 对照

R 中很多字符串操作：

str_replace(...)str_detect(...)

Python 原生 string 则大量使用：

seq.replace(...)motif in seq

44. 一个重要思维差异：函数 vs 方法

Python：

len(seq)

这是：

function

而：

seq.upper()

这是：

method

看语法区别：

function:

function(object)

例如：

len(seq)

method：

object.method()

例如：

seq.upper()

后面你会大量遇到这种模式：

df.head()array.mean()path.exists()

理解“对象的方法”是 Python 面向对象语法的第一步。

45. 为什么是 len(seq) 而不是 seq.len()？

这是 Python 自己的设计。

你不需要现在纠结为什么。

先熟悉：

len(seq)

但：

seq.upper()seq.count("A")seq.replace("T", "U")

这些是字符串自己的 method。

46. String interpolation：f-string 再练一次

seq = ”ATGCGT”print(f”Sequence: {seq}”)print(f”Length: {len(seq)}”)print(f”GC count: {seq.count('G') + seq.count('C')}”)

所以科研代码推荐：

gc_count = seq.count(”G”) + seq.count(”C”)gc_content = gc_count / len(seq)print(f”GC content: {gc_content:.2%}”)

48. 判断 DNA 是否只含合法碱基

seq = ”ATGCGT””A” in seq

下一课学：

for

后我们就能自己检查每个 base。

再之后学：

set

会得到非常漂亮的写法：

set(seq) <= {”A”, ”C”, ”G”, ”T”}

49. Rosalind 题 1：DNA nucleotide counting

给定：

seq = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC"

要求输出：

A C G T

数量。

seq = ”AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC”a = seq.count(”A”)c = seq.count(”C”)g = seq.count(”G”)t = seq.count(”T”)print(a, c, g, t)

50. Rosalind 题 2：Transcription

dna = "GATGGAACTTGACTACGTAAATT"

dna = ”GATGGAACTTGACTACGTAAATT”rna = dna.replace(”T”,”U”)print(rna)

51. Rosalind 题 3：Reverse complement

dna = "AAAACCCGGT"

dna = ”AAAACCCGGT”table = str.maketrans(”ACGT”,”TGCA”)rev_comp = dna.translate(table)[::-1]print(rev_comp)

dna = ”AAAACCCGGT”rev = dna[::-1]revrev_comp = rev.replace(”A”,”t”).replace(”T”,”A”).replace(”t”,”T”)rev_comp = rev_comp.replace(”C”,”g”).replace(”G”,”C”).replace(”g”,”G”)rev_comp

dna = ”AAAACCCGGT”complement = {    "A": "T",    "T": "A",    "G": "C",    "C": "G",}rev_comp = ""for base in dna[::-1]:    rev_comp += complement[base]dna

52. Motif：先学最简单情况

seq = ”GATATATGCATATACTT”motif = ”ATAT”seq.find(motif) + 1

53. 但是 .find() 只能直接找到第一处

里面其实出现不止一次。

.find()：

seq.find(motif)

只会告诉你第一处。

如果我们要找：

所有 motif positions

就需要：

for

所以这正好是下一课要解决的问题。

54. slicing 和 k-mer 的关系

seq = ”ATGCGT”# k = 3print(seq[0:3])print(seq[1:4])print(seq[2:5])print(seq[3:6])# seq[i:i+k]

以后 genome assembly 中的核心操作：

kmer = seq[i:i+k]

55. 再看 de Bruijn graph 的一个伏笔

比如

kmer = ”ATGC”

prefix：

kmer[:-1]

suffix：

kmer[1:]

于是一个 k-mer：

ATGC

可以形成一条 edge：

ATG → TGC

这就是 de Bruijn graph 的最基本结构。

所以今天的：

[:-1][1:]

以后会直接出现在 genome assembly 代码里。

56. 请特别记住几个 slice

sqe = ”ATGCGT”

第一个字符

seq[0]

最后一个字符

seq[-1]

前三个

seq[:3]

后三个

seq[-3:]

去掉第一个

seq[1:]

去掉最后一个

seq[:-1]

反转

seq[::-1]

57. 常见错误 1：忘记 0-based

seq[1]# 是第二个

58. 常见错误 2：以为 slice 包含终点

seq[2:5]# 是第2个到第4个

59. 常见错误 3：直接修改字符串

seq[0] = ”A”# 不改变seq

60. 常见错误 4：忘记保存方法返回值

seq = ”atgc”seq.upper()print(seq)seq = seq.upper()print(seq)

61. 常见错误 5：文件中的换行符

seq.strip()

62. Lesson 02 速查表

63. 本课练习

这次建议你真正动手写。

Exercise 1：Indexing

seq = "ATGCCGTAGC"

请分别得到：

第一个碱基 第三个碱基 最后一个碱基 倒数第二个碱基

seq = ”ATGCCGTAGC”print(seq[0])print(seq[2])print(seq[-1])print(seq[-2])

Exercise 2：Slicing

还是：

seq = "ATGCCGTAGC"

请提取：

前 3 个碱基 后 3 个碱基 index 2 到 index 6 之前 去掉第一个碱基 去掉最后一个碱基

seq = ”ATGCCGTAGC”print(seq[:3])print(seq[-3:])print(seq[2:7])print(seq[1:])print(seq[:-1])

Exercise 3：Reverse

seq = "ATGCCGTA"

得到 reversed sequence。

seq = ”ATGCCGTA”seq[::-1]

Exercise 4：DNA → RNA

dna = "ATGCGTACCTTAG"

转为 RNA。

dna = ”ATGCGTACCTTAG”dna.replace(”T”,”U”)

Exercise 5：GC content

seq = "ATGCGCGTAACCGT"

输出：

Length: ... A: ... C: ... G: ... T: ... GC content: ...%

seq = ”ATGCGCGTAACCGT”seq_len = len(seq)a = seq.count(”A”)c = seq.count(”C”)g = seq.count(”G”)T = seq.count(”T”)gc_content = (c + g) / seq_lenprint(f”Length: {seq_len}”)print(f”A: {a}”)print(f”C: {c}”)print(f”G: {g}”)print(f”T: {t}”)print(f”GC content: {gc_content:.2%}”)

64. Exercise 6：Reverse complement

给定：

seq = "AGTCCTAG"

请用：

str.maketrans() translate() [::-1]

得到 reverse complement。

seq = ”AGTCCTAG”table = str.maketrans(”ATCG”,”TAGC”)rev_comp = seq.translate(table)[::-1]print(rev_comp)

65. Exercise 7：motif

seq = "AAGCTTAGCTAGCTA" motif = "GCTA"

请回答：

motif 是否存在？

第一处 Python index 是多少？

第一处按 1-based position 是多少？

seq = ”AAGCTTAGCTAGCTA”motif = ”GCTA”

motif in seq

seq.find(motif)

seq.find(motif) + 1

66. Challenge：用 slicing 理解 codon

给定：

cds = "ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG"

请取出：

第一个 codon 第二个 codon 最后一个 codon

cds = ”ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG”print(cds[0:3])print(cds[3:6])print(cds[-3:])

67. Genome assembly 预热题

给定：

kmer = "ATGCG"

请分别得到：

prefix = ATGC suffix = TGCG

kmer = ”ATGCG”prefix = kmer[:-1]suffix = kmer[1:]print(prefix)print(suffix)

如果这一题你能非常自然地写出来，那么以后看到：

graph[kmer[:-1]].append(kmer[1:])

就不会觉得神秘。