Python：从 fasta 文件到字典丢失核苷酸

Question

我正在尝试编写代码来提取 fasta 文件中最长的 ORF。它来自 Coursera Genomics 数据科学课程。

该文件为练习文件："dna.example.fasta"

数据在这里：https://d396qusza40orc.cloudfront.net/genpython/data_sets/dna.example.fasta

我的部分代码在下面提取阅读框2（从序列的第二个位置开始。例如：seq：ATTGGG，以获得阅读框2：TTGGG） :

#!/usr/bin/python


import sys
import getopt

o, a = getopt.getopt(sys.argv[1:], 'h')

opts = dict()
for k,v in o:
        opts[k] = v
        if '-h' in k:
                print "--help\n"

if len(a) < 0:
        print "missing fasta file\n"



f = open(a[0], "r")


seq = dict()
for line in f:

        line = line.strip()

        if line.startswith(">"):
                name = line.split()[0]

                seq[name] = ''

        else:
                seq[name] = seq[name] + line[1:]


k = seq[">gi|142022655|gb|EQ086233.1|323"]
print len(k)

这个特定序列的长度应该是 4804 bp。因此，通过单独使用这个序列，我可以获得正确的答案。然而，使用代码，在字典中，这个特定序列变成只有 4736 bp。

我是 python 的新手，所以我不知道那 100 个基点去了哪里？

谢谢，

晓

Answer 1

再看看你的数据文件

一些行的例子：

>gi|142022655|gb|EQ086233.1|43 marine metagenome JCVI_SCAF_1096627390048 genomic scaffold, whole genome shotgun sequence
TCGGGCGAAGGCGGCAGCAAGTCGTCCACGCGCAGCGCGGCACCGCGGGCCTCTGCCGTGCGCTGCTTGG
CCATGGCCTCCAGCGCACCGATCGGATCAAAGCCGCTGAAGCCTTCGCGCATCAGGCGGCCATAGTTGGC

注意序列是如何从每行的第一个值开始的。您的添加行 seq[name] = seq[name] + line[1:] 是在第一个字符之后添加该行的所有内容，不包括第一个字符（Python 2 索引从零开始）。事实证明，您丢失的核苷酸数量是构建该基因组所需的行数，因为您每次都丢失了第一个字符。

修改后的方式是seq[name] = seq[name] + line，它只是添加了一行而不丢失第一个字符。

找到这类调试错误的最快方法是使用正式的调试器，或者在您的代码中添加一堆打印语句并使用文件的一小部分进行测试——您可以看到输出并检查自己是否正确。一个可能包含 50 个核苷酸而不是 5000 个核苷酸的短文件更容易手动评估并确保代码正在执行您想要的操作。这就是我在大约 5 分钟内得出问题答案的方法。

另外，为了将来参考，请事先说明您使用的 python 版本。 python 2（您正在使用的那个）和 python 3.

之间有很多区别

我对你的代码做了一些额外的测试，如果你在最后得到任何额外的字符，它们可能是空格。确保在将每一行添加到字符串之前使用 .strip() 方法，这会清除空格。

处理您的评论，

要仅从序列第一行的第二个位置开始，然后使用整行直到下一个核苷酸，您可以利用文件的线性格式，只需在 if 语句中再添加一个子句，一个小精灵。这将测试我们是否在序列的第一行，如果是，则使用从第二行开始的字符，如果我们在任何其他行，则使用整行。

        if line.startswith(">"):
                name = line.split()[0]

                seq[name] = ''

             #If it's the first line in the series, then the dict's value
             # will be an empty string, so this elif means "If we're at the
             # start of the series..."
        elif seq[name] == '':
                seq[name] = seq[name] + line[1:]
        else:
                seq[name] = seq[name]

这种适应将从基因组中的第二个核苷酸开始，而不会丢失其余核苷酸中每一行的第一个。

Python：从 fasta 文件到字典丢失核苷酸

Python: losing nucleotides from fasta file to dictionary

python

fasta