目录
4.1序列数据的表征
4.2 存储 DNA 序列的程序
4.3 连接DNA片段
4.4 转录:从 DNA 到 RNA
4.5 使用Python文档
4.6 在Python中计算反向互补
4.7 蛋白质、文件和数组
4.8 从文件中读取蛋白质序列数据
4.9 列表
4.10 练习题
在本章中,我们将开始编写Python程序,来处理DNA和蛋白质的生物序列数据。在你的计算机上,一旦有了这样的序列,你就可以编写程序对序列数据进行下列的处理了:
- •把DNA转录成RNA
- •把序列连接起来
- •获取反向互补序列
- •从文件中读取序列
此外,你还将编写获取序列信息的程序。DNA的GC含量如何?蛋白质的疏水性如何?你将学习相关的编程技术,利用它们就可以来回答这些类似的问题了。
在本章中你将学习到的技能涉及Python语言的基础,此处列出了其中的一部分:
- •标量变量
- •数组变量
- •替换、翻译等字符串的操作
- •从文件中读取数据
4.1序列数据的表征
本书的大部分内容都是处理表征DNA和蛋白质生物学序列的符号。在生物信息学中,使用这些符号来表征序列,这些符号和生物学家日常工作中用来表征序列的符号是完全一样的。
如前所述,DNA是由四种基础分子组成的,它们就是核酸,也叫做核苷酸或碱基;而蛋白质则是由20中基础分子组成的,它们就是氨基酸,也叫做残基。蛋白质的片段叫做肽(即缩氨酸)。不管是DNA还是蛋白质,本质上都是多聚物,由构成它们的基础分子首尾相连而形成。所以,仅仅通过碱基或氨基酸序列就可以表征DNA或蛋白质的基本结构。
这些都是最基本的定义。我假设你对这些内容都比较熟悉,或者打算查阅一本入⻔性的参考书来学习更加详细的内容。表4.1中列出的是碱基;在碱基上加上一个糖基,就可以得到核苷:腺苷、⻦苷、胞苷、胸苷和尿苷;在核苷上再加一个磷酸基团,就可以得到核苷酸:腺苷酸、⻦苷酸、胞苷酸、胸苷酸和尿苷酸。核酸就是由核苷酸通过化学键相连形成的序列。肽是由数个氨基酸相连而成的,更⻓的话就叫做多肽了。蛋白质是由一个或多个多肽组成的生物学功能基团。残基指的就是多肽链中的氨基酸。为了方便,如表4.1和表4.2所示,常用一个字⺟或三个字⺟的代码来表示核酸和氨基酸。(本书中主要用单字⺟代码来表示氨基酸。)
代码 |
核酸 |
A |
腺嘌呤 |
C |
胞嘧啶 |
G |
⻦嘌呤 |
T |
胸腺嘧啶 |
U |
尿嘧啶 |
M |
A 或 C (aMino, 氨基) |
R |
A 或 G (puRine, 嘌呤) |
W |
A 或 T (Weak, 作用力弱) |
S |
C 或 G (Strong, 作用力强) |
Y |
C 或 T (pYrimidine, 嘧啶) |
K |
G 或 T (keto, 酮基) |
V |
A 或 C 或 G(非 T) |
H |
A 或 C 或 T(非 G) |
D |
A 或 G 或 T(非 C) |
B |
C 或 G 或 T(非 A) |
N |
A 或 G 或 C 或 T (aNy, 任何一个碱基) |
表 4.1中的核酸代码不仅包括四种基本核酸的字⺟缩写,还定义了二个核酸、三个核酸或四个核酸所有可能组合的单字⺟缩写。在本书的大多数例子中,我仅使用 A、C、G、T、U 和 N。A、C、G 和 T 字⺟代表了 DNA 的核酸,而当 DNA(脱氧核糖核酸)转录成RNA(核糖核酸)时 U 将替换 T。当测序仪无法确定碱基时,常用 N 来表示“未知碱基”。稍后,在第 9 章中,在编程处理限制性酶切图谱时,我们还需要用到表示各种核酸组合的其他代码。有时,也会使用这些单字⺟代码的小写形式,这对 DNA 来说比较常⻅,但在蛋白质中则很少这样使用。
单字⺟代码 |
氨基酸 |
三字⺟代码 |
A |
Alanine, 丙氨酸 |
Ala |
B |
Aspartic acid or Asparagine, 天冬氨酸 或天冬酰胺 |
Asx |
C |
Cysteine, 半胱氨酸 |
Cys |
D |
Aspartic acid, 天冬氨酸 |
Asp |
E |
Glutamic acid, 谷氨酸 |
Glu |
F |
Phenylalanine, 苯丙氨酸 |
Phe |
G |
Glycine, 甘氨酸 |
Gly |
H |
Histidine, 组氨酸 |
His |
I |
Isoleucine, 异亮氨酸 |
Ile |
K |
Lysine, 赖氨酸 |
Lys |
L |
Leucine, 亮氨酸 |
Leu |
M |
Methionine, 甲硫氨酸 |
Met |
N |
Asparagine, 天冬酰胺 |
Asn |
P |
Proline, 脯氨酸 |
Pro |
Q |
Glutamine, 谷氨酰胺 |
Gln |
R |
Arginine, 精氨酸 |
Arg |
S |
Serine, 丝氨酸 |
Ser |
T |
Threonine, 苏氨酸 |
Thr |
V |
Valine, 缬氨酸 |
Val |
W |
Tryptophan, 色氨酸 |
Trp |
X |
Unknown, 未知氨基酸 |
Xxx |
Y |
Tyrosine, 酪氨酸 |
Tyr |
Z |
Glutamic acid or Glutamine, 谷氨酸或谷氨酰胺 |
Glx |
对于表 4.1和表 4.2中的代码来说,计算机科学中的术语和生物学中的术语还是有一定差别的。从计算机科学的⻆度来看,这两个表定义了两个按字⺟顺序排列的有限的符号集合,使用它们可以来构建字符串。字符串指的就是符号序列。比如,这句话本身就是一个字符串(this sentence is a string)。语言就是一个(有限或无限)字符串的集合。在本书中,语言主要就是 DNA 和蛋白质的序列数据。和在序列数据中的具有生物学含义的表征不同,生物信息学家常常会把真正的 DNA 或蛋白质序列称作“字符串”。两个不同学科中的术语会交叉使用,这就是一个例子。
如同你在表格中看到的那样,我们会使用简单的字⺟来表征数据,这和在纸张上手写时使用的方式是一样的。计算机实际上会用另外的代码来表征简单的字⺟,但你不用担心这些,只要记住在使用文本编辑器时将它们保存为 ASCII 或者纯文本即可。
ASCII 是计算机在内存中存储文本(和控制信息)数据的一种方式。当文本编辑器等程序读取数据时,计算机知道它正在读取 ASCII,因为计算机知道每个代码代表什么,所以它就会在屏幕上以一种容易识别的方式把相应的⺟显示出来。总而言之,知道 ASCII 是计算机表征文本的一种代码就足够了。
4.2 存储 DNA 序列的程序
让我们来写一个小程序吧,它把 DNA 存储在变量中,然后把它打印输出到屏幕上。
我们会用最常用的方式——由 A、C、G 和 T 四个字⺟组成的字符串——来书写 DNA,并且把存储 DNA 的变量叫做 DNA。换言之,DNA 就是程序中 DNA 序列数据的名称。注意这一点,在 Python 中,变量就是你打算处理的数据的名称,使用该名称,你可以对数据进行完全的访问。例 4.1是完整的程序。
#!/usr/bin/env python # Storing DNA in a variable, and printing it out # First we store the DNA in a variable called DNA DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC' # Next, we print the DNA onto the screen print(DNA) # Finally, we'll specifically tell the program to exit. exit()
在第 2 章中,我们已经学习了文本编辑器和运行 Python 程序的知识,运用这些知识,输入例子中的代码(或者从书籍官网上把它复制下来)并保存到一个文件中。牢记一定要把程序保存成 ASCII 或者纯文本格式,否则 Python 在读取该文件时可能会出现问题。
接下来就是运行程序了。运行程序的具体步骤取决于你使用的计算机(参看第 2 章)。
我们假定程序是你计算机中一个叫做 example4-1 的文件。回顾第 2 章中相关的知识,如果你想在 Unix、Linux 或者 Mac 中运行这个程序,需要在 shell 窗口中键入:
python example4-1
在 Windows 中,在 MS-DOS 命令窗口中键入:
python example4-1
如果你成功运行了该程序,在你的计算机屏幕上你将看到相应的输出。
4.2.1 控制流
例 4.1展示了所有 Python 程序都将用到的许多理念,其中一个便是控制流的理念,即计算机是以什么顺序来执行程序中的语句的。
所有的程序都从第一行开始,除非明确指明了其他的运行顺序,否则它将一条一条地执行语句,直到程序的最后一行。例 4.1只是简单的从头到尾执行程序,并没有其他的运行支路。
在后续的章节中,你将学习到如何编程控制程序的执行顺序。
4.2.2 再说注释
现在让我们看一下例 4.1程序中的细节。你会发现其中有许多空行,它们的存在使得程序更容易被人阅读。另外,注意以 # 起始的注释。在第 3 章中提到过,当 Python 运行时,它会把这些注释和空行全部忽略掉。事实上,对于 Python 来说,下面的程序和例 4.1中的程序是完全一样的:
#!/usr/bin/python DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC'; print(DNA); exit();
在例 4.1中,我使用了大量的注释。代码开始的注释指明了程序的用途、作者以及其他信息,当其他人需要理解代码时,这些信息会对他有很大的帮助。注释还解释了代码每一部分的作用,有时还对代码的工作原理进行了阐释。
了解注释的重要性是非常必要的。在大多数高校的计算机科学课程的课堂作业中,没有注释的程序通常会得到很低的分数甚至不及格;而在工作中不对代码进行注释的程序员,他的职业生涯将是短暂而失败的。
4.2.3 命令解释
程序的第一行以 # 起始,这使得它看上去像是一行注释,但又不像是有什么信息含量的注释:
#!/usr/bin/python
这是比较特殊的一行,叫做命令解释,它告诉运行 Unix 或者 Linux 的计算机,这是一个 Python 程序。在不同的计算机中,这一行可能会有些许的差异。在某些计算机中,这一行并不是必需的,因为计算机可以通过其他信息识别出 Python。在 Windows 计算机中,通常会配置成把以.py 结尾的任何程序都假定成 Python 程序。在 Unix 或 Linux 中、Windows 的命令窗口中、或者 MacOS X 的 shell 中,你可以键入 python my_program,这样你的 Python 程序 my_program 中就不需要这样特殊的一行了。然而,通常都会写上这一行,所以在我们所有程序的开头也都会有它。
4.2.4 语句
例 4.1中的下一行代码把 DNA 存储到了一个变量中:
DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC'
这在计算机语言中是非常常⻅、非常重要的,所以我们将对它进行详细的解释。总的来说,你将会看到 Python 和编程语言的一些基本特性,所以不要跳读了,停下来慢慢阅读学习吧。
这行代码叫做语句。在 Python 中,语句以换行结尾。
更准确的来说,这一行是一个赋值语句。在该程序中,它作用就是把 DNA 存储到DNA 变量中。正如在接下来的小节中你将看到的那样,此处有许多基本的事件发生。
变量
首先,让我们来看看 DNA 变量。它的名字有些随意,你可以给它起另外一个名字,程序仍会正常运行。举个例子,如果你把下面这两行:
DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC' print(DNA)
换成这样的两行:
A_poem_by_Seamus_Heaney = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC' print(A_poem_by_Seamus_Heaney)
程序仍将像原来那样正常运行,把 DNA 打印输出到计算机屏幕上。不管怎样,计算机程序中变量的名字都是由你来起的。(要满足特定的限制:在 Python 中,变量的名字只能由大小写字⺟、数字和下划线 _ 组成,而且第一个字符不能是数字。)
前面已经强调过,使用空行和注释可以使代码更加易读,变量的命名也存在同样的问题。不管变量名是 DNA 还是 A_poem_by_Seamus_Heaney,对于计算机来说都没有什么特殊的含义,但对于阅读程序的人来说就不一样了。有意义的变量名,可以清晰地表明程序中变量的作用,使得程序更加容易理解。其他的名字可能会使得程序的功能和变量的作用不甚明朗。使用精心选择的变量名是自文档化代码的一部分(self-documenting code)。精心选择变量名的话,你仍然需要注释,但就不需要那么多的注释了。
你会注意到 DNA 变量名以一个美元符号起始。在 Python 中,这样的变量叫做标量变量,它是存储单个数据项目的变量。在存储字符串或者各种各样的数字(如,字符串 hello,或者 25、6.234、3.5E10、-0.8373 这样的数字)时,需要使用标量变量。一个标量变量一次只能存储数据中的一个项目。
字符串
在例 4.1中,标量变量 $DNA 存储着用 A、C、G 和 T 表征的 DNA。如前所述,在计算机科学中,字⺟序列叫做字符串。在 Python 中,你需要把它放在引号中来表明它是字符串。可以使用单引号,就像例 4.1中那样,也可以使用双引号。(稍后你会看到两者的区别。)因此,DNA 就被表征成:
'ACGGGAGGACGGGAAAATTACTACGGCATTAGC'
赋值
在 Python 中,要把一个变量设成特定的值,需要使用 = 符号。= 符号因此被叫做赋值操作符。在例 4.1中,值
'ACGGGAGGACGGGAAAATTACTACGGCATTAGC'
被赋给了 DNA 变量。赋值后,你可以通过变量名来获取它的值,就像例 4.1中的 print 语句那样。
在赋值语句中,每个部分的顺序是非常重要的。要赋给变量的值在赋值操作符的右边,而需要赋值的变量总在赋值操作符的左边。在编程手册中,有时你可能会遇到 lvalue 和 rvalue 这样的术语,它们分别指代赋值操作符左边和右边的项目。在编程语言中,使用 = 符号进行赋值有一段很⻓的历史。然而,这也导致了某种形式的混乱:比如说,在数学中,使用 = 时表示等号两边的数是相等的。所以,一定要牢记,在 Python 中,= 符号并不表示相等,而是把值赋给一个变量。(稍后我们会看到如何表示相等)。
来总结一下,对于这个语句,我们都学习到了那些知识:
DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC'
这是一个赋值语句,它把表征 DNA 的字符串赋给了变量 DNA,作为这个变量的值。
打印输出
该语句:
print(DNA)
把 ACGGGAGGACGGGAAAATTACTACGGCATTAGC 打印输出到计算机屏幕上。注意,print 语句处理的是标量变量,它把它的值打印输出出来,此处就是变量 DNA 包含的字符串。
稍后,你将看到更多关于打印输出的内容。
退出
最后,exit() 语句告诉计算机退出程序。在 Python 语言中,在程序的最后并不需要 exit() 语句,一旦运行到结尾,程序就会自动退出。但放上这么一个语句也没什么坏处,还明确表示了程序的结束。你会看到,在正常结束之前,程序会因为某些错误而退出,所以 exit() 语句还是非常有用的。
4.3 连接DNA片段
现在,我们对例 4.1进行简单的修改,来演示一下如何把 DNA 片段连接起来。所谓连接指的就是把一个东西附加在另一个东西的结尾上。生物学家都知道,在生物学实验室中把 DNA 序列连接起来是非常常⻅的工作,比如把克隆插入到细胞载体中,或者在基因表达过程中把剪切的外显子连接起来。许多生物信息学的软件包都能够进行这样的工作,因此这里只是把它作为一个实例来讲解。例 4.2演示了关于字符串、变量和打印输出语句的更多内容。
#!/usr/bin/env python # Concatenating DNA # Store two DNA fragments into two variables called DNA1 and DNA2 DNA1 = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC' DNA2 = 'ATAGTGCCGTGAGAGTGATGTAGTA' # Print the DNA onto the screen print("Here are the original two DNA fragments: ") print(DNA1, " ") print(DNA2, " ") # Concatenate the DNA fragments into a third variable and print them # Using "% format" DNA3 = "%s%s" % (DNA1, DNA2) print("Here is the concatenation of the first two fragments (version 1): ") print("%s " % DNA3) # An alternative way using the "add operator": # Concatenate the DNA fragments into a third variable and print them DNA3 = DNA1 + DNA2 print("Here is the concatenation of the first two fragments (version 2): ") print("%s " % DNA3) # Print the same thing without using the variable DNA3 print("Here is the concatenation of the first two fragments (version 3): ") print(DNA1, DNA2, " ", sep='') exit()
就像你看到的那样,这里有三个变量:DNA1、DNA2 和 DNA3。为了对运行过程进行说明,我添加了一些 print 语句,这样打印到计算机屏幕上的输出就不仅仅是一个接一个 DNA 片段了,看起来会更加明了一些。下面是例 4.2的输出:
Here are the original two DNA fragments: ACGGGAGGACGGGAAAATTACTACGGCATTAGC ATAGTGCCGTGAGAGTGATGTAGTA Here is the concatenation of the first two fragments (version 1): ACGGGAGGACGGGAAAATTACTACGGCATTAGCATAGTGCCGTGAGAGTGATGTAGTA Here is the concatenation of the first two fragments (version 2): ACGGGAGGACGGGAAAATTACTACGGCATTAGCATAGTGCCGTGAGAGTGATGTAGTA Here is the concatenation of the first two fragments (version 3): ACGGGAGGACGGGAAAATTACTACGGCATTAGCATAGTGCCGTGAGAGTGATGTAGTA
例 4.2和例 4.1有许多相似之处。让我们来看一下两者的不同吧。作为开始,我们会发现 print 函数语句中多了一些看起来并不直观的东西:
print(DNA1) print(DNA1, DNA2, " ", sep='')
像前述一样,print 语句中有包含 DNA 的变量,但现在后面又多了逗号和 " ",关键字参数“sep=''”。 ”是换行符,告诉计算机继续下一行的开头进行后续打印。关键字参数“sep=''”告诉每个打印对象(DNA1、DNA2)之间使用无缝连接,默认是一个空格。
看一下例 4.2的代码,确保已经明白这些换行符是如何决定输出效果的。空行指的就是没有任何打印输出的一行。取决于你的操作系统,它可能仅仅是一个换行符,也可能是换⻚符和回⻋符的组合(在这种情况下,它也被叫做空白行),还有可能包含空格和制表符等非打印的空白字符。注意在 print 语句中还有逗号。逗号分隔列表中的项目,print 语句会打印输出列表中的所有项目。仅此而已。现在,让我们看一下把 DNA1 和 DNA2 两个 DNA 片段连接到 DNA3 变量中的语句吧:
DNA3 = "%s%s" % (DNA1, DNA2)
对 DNA3 进行的赋值是一个典型的赋值操作,就像你在例 4.1看到的那样,变量名后面跟着 = 符号,= 后则是要被赋予的值。
赋值语句中右边的值是包裹在双引号中的字符串。双引号会把字符串中的%s替换为变量的值,这叫做字符串格式化。所以事实上,此处的字符串就是 DNA1 变量中的 DNA,后面紧跟着 DNA2 变量中的 DNA。两个 DNA 片段连接后就被赋值给了变量 DNA3。
在把连接结果赋值给 DNA3 变量后,你把它打印出来,后面跟着一个空行:
print("%s " % DNA3)
程序的下面一部分就演示了连接两个字符串的另外一种办法——使用加法操作符。当把加法操作符放在两个字符串中间时,它会把原来的两个字符串连接起来,产生一个新的字符串。所以这一行,演示了加法操作符的使用。
DNA3 = DNA1 + DNA2
最后,来练习一下 Python 的另外一种方法,仅仅使用 print 语句就可以完成同样的连接任务:
print(DNA1, DNA2, " ", sep='')
此处的 print函数语句有用逗号分隔开的四部分:存储在两个变量中的两个 DNA 片段、一个换行符和连接符参数。
在结束这一小节之前,让我们来看看 Python变量的其他用法吧。你已经看到,使用变量可以存储 DNA 序列数据的字符串。还有其他类型的数据,编程语言也需要变量来存储它们。在 Python 中,一个像 DNA 这样的标量变量可以存储字符串、整数、浮点数(有小数点的数字)、布尔值(True 或 False)等。现在,试着在例 4.1或例 4.2中添加如下几行,在标量变量中存储一个数字并把它打印出来:
number = 17 print(number, " ")
4.4 转录:从 DNA 到 RNA
作为生物信息学中的 Python 程序员,你很大一部分时间都是在做类似于例 4.1和例 4.2那样的事情:你获取一些数据,可能是 DNA、蛋白质、GenBank 条目或者其他数据,处理这些数据,并把处理结果输出打印出来。例 4.3是处理 DNA 的另一个程序:它把 DNA 转录成 RNA。在细胞中,把 DNA 转录成 RNA 是由精巧、复杂且有勘误功能的分子机器完成的。此处仅是简单的替换而已。当 DNA 转录成 RNA 时,所有的 T 都会被替换成 U,这也正是我们的程序所要做的全部工作。
#!/usr/bin/env python # Transcribing DNA into RNA # The DNA DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC' # Print the DNA onto the screen print("Here is the starting DNA: ") print("%s " % DNA) # Transcribe the DNA to RNA by substituting all T's with U's. RNA = DNA.replace('T', 'U') # Print the RNA onto the screen print("Here is the result of transcribing the DNA to RNA: ") print("%s " % RNA) # Exit the program. exit()
这是例 4.3的输出:
Here is the starting DNA: ACGGGAGGACGGGAAAATTACTACGGCATTAGC Here is the result of transcribing the DNA to RNA: ACGGGAGGACGGGAAAAUUACUACGGCAUUAGC
这个简短的程序展示了 Python 重要的特性:轻松处理 DNA 字符串等文本数据的能力。
类似的处理可能会有多种:翻译、反转、替换、删除和重排序等等。总的来说,Python 在此类任务中的便利是其能够在生物信息学领域成功及在程序员中广泛应用的主要原因。
首先,程序生成了一个 DNA 的拷⻉,并把它存储在叫做 RNA 的变量中:
RNA = DNA
注意在执行这个语句后,叫做 RNA 的这个变量存储的就是 DNA。你可以给变量起任何你喜欢的名字,这是完全合法的,但不准确的变量名可能会导致一些混乱。在这个例子中,拷⻉完变量值后,紧跟着的是内容丰富的注释,之后便是让 RNA 变量真正包含 RNA 的语句,所以此处给它起名为 RNA 完全没有问题。有一个让 RNA 只包行 RNA 而不包含其他内容的方法:
RNA = DNA.replace('T', 'U')
在例 4.3中,转录过程发生在这个语句中:
RNA = DNA.replace('T', 'U')
在这个语句中有两个新的项目:赋值操作符(=)和替换函数 replace。很明显,赋值操作符 = 用于包含字符串的变量,如此处的 RNA 变量储存 DNA 序列转录后的RNA序列数据。replace函数就是“把 DNA 变量存储的字符串数据中的所有 T 都替换成 U”。
字符串操作对于文本处理来说至关重要,在后续的章节中你将看到,字符串操作是 Python 最为强大的特性之一。
4.5 使用Python文档
对于 Python 程序员来说,最重要的资源便是 Python 的文档。它应该已经安装在了你的电脑上,另外通过因特网在 Python 网站上也可以找到它。对于不同的计算机系统来说,Python 文档的格式可能有少许的差别,但网络版对任何一个人来说都是一样的,这也正式我在本书中参考的版本。参看附录 A中的参考资料,你会找到对于 Python 文档不同资源的讨论。
来试一下,让我们找找 print 操作符的文档吧。首先,打开你的网⻚浏览器,进入 https://www.python.org/ 网站,然后点击python3.x文档链接,依次选择“一般索引”(”General Index”)。你会看到一个把 Python 对象按字⺟顺序进行排列的一个冗⻓的列表。一旦你找到这个⻚面,你可能需要把它收藏到浏览器的书签中,因为你会频繁地访问该⻚面。
现在点击 print 来查看 print 函数的文档吧。
看一下文档中的例子,看看 Python 语言的特性是如何被运用的。这往往是你找到所需内容的最快方法。
一旦在你的屏幕中打开了文档,你会发现通过阅读它会找到一些答案,但也会产生其他一些疑问。文档试图以一种简洁的形式把所有的内容都包含进去,但这却会让初学者们心生胆怯。比如,print 函数文档的开头还比较简单:“将对象打印到文本流文件,以sep分隔,然后结束。 sep,end,file和flush(如果存在)必须作为关键字参数给出”。但之后就是一堆的胡言乱语了(在你学习的现阶段它看起来确实是这样):sep?end?file?flush?
文档中的所有信息都是必需的,毕竟,你需要在某个地方找到这样的全部内容!通常情况下,你可以忽略掉那些对你来说毫无意义的内容。
Python 文档中也包含了一些对学习 Python 有很大帮助的教程。有时,它们会假定你掌握了比一个编程语言初学者应当掌握的知识更多的知识,但你会发现它们仍然非常有用。翻阅文档是在学习 Python 语言过程中快速成⻓的绝佳途径。
4.6 在Python中计算反向互补
第 1 章中已经提到了,DNA 聚合物是由核苷酸构成的。考虑到 DNA 双螺旋中两条链的亲密关系,最好能编写这样一个程序:给出一条链,输出另一条链。这样的工作对许多生物信息学应用程序来说都是非常重要的一部分。比如,当在数据库中查询某条 DNA 时,常常也需要自动查询该 DNA 的反向互补序列,因为你手上的序列有可能是某个已知基因的负链。回到正题,这是例 4.4,它使用了一些新的 Python 特性。就像你将看到的那样,它首先尝试一种方法,失败了,然后尝试另外一种方法,最终取得了成功。
#!/usr/bin/perl -w # Calculating the reverse complement of a strand of DNA # The DNA DNA = 'ACGGGAGGACGGGAAAATTACTACGGCATTAGC' # Print the DNA onto the screen print("Here is the starting DNA: ") print("%s " % DNA) # Calculate the reverse complement # Warning: this attempt will fail! # # First, copy the DNA into new variable revcom # (short for REVerse COMplement) # Notice that variable names can use lowercase letters like # "revcom" as well as uppercase like "DNA". In fact, # lowercase is more common. # # It doesn't matter if we first reverse the string and then # do the complementation; or if we first do the complementation # and then reverse the string. Same result each time. # So when we make the copy we'll do the reverse in the same statement. # revcom = DNA[::-1] # # Next substitute all bases by their complements, # A->T, T->A, G->C, C->G # revcom = revcom.replace('A', 'T') revcom = revcom.replace('T', 'A') revcom = revcom.replace('G', 'C') revcom = revcom.replace('C', 'G') # Print the reverse complement DNA onto the screen print("Here is the reverse complement DNA: ") print("%s " % revcom) # # Oh-oh, that didn't work right! # Our reverse complement should have all the bases in it, since the # original DNA had all the bases--but ours only has A and G! # # Do you see why? # # The problem is that the first two substitute commands above change # all the A's to T's (so there are no A's) and then all the # T's to A's (so all the original A's and T's are all now A's). # Same thing happens to the G's and C's all turning into G's. # print(" That was a bad algorithm, and the reverse complement was wrong! ") print("Try again ... ") # Make a new copy of the DNA (see why we saved the original?) revcom = DNA[::-1] # See the text for a discussion of str.maketrans('ATCGatcg', 'TAGCtagc') revcom =revcom.translate(str.maketrans('ATCGatcg', 'TAGCtagc')) # Print the reverse complement DNA onto the screen print("Here is the reverse complement DNA: ") print("%s " % revcom) print(" This time it worked! ") exit()
这是你屏幕上例 4.4的输出:
Here is the starting DNA: ACGGGAGGACGGGAAAATTACTACGGCATTAGC Here is the reverse complement DNA: GGAAAAGGGGAAGAAAAAAAGGGGAGGAGGGGA That was a bad algorithm, and the reverse complement was wrong! Try again ... Here is the reverse complement DNA: GCTAATGCCGTAGTAATTTTCCCGTCCTCCCGT This time it worked!
通过从不同的末端开始读起,也就是说一条链从左向右读,另一条链从右向左读,你可以检查一下 DNA 的两条链是不是反向互补的。当你读这两条链的时候,比较它们对应的碱基,应该总是 C 和 G 对应、A 和 T 对应。
在第一次尝试中,试着从原始的 DNA 和反向互补后的 DNA 中读几个字符,你会发现反向互补的第一次尝试失败了。这是一个错误的算法。
这是在你程序中经常要经历的一种体验。你写一个程序来完成某项任务,但却发现它并没有像你期望的那样工作。在这种情况下,我们需要运用已掌握的知识来尝试解决全新的问题。它并没有如期望那样完成任务,哪里出错了呢?
你会发现这样的经历会非常相似:你编写代码,但它并不工作!然后你就修正语法(这通常是最简单的,而且利用错误信息提供的线索就可以轻松修正语法),或者对问题进行思考,找到问题所在,并尝试设计一种新的可以成功的方法。通常情况下,这需要你浏览语言的文档,查找语言工作过程的细节内容,同时期望能找到一个解决问题的特性。
如果这个问题在计算机上能够解决,那么你用 Python 也可以将它解决。问题是,能够精确到什么程度?
在例 4.4中,计算反向互补的第一次尝试失败了。使用四个全局替换,序列中的每一个碱基都作为一个整体进行了处理。还需要另外的方法。你可以从左向右查阅 DNA,每次只查找一种碱基,把它替换成互补的碱基,然后再在 DNA 中查找另外一种碱基,一直到字符串的结尾。然后把字符串反转过来,任务就完成了。事实上,这是一个非常好的方法,在 Python 中也不难实现。但还需要学习语言的其他内容才行,这将在第 5 章中进行讲解。
然而,在这个例子中,maketrans先构建了ATCG对应的转换关系,然后translate根据转换关系翻译序列。
序列切片也是我们需要的,虽然有点大材小用了。它被设计用来反转元素的顺序,包括字符串,就像在例 4.4中看到的那样。
4.7 蛋白质、文件和数组
到现在为止,我们已经编写了处理 DNA 序列数据的程序,现在,我们也将要处理同样重要的蛋白质序列数据。接下来的几个小节将要学习一些新的知识,这里是一个简单的概述:
- 在 Python程序中如何使用蛋白质序列数据
- 如何从一个文件中读入蛋白质序列数据
- Python 语言中的列表在本章的剩余部分中,蛋白质和 DNA 序列数据都将使用到。
4.8 从文件中读取蛋白质序列数据
程序要和计算机磁盘中的文件进行交互。这些文件可以存储在任何永久性存储介质中 ——硬盘、光盘、软盘、Zip 磁盘、磁带等。
让我们看看如何从文件中读取蛋白质序列数据吧。首先,在你的计算机上创建一个文件(使用文本编辑器),并在其中存储一些蛋白质序列数据。给这个文件起名为
NM_021964fragment.pep(你可以在书籍网站上下载到该文件)。你将使用下列数据(人类锌指蛋白 NM_021964 的一部分):
MNIDDKLEGLFLKCGGIDEMQSSRTMVVMGGVSGQSTVSGELQD SVLQDRSMPHQEILAADEVLQESEMRQQDMISHDELMVHEETVKNDEEQMETHERLPQ GLQYALNVPISVKQEITFTDVSEQLMRDKKQIR
你可以给它起任意一个名字,只要和文件夹中已有的文件不重名即可。
就像好的变量名对于理解程序至关重要一样,好的文件名和文件夹名也是非常重要的。如果你有一个项目,这个项目会生成大量的文件,那你就需要认真考虑如何对这些文件和文件夹命名和组织了。不管是对于某一个研究人员还是对于一个庞大的多国团队来说,这都是非常现实的问题。花一些功夫来给文件起一个富含信息量的名字是非常重要的。文件名 NM_021964fragment.pep 来源于蛋白质的 GenBank ID。同时,它还表明了这只是数据的一部分,而文件的后缀.pep 则提醒你文件中保存的是肽或蛋白质序列数据。当然,其他的命名方案可能会更加适合于你。不管怎样,关键的一点就是不需要打开文件,仅仅通过文件名就可以对文件保存的数据有所了解。你已经创建或者下载了保存有蛋白质序列数据的文件,那我们就来编写一个程序吧,这个程序从文件中读入蛋白质序列数据并把它保存到变量中。例 4.5是第一次尝试,随着学习我们会逐步对它进行扩展。
#!/usr/bin/env python # Reading protein sequence data from a file # The filename of the file containing the protein sequence data proteinfilename = 'NM_021964fragment.pep' # First we have to "open" the file, and associate # a "filehandle" with it. We choose the filehandle # PROTEINFILE for readability. PROTEINFILE = open(proteinfilename) # Now we do the actual reading of the protein sequence data from the file, # by using the read function to get the input from the # filehandle. We store the data into our variable protein. protein = PROTEINFILE.readline() # Now that we've got our data, we can close the file. PROTEINFILE.close() # Print the protein onto the screen print("Here is the protein: ") print(protein) exit()
下面是例 4.5的输出:
Here is the protein: MNIDDKLEGLFLKCGGIDEMQSSRTMVVMGGVSGQSTVSGELQD
注意只有文件的第一行被打印输出出来,稍后我会解释为什么会这样。
让我们仔细研究一下例 4.5吧。在把文件名保存到 proteinfilename 变量中后,下面这个语句会打开文件:
PROTEINFILE = open(proteinfilename)
打开文件后,你就可以对它进行各种操作了,比如读取、写入、查找、定位到文件的特定位置、清除文件中的所有数据,等等。注意,程序假设保存在 proteinfilename 变量中的文件是存在的,而且可以打开。你将会看到如何来确认这一点,现在你可以进行一下这样的尝试:更改 proteinfilename 变量中文件的名字,这样计算机上就没有叫原来那个名字的文件了,之后运行一下程序。如果文件并不存在,你会看到一些错误信息。
如果你查阅 open 函数的文档,你会看到好多选项,大多数选项都是在打开文件后让你精确指定如何使用文件的。
让我们来看一下 PROTEINFILE 这个数据,它叫做文件句柄。没有必要理解文件句柄真正是什么,它们就是你处理文件时使用的东西。对于文件句柄,不一定必须使用大写字⺟,但这是一个广为接受的惯例。在使用 open 语句对文件句柄赋值后,对文件进行的任何交互操作都将通过使用文件句柄来进行。
使用这个语句,就可以把数据读入到程序中了:
protein = PROTEINFILE.readline()
通过调用文件句柄的readline函数,你就可以从程序外部的来源中读入数据了。在这个例子中,我们读入了 NM_021964fragment.pep 文件中的数据,该文件的名字保存在 proteinfilename 变量中,而 open 语句则把它和文件句柄关联了起来。数据现在就保存到了 protein 变量中,之后可以把它打印输出出来。
然而,就像我们前面已经注意到的那样,只有这个多行文件的第一行被打印了出来。这是问什么呢?因为对于读取文件还有一些知识需要学习。有许多方法可以读取整个文件,例 4.6演示了其中的一种。
#!/usr/bin/envl python # Reading protein sequence data from a file, take 2 # The filename of the file containing the protein sequence data proteinfilename = 'NM_021964fragment.pep' # First we have to "open" the file, and associate # a "filehandle" with it. We choose the filehandle # PROTEINFILE for readability. PROTEINFILE = open(proteinfilename) # Now we do the actual reading of the protein sequence data from the file, # by using the function readline to get the input from the # filehandle. We store the data into our variable $protein. # # Since the file has three lines, and since the read only is # returning one line, we'll read a line and print it, three times. # First line protein = PROTEINFILE.readline() # Print the protein onto the screen print(" Here is the first line of the protein file: ") print($protein) # Second line protein = PROTEINFILE.readline() # Print the protein onto the screen print(" Here is the second line of the protein file: ") print(protein) # Third line protein = PROTEINFILE.readline() # Print the protein onto the screen print(" Here is the third line of the protein file: ") print(protein) # Now that we've got our data, we can close the file. PROTEINFILE.close() exit()
下面是例 4.6的输出:
Here is the first line of the protein file: MNIDDKLEGLFLKCGGIDEMQSSRTMVVMGGVSGQSTVSGELQD Here is the second line of the protein file: SVLQDRSMPHQEILAADEVLQESEMRQQDMISHDELMVHEETVKNDEEQMETHERLPQ Here is the third line of the protein file: GLQYALNVPISVKQEITFTDVSEQLMRDKKQIR
这个程序中比较有趣的一点就是,它演示了如何成功得从文件中读取数据。每当把数据读取到 protein 这样的变量中后,都会对文件的下一行进行读取。上一次读取到哪儿了,现在需要读取哪一行,程序对这些信息都有所记录。
另一方面,程序的缺陷也是显而易⻅的。对于输入文件的每一个行,都需要编写一行代码来读入,这样太不方便了。但是,在 Python 中有两个特性可以很好的解决这个问题:列表(下一小节进行介绍)和循环(参看第 5 章)。
4.9 列表
在计算机语言中,列表是存储多个标量值的变量。这些标量的值可以是数字、字符串,或者,像此处的例子一样,也可以是蛋白质序列数据输入文件中的每一行。让我们看看如何来使用列表吧。例 4.7演示了如何使用列表来把输入文件中的所有行都读入进来。
#!/usr/bin/env python # Reading protein sequence data from a file, take 3 # The filename of the file containing the protein sequence data proteinfilename = 'NM_021964fragment.pep' # First we have to "open" the file PROTEINFILE = open(proteinfilename) # Read the protein sequence data from the file, and store it # into the array variable proteins proteins = PROTEINFILE.readlines() # Print the protein onto the screen print(*proteins, sep='') # Close the file. PROTEINFILE.close() exit()
下面是例 4.7的输出:
MNIDDKLEGLFLKCGGIDEMQSSRTMVVMGGVSGQSTVSGELQD SVLQDRSMPHQEILAADEVLQESEMRQQDMISHDELMVHEETVKNDEEQMETHERLPQ GLQYALNVPISVKQEITFTDVSEQLMRDKKQIR
就像你看到的,这就是文件中的全部数据。我们成功了!使用数组的方便性显而易⻅——只需要一行代码就可以把所有的数据都读入到程序中了。
请注意,Python中的变量可以赋值为任何类型,另外,“*”将多个元素的列表解包为print函数的多个对象,然后print函数打印每个对象。
列表是存储多个标量值的变量。列表中的每一个项目或者元素都是标量值,可以通过在列表中的位置(它的下标或者偏移量)对它进行引用。让我们来看几个列表和列表常用操作的实例吧。我们定义一个叫做 bases 的输出,其中存储着 A、C、G 和 T 四个碱基。现在我们对它应用几个常⻅的列表操作。
这是一个代码片段,它演示了如何初始化列表,以及如何使用下标来访问列表中的单个元素:
# Here's one way to declare an array, initialized with a list of four scalar values. bases = ['A', 'C', 'G', 'T'] # Now we'll print each element of the array print("Here are the array elements:") print(" First element: ") print(bases[0]) print(" Second element: ") print(bases[1]) print(" Third element: ") print(bases[2]) print(" Fourth element: ") print(bases[3])
这个代码片段将会输出:
First element: A Second element: C Third element: G Fourth element: T
你可以像这样把元素一个接一个的输出出来:
bases = ['A', 'C', 'G', 'T'] print(" Here are the array elements: ") print(*bases, sep='')
它会产生这样的输出:
Here are the array elements: ACGT
你也可以输出用空格分隔的元素(注意 print 语句中的sep参数):
bases = ['A', 'C', 'G', 'T'] print(" Here are the array elements: ") print(*bases)
它会产生这样的输出:
Here are the array elements: A C G T
使用 pop,你可以从列表的末尾拿掉一个元素:
bases = ['A', 'C', 'G', 'T'] base1 = bases.pop() print("Here's the element removed from the end: ") print(base1, " ") print("Here's the remaining array of bases: ") print(*bases)
它会产生这样的输出:
Here's an element removed from the beginning: A Here's the remaining array of bases: C G T
使用 pop,你也可以从列表的开头拿掉一个碱基:
bases = ['A', 'C', 'G', 'T'] base2 = bases.pop(0) print("Here's an element removed from the beginning: ") print(base2, " ") print("Here's the remaining array of bases: ") print(*bases)
它会产生这样的输出:
Here's an element removed from the beginning: A Here's the remaining array of bases: C G T
使用insert,你可以把一个元素添加到列表的开头:
bases = ['A', 'C', 'G', 'T'] base1 = bases.pop() bases.insert(0, base1) print("Here's the element from the end put on the beginning: ") print(*bases, " ")
它会产生这样的输出:
Here's the element from the end put on the beginning: T A C G
使用append,你可以把一个元素添加到列表的末尾:
bases = ['A', 'C', 'G', 'T'] base2 = bases.pop(0) bases.append(base2) print("Here's the element from the beginning put on the end: ") print(*bases, " ")
它会产生这样的输出:
Here's the element from the beginning put on the end: C G T A
使用reverse函数,反转列表:
bases = ['A', 'C', 'G', 'T'] reverse = bases.reverse() print("Here's the array in reverse: ") print(reverse, " ")
它会产生这样的输出:
Here's the array in reverse: T G C A
使用len函数,计算列表长度:
bases = ['A', 'C', 'G', 'T'] print("Here's the length of the array: ", len(bases))
它会产生这样的输出:
Here's the length of the array: 4
使用insert函数,你可以在数组的任意一个位置插入一个元素:
bases = ['A', 'C', 'G', 'T'] bases.insert(2, 'X') print("Here's the array with an element inserted after the 2nd element: ", *bases)
它会产生这样的输出:
Here's the array with an element inserted after the 2nd element: A C X G T
4.10 练习题
习题 4.1
探索编程语言对语法错误的敏感性。对于一个可以运行的程序,试着把其中任意一个语句末尾的分号删除掉,看看会产生什么样的错误信息。试着改变其他的语法项:添加一个小括号或者大括号,把“print”或其他保留字拼错,键入或者删除任何一些内容。程序员对这样的错误习以为常。即使精通语言后,在你一点点编写代码时,仍然会常常出现这样的语法错误。注意一个错误是如何导致多行错误报告的。Python 报告的出现错误的行是不是完全准确?
习题 4.2
编写一个程序,把一个整数保存到变量中并把它打印出来。
习题 4.3
编写一个把 DNA(原始序列可以是大写或者小写)以小写形式(acgt)输出的程序;再编写一个把 DNA 以大写形式(ACGT)输出的程序。
习题 4.4
任务和练习 4.3 中的一样,但这次使用maketrans和translate函数。
习题 4.5
有时,信息也可以从 RNA 流向 DNA。编写一个把 RNA 反转录成 DNA 的程序。
习题 4.6
读取两个文件的数据,在输出第一个文件的内容后紧接着输出第一个文件的内容。
习题 4.7
这是一个更有难度的练习题。编写一个程序,读去一个文件,然后逆序输出它的每一行,也就是首先输出它的最后一行。你可能需要看一下pop、insert 和 reverse 这三个函数,从中选择一个或多个来完成此练习。你可能需要预习一下第 5 章,这样就可以在程序中使用循环了。但取决于你采取的方案,这并不是必需的。或者,你可以对包含所有行的列表使用 reverse 函数。