Python：从fasta文件读取和打印序列

Question

我试图在输入文件中仅打印序列id，而不是整个描述行，以及fasta文件中它旁边的GC内容，例如：

Seq1   40%
Seq2   37%
Seq3   12%

当我运行这段代码时，什么都没有发生

def main():
    calcGC()

def calcGC():
    fileReader = open("Sequences.fasta",'r')
    for line in fileReader:
        seqID = line.startswith (">")
        seq = line[0:]

    gc_count = float((seq.count("G") + seq.count("C"))) / len(seq)*100
    print(seqID+"   "+ gc_count)

    fileReader.close

main()

Answer 1

这将打印某些内容，而不是您想要的输出。您的代码中有几个错误：

压痕是错误的

SeqID是一个布尔值，它只检查行是否以<字符开头。所以我猜你想测试你是否应该打印这行。无论如何，如果您想要打印的是行号，如果它以<开头，它将是这样的：

def calcGC():
    fileReader = open("Sequences.fasta", 'r')
    for seqID, line in enumerate(fileReader):
        if line.startswith(">"):
            seq = line[0:]

            gc_count = float((seq.count("G") + seq.count("C"))) / len(seq) * 100
            print(seqID, "   ", gc_count)

    fileReader.close()

if __name__ == "__main__":
    calcGC()

Answer 2

首先，这里的人似乎不熟悉.fasta格式(Wikipedia article) header/seq_ID/info行始终以“>；”字符开头，后跟作者放在那里的任何信息，因此检查行是否以“>；”开头

所以如果

if line.startswith(">"):

如果是真的，那么我们在标题行。要获取标题，可以使用

seq_ID = line[1:] # not line[0:], as we don't need the ">"

接下来的所有行都是属于该标题的序列数据，直到下一行以“>；”开头或文件结束

收集完所有序列数据后（在标题后面），然后计算序列数据的GC内容，而不是标题（因为问题和两个答案都是相反的，这在生物学上没有意义）

如果可以，只需使用biopython软件包：

from Bio import SeqIO
from Bio.SeqUtils import GC

input_file = "Sequences.fasta"

for seq_record in SeqIO.parse(input_file, "fasta"):
    print(f'{seq_record.id} {round(GC(seq_record.seq))}%')

（代码只是稍微改编自biopython tutorial）

Answer 3

我首先看到的是范围问题。 在lobal作用域中调用line，但line只能在calcGC函数作用域的循环中访问。在python中，范围由标识定义

现在我在这里看到的第二件事，我不理解的是.startswith()方法的使用。它将返回一个Bollean，而不是seqID…也许在这里添加一个if语句

最后一件事是：您应该使用with语句打开一个文件，它将为您关闭该文件并为您提供一个生成器。 顺便说一下，这里不需要行：seq = line[0:]，您可以在循环中直接使用“seq”，如下所示：

def main():
    calcGC()

def calcGC():
    with open("Sequences.fasta", 'r') as fp :
        # using the enumerate here will give you both the index and the line itself. I assume here that the seqID you wanted to use is the Line index....
        for seqID, seq in enumerate(fp):
            if seq.startswith(">"):
                gc_count = float((seq.count("G") + seq.count("C"))) / len(seq) * 100
                print("{} {}".format(seqID,gc_count))


main()