从FASTA文件中批量提取指定序列【Python脚本】

前言

背景：学测序流程的时候，做到mapping的时牛的基因组有两个多G，因为是在个人PC上初步学习，建立index实在太慢了，而且临时也没有现成的index。于是决定只挑基因组前十条染色体拿来练习（所以需要从基因组文件里选取序列，尝试自己用python写脚本处理）。自己的python学习之前因为脑子实在不聪明也就刚学到字典部分，借鉴了网上的一些内容整理了一下，写了几种情况关于利用python脚本提取fasta基因序列。

一：读取含特定字符的序列并输出

比如现在有一个从NCBI上下载的fasta格式的文件，里面有很多条序列，我现在需要从中选取序列名中含有“XXX”的这些序列，那就可以用我在网上找到并改进的这个脚本

举个例子（用法中的三个参数）：
[fasta_file]下载的数据源fasta文件

[namelist_file] 把我们要找的名字另建一个txt文件，一行一个

[outfile_name] 输出的结果文件，自己命名

# -*- coding: utf-8 -*-

import sys
import time
start=time.clock()

def usage():
    print('Usage: python script.py [fasta_file] [namelist_file] [outfile_name]')


def main():
    outf = open(sys.argv[3],'w')
    dict = {}
    with open(sys.argv[1], 'r') as fastaf:
        for line in fastaf:
            if line.startswith('>'):
                name = line
                dict[name] = ''
            else:
                dict[name] += line.replace('\n','') #读取整个fasta文件构成字典
         

    with open(sys.argv[2],'r') as listf:
        for row in listf:
            row = row.strip()
            for key in dict.keys(): #选取包含所需名称的基因名和序列
                if row in key:
                    outf.write(key)
                    outf.write(dict[key] + '\n')
    outf.close()

try:
    main()
except IndexError:
    usage()

end=time.clock()
print 'Running time: %s Seconds'%(end-start) #运行时间

这个脚本不能用来处理大文件！！！

演示

[fasta_file]

这是需要被筛选读取的源fasta文件示例

[namelist_file]

这是事先准备的包含要选出来的序列名的txt文件（每个名字独立一行），比如这里我要选取在上面那个test文件中（也就是[fasta_file]）基因名包含05/01/04的所有序列

然后，运行该脚本，我是在Linux下（也可以直接在Windows下的python环境，我这里Linux用的python 2.7）

被读取的[fasta_file]也就是我测试的test.fasta大小为4.5M（这个脚本不能用来处理大文件！！！），这里实测耗时0.047秒，然后会生成下面的结果文件

[outfile_name]

这是我输出的结果（即[outfile_name]）

二：读到某一个字符之前的全部输出

背景：找到上面那个方法之后用了一下，发现对于这种几个G的大文件根本行不通，先要读取变成字典，再从字典里去找，太耗时了（我花了几个小时也没运行完得到结果）。所以说不能被别人的脚本限制了思维，想了想我反正只要前十条，那读到第十一条染色体是时候停止读写不就行了？于是自己重写了个

同样，这个也适用于**“输出某个字符之前的所有内容”**这一场景，而在这个例子中：
'result.fasta' 就是我最后的结果文件
'fasta.fasta' 这里就是我的基因组文件，即源数据fasta文件
上面两个参数如果和脚本不在同一个目录下则可以自己在脚本中添加地址，文件名自行替换

脚本中的'chromosome 11'可以根据自己的需要替换

# -*- coding: utf-8 -*-

import sys
import time
start=time.clock()


outf=open('result.fasta','w')
inf=open('fasta.fasta', 'r')
for line in inf:
    if 'chromosome 11' not in line: #没遇到11号染色体就一直读并输出
        outf.write(line)
    else: 
        break  #读到11号就跳出循环
outf.close()
inf.close()


end=time.clock()
print 'Running time: %s Seconds'%(end-start)  

使用方法

直接把该脚本放到需要处理的fasta文件目录下运行就可以（比如我直接IDLE），也可以手动修改成 一：读取含特定字符的序列并输出 的版本，同样结果也生成在该目录下

这个脚本可以用来处理大文件，因为是读一行写一行。

三：输出前n条序列

背景：朋友给了我这个解法我才知道基因组文件是多少条染色体就有多少条序列，比如牛是30条染色体，基因组fna文件中就是30条序列，从1到30按顺序排列，那我要前十条染色体的序列的话，逐条读取输出，并在第11条处终止就行了，如下

if chr_num<=10: 这条条件语句可以根据需要序列数自行更改

fa_file = open("fasta.fasta","r")
out_file=open("result.fasta",'w')

chr_num=0
for each_line in fa_file:
    each_line = each_line.strip()
    if each_line.startswith(">"):
        chr_num +=1 #利用序列名判断
    if chr_num<=10: 
        out_file.write(each_line+"\n") #输出前10条
    else:
        break

fa_file.close()
out_file.close()

使用方法

同前第二种，不赘述了

可以处理大fasta文件

---

总结

上述这三种都可以根据需要自行调整那几个关键参数，第一种是我借鉴的其他博主的，但发现不适用于大的fasta序列文件，后两种自己写的，实测对将近3G的基因组也就一分钟以内。但第一种可以挑选多个条件（不同基因名）的序列，而后两种的限定条件比较窄，特定字符或者特定数目。
边学边练习（实际上我就刚刚学完了字典……所以这三个脚本基本都是一定看得懂的吧），肯定有很多不对的地方，欢迎交流，请多指教~

参考资料：

[根据ID从FASTA文件中批量提取序列【Python脚本】