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

>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]
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

• So, how long is a "very very long time"? Multiple days per file? Minutes? Seconds, but you have a million files? – Graipher Oct 29 '18 at 16:59
• Also, are there any other letters in the genome except "GALMFWKQESPVICYHRNDT"? Because if not, genomesize = len(fasta.seq) would probably work and be a lot faster. – Graipher Oct 29 '18 at 17:02
• Finally, I get 167.77777777777777 with the supplied file (both with your and my code). – Graipher Oct 30 '18 at 7:39
• I had 100,000+ files and last time I let it run it hadn't finished after a week, yes sorry the example input wasn't the whole input, 167.7 is right for that subset! I think len(fasta.seq) would still include stop codons in the final length! (marked by * in the files I'm looking at, apologies again, forgot to include in my input file) – Biomage Oct 30 '18 at 10:20
• What are stop codons (I have bascially no clue of biology)? And I would have assumed that the implementation in the biopython package would be able to deal with them in some way if they have a bioinformatics meaning. Maybe you could add an input where they are included? – Graipher Oct 30 '18 at 10:25

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


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.