| 如何批量更改Fasta文件的序列名？ | 您所在的位置：网站首页 › python批量替换文件名称 › 如何批量更改Fasta文件的序列名？ |
如何批量更改Fasta文件的序列名？
|
在对fasta文件处理的时候，许多小伙伴经常会遇到要修改fasta文件序列名的情况。当要修改的序列只有几条的时候，手动就好了。但当序列成百上千的时候，手动修改就显得很笨了。 笔者最近也遇到了要大量修改序列名的情况，所以写了一个简单的python脚本，在这里和大家共享： 文件准备首先，我们要准备输入文件。 输入文件一：序列文件输入文件一是我们要处理的fasta文件，可以是核苷酸序列，也可以是氨基酸序列，两个都可以，取决于你要处理的序列。 fasta文件的格式很简单，由两部分组成： 一部分是序列名，该部分以>开头，后边跟着该条序列的名字。有时候的时候名字由多个由空格隔开的字符串组成，在分析的时候，我们通常会简化名字，只保留重要的那一个字符串。这是因为很多生物信息软件在处理复杂的序列名时，会出现意想不到的错误，比如基因结构注释软件Braker3。「注意」：序列名只能是一行。 另一部分是序列，可以由多行组成，也可以是一行。 在我们的脚本里，支持对含有空格的序列名处理，但不支持以\t（制表符）分割的，下边举几个例子： >seq1 MAASTMALSSPAFAGKAVNLSPAASEVLGSGRVTMRKTVAKPKGPSGSPWYGSDRVKYLG PFSGESPSYLTGEFPGDYGWDTAGLSADPE >seq2 LHCB1 chr1:123124-124078 MAAATAMHSTTLAGQSLLKPVNELSRKVGSSEARVTMRRTVKSTSDSIWYGADRPKYLGP FSGETPSYLTGEFAGDYGWDTAGLSSDPETFARNRELEVIHARWAMLGALGCLTPELLAK SGVKFGEAVWFKAGAQIFSEGGLDY >seq3 ATGCGCGTGATCGCGCTATCTGCCGGGGTTACCCCGTCCCTTAGGATCGGCTTACCCATGTGCAAGTGCC GTTCACATGGAACCTTTCTCCTCTTCGGCCTTCAAAGTTCTCATTTGAATATTTGCTACTACCACCAAGA TCTGCCGCGACGGCCGCTCTGCCCGGGCTCGCGCCCCGGGTTTTGCAGCGGCCGCCGCGCCCTCCTACTC ATCGGGGCATGGCGCTCGTCCAGATGGCCGGGTGTGGGTCGCGCGCTTTAG >seq4 LHCB2 MAASTMALSSPAFAGKAVNLSPAASEVLGSGRVTMRKTVAKPKGPSGSPWYGSDRVKYLG PFSGESPSYLTGEFPGDYGWDTAGLSADPE猜猜以上的序列哪个没办法使用我们的脚本，没错第四个不行，因为第四个的seq4和LHCB2之间是制表符分割的，这会在后边的处理中引起报错。 除了第四个序列外，均可以使用我们的脚本进行处理，不管是氨基酸序列还是核苷酸序列，只要符合fasta格式，并且序列名中无制表符，都可以使用。 输入文件二：替换集输入文件二是替换集，该文件有两列组成，默认以\t（制表符）分割。替换集的行数取决于我们要更改的基因数，每个要更改的基因一行。「备注」：分隔符可以更改，但更改的时候，需要同步更改代码处理这部分的代码。 下边是一个示例替换集： seq1 gene1 seq2 LHCB1 chr1:123124-124078 gene2 seq3 gene3seq4不符合规范，故没有放进来。 不知道自己的序列名里有没有制表符的可以使用cat -A [filename]查看，只要序列中间的空袭不是^I填充的就没有问题，^I 是一个由 "caret" 符号（^）和大写字母 "I" 组成的字符串，表示制表符。 脚本代码详解以下是我们的脚本（batch_change_FastaSeqName-command_line.py）： #加载包 from sys import argv import sys # 输入文件名 input_file = argv[1] # 更改为您的FA文件名 replace_file = argv[2] # 更改为您的替换文件名 output_file = argv[3] # 更改为新的输出文件名 # 读取fasta文件 with open(input_file, 'r') as f: fasta_lines = f.readlines() # 读取替换文件 replacements = {} with open(replace_file, 'r') as f: for line in f: (old_str, new_str) = line.strip().split('\t')#更改此处，可以更换替换集的分隔符 replacements[old_str] = new_str # 对fasta文件进行替换并保存 with open(output_file, 'w') as f: for line in fasta_lines: # 判断当前行是否以>开头 if line.startswith('>'): # 处理注释行并进行替换 seq_id = line.split('>')[1].strip() # 获取序列ID if seq_id in replacements: # 如果序列ID存在于替换文件中 new_seq_id = replacements[seq_id] line = line.replace(seq_id, new_seq_id) f.write(line) else: # 当seq_id不在replacements中时输出日志并退出程序 print(seq_id + '不在替换集中，程序退出运行，请检查您的替换集') sys.exit() else: # 不以>开头的行不做修改 f.write(line)保存该文件为*.py的脚本就可以了，随后运行代码： python *.py input_file replace_file output_file就可以得到你想要的结果啦~ 简洁代码 from sys import argv import sys input_file = argv[1] replace_file = argv[2] output_file = argv[3] with open(input_file, 'r') as f: fasta_lines = f.readlines() replacements = {} with open(replace_file, 'r') as f: for line in f: (old_str, new_str) = line.strip().split('\t') replacements[old_str] = new_str with open(output_file, 'w') as f: for line in fasta_lines: if line.startswith('>'): seq_id = line.split('>')[1].strip() if seq_id in replacements: new_seq_id = replacements[seq_id] line = line.replace(seq_id, new_seq_id) f.write(line) else: print(seq_id + ' is not in the replacement set, the program exits running, please check your replacement set!') sys.exit() else: f.write(line) 结束语本篇文章就说这么多，如果你有什么问题欢迎在评论区留言或者私信我~ 「相关代码和示例文件我放在github上了，点击阅读原文即可下载！」
|
【本文地址】