提取序列的ID
本人是一枚生物学的学生，由于对生物信息学特别感兴趣，于是想自学生物信息学（新手莫怪）。了解到生物信息学要有编程基础，尤其是要会一门编程语言，例如：R语言、Python、Perl等，还要熟悉Linux系统，作为生信小白，听说Python挺简单的，于是就自学了Python，花了两天时间了解了Python的基础语法后，今天想做个练习题试试手（实践是检验真理的唯一标准），下面是练习题：(试题来源：《PYTHON生物信息学数据管理》)
**1 读取 FASTA 格式序列文件，并只将序列标题写到一个新文件中 **
FASTA 文件中的每个序列记录由两部分组成: 一个标题行和多个 64 个字符长的序列 行(这个数字取决于序列长度〉。 标题行在行首位置用大于符号(> )标记，用户可以用它来区 分标题行和序列行。
由于这部分比较简单，就直接上注释完的脚本：（这个脚本可用于提取序列的ID）
fasta_file = open('C:/GooF.fasta','r') #打开要读取的文件
out_file = open('C:/GooF_ID.txt','w') #要输出的文件
for line in fasta_file: #分行读取
if line[0:1] == '>': #读取以>开头的行，并匹配
out_file.write(line) #将匹配的结果按行输出
print(line) #输出结果
out_file.close() #关闭文件夹
fasta_file.close()
如何从多序列 FASTA 文件中提取登记码的列表
Python 脚本所用的输入数据文件是如何产生的呢?考虑 SwissProtSeq.fasta，这是 FASTA 文件的一种形式 。 登记码 (Accession Code)可以在类似下面的标题行中获取，两边附上了管道符号" I " :
下面这个简单的脚本从文件的每个标题行提取登记码，将其添加到列表中，并打印列表:（我使用的是别的数据库中的fasta文件，类似 SwissProtSeq.fasta文件）
input_file = open("C:/Sorbi1_aa.fasta","r" )
out_file = open('C:/SB_ID.txt','w') #要输出的文件
ac_list = []
for line in input_file:
if line [0] == '>':
fields = line.split('|')
out_file.write(fields[2]+'\n')
ac_list.append(fields[2])
print(ac_list)
out_file.close()
输出结果为：
745699', '4745700', '4745701', '4745702', '4745703', '4745704', '4745705',
.........
Process finished with exit code 0
注意，列表使用 appendO方法填充之前，已被初始化为空列表。为了从每个标题行中提取登 记码(AC) ，先要通过一个 if 条件(此处使用了 FASTA 文件标题行由一个初始"> "符号来标记的事实) ，然后用字符串对象的 splitO方法将字符串切分，它返回的是一个 Python 列表， 由括号中的参数(在这里是"1"，作为分隔符)所分隔的子串构成。 像通常一样，这一方法与 它的对象通过点连接。
这个例子展示了一件重要的事情:如果想程式化地解析记录文件(任何记录文件) ，就必 须首先分析文件格式和结构，并找出一个"窍门"(如，一个反复出现的元素，它标志着专用 行或特定的记录，或者一个与特定文字相关的符号，一个合适的能将一行拆成几列的分隔符 等等)来选择性地提取相关的信息。
日常结尾：
虽然这是个小小的计算程序，但对于初学者的我来说每一次对原代码的升级改造，哪怕是读懂后的注释都感觉是一次进步提升，总之代码虽小，动手最重要！希望更多学习Python的爱好者不要像我一样眼高手低，学习编程就是要，思考，敲码，思考，敲码，敲码，再敲码！
网友评论