提取FASTA文件描述行中含有特定字符记录的python代码实现

这个代码是提取FASTA文件描述行中含有特定字符的记录。

这个脚本代码量很少,但是需要动一下脑子,难点在于序列部分的行数是不固定的。

使用软件

  • python 3.8 及以上版本

输入输出文件格式

输入文件为一个FASTA文件,示例数据见 1_fa.txt。输出数据为提取特定记录后生成的FASTA文件,示例结果文件见 out_fa.txt

FASTA 可以存储任何的序列数据,比如参考基因组文件,蛋白质序列,coding DNA sequence (CDS) , 转录序列等。

FASTA 文件的每条记录包含两个部分:一个是描述部分(> 开始,包含序列名称和其他信息),一个是序列数据部分(第二行直至到描述行,可占据多行)。

FASTA 文件的灵活性导致它的结构定义很宽泛,结构不严谨,运行程序可能会遭遇未知的错误。

举个例子,示例数据的第一条记录见下,第一行是描述行,第二行和第三行是序列(应该是RNA序列)

1
2
3
>cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA

运行代码

代码文件见:tiao_Sus.py

将输入文件和本程序放在同一文件夹下,运行命令示范如下:

1
python tiao_Sus.py 1_fa.txt Sus out_fa.txt

参数说明:

1_fa.txt: 第一个参数为输入 FASTA 文件名称

Sus: 第二个参数为提取的描述部分字段

out_fa.txt: 第三个参数为输出文件名称

代码说明

输入参数,打开文件

导入 sys 模块,输入参数

1
2
3
4
5
6
7
#coding=utf-8
#挑不规则行数据

import sys
in_file = open(sys.argv[1],'r')
field = sys.argv[2]
out_file = open(sys.argv[3],'w')

遍历输入文件

先创建一个空字典,然后遍历输入FASTA文件,字典 dick 的键是行索引 (从0开始),值为该行的内容

1
2
3
4
5
6
7
dick = {}

j=0
for i in in_file:
dick[j] = i #行:内容 从0开始
j+=1

遍历字典

这里是重点。遍历字典 dick ,这里其实是 j 的递增来实现文件内容的换行;j 每增加1,相当于文件跳到下一行。通过第一个while 循环 (while j < len(dick)) 保证迭代到文件的最后一行。之后就是一个判断,需要同时满足两个条件:> 在这一行中并且需要查找的特定字符也在这一行中。如果满足这个条件,将这一行内容写入到输出文件中,跳到下一行,再通过 while 循环将后面跟着的序列行也写入到输出文件中(这里用 while 循环是解决序列行数不固定的问题)。如果不满足上面的判断,那么直接跳到下一行。

上面说了一堆逻辑,有点啰嗦,其实你看着输入文件和脚本,在脑子里运行一遍就清楚了。我写这个脚本的时候,也是在脑子里运行一遍才明确没有问题的。

1
2
3
4
5
6
7
8
9
10
11
j=0 
while j < len(dick): #不是<=
if ('>' in dick[j]) and (field in dick[j]):
out_file.write(dick[j])
j+=1
while ('>' not in dick[j]):
out_file.write(dick[j])
j+=1
else:
j+=1

关闭文件

1
2
in_file.close()
out_file.close()

闲话

这是几年前给涛哥写的小代码,虽然代码很短,但是确实花了心思。

  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2019-2022 Vincere Zhou
  • 访问人数: | 浏览次数:

请我喝杯茶吧~

支付宝
微信