Counting average protein length from a proteome fasta file

Question

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):

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

Output of code:

302.7408088235294   Methanococcus voltae

Answer 1

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 >>proteinlength.csv

Answer 2

In Python there is no need to initialize variables before using them. So e.g. avgproteinlength = 0 is not needed.

Python also has an official style-guide, PEP8, which recommend using lower_case_with_underscores for variable and function names, so that should probably be named average_protein_length.

While we are renaming things, species_name_function is probably better named get_species_name.

Now, Bio.Seq objects actually support getting their length, so your main block can be greatly reduced in complexity (which should also cut down on the runtime). You can also pull out getting the species name from the loop (it should stay the same, and if not you are not taking that into account at the moment either). You are also not using the frequency variable, nor the listofbases.

from statistics import mean

if __name__ == '__main__':
    infile = sys.argv[1]
    species_name = get_species_name(infile)
    genomes = (fasta.seq for fasta in SeqIO.parse(open(infile), "fasta")) # a generator
    average_protein_length = mean(len(seq) for seq in genomes)

    with open("proteinlength.csv", "a") as myfile:
        myfile.write(f"{average_protein_length}\t{species_name}\t\n")

I used the Python 3.4+ statistics.mean function and the Python 3.6+ f-strings.

This of course assumes that the length of the genome which you calculate by hand by summing all frequencies, is the same as the length of the genome sequence.