FASTA序列字典

Question

我正试图从FASTA序列中创建一本字典。但是，我一直收到同样的错误，即变量是没有定义的，即使它是定义的。有人能帮我理解一下这段代码是否正确吗？

#!/usr/bin/pyhton
# Build a dictionary containing all sequences from a FASTA file from typing import List



try:
    f= open("FASTA.rtf",'r')
except IOError:
    print('The file "FASTA" does not exist')

seqs = {}
for line in f:
    # Let's discard the newline at the end (if any)
    line = line.rstrip()
    if line[0]==">": #or line.startswith('>')
        words=line.split()
        name= words[0][1:]
        seqs[name]=""
    else: #sequence, not header
        seqs[name]=seqs[name]+line

for name,seq in seqs.items():
    print(name,seq)

Answer 1

首先，为什么在.rtf (富文本格式)中使用FASTA，而不是使用.fa/.fasta或.fq/.fastq？这可能导致了你的一些问题。下面是我用来编辑FASTA文件基因标签的代码片段。

def remove_trailing_inplace(filename):

    with fileinput.input(files=str(filename), inplace=True, backup=".bak") as f
        for line in f:
            if line[0] == ">": # only engages with lines that have genes
                a = line.split(" ") # series of splits to properly modify each
                idx = a[0].split(".")[0]
                gene = a[3].split(".")[0]
                rest = " ".join(a[4:])
                print(idx + " " + a[1] + " " + a[2] + " " + gene + " " + rest.str)
            else:
                print(line, end="")

这个代码片段工作得很好，捕捉到了我需要的所有基因。最大的区别是它直接获取FASTA文件，而不是fasta格式的rtf。此外，我不使用任何形式的剥离: FASTA格式不需要它，据我所知。

此外，正如@Frank所评论的，您的if/else语句在逻辑上是不正确的，因为在else部分中，您试图使用仅在if部分中定义的name。if和else语句是相互排斥的，因此不能依赖在另一个变量中声明的变量。

幸运的是，修复并不太困难。我建议的是在循环开始时(在name之外)声明if/else。

seqs = {}
names = []
for line in f:
    
    if line[0]==">":
        words=line.split()
        name= words[0][1:]
        names.append(name)
        seqs[name]=""
    else: #sequence, not header
        seqs[names[-1]]=seqs[name[-1]]+line

在编写劣质代码的同时，应该避免名称错误。此外，由于FASTA是这样的，我们可以假设我们总是有一个'>‘的开始.

编辑:对于那些不熟悉FASTA格式的人来说，下面是维基百科的一个例子：

>SEQUENCE_1
SATVSEINSETDFVAKNDQFIALTKDTTAHIQSNSLQSVEELHSSTINGVKFEEYLKSQI
ATIGENLVVRRFATLKAGANGVVNGYIHTNGRVGVVIAAACDSAEVASKSRDLLRQICMH

其中带有>的行包含序列信息，而后面的所有行(直到下一个lines‘> '>')都表示序列本身。因此，在完成第一个序列之前，我们总是可以假设我们永远不会碰到另一个序列。这就是为什么上面的代码可以工作，以及为什么使用names[-1]是可以的。