用Python解析公共数据库中生物序列的脚本

Question

大家好，

我现在正在学习生物信息学课程，这是我生物医学学位的一部分（我基本上是个Python新手），这次编程作业需要完成以下任务：

提取“基序序列”（氨基酸序列，简单来说就是程序中的字符串），这些序列是从实现多序列比对的算法中提取出来的，之后又通过反复扫描数据库来生成最佳保守序列。最终的目的是从这些“基序”中推断出它们的功能意义。

这些基序存储在一个公共数据库中，文件里有多个数据字段对应每个蛋白质（比如uniprot ID、登录号，还有存储在超链接 .seq 文件中的比对结果），目前我关注的是其中一个字段。这个数据字段叫做“提取的基序集合”。

我想知道如何编写一个脚本，基本上就是解析这些“基序字符串”，并将它们输出到一个文件中。我现在写的脚本大致如下（我还没有将结果写入文件）：

import os, re, sys, string 

printsdb = open('/users/spyros/folder1/python/PRINTSmotifs/prints41_1.kdat', 'r')

protname = None  
final_motifs = []

for line in printsdb.readlines():
 if line.startswith('gc;'):
        protname = line.lstrip()    
        #string.lower(name)  # convert to lowercase
        break

def extract_final_motifs(protname):

"""Extracts the sequences of the 'final motifs sets' for a PRINTS entry.
Sequences are on lines starting 'fd;' A simple regex is used for retrieval"""

for line in printsdb.readlines():
        if line.startswith('fd;'):
                final_motifs = re.compile('^\s+([A-Z]+)\s+<')
                final_motifs = final_motifs.match(line)
                #print(final_motifs.groups()[0])
                motif_dict = {protname : final_motifs}
                break 
return 

motif_dict = extract_final_motifs('ADENOSINER')
print(motif_dict)  

现在的问题是，我的代码在处理一个原始数据库文件（prints41_!.kdat）时，输出的结果在Python命令行上只是“none”，而不是应该生成像[AAYIGIEVLI, AAYIGIEVLI, AAYIGIEVLI, 等等..]这样的列表。

有没有人知道逻辑错误出在哪里？任何建议都非常感谢！我为文字有点多而感到抱歉，我只是希望能尽量清楚。提前感谢任何帮助！

Answer 1

你有没有注意到你的函数没有返回任何值? 这就是为什么会出现None的原因。

Answer 2

首先，你的做法基本上是对的，但你需要把第2行的 "extracted motif sets" 改成一个变量，比如叫 line。这个 for 循环的作用是逐行读取文件中的数据，每次读取的内容会存储在你指定的变量里，这里就是 line。接下来要弄清楚的是 lysozyme.seq 文件的格式。听起来这个文件里的数据字段之间没有空格。如果是这样的话，你可以用 line.split(" ") 或者 line.split("\t") 来分割数据，其中 \t 代表制表符。这个分割操作会按照你在程序里写的内容，把字符串每次遇到 " " 或 "\t" 的地方切开。

找到文件的目录应该不难，可能这里有一些相关的问题可以参考。

如果你能把数据或者文件中的一部分发给我们看看，我们可能能帮你解析它 :).