FASTA转换脚本

Question

我编写了一个脚本，该脚本以多序列.fasta文件作为输入，并以.tsv格式输出每个氨基酸/序列的数量。

以下是.fasta文件的外观：

>sp|P21515|ACPH_ECOLI Acyl carrier protein phosphodiesterase OS=Escherichia coli (strain K12) OX=83333 GN=acpH PE=1 SV=2 
MNFLAHLHLAHLAESSLSGNLLADFVRGNPEESFPPDVVAGIHMHRRIDVLTDNLPEVREAREWFRSETRRVAPITLDVMWDHFLSRHWSQLSPDFPLQEFVCYAREQVMTILPDSPPRFINLNNYLWSEQWLVRYRDMDFIQNVLNGMASRRPRLDALRDSWYDLDAHYDALETRFWQFYPRMMAQASRKAL

我知道有些库会做同样的事情，但是，我不得不把这个脚本写成家庭作业，而不需要BioPython或其他类似库的帮助。

我的问题是:如何提高代码的效率？

import sys
import os


def readFasta( path ):
    
    seq = { }
    aCount = { }
    alphabet = "ABCDEFGHIJKLMNOPQRSTUVWYZX*-"  

    for fileName in path:
        
        with open( fileName ) as file:
            fastaLine = file.readlines( )
            
            # sequence storage
            
            for line in fastaLine:
                if line[ 0 ] == ">" :
                    headerTitle = line.split( "|" )[ 1 ]
                    seq[ headerTitle ] = ""
                else:
                    seq[ headerTitle ] += line.strip( )

            # aminoacid count

            for id in seq:
                aCount[ id ] = { }
                for sequence in seq[ id ]:
                    if sequence not in alphabet:
                        print( "Check your sequence! Character: " + str(sequence) + " from sequence: " + id )
                    for letter in alphabet:
                        aCount[ id ][ letter ] = 0
            
            # count aminoacid occurence/sequence
            
            for name in seq:
                for char in seq[ name ]:
                    if char in aCount[ name ]:
                        aCount[ name ][ char ] += 1
    
    # write to .csv file
    
    outFile = open( "output.csv" , "w" )
    
    
    outFile.write( "Name" + '\t' )
    
    header = ""
    for letter in alphabet:
        header += letter + '\t'

    outFile.write( header + "Total" )
    outFile.write( '\n' )          
    
    for name in aCount:
        entry = ""
        counter = 0
        for letter in aCount[ name ]:
            entry += str(aCount[ name ][ letter ]).strip( ) + '\t'
            counter += aCount[ name ][ letter ]
        outFile.write( name + '\t' + entry + '\t' + str( counter ) )
        outFile.write( '\n' )

# files to list

entries = os.listdir( path )

# filter for .fasta/.fa files

new_list = [ ]
for entry in entries:
    if ".fasta" in entry:
        new_list.append( entry )
    elif ".fa" in entry:
        new_list.append( entry )
        
# execute func

readFasta( new_list )

Answer 1

除了在其他答复中提出的问题外：

外来进口

第一行是import sys，但我没有看到任何地方都使用sys；这是可以实现的。

显式命名

通常使用显式名称更好。因此，不只是seq (例如，sequence或sequences的缩写)，而是把它拼出来。在这种情况下，您应该使用sequences。类似地，aCounts没有告诉我你在计算什么，除了涉及到一个a之外，只需将其称为amino_acid_counts。

一致命名

当您在代码中的两个不同位置谈论同一个对象时，您应该更喜欢使用一致的名称。例如，在读取文件时使用headerTitle作为seq键，但在下一个块中使用id表示相同的键，在该块中使用name。同样，我将使用一个更显式的名称，比如sequence_id。

阴影内置对象

说到id，那实际上是一个内建函数。重新命名语言中内置的东西通常是个坏主意。这被称为“阴影”，通过这样做就有可能破坏代码--而且通常认为是不好的形式。

使文件打开的时间比它需要的时间长，

一旦您运行了fastaLine = file.readlines( )，您就完成了该文件，并且您可以关闭它--这意味着您已经完成了with open块，并且您可以在它之后取消缩进块。当然，正如Reinderien所指出的，在迭代每一行时最好只处理它们。(这对于大型文件尤为重要，您可能希望一次只将一行加载到内存中，而不是将整个文件加载到内存中。)因此，您希望with open块看起来更像这样：

        # sequence storage
        with open(path, "r") as file:
            for line in file:
                if line[0] == ">" :
                    sequence_id = line.split("|")[1]
                    sequences[sequence_id] = ""
                else:
                    sequences[sequence_id] += line.strip()

        # amino acid count

请注意，我没有缩进氨基酸计数注释(以及后面的所有内容)。

可能过度循环

我猜您只是想警告用户在给定行中的所有未识别字符，而不是每一次这样的字符的出现。在这种情况下，最好对序列中的所有字符创建一个set，并检查是否有不在alphabet中的字符：

sequence_characters = set(seq[id])
if not sequence_characters.issubset(alphabet):
    print(
        "Check your sequence! Extra characters: '" + str(sequence_characters - alphabet)
        + "' from sequence: " + str(sequence_id)
    )

绝对过度循环

在执行for sequence in seq[ id ]时，您将遍历序列中的每个字符。虽然我怀疑这是您在前面的行中想要的，但这肯定不是您在执行for letter in alphabet时想要的；这个循环应该是不缩进的，因此它只对每个id执行一次。

使用内置的str.count方法来计算

。

Python已经知道如何计算一个字符在字符串中出现的次数:这是这个str.count方法。因此，与其循环遍历字符串中的每个字符并将1添加到适当的字符中，不如在字母表上循环。您可以删除将计数初始化为0的循环，然后计算出计数。将这些内容与前面的项目(包括重命名一些变量)放在一起，您将得到如下内容：

        # amino acid count
        for sequence_id, sequence in sequences.items():
            amino_acid_counts[sequence_id] = { }
            sequence_characters = set(sequence)
            if not sequence_characters.issubset(alphabet):
                print(
                    "Check your sequence! Extra characters: '" + str(sequence_characters - alphabet)
                    + "' from sequence: " + str(sequence_id)
                )
            for letter in alphabet:
                amino_acid_counts[sequence_id][letter] = sequence.count(letter)

而# count aminoacid occurence/sequence下的整个循环可以被删除。这样会快得多。

使用with open块编写文件，也使用

正确地使用with open块读取文件；在编写CSV文件时也要使用它：

    with open("output.csv", "w") as outFile:

(以及缩进任何跟随和使用outFile的内容)。

不需要添加显式字符串

当您编写"Name" + '\t'时，您还可以编写"Name\t"。同样的，当你写

    outFile.write( header + "Total" )
    outFile.write( '\n' )

你也可以写

    outFile.write(header + "Total\n")

Answer 2

欢迎来到代码评审！

我再加上莱因德林的另一个答案。

PEP-8

在python中，通常(并建议)遵循PEP-8风格的指南来编写干净、可维护和一致的代码。

函数和变量应该在lower_snake_case中命名，类命名为UpperCamelCase，常量命名为UPPER_SNAKE_CASE。

在定义或调用函数时，参数周围不需要不必要的空格。这也适用于引用数组/列表元素。

函数

将您的代码拆分成单独的较小的函数，执行奇异的任务。几个例子是:读取/解析fasta文件，验证基因序列，计数出现的次数(尝试collections.Counter包来实现这一点)。

CSV/TSV

Python附带内置包csv，它也可以用来编写tsv文件。您不需要自己修改/写每一行。

Glob搜索文件

Python 3+还有另一个内置的包pathlib，它支持使用glob模式获取所有文件。同样，您不需要手动(比方说)聚合所有具有.fa或.fasta扩展名的文件。

if __name__块

对于脚本，将可执行特性放入if __name__ == "__main__"子句是一种很好的做法。

Answer 3

路径？

当然，path并不是一条单一的路径，因为您可以遍历它。因此，至少，这是一个很糟糕的名字，应该是paths。迭代它的变量fileName应该是file_name；与其他函数和变量名称类似。

线迭代

而不是调用readlines()，只需迭代file本身，这将产生同样的效果。

字母表

为什么它坏了- WYZX？

冗余if

if ".fasta" in entry:
    new_list.append( entry )
elif ".fa" in entry:
    new_list.append( entry )

可以将if折叠为

if '.fa' in entry:
    new_list.append(entry)

如果您查看子字符串模式，则第二个包含第一个。另外，不要在父母的内部添加空格。