-9

我需要一个 python 脚本来将 fasta 格式解析为基于行的序列。

我需要改变这个：

>GeneID12345
ATTACATATACCATACC
CCATATTAATCCGAGGG
TTACCTATAGGTATACC
>GeneID12346
TTGATACCATATATCCC
ATATGCCCTATATTCCT
TTACCTATC

对此：

GeneID12345 ATTACATATACCATACCCCATATTAATCCGAGGGTTACCTATAGGTATACC
GeneID12346 TTGATACCATATATCCCATATGCCCTATATTCCTTTACCTATC
4

4 回答 4

3

这将完成这项工作：

import sys

sep=''
with open(sys.argv[1]) as f:
    for line in f:
        if line.startswith(">GeneID"):
            sys.stdout.write(sep+line.strip()+' ')
        else:
            sys.stdout.write(line.strip())
            sep='\n'

像运行python parse_fasta.py input.fasta。

于 2013-04-20T23:05:21.803 回答
3

如果您不限于python以下简洁的解决方案awk：

$ awk 'NR>1&&/^>/{printf "\n"}{printf "%s",$0}/^>/{printf " "}' file
>GeneID12345 ATTACATATACCATACCCCATATTAATCCGAGGGTTACCTATAGGTATACC
>GeneID12346 TTGATACCATATATCCCATATGCCCTATATTCCTTTACCTATC
于 2013-04-20T23:18:01.493 回答
1 
data = []
last = None
for line in f:
    if line.strip('ATCG') != '': # contains non-nucleobases, i.e. is an id
        if last: # save previous data
            data.append(last)
        last = line + ' '
    else:
        last += line
if last:
    data.append(last)

# now pretty-print
for gene in data:
    print(gene)

结果是：

GeneID12345 ATTACATATACCATACCCCATATTAATCCGAGGGTTACCTATAGGTATACC
GeneID12346 TTGATACCATATATCCCATATGCCCTATATTCCTTTACCTATC

这假设f是包含原始数据的文件；即循环遍历文件的行。

于 2013-04-20T23:13:20.407 回答
1

步骤#1：安装BioPython。不要浪费时间重新发明轮子。如果您正在使用 FASTA 文件，您可以编写自己的半功能解析器，或者您可以使用已经存在的解析器，以及教程和所有内容。

步骤#2：嗯，实际上，步骤#1 主要是它。代码：

from Bio import SeqIO

with open("example.fasta") as fp_in, open("newformat.txt", "w") as fp_out:
    for record in SeqIO.parse(fp_in, "fasta"):
        fp_out.write("{} {}\n".format(record.id, record.seq))

生产

~/coding$ cat newformat.txt 
GeneID12345 ATTACATATACCATACCCCATATTAATCCGAGGGTTACCTATAGGTATACC
GeneID12346 TTGATACCATATATCCCATATGCCCTATATTCCTTTACCTATC

更重要的是，您会获得SeqRecord易于操作的实例，并且可以正确获取可能存储在更复杂示例中的其他信息：

>>> record
SeqRecord(seq=Seq('TTGATACCATATATCCCATATGCCCTATATTCCTTTACCTATC',
SingleLetterAlphabet()), id='GeneID12346', name='GeneID12346', 
description='GeneID12346', dbxrefs=[])
于 2013-04-20T23:25:28.783 回答