Counting average protein length from a proteome fasta file this code is counting the average protein length (number of amino acids) across the proteome of some microbes.

As it stands the code takes a very very long time to run, I think I'm being very inefficient somewhere but not 100% sure where.

The way I run the code is using a bash-loop (using find for each example of particular file type, like *.faa) like this:

for FAA in $(find . -name "*.faa")
do
    python proteinlengthgen.py $FAA
done


My code:

from Bio import SeqIO
import sys
from collections import Counter

def species_name_function(infile):
    for record in SeqIO.parse(infile, "fasta"):
        speciesname = record.description.split('[', 1)[1].split(']', 1)[0]
    return speciesname

if __name__ == '__main__':
    frequency = []
    infile = sys.argv[1]
    avgproteinlength = 0
    totalproteinlength = 0
    fastacounter = 0

    for fasta in SeqIO.parse(open(infile), "fasta"):
        dna = str(fasta.seq)
        freq = Counter(dna)
        fastacounter = fastacounter + 1
        frequency.append(freq)
        species_name = species_name_function(infile)

        genomesize = freq['G'] + freq['A'] + freq['L'] + freq['M'] + freq['F'] + freq['W'] + freq['K'] + freq['Q'] + freq['E'] + freq['S'] + freq['P'] + freq['V'] + freq['I'] + freq['C'] + freq['Y'] + freq['H'] + freq['R'] + freq['N'] + freq['D'] + freq['T']

        listofbases = ['G', 'A', 'L', 'M', 'F', 'W', 'K', 'Q', 'E', 'S', 'P', 'V', 'I', 'C', 'Y', 'H', 'R', 'N', 'D', 'T']

        totalproteinlength = genomesize + totalproteinlength

    avgproteinlength = totalproteinlength / fastacounter

    towrite = str(avgproteinlength) + '\t' + species_name + '\t' + '\n'

    with open("proteinlength.csv", "a") as myfile:
        myfile.write(towrite)


Example of a subset of input file (I cannot use a whole one as they are very big):

&gt;WP_013179448.1 DNA-directed RNA polymerase [Methanococcus voltae]
MYKILTIEDTIRIPPKMFGNPLKDNVQKVLMEKYEGILDKDLGFILAIEDIDQISEGDIIYGDGAAYHDTTFNILTYEPE
VHEMIEGEIVDIVEFGAFIRLGPLDGLIHISQVMDDYVAFDPQREAIIGKETGKVLEKGDKVRARIVAVSLKEDRKRGSK
IALTMRQPALGKLEWLEDEKLETMENAEF
&gt;WP_013179449.1 DNA-directed RNA polymerase subunit E'' [Methanococcus voltae]
MARKGLKACTKCNYITHDDFCPICQHETSENIRGLLIILDPVNSEVAKIAQKDIKGKYALSVK
&gt;WP_013179451.1 30S ribosomal protein S24e [Methanococcus voltae]
MDIKVVSEKNNPLLGRKEVKFALKYEGATPAVKDVKMKLVAILNANKELLVIDELAQEFGKMEANGYAKIYESEEAMNSI
EKKSIIEKNKIVEEAEEAQE
&gt;WP_013179452.1 30S ribosomal protein S27ae [Methanococcus voltae]
MAQKTKKSDYYKIDGDKVERLKKSCPKCGEGVFMAEHLNRFACGKCGYMEYKKNEKAEKEE
&gt;WP_013179453.1 hypothetical protein [Methanococcus voltae]
MNELNELKNPDKIDGKNNNTKNNNNNNNKDSNTENSITEIIKADNETQDNLSDLCVLEDIKTLKSKYKVYKTSKYLTKKD
INNIIEKDYDEIIMPQSIYKLLNEKNKSSMEKLRLCGIIVKTTDNVGRPKKITKYDKDKIKELLVDGKSVRKTAEIMDMK
KTTVWENIKDCMNEIKIEKFRKMIYEYKELLIMQERYGSYVESLFLELDIYINNEDMENALEILNKIIIYVKSEDKKD
&gt;WP_013179454.1 integrase [Methanococcus voltae]
MKNKRINNNQKSKWETMRTDVINTQRNQNINSKNKQYRVKKHYCKEWLTKEELKVFIETIEYSEHKLFFKMLYGMALRVS
ELLKIKVQDLQLKEGVCKLWDTKTDYFQVCLIPDWLINDIKEYIALKSLDSSQELFKFNNRKYVWELAKKYSKMAELDKD
ISTHTFRRSRALHLLNDGVPLEKVSKYLRHKSIGTTMSYIRITVVDLKQELDKIDDWYEL
&gt;WP_013179455.1 hypothetical protein [Methanococcus voltae]
MNTQNAIKKTLKTSKVNKNISNVIIGYSAILLDTYSNNKNLLLVKYDKLFKGFLNSSSITEKQYNKLYDTVLNSLF
&gt;WP_013179456.1 hypothetical protein [Methanococcus voltae]
MVVKLVKISNGGYVSSLELKRINDIILSQLTNEFTIKDIVNMYSNKYDDCNNNAIAQKTRRLLNNHIESGVFTVRNALKN
KKIYKFKDVFVPASAGDTNTSLLFYSTSMKNSNHIEKQKKNNNKYNTNVNKPTITPDQIRVMAGIVNNPQIKSLKKERFK
SILHLNCKHMLNEEDRTELLENFKEYIIKASSQNLVLERTRYHKNKPKYITFPYLTRFTNSKQLKRQLAQYNCIFEQKAI
KYNRGVHLTLTTDPKLFRNIYEANKHFSKAFNRFMSFLSKRNKDVNGKSRRPVYLAAYEYTKSGLCHAHILIFGKKYLLD
QRVITQEWSRCGQGQIAHIYSIRRDGINWIYGRARPEEIQTGEYKNANEYLKKYLVKSLYSELDGSMYWAMNKRFFTFSK
SLKPAMMKRPKPEKYYKFVGIWTDEEFTDEINQAFKTGIPYDEIKKIHWKNLNNGLSCG
&gt;WP_013179457.1 hypothetical protein [Methanococcus voltae]
MVRGRYPVFSGFKKFNKINLGKEKRNEGVYKYYNQDKTLLYVGVSNEVKLRLLSAYYGRSDYAVLENKKKLRQNIAYYKV
KYCGLDQARKIEHRIKKQCKYNLN


Output of code:

302.7408088235294   Methanococcus voltae
The first thing to do to accelerate processing is to read the file only once. species_name_function reads the entire file, and is called multiple times. Extract the line which parses it and inline into the main loop:

    for fasta in SeqIO.parse(open(infile), "fasta"):
        species_name = fasta.description.split('[', 1)[1].split(']', 1)[0]




Now,


        dna = str(fasta.seq)
        freq = Counter(dna)



does unnecessary work. There's no need to convert to a string before counting.



listofbases is unused, and genomesize is calculated by a large explicit sum. My best guess is that you were experimenting, trying to come up with

        genomesize = sum(freq[base] for base in "GALMFWKQESPVICYHRNDT")




Using those changes and applying some further simplifications which I don't think need individual explanation, I get:

from Bio import SeqIO
import sys
from collections import Counter

if __name__ == '__main__':
    infile = sys.argv[1]
    genomesizes = []

    for fasta in SeqIO.parse(open(infile), "fasta"):
        species_name = fasta.description.split('[', 1)[1].split(']', 1)[0]
        freq = Counter(fasta.seq)
        genomesizes.append(sum(freq[base] for base in "GALMFWKQESPVICYHRNDT"))

    avgproteinlength = sum(genomesizes) / len(genomesizes)
    towrite = str(avgproteinlength) + '\t' + species_name + '\t' + '\n'

    with open("proteinlength.csv", "a") as myfile:
        myfile.write(towrite)




However, I'm not very convinced by the explicit output filename. Also, I would be inclined to simplify the usage by supporting multiple files:

from Bio import SeqIO
import sys
from collections import Counter

if __name__ == '__main__':
    for infile in sys.argv[1:]:
        genomesizes = []

        for fasta in SeqIO.parse(open(infile), "fasta"):
            species_name = fasta.description.split('[', 1)[1].split(']', 1)[0]
            freq = Counter(fasta.seq)
            genomesizes.append(sum(freq[base] for base in "GALMFWKQESPVICYHRNDT"))

        avgproteinlength = sum(genomesizes) / len(genomesizes)
        print(str(avgproteinlength) + '\t' + species_name)


called as

find . -name "*.faa" | xargs python proteinlengthgen.py &gt;&gt;proteinlength.csv