第 5 章基序和循环

本章将在第 4 章的基础上进一步介绍 Python 语言的基础知识。到本章结束的时候，你将学会：

在 DNA 或蛋白质中查找基序
通过键盘与用户进行交互
把数据写入文件
使用循环
使用基本的正则表达式
根据条件测试的结果采取不同的行动
通过字符串和列表的操作对序列数据进行细致的处理所有这些主题，加上在第 4 章学习的知识，将使你能够编写出实用的生物信息学程序。在本章中，你将学习编写一个在序列数据中查找基序的程序。

5.1 流程控制

流程控制指的就是程序中语句执行的顺序。除非明确指明不按顺序执行，否则程序将从最顶端的第一个语句开始，顺序执行到最底端的最后一个语句。有两种方式可以告诉程序不按照顺序执行：条件语句和循环。条件语句只会在条件测试成功的前提下执行相应的语句，否则就会直接跳过这些语句。循环会一直重复那些语句，直到相应的测试失败为止。

5.1.1 条件语句

让我们再看一下 open 语句吧。回忆一下，当你试图打开一个并不存在的文件时，你会看到错误信息。在你尝试打开文件之前，可以对文件是否存在进行明确的测试。事实上，类似的测试都是计算机语言最强大的特性之一。if和if-else 条件语句就是 Python 语言中用来进行类似测试的两个语句。

这类语句的最主要特性就是可以对条件进行测试。条件测试的结果可以是 True 或者是 False 。如果测试结果为True ，那么后面的语句就会被执行；与之相反，如果测试结果为False ，这些语句就会被跳过而不执行。

那么，“什么是真呢”？对于这个问题，不同计算机语言的回答会稍有不同。

本节包含了一些演示 Python 条件语句的实例，其中的真-假条件测试就是看看两个数字是否相等。注意，因为一个等号 = 是用来对变量进行赋值的，所以两个数字的相等要用两个等号 == 来表示。

下面这个例子演示了条件测试的结果是 True 还是 False 。平时你很少会用到这么简单的测试；通常，你会测试那些预先并不知道结果的值，比如读入变量的值，或者函数调用的返回值。

条件测试为True的 if 语句：

if 1 == 1：
　　print（"1 equals 1\n\n"）

它会输出：

1 equals 1

我们进行的测试是 1 == 1，用口语来说就是，“1 是否等于 1”。1 确实等于 1，所以条件测试的结果为 True ，因此和 if 语句相关联的语句就会被执行，信息被打印了出来。

你也可以这样写：

if 1:
    print("1 evaluates to true\n\n")

它会输出：

11evaluates to true

条件测试为 False 的 if 语句：

if 1 == 0 :
    print("1 equals 0\n\n")

它没有任何输出！我们进行的测试是 1 == 0，用口语来说就是，“1 是否等于 0”。 1 当然不等于 0 了，所以条件测试的结果为 False ，因此和 if 语句相关联的语句就不会被执行，不会输出任何信息。你也可写成：

if 0:
    print("0 evaluates to true\n\n")

现在，让我们看一下条件测试为的 if-else 语句：

if 1 == 1:
    print("1 equals 1\n\n")
else:
    print("1 does not equal 1\n\n")

它会输出：

1 equals 1

当测试结果为True时，if-else 会做某件事情，当测试结果为False时，它则会另一件事情。这是条件测试为True的 if-else 语句：

if 1 == 0:
    print("1 equals 0\n\n")
else:
    print("1 does not equal 0\n\n")

它会输出：

1 does not equal 0

条件测试和括号成对

对这些语句和它们的条件测试，还有两点注释：

第一点，在这些语句的条件部分，有许多测试可以使用。除了在前述实例中使用的数字相等 == 外，你还可以对不等 !=、大于 > 和小于 < 等进行测试。

也有一些文件测试操作符，允许你对文件进行相应的测试：是否存在，是否为空，是否设置了特定的权限，等等（参看附录 B）。另一个比较常用的是对变量名进行测试：如果变量存储的是 0，那它就被认为是 False 的；任何其他的数字都被认为是的。如果变量中有非空的字符串，那它就是 True 的；空的字符串用 "" 或者 '' 来表示，它是 False 的。

第二点，注意条件测试之后的语句都缩进在后面，叫做块，这在 Python 中非常常⻅。每一行代码做的事情不要太多，使用缩进让代码块更加突出明显一些（参看第 5.2 节）。

回过头来再看看条件语句。就像例 5.1中演示的那样，if-else 还有 if-elsif-else 这样的形式。第一个 if 和 elsifs 中的条件被交替测试，一旦测试结果为True ，相应的代码块就会

被执行，并且剩余的条件测试也会被忽略掉。假设没有一个条件的测试结果为True，如果存在 else 代码块，那它就会被执行。else 代码块是可有可无的。

例 5.1：if-elsif-else

#!/usr/bin/env python
# if-elif-else
 
word = 'MNIDDKL'
 
# if-elif-else conditionals
if word == 'QSTVSGE':
    print("QSTVSGE\n")
elif word == 'MRQQDMISHDEL':
    print("MRQQDMISHDEL\n")
elif word == 'MNIDDKL':
    print("MNIDDKL--the magic word!\n")
else:
    print('Is \"%s\" a peptide? This program is not sure.\n' % word)
 
exit()

注意 else 代码块的 print 语句中的 \"，它可以在双引号包裹起来的字符串中打印出双引号（"）。反斜线告诉 Python，把后面紧跟的” 当成引号本身，而不是把它当做字符串结束的标志。此外还请注意使用了 == 来检测字符串的相等。例 5.1的输出：

MNIDDKL--the magic word!

5.1.2 循环

循环允许你重复执行被成对大括号包裹起来的语句块。在 Python 中有多种方法实现循环：while 循环、for 循环等等。例 5.2（来源于第 4 章）演示了如何使用 while 循环从一个文件读取蛋白质序列数据。

例 5.2：从文件中读取蛋白质序列数据，第四次尝试

#!/usr/bin/env python
import os
# Reading protein sequence data from a file, take 4
 
# The filename of the file containing the protein sequence data
proteinfilename = 'NM_021964fragment.pep'
 
# First we have to "open" the file, and in case the
# open fails, print an error message and exit the program.
if os.path.exists(proteinfilename):
    PROTEINFILE = open(proteinfilename)
else:
    print("Could not open file %s!\n" % proteinfilename)
    exit()
     
 
# Read the protein sequence data from the file in a "while" loop,
# printing each line as it is read.
protein = PROTEINFILE .readline()
while protein:
    print("  ######  Here is the next line of the file:\n")
    print(protein)
    protein = PROTEINFILE .readline()
 
 
# Close the file.
PROTEINFILE.close()
 
exit()

下面是例 5.2的输出：

  ######  Here is the next line of the file:
MNIDDKLEGLFLKCGGIDEMQSSRTMVVMGGVSGQSTVSGELQD
  ######  Here is the next line of the file:
SVLQDRSMPHQEILAADEVLQESEMRQQDMISHDELMVHEETVKNDEEQMETHERLPQ
  ######  Here is the next line of the file:
GLQYALNVPISVKQEITFTDVSEQLMRDKKQIR

在 while 循环中，注意在循环至文件的下一行时 protein 变量是如何一次次被赋值的。在 Python 中，赋值语句会返回赋予的值。此处，条件测试检测的是赋值语句是否成功得读取了另一行。如果有另一行被读入，就会发生赋值行为，条件测试就为 True，新的一行会被存储在 protein 变量中，而包含两个 print 语句的代码块也会被执行。如果没有更多行了，赋值就会是未定义，条件测试就为 False，而程序也会跳过包含两个 print 语句的代码块，退出 while 循环，继续执行程序的后续部分（在这个例子中，就是 close 和 exit 函数）。

open 和 os

open 是一个系统调用，因为要打开一个文件，Python 必须向操作系统提出请求。操作系统可以是 Unix 或 Linux、Microsoft Windows 或 Apple Maintosh 等的某个版本。文件被操作系统管理者，因此也只有操作系统能够访问它。

检查系统调用的成功与否是一个好习惯，尤其是当打开文件时更是如此。如果系统调用失败了，而你也没有对它进行检查，那么程序将继续运行，可能会尝试读取或写入一个你并没有打开过的文件。你应该总是检查这种失败，当文件无法打开时，要立即告知程序的使用者。当打开文件失败时，通常情况下你会希望立即退出程序，并尝试打开另一个文件。在例 5.2中，使用os模块来判断文件是否存在。

os模块是python内置模块，包含普遍的操作系统功能，与具体的平台无关。例5.2中，import os语句导入os模块后，可以使用os模块的文件判断函数os.path.exists，如果存在，就调用open函数打开文件；反之，如果不存在，程序会输出错误信息，随后退出。

总结一下，在 Python 中，条件和循环都是比较简单的概念，学起来并不困难。它们是编程语言最强大的特性之一。条件可以使你的程序适应不同的状况，使用这种方式，就可以针对不同的输入采取不同的方案了。在人工智能的计算机程序中，它占了相当大的比重。循环充分利用了计算机的速度，利用循环，仅仅使用几行代码，你就可以处理大量的输入，或者对计算进行重复迭代与提炼。

5.2 代码风格

1. 命名风格

总体原则，新编代码必须按下面命名风格进行，现有库的编码尽量保持风格。
尽量以免单独使用小写字母'l'，大写字母'O'，以及大写字母'I'等容易混淆的字母。
模块命名尽量短小，使用全部小写的方式，可以使用下划线。
包命名尽量短小，使用全部小写的方式，不可以使用下划线。
类的命名使用CapWords的方式，模块内部使用的类采用_CapWords的方式。
异常命名使用CapWords+Error后缀的方式。
全局变量尽量只在模块内有效，类似C语言中的static。实现方法有两种，一是__all__机制;二是前缀一个下划线。对于不会发生改变的全局变量，使用大写加下划线。
函数命名使用全部小写的方式，可以使用下划线。
常量命名使用全部大写的方式，可以使用下划线。
使用 has 或 is 前缀命名布尔元素，如: is_connect = True; has_member = False
用复数形式命名序列。如：members = ['user_1', 'user_2']
用显式名称命名字典，如：person_address = {'user_1':'10 road WD', 'user_2' : '20 street huafu'}
避免通用名称。诸如 list, dict, sequence 或者 element 这样的名称应该避免。又如os, sys 这种系统已经存在的名称应该避免。
类的属性（方法和变量）命名使用全部小写的方式，可以使用下划线。
对于基类而言，可以使用一个 Base 或者 Abstract 前缀。如BaseCookie、AbstractGroup
内部使用的类、方法或变量前，需加前缀'_'表明此为内部使用的。虽然如此，但这只是程序员之间的约定而非语法规定，用于警告说明这是一个私有变量，外部类不要去访问它。但实际上，外部类还是可以访问到这个变量。import不会导入以下划线开头的对象。
类的属性若与关键字名字冲突，后缀一下划线，尽量不要使用缩略等其他方式。
双前导下划线用于命名class属性时，会触发名字重整；双前导和后置下划线存在于用户控制的名字空间的"magic"对象或属性。
为避免与子类属性命名冲突，在类的一些属性前，前缀两条下划线。比如：类Foo中声明__a,访问时，只能通过Foo._Foo__a，避免歧义。如果子类也叫Foo，那就无能为力了。
类的方法第一个参数必须是self，而静态方法第一个参数必须是cls。
一般的方法、函数、变量需注意，如非必要，不要连用两个前导和后置的下线线。两个前导下划线会导致变量在解释期间被更名。两个前导下划线会导致函数被理解为特殊函数，比如操作符重载等。

2 关于参数

要用断言来实现静态类型检测。断言可以用于检查参数，但不应仅仅是进行静态类型检测。 Python 是动态类型语言，静态类型检测违背了其设计思想。断言应该用于避免函数不被毫无意义的调用。
不要滥用 *args 和 **kwargs。*args 和 **kwargs 参数可能会破坏函数的健壮性。它们使签名变得模糊，而且代码常常开始在不应该的地方构建小的参数解析器

3 代码编排

缩进。优先使用4个空格的缩进（编辑器都可以完成此功能），其次可使用Tap，但坚决不能混合使用Tap和空格。
每行最大长度79，换行可以使用反斜杠，最好使用圆括号。换行点要在操作符的后边敲回车。
类和top-level函数定义之间空两行；类中的方法定义之间空一行；函数内逻辑无关段落之间空一行；其他地方尽量不要再空行。
一个函数 : 不要超过 30 行代码, 即可显示在一个屏幕类，可以不使用垂直游标即可看到整个函数；一个类 : 不要超过 200 行代码，不要有超过 10 个方法；一个模块不要超过 500 行。

4. 文档编排

模块内容的顺序：模块说明和docstring—import—globals&constants—其他定义。其中import部分，又按标准、三方和自己编写顺序依次排放，之间空一行。
不要在一句import中多个库，比如import os, sys不推荐。
如果采用from XX import XX引用库，可以省略‘module.’。若是可能出现命名冲突，这时就要采用import XX。

5. 空格的使用

总体原则，避免不必要的空格。
各种右括号前不要加空格。
函数的左括号前不要加空格。如Func(1)。
序列的左括号前不要加空格。如list[2]。
逗号、冒号、分号前不要加空格。
操作符（=/+=/-+/==/</>/!=/<>/<=/>=/in/not in/is/is not/and/or/not)左右各加一个空格，不要为了对齐增加空格。如果操作符有优先级的区别，可考虑在低优先级的操作符两边添加空格。如：hypot2 = x*x + y*y; c = (a+b) * (a-b)
函数默认参数使用的赋值符左右省略空格。
不要将多句语句写在同一行，尽管使用‘；’允许。
if/for/while语句中，即使执行语句只有一句，也必须另起一行。

6. 注释

总体原则，错误的注释不如没有注释。所以当一段代码发生变化时，第一件事就是要修改注释！避免无谓的注释
注释必须使用英文，最好是完整的句子，首字母大写，句后要有结束符，结束符后跟两个空格，开始下一句。如果是短语，可以省略结束符。
行注释：在一句代码后加注释，但是这种方式尽量少使用。。比如：x = x + 1 # Increment
块注释：在一段代码前增加的注释。在‘#’后加一空格。段落之间以只有‘#’的行间隔。比如：

# Description : Module config.
#
# Input : None
#
# Output : None

7. 文档描述

为所有的共有模块、函数、类、方法写docstrings；非共有的没有必要，但是可以写注释（在def的下一行）。
如果docstring要换行，参考如下例子,详见PEP 257

"""Return a foobang
 
Optional plotz says to frobnicate the bizbaz first.
 
"""

8. 编码建议

编码中考虑到其他python实现的效率等问题，比如运算符‘+’在CPython（Python）中效率很高，都是Jython中却非常低，所以应该采用.join()的方式。
与None之类的单件比较，尽可能使用‘is’‘is not’，绝对不要使用‘==’，比如if x is not None 要优于if x。
使用startswith() and endswith()代替切片进行序列前缀或后缀的检查。比如：建议使用if foo.startswith('bar'): 而非if foo[:3] == 'bar':
使用isinstance()比较对象的类型。比如：建议使用if isinstance(obj, int): 而非if type(obj) is type(1):
判断序列空或不空，有如下规则：建议使用if [not] seq: 而非if [not] len(seq)
字符串不要以空格收尾。
二进制数据判断使用 if boolvalue的方式。
使用基于类的异常，每个模块或包都有自己的异常类，此异常类继承自Exception。错误型的异常类应添加"Error"后缀，非错误型的异常类无需添加。
异常中不要使用裸露的except，except后跟具体的exceptions。
异常中try的代码尽可能少。比如：

try:
    value = collection[key]
except KeyError:
    return key_not_found(key)
else:
    return handle_value(value)

5.3 查找基序

在生物信息学中我们最常做的事情之一就是查找基序，即特定的 DNA 或蛋白质片段。这些基序可能是 DNA 的调控元件，也可能是已知在多个物种中保守的蛋白质片段。（PROSITE 网站——http://www.expasy.ch/prosite/——有关于蛋白质基序的更多信息。）

在生物学序列中你要查找的基序往往不是特定的一个序列，它们可能有一些变体—— 比如，某个位置上的碱基或者残基是什么并不重要。它们也有可能有不同的⻓度。但通常它们都可以用正则表达式来进行表征，在接下来的例5.3、第 9 章以及本书中的许多地方，你都可以看到更多关于正则表达式的讨论。

Python 有一系列在字符串中进行查找的便利的特性，这同样使得 Python 成为生物信息学领域中流行的语言。例 5.3对这种字符串搜索的能力进行了演示；它做的事情真的非常实用，许多类似的程序在生物学研究中一直被使用着。它做的事情包括：

从文件中读入蛋白质序列数据
为了便于搜索，把所有的序列数据都放到一个字符串中
查找用户通过键盘键入的基序

例 5.3：查找基序

#!/usr/bin/env python
import os
# Searching for motifs
 
# Ask the user for the filename of the file containing
# the protein sequence data, and collect it from the keyboard
print "Please type the filename of the protein sequence data: ";
 
proteinfilename = input().rstrip()
 
# open the file, or exit
if os.path.exists(proteinfilename):
    PROTEINFILE = open(proteinfilename)
else:
    print("Could not open file %s!\n" % proteinfilename)
    exit()
 
# Read the protein sequence data from the file, and store it
# into the array variable proteins
proteins = PROTEINFILE.readlines()
 
# Close the file - we've read all the data into @protein now.
PROTEINFILE.close()
 
# Put the protein sequence data into a single string, as it's easier
# to search for a motif in a string than in an array of
# lines (what if the motif occurs over a line break?)
protein = ''.join(proteins)
 
# Remove whitespace
protein = protein.replace('\n', '')
 
# In a loop, ask the user for a motif, search for the motif,
# and report if it was found.
# Exit if no motif is entered.
 
 
while True:
    print("Enter a motif to search for: ")
     
    motif = input().rstrip()
     
    # exit on an empty user input
    if not motif:
        break
 
    # Look for the motif
 
    if protein.find(motif) != -1:
        print("I found it!\n\n")
    else:
        print("I couldn\'t find it.\n\n")
 
# exit the program
exit()

这是例 5.3典型的一个输出：

Please type the filename of the protein sequence data: NM_021964fragment.pep
Enter a motif to search for: SVLQ
I found it!
 
Enter a motif to search for: jkl
I couldn't find it.
 
Enter a motif to search for: QDSV
I found it!
 
Enter a motif to search for: HERLPQGLQ
I found it!
 
Enter a motif to search for:
I couldn't find it.

就像你在结果中看到的那样，这个程序查找用户通过键盘键入的基序。利用这样一个程序，你就不用再在海量的数据中手工进行查找了。让计算机来做这样的工作，它比人工做的更快、更准确。

如果程序不仅报告它找到了基序，还能告诉我们在哪里找到了它，那就更好了。在第 9 章中你将看到如何来实现这个目标，那一章的一个练习要求你对这个程序进行修改，使得它能够报告基序的位置。

接下来的小节将对例 5.3中这些新的内容进行介绍与讨论：

获取用户的键盘输入
把文件的所有行连接成一个字符串
find和rstrip函数
not

5.3.1 获取用户的键盘输入

在例 4.5中你第一次遇到文件句柄。在例 5.3中（在例 4.3中也是这样），使用文件句柄和readlines函数从打开的文件中读入数据并保存到列表中，就像这样：

proteins = PROTEINFILE.readlines()

Python 使用同样的语法获取用户的键盘输入。在例 5.3中，为了获取用户的键盘输入，使用了一个名为input的函数接受用户输入数据，返回为 string 类型。

proteinfilename = input()

文件句柄可以和一个文件相关联，也可以和用户回答程序问题时的键盘输入相关联。

在例 5.3中，用户被要求输入一个包含蛋白质序列数据的文件的文件名。在通过这种方式得到文件名后，打开文件之前还需要一步操作。当用户键入文件名并通过 Enter 键

（也叫做 Return 键）换行时，文件名和其末尾的换行符一起被保存进了变量中。这个换行符并不是文件名的一部分，所以在使用 open 系统调用之前应该把它去除掉。Python函数 rstrip() 会去除掉字符串末尾的换行符 newlines（或者它的表亲换行符 linefeeds 和回⻋符 carriage returns）。所以这一个部分还有一点附带的知识：删除获取的用户键盘输入中的换行符。试着不调用 rstrip 函数，你会看到 open 调用失败，因为并没有末尾包含换行符的文件名。（对于文件名中可以使用的字符，操作系统有自己的规定。）

5.3.2 使用 join 把列表合并成标量

常常可以发现，蛋白质序列数据被打断成了短的片段，每个片段有大约 80 个字符。

原因很简单：当要把数据打印到纸张上或展示在屏幕上时，根据纸张或屏幕的空间需要把它打断成数行。你的数据被打断成了片段，但对你的 Python 程序来说却多有不便。当你要查找一个基序，而它却被换行符分割开了，这会怎样呢？你的程序将无法找到这个基序。

事实上，例 5.3中查找的一些基序确实被换行符分割开了。在 Python 中，你可以使用 join 函数来处理类似的片段数据。在例 5.3中，join 把存储在 protein 列表中的数据的所有行合并成一个单独的字符串，并把它保存在了新的标量变量 protein 中：

protein = ''.join(proteins)

当合并列表时，你需要指定一个字符串，这个字符串会放在列表的各个元素之间。在这个例子中，你指定的是空字符串，它会放在输入文件的各行之间。空字符串就夹在成对的单引号 '' 中（当然也可以使用双引号 ""）。

回忆一下例 4.2，在那个例子中，我介绍了好几种连接两个 DNA 片段的方法。join 函数的作用与之类似，它取出列表元素的变量值并把它们连接成一个单独的标量值。回忆一下例 4.2中的下列语句，它是连接两个字符串的方法之一：

DNA3 = DNA1 + DNA2

另一种实现同样连接的方法是使用 join 函数：

DNA3 = ''.join([DNA1, DNA2])

在这个例子中，我没有给出列表名，而是指定了一个标量元素的列表：

[DNA1, DNA2]

5.3.3python中实现do-until类似循环例子5-3中，使用while和if判断实现了类似do-until功能的循环。它首先执行一个块，然后进行if判断，如果测试成功，就用break跳出循环。例子5-3首先打印用户提示信息，获取用户输入，调用find函数搜索模体并报告结果（是否为-1，-1意味着没有找到）。在执行查找操作之前，用if语句判断用户是否输入的是空行，空行意味着用户用户没有更多的模体序列要查找，退出循环。
5.3.4字符串函数find和索引Python可以轻松操作各种字符串，例如DNA和蛋白质序列数据。字符串内含函数find用来搜索子字符串（模体序列）是否出现在字符串（蛋白质序列）中，如果找到则返回子字符串的起始索引，否则返回-1。
5.4 计数核苷酸对于一段 DNA 序列，有许多信息你可能想知道。它是编码的还是非编码的？它是否含有调控元件？它是否和其他已知的 DNA 相关，如果相关，是怎么相关呢？DNA 中四种核苷酸的数目各是多少？事实上，在许多物种中，编码区域都有特定的核苷酸偏向性，所以，最后这个问题对于基因识别来说是非常重要的。此外，不同的物种有不同的核苷酸使用模式。所以对核苷酸进行计数不但有趣而且非常有用。
接下来的小节中有两个程序，例 5.4和例 5.6，它们对 DNA 序列中每种核苷酸都进行计数。这两个程序展示了 Python 的一些新的知识：
“拆解”字符串
查看字符串的特定位置
对列表进行迭代
对字符串的⻓度进行迭代
为了对 DNA 中的每一种核苷酸进行计数，你需要查看每一个碱基，看看它是什么，并且要为每一种碱基记录一个计数，一共四个计数。我们将通过两种方法来实现：
把 DNA 拆解成单个碱基，存储到列表中，然后对列表进行迭代（换言之，一个接一个的对列表元素进行处理）
首先，让我们从该任务的伪代码开始。之后，我们会给出更加详细的伪代码，最终，编写出这两种方案的 Python 程序。
下面的伪代码描述了基本的要求：

for each base in the DNA
    if base is A
        count_of_A = count_of_A + 1
    if base is C
        count_of_C = count_of_C + 1
    if base is G
        count_of_G = count_of_G + 1
    if base is T
        count_of_T = count_of_T + 1
done
 
print count_of_A, count_of_C, count_of_G, count_of_T

就像你看到的一样，这是一个非常简单的想法，它模拟了你手工计数时的工作。（如果你想计算所有人类基因中碱基的相对频率，手工计数就不现实了，因为数量太过巨大，你将不得不使用类似的程序。这就是生物信息学。）现在，让我们看看如何用 Python 来实现它吧。
5.5把字符串拆解成列表假设你打算把 DNA 字符串拆解成列表。所谓拆解，我指的是字符串中的每一个字⺟ 分离开来，就像把字符串分离成 bit 一样。换句话说，DNA 字符串中表示碱基的字⺟都被分离开，每一个字⺟都将成为数据中的标量值。然后一个接一个的查看列表元素（每一个都是一个单独的字符），这样你就可以进行计数了。这与第 5.3.2 小节中的 join 函数是相反的，它把列表中的字符串合并成单个的标量值。（在把字符串拆解成列表后，如果你愿意，可以使用 join 再把这个列表合并成和原来一模一样的字符串。）
在刚才的伪代码中，我将加上下面的这些操作指令：从文件中获取 DNA，操作文件数据使 DNA 序列成为单一的字符串。所以，首先，你把包含原始文件数据行的列表合并起来，通过删除空白把它清理干净，直到只有序列保留下来为止，然后，把它拆解成列表。
当然，关键的一点就是最后的这个列表就是我们所需要的，使用它可以在循环中方便地进行计数。与包含行、换行符和其他可能的无用的字符的列表不同，这个列表就是包含单个碱基的列表。

read in the DNA from a file
 
join the lines of the file into a single string DNA
 
# make an array out of the bases of DNA
DNA = list(DNA)
 
# initialize the counts
count_of_A = 0
count_of_C = 0
count_of_G = 0
count_of_T = 0
 
for each base in DNA
 
    if base is A
        count_of_A = count_of_A + 1
    if base is C
        count_of_C = count_of_C + 1
    if base is G
        count_of_G = count_of_G + 1
    if base is T
        count_of_T = count_of_T + 1
done
 
print count_of_A, count_of_C, count_of_G, count_of_T

如前所述，这里的伪代码更加详细一些。通过把 DNA 字符串拆解成包含单个字符的数组，它提供了一种查看每个碱基的办法。为了保证计数正确，它还把所有的计数都初始化为 0。如果你理解了程序中的初始化，就很容易看出具体发生了什么，这可以防止代码中出现某些错误。
我们已经设计好了程序，现在就把它转化成 Python 代码。例 5.4是一个可以工作的程序。随着本章的学习，你将看到其他完成该任务的方法，它的速度会更快一些，但此处速度并不是我们关注的焦点。
例 5.4：确定核苷酸的频率

#!/usr/bin/env python<br>import os
# Determining frequency of nucleotides
 
# Get the name of the file with the DNA sequence data
print("Please type the filename of the DNA sequence data: ")
 
dna_filename = input()
 
 
# open the file, or exit
if os.path.exists(dna_filename):
    DNAFILE = open(dna_filename)
else:
    print("Could not open file %s!\n" % dna_filename)
    exit()

# Read the DNA sequence data from the file, and store it
# into the array variable DNAs
DNAs = DNAFILE.readlines()
 
# Close the file
DNAFILE.close()
 
# From the lines of the DNA file,
# put the DNA sequence data into a single string.
DNA = ''.join(DNAs)
 
# Remove whitespace
DNA = DNA.replace('\n', '')
 
# Now explode the DNA into an array where each letter of the
# original string is now an element in the array.
# This will make it easy to look at each position.
# Notice that we're reusing the variable DNA for this purpose.
DNA = list(DNA)
 
# Initialize the counts.
# Notice that we can use scalar variables to hold numbers.
count_of_A = 0
count_of_C = 0
count_of_G = 0
count_of_T = 0
errors     = 0
 
# In a loop, look at each base in turn, determine which of the
# four types of nucleotides it is, and increment the
# appropriate count.
for base in DNA:
 
    if base == 'A':
        ++count_of_A
    elif base == 'C':
        ++count_of_C
    elif base == 'G':
        ++count_of_G
    elif base == 'T':
        ++count_of_T
    else:
        print("!!!!!!!! Error - I don\'t recognize this base: %s\n" % base)
        ++errors
 
# print the results
print("A = %s\n" % count_of_A)
print("C = %s\n" % count_of_C)
print("G = %s\n" % count_of_G)
print("T = %s\n" % count_of_T)
print("errors = %s\n" % errors)
 
# exit the program
exit()

为了演示例 5.4，我创建了下面这个小的 DNA 文件，把它命名为 small.dna：

AAAAAAAAAAAAAAGGGGGGGTTTTCCCCCCCC
CCCCCGTCGTAGTAAAGTATGCAGTAGCVG
CCCCCCCCCCGGGGGGGGAAAAAAAAAAAAAAATTTTTTAT
AAACG

你可以使用钟爱的文本编辑器在计算机上键入 small.dna 文件，也可以从本书的网站上下载它。注意文件中有一个 V，这是一个错误。下面是例 5.4的输出：

Please type the filename of the DNA sequence data: small.dna
!!!!!!!! Error - I don't recognize this base: V
 
A = 40
C = 27
G = 24
T = 17

现在我们来看看程序中出现的新的知识。打开和读取序列数据与前面的程序是一样的。第一个新知识出现在这一行中：

DNA = list(DNA)

list函数将DNA字符串分割成了单个字母组成的列表。在交互式环境输入help（list）或使用文档查看list函数使用方法，list函数接受一个可迭代的对象，将其转换为列表。
在Python中，字符串就是一个可迭代对象，并且可迭代对象是可以使用for循环。因此，上述程序可去掉“DNA = list(DNA)”，写成如下方式：

#!/usr/bin/env python
import os
# Determining frequency of nucleotides
 
# Get the name of the file with the DNA sequence data
print("Please type the filename of the DNA sequence data: ")
 
dna_filename = input()
 
 
# open the file, or exit
if os.path.exists(dna_filename):
　　DNAFILE = open(dna_filename)
else:
　　print("Could not open file %s!\n" % dna_filename)
　　exit()
 
 
# Read the DNA sequence data from the file, and store it
# into the array variable DNAs
DNAs = DNAFILE.readlines()
 
# Close the file
DNAFILE.close()
 
# From the lines of the DNA file,
# put the DNA sequence data into a single string.
DNA = ''.join(DNAs)
 
# Remove whitespace
DNA = DNA.replace('\n', '')
 
 
# Initialize the counts.
# Notice that we can use scalar variables to hold numbers.
count_of_A = 0
count_of_C = 0
count_of_G = 0
count_of_T = 0
errors     = 0
 
# In a loop, look at each base in turn, determine which of the
# four types of nucleotides it is, and increment the
# appropriate count.
for base in DNA:
 
    if base == 'A':
        ++count_of_A
    elif base == 'C':
        ++count_of_C
    elif base == 'G':
        ++count_of_G
    elif base == 'T':
        ++count_of_T
    else:
        print("!!!!!!!! Error - I don\'t recognize this base: %s\n" % base)
        ++errors
 
# print the results
print("A = %s\n" % count_of_A)
print("C = %s\n" % count_of_C)
print("G = %s\n" % count_of_G)
print("T = %s\n" % count_of_T)
print("errors = %s\n" % errors)
 
# exit the program
exit()

5.6 操作字符串要查看每一个字符，其实也没必要把字符串拆解成数组。事实上，有时你会想避免那种做法。一个大的字符串会占用你计算机的大量内存，对于大的数组也是如此。当你把字符串拆解成数组时，原始的字符串仍然存在着，而你还要对每一个字符制作一个拷⻉，作为元素存储到创建的数组中。如果你有一个大的字符串，它已经占用了一大部分可用内存，创建另一个数组可能会导致内存不足。当内存不足时，计算机的性能会大打折扣，它会慢的像乌⻳一样，甚至崩溃或冻结（“挂起”）。到现在为止，这些都还不是什么令人担忧的因素，但当你处理大的数据集（比如人类基因组）时，你将不得不考虑这些问题。
那么假设你不想制作一个 DNA 序列数据的拷⻉存进另一个变量。有没有什么方法可以直接查看 DNA 字符串并对碱基进行计数吗？答案是肯定的。下面是伪代码，紧随其后的是 Python 程序：

read in the DNA from a file
 
join the lines of the file into a single string of $DNA
 
# initialize the counts
count_of_A = 0
count_of_C = 0
count_of_G = 0
count_of_T = 0
 
for each base at each position in DNA
 
    if base is A
        count_of_A = count_of_A + 1
    if base is C
        count_of_C = count_of_C + 1
    if base is G
        count_of_G = count_of_G + 1
    if base is T
        count_of_T = count_of_T + 1
done
 
print count_of_A, count_of_C, count_of_G, count_of_T

例5.6是查看 DNA 字符串中每个碱基的程序。
例 5.6：确定核苷酸频率，第二次尝试

#!/usr/bin/env python
import os
# Determining frequency of nucleotides
 
# Get the name of the file with the DNA sequence data
print("Please type the filename of the DNA sequence data: ")
 
dna_filename = input()
 
 
# open the file, or exit
if os.path.exists(dna_filename):
　　DNAFILE = open(dna_filename)
else:
　　print("Could not open file %s!\n" % dna_filename)
　　exit()
 
 
# Read the DNA sequence data from the file, and store it
# into the array variable DNAs
DNAs = DNAFILE.readlines()
 
# Close the file
DNAFILE.close()
 
# From the lines of the DNA file,
# put the DNA sequence data into a single string.
DNA = ''.join(DNAs)
 
# Remove whitespace
DNA = DNA.replace('\n', '')
 
 
# Initialize the counts.
# Notice that we can use scalar variables to hold numbers.
count_of_A = 0
count_of_C = 0
count_of_G = 0
count_of_T = 0
errors     = 0
 
# In a loop, look at each base in turn, determine which of the
# four types of nucleotides it is, and increment the
# appropriate count.
for position in range(len(DNAs)):
    base = DNAs[position]
 
    if base == 'A':
        ++count_of_A
    elif base == 'C':
        ++count_of_C
    elif base == 'G':
        ++count_of_G
    elif base == 'T':
        ++count_of_T
    else:
        print("!!!!!!!! Error - I don\'t recognize this base: %s\n" % base)
        ++errors
 
# print the results
print("A = %s\n" % count_of_A)
print("C = %s\n" % count_of_C)
print("G = %s\n" % count_of_G)
print("T = %s\n" % count_of_T)
print("errors = %s\n" % errors)
 
# exit the program
exit()

下面是例 5.6的输出：

Please type the filename of the DNA sequence data: small.dna
!!!!!!!! Error - I don't recognize this vase: V
A = 40
C = 27
G = 24
T = 17
errors = 1

在例 5.6中，我添加了一行代码，来看看文件是否存在：

if os.path.exists(dna_filename)

注意文件有许多属性，比如大小、权限、文件系统中的位置、文件类型等，对于许多这样的属性都可以调用os模块中的函数来实现。
此外，注意我保留了关于for循环的详细注释，此处的少许注释可以帮助读者跳过代码。

for position in range(len(DNAs)):

这里的 for 循环和下面的 while 循环是完全等价的：

position = 0
 
while position < len(DNAs):

    ++position

比较这两个循环，不难发现for循环只是while循环的简写，循环判断条件是position是否小于字符串DNA的长度，这里使用了字符串索引方式。默认情况下，Python假定字符串从0索引开始，其最后一个字符编号为字符串长度减一。
在上述程序中，我们还使用了python内置函数range，这是一个生成器函数（生成迭代器或列表），例如，range（5）等效生成[0, 1, 2, 3, 4]五个元素的列表。
5.7 写入文件
例 5.7演示了对 DNA 字符串中的核苷酸进行计数的另一种方法。它使用了一个 Python 的技巧，python字符串内置计数函数count。在 while 循环的测试中，它进行了一个全局的碱基计数，就像你将看到的，这是一种对字符串中字符进行计数的简洁的方法。
Python 比较好的一点就是，对于那些相当常⻅的事情，它很可能提供了一种相对简洁的处理方法。（而它的弊端就是对于 Python，有大量知识需要你去学习。）
例 5.7的结果，不仅会输出打印到屏幕，还会写入到文件中。实现写入文件的代码如下：

# Also write the results to a file called "countbase"

outputfile = "countbase";

COUNTBASE = open(outputfile, "w")

print("Cannot open file \"%s\" to write to!!\n\n" % outputfile)

print("A=%s C=%s G=%s T=%s errors=%s\n" % (a, c, g, t, e), file=COUNTBASE)

COUNTBASE.close()

exit()

正如你所⻅，要写入到文件，你要像读取文件那样调用 open，但有一点不同：在函数中使mode参数为"w"。文件句柄成了 print 语句的file参数，这使得 print 语句把它的输出定向到了文件。
例 5.7检查 DNA 字符串中每个碱基的第三个版本的 Python 程序。
例 5.7：确定核苷酸频率，第三次尝试

#!/usr/bin/env python
import os
# Determining frequency of nucleotides
 
# Get the name of the file with the DNA sequence data
print("Please type the filename of the DNA sequence data: ")
 
dna_filename = input()
 
 
# open the file, or exit
if os.path.exists(dna_filename):
　　DNAFILE = open(dna_filename)
else:
　　print("Could not open file %s!\n" % dna_filename)
　　exit()
 
 
# Read the DNA sequence data from the file, and store it
# into the array variable DNAs
DNAs = DNAFILE.readlines()
 
# Close the file
DNAFILE.close()
 
# From the lines of the DNA file,
# put the DNA sequence data into a single string.
DNA = ''.join(DNAs)
 
# Remove whitespace
DNA = DNA.replace('\n', '')
 
 
# In a loop, look at each base in turn, determine which of the
# four types of nucleotides it is, and increment the
# appropriate count.
count_of_A  = DNA.count('A') + DNA.count('a')
count_of_C  = DNA.count('C') + DNA.count('c')
count_of_G  = DNA.count('G') + DNA.count('g')
count_of_T  = DNA.count('T') + DNA.count('t')
errors     = len(DNA) - count_of_A  -count_of_C - count_of_G - count_of_T
 
# print the results
print("A=%d C=%d G=%d T=%d errors=%d\n" % (count_of_A, count_of_C, count_of_G, count_of_T, errors))
 
# Also write the results to a file called "countbase"
outputfile = "countbase"
 
COUNTBASE  = open(outputfile, 'w')
 
COUNTBASE.write("A=%d C=%d G=%d T=%d errors=%d\n" % (count_of_A, count_of_C, count_of_G, count_of_T, errors))
 
COUNTBASE.close()
# exit the program
exit()

当你运行例 5.7时，会看到类似的输出：

Please type the filename of the DNA sequence data: small.dna
A=40 C=27 G=24 T=17 errors=1

在你运行例 5.7后，计数碱基的输出文件中包含以下内容：

A=40 C=27 G=24 T=17 errors=1

5.8   练习题习题 5.1
使用循环写一个永不挂起的程序。每次迭代时，条件测试都应该永远为 True。注意，有的系统会注意到你在一个无穷的循环中，从而自动结束程序。依据你使用的操作系统的不同，结束程序的方法也会有所不同。在 Unix 和 Linux、Windows MS-DOS 命令窗口或者 MacOS X shell 窗口中，使用 Ctrl-C 就可以结束程序。
习题 5.2
提示用户输入两个（短的）DNA 字符串。通过使用 .= 赋值操作符，把第二个附加到第一个 DNA 字符串的后面，从而把两者连接起来。先打印输出连接后的两个字符串，然后打印输出第二个字符串，但要与连接后字符串末尾对应的字符串对⻬。举个例子，如果输入的字符串是 AAAA 和 TTTT，则要打印输出：
AAAATTTT
         TTTT
习题 5.3
编写一个程序，输出从 1 到 100 的所有数字。当然，你的程序代码应该远远小于100 行。
习题 5.4
编写一个程序，计算 DNA 链的反向互补链。不要使用 s///或者 tr 函数。使用 substr 函数，当你计算反向互补链时，要一次检查原始 DNA 中的一个碱基。（提示：你可能会发现检查原始 DNA 字符串时，尽管从左到右也是可以的，但从右向左要更加方便一些)
习题 5.5
编写一个程序，报告蛋白质序列中疏水氨基酸的百分比。（想要知道哪些氨基酸是疏水性的，可以参看任何关于蛋白质、分子生物学或细胞生物学导论性的介绍。在附录 A中你会找到一些有用的资源。）
习题 5.6
编写一个程序，检查一下作为参数提供的两个 DNA 字符串是不是反向互补的。
使用 Python 的内置函数 insert、pop、reverse和 ==（==实际上是一个操作符）。
习题 5.7
编写一个程序，报告序列的 GC 含量。（换言之，就是给出 DNA 中 G 和 C 的百分比）。
习题 5.8
修改例 5.3，使它不仅能通过正则表达式找到基序，还能打印输出它找到的基序。
习题 5.9
编写一个程序，交换 DNA 字符串中特定位置的两个碱基。
习题 5.10                                                                                                ·           ·
编写一个程序，写入一个临时文件然后把它删除掉。unlink 函数可以删除一个文件，举个例子，这样来使用：
1 unlink "tmpfile";
但要检查一下看看 unlink 是否成功了。

posted @ 2018-08-12 20:37 青蛙快飞阅读(526) 评论(0) 编辑收藏举报

刷新页面返回顶部

加载中......

生信小栈

亚恒的生信学习笔记

第 5 章基序和循环

5.1 流程控制

5.1.1 条件语句

5.1.2 循环

5.2 代码风格

5.3 查找基序

5.3.1 获取用户的键盘输入

5.3.2 使用 join 把列表合并成标量

5.3.3python中实现do-until类似循环

5.3.4字符串函数find和索引

5.4 计数核苷酸

5.5把字符串拆解成列表

5.6 操作字符串

5.8 练习题

公告

加载中......

生信小栈

亚恒的生信学习笔记

第 5 章 基序和循环

5.1 流程控制

5.1.1 条件语句

5.1.2 循环

5.2 代码风格

5.3 查找基序

5.3.1 获取用户的键盘输入

5.3.2 使用 join 把列表合并成标量

5.3.3python中实现do-until类似循环

5.3.4字符串函数find和索引

5.4 计数核苷酸

5.5把字符串拆解成列表

5.6 操作字符串

5.8 练习题

公告

第 5 章基序和循环