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

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

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

使用软件

• python 3.8 及以上版本

输入输出文件格式

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

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

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

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

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

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

运行代码

代码文件见：tiao_Sus.py

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

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

参数说明：

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

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

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

代码说明

输入参数，打开文件

导入 sys 模块，输入参数

#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开始)，值为该行的内容

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 循环是解决序列行数不固定的问题）。如果不满足上面的判断，那么直接跳到下一行。

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

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

关闭文件

in_file.close()
out_file.close()

闲话

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