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

生信问题记录

我的需求

input：

FASTA文件，含六千余个蛋白序列。命名为FA.fasta
txt文件，经过interpro注释后，筛选出五千余个蛋白，将五千余个蛋白ID导出到txt文件中，每行一个。命名为ID.txt

output：

FASTA文件，根据ID.txt里的ID从FA.fasta里提取出对应的五千余个序列。提取结果也以文件的格式保存为out_file.fasta

环境

Windows 10 64x
命令行运行py

解决流程

如果编程娴熟，可以直接用perl或Python写脚本；碍于个人技艺生疏，直接面向百度编程。搜索发现，解决方案非常多，最终某位网友的Python脚本成功运行

提示如果你是用ID行，形如 >RLS74643.1 dihydrofolate reductase [Planctomycetes bacterium] 来提取序列，那么更为简单，可以使用联川生物的免费在线小工具

参考

1c43f522e1c3的帖子

https://www.jianshu.com/p/b7032dfae44e

步骤

直接复制网友的代码会出很多错，要先安装click模块和biopython的包

安装包

pip install click
pip install biopython

复制粘贴，另存为get_seqs_by_id.py

#coding:utf-8import clickfrom Bio import SeqIO@click.command()@click.option('-f', '--fastafile', help='Input a fasta file', required=True)@click.option('-i', '--idfile', help='Input an idlist', required=True)@click.option('-o', '--outfile', help='Input the name of result file', default='result_out.fa')def main(fastafile="FA.fa",idfile="ID.txt",outfile= "result_out.fa"):with open(idfile) as id_handle:wanted = set(line.rstrip("\n").split(None,1)[0] for line in id_handle)print("Found %i unique identifiers in %s" % (len(wanted), idfile))records = (r for r in SeqIO.parse(fastafile, "fasta") if r.id in wanted)count = SeqIO.write(records, outfile, "fasta")print("Saved %i records from %s to %s" % (count, fastafile, outfile))if count < len(wanted):print("Warning %i IDs not found in %s" % (len(wanted) - count, fastafile))if __name__ == '__main__':main()

将输入文件和py文件放在同一目录下，命令行输入

python get_seqs_by_id.py -f **.fasta -i ID.txt -o out_file.fasta

如果顺利出结果再好不过了，然而我报了错

报错不要慌，只有最后一行是关键信息

debug

错误的意思是：Unicode的解码（Decode）出现错误（Error）了，以gbk编码的方式去解码（该字符串变成Unicode），但是此处通过gbk的方式，却无法解码（can’t decode ）。“illegal multibyte sequence”意思是非法的多字节序列，即没法（解码）了。

于是我在open函数这一行加入encoding=‘UTF-8’

with open(idfile,encoding='UTF-8') as id_handle:

结果又双叒报错了！

耐心搜索，发现可能是用Python读文件（txt或者csv），出现编码错误
打开我的ID.txt一看，从Excel复制的信息居然是UTF-16！怪不得一直报错。把txt文件用UTF-8另存就好了……

运行成功

运行成功是这样的，输出一个fasta文件，正是我想要的。

根据ID从FASTA文件中批量提取序列【Python】相关推荐

从FASTA文件中批量提取指定序列【Python脚本】
文章目录前言一:读取含特定字符的序列并输出演示二:读到某一个字符之前的全部输出使用方法三:输出前n条序列使用方法总结前言背景:学测序流程的时候,做到mapping的时牛的基因组有两 ...
（连载1.1）从招股说明书pdf文件中批量提取财务报表数据
目录标题背景说明阅读代码时注意点源代码背景说明本文选取的是当前日期上交所科创板所有的上市公司样本. 用八爪鱼从上交所公告页面爬取公告下载链接,使用迅雷批量下载. 阅读代码时注意点流程思路: ...
利用NCO或者CDO从nc文件中批量提取数据
nco提取数据方法: ncrcat -v bsf MMEAN*.nc -o bsf.nc cdo提取数据方法: cdo select,name=PRECT NBF1850_f19_tn11*.nc P ...
linux提取fasta文件的id,从大的fasta文件中提取特定的fasta序列
我想使用以下脚本从大的fasta文件中提取特定的fasta序列,但输出为空.从大的fasta文件中提取特定的fasta序列 transcripts.txt文件包含我想从assembly.fasta到s ...
实战：从Mysql数据库frm文件中，提取表结构创建SQL语句
需求在某些特殊的场景下,例如你的mysql数据库无法启动,需要你将表的ibd文件拷贝到另一个数据库中,恢复业务数据库,恢复业务数据的前提,是你需要在另一个数据库中,创建好一模一样的表结构.这时你就需 ...
fasta文件中序列的排序
同样的名为read_1.fa 的fasta文件,里面有若干序列,如: >@r1 TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAA ...
怎么将文件名称批量提取？
如何将文件名称批量提取?有些小伙伴是从事文件管理工作的,有时候为了便于对电脑中的文件进行管理,我们就需要将某些重要文件的名称提取出来,然后整理归档.如果我们将所有文件的名称一个一个的复制,然后再粘贴到 ...
Excel从详细地址中批量提取分列出省、市、县区
Excel快速从详细地址中批量提取分列出省.市.县区! 这个操作必修给跪啊,从终端收集来一堆手机网点的地址信息,现在领导要求汇总出某个省共计有多少网点要开展活动,可是员工小明一打开表懵圈了, 怎么回事 ...
fasta文件中DNA to RNA
同样的名为read_1.fa 的fasta文件,里面有若干序列,如: >@r1 TGAATGCGAACTCCGGGACGCTCAGTAATGTGACGATAGCTGAAAACTGTACGATAA ...