6
\$\begingroup\$

I wrote this code to count the number of amino acids that appear in a proteome. However, I was recently told by someone I could improve it quite a lot to make it more readable and efficient by using dictionaries.

This is the code:

import Bio
from Bio import SeqIO
from Bio import AlignIO
import itertools
import numpy as np
import glob
import pandas as pd
import csv
import re
from collections import defaultdict
import time
#Packages being imported

GAP = 0
SEQCOUNTERPERALIGNMENT = 0
AMINOACIDSPERSTRAIN = 0
NONE = 0

G = 0
A = 0
L = 0
M = 0
F = 0
W = 0
K = 0
Q = 0
E = 0
S = 0
P = 0
V = 0
I = 0
C = 0
Y = 0
H = 0
R = 0
N = 0
D = 0
T = 0

#These variables are all used for counting amino acid proportions later on, here I am establishing them as zero,


#This is a list of all amino acids, these will evenetually become percentages, but are being created here for use later on
current_time = time.strftime("%d.%m.%y %H:%M", time.localtime())
output_name = 'test#%s.txt' % current_time
file = open(output_name, "w+")
#this opens a results file for writing, with a time stamp

files = glob.glob('*.faa')
#this collects all the files ending in .faa, within the current directory as a list
for x in files:
    #for each file
    SEQCOUNTERPERALIGNMENT = 0
    #the sequence counter is reset
    FILE = x
    #we are establishing the file variable within the loop from x
    for record in SeqIO.parse (FILE, "fasta"):
        #using the parse function of biopyton, we read the fasta file
        sequence=record.seq
        #the sequence is established from each record of the fasta file
        for character in sequence:
            #for each variable within the record
            if character in ['G', 'A', 'L', 'M', 'F', 'W', 'K', 'Q', 'E', 'S', 'P', 'V', 'I', 'C', 'Y', 'H', 'R', 'N', 'D', 'T']:
                #if the variable is one of the 20 amino acids, we increase the count by 1, if not we either describe it as a gap or increase none by 1, none can be used to check if there is any issue with the code
                SEQCOUNTERPERALIGNMENT = SEQCOUNTERPERALIGNMENT + 1
                AMINOACIDSPERSTRAIN = AMINOACIDSPERSTRAIN + 1
                if character in 'G':
                    G = G + 1
                elif character in 'A':
                    A = A + 1
                elif character in 'L':
                    L = L + 1
                elif character in 'M':
                    M = M + 1
                elif character in 'F':
                    F = F + 1
                elif character in 'W':
                    W = W + 1
                elif character in 'K':
                    K = K + 1
                elif character in 'Q':
                    Q = Q + 1
                elif character in 'E':
                    E = E + 1
                elif character in 'S':
                    S = S + 1
                elif character in 'P':
                    P = P + 1
                elif character in 'V':
                    V = V + 1
                elif character in 'I':
                    I = I + 1
                elif character in 'C':
                    C = C + 1
                elif character in 'H':
                    H = H + 1
                elif character in 'R':
                    R = R + 1
                elif character in 'N':
                    N = N + 1
                elif character in 'D':
                    D = D + 1
                elif character in 'T':
                    T = T + 1
                elif character in 'Y':
                    Y = Y + 1
            elif character in '-':
                GAP = GAP + 1
            else:
                NONE = NONE + 1

        pG = 100.* G / AMINOACIDSPERSTRAIN
        pA = 100.* A / AMINOACIDSPERSTRAIN
        pL = 100.* L / AMINOACIDSPERSTRAIN
        pM = 100.* M / AMINOACIDSPERSTRAIN
        pF = 100.* F / AMINOACIDSPERSTRAIN
        pW = 100.* W / AMINOACIDSPERSTRAIN
        pK = 100.* K / AMINOACIDSPERSTRAIN
        pQ = 100.* Q / AMINOACIDSPERSTRAIN
        pE = 100.* E / AMINOACIDSPERSTRAIN
        pS = 100.* S / AMINOACIDSPERSTRAIN
        pP = 100.* P / AMINOACIDSPERSTRAIN
        pV = 100.* V / AMINOACIDSPERSTRAIN
        pI = 100.* I / AMINOACIDSPERSTRAIN
        pC = 100.* C / AMINOACIDSPERSTRAIN
        pH = 100.* H / AMINOACIDSPERSTRAIN
        pR = 100.* R / AMINOACIDSPERSTRAIN
        pN = 100.* N / AMINOACIDSPERSTRAIN
        pD = 100.* D / AMINOACIDSPERSTRAIN
        pT = 100.* T / AMINOACIDSPERSTRAIN
        pY = 100.* Y / AMINOACIDSPERSTRAIN


    Speciesname = record.description.split('[', 1)[1].split(']', 1)[0]



    file.write ('\nG,' + str(pG) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nA,' + str(pA) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nL,' + str(pL) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nM,' + str(pM) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nF,' + str(pF) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nW,' + str(pW) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nK,' + str(pK) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nQ,' + str(pQ) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nE,' + str(pE) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nS,' + str(pS) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nP,' + str(pP) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nV,' + str(pV) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nI,' + str(pI) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nC,' + str(pC) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nH,' + str(pH) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nN,' + str(pN) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nR,' + str(pR) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nD,' + str(pD) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nT,' + str(pT) +  ',' + str(FILE) + ',' + str(Speciesname))
    file.write ('\nY,' + str(pY) +  ',' + str(FILE) + ',' + str(Speciesname))

    #writes the seperate amino acid proportions to the file
    AMINOACIDSPERSTRAIN = 0
    GAP = 0
    NONE = 0

    G = 0
    A = 0
    L = 0
    M = 0
    F = 0
    W = 0
    K = 0
    Q = 0
    E = 0
    S = 0
    P = 0
    V = 0
    I = 0
    C = 0
    Y = 0
    H = 0
    R = 0
    N = 0
    D = 0
    T = 0

#this resets the variables for the next loop over the next file

The python script establishes a bunch of variables at the beginning to be zero, and then on a loop searches the file for certain characters, it then uses the count of the characters, to calculate the percentage of each of these characters, over the whole file, and then write it to a next text file. After this is complete, it then resets the variables to 0, and then does it again on the loop.

Ideally, I'd like to use dictionaries to streamline this code. but I'm not sure where to start.

Example input (a snippet from a .faa file):

>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
>WP_013179458.1 hypothetical protein [Methanococcus voltae]
MVLNLEDLDKLDSIFSDGGIDKIENKTKNYNNDSDSFNVLDALKEVNKIFENWRSIRGIPKAQNIQPLKEYQVSKEKQTE
VKKDSNEITNVSNTNNINKNISAQDIYDNFLQALEFFKSSYGDMPVSEMVSTLKENKEDILSVINLSMGDVNGA
>WP_013179459.1 hypothetical protein [Methanococcus voltae]
MGRKPLDPKAIKKKLEEHEAGTIKLPYSTLQRYKNTLDKQDLKEKDEEYKQNIDLDDELNNIDLNSEYVNYYDIIDFKNP
FSMCVFGIKRQGKTTLLKHMAYSNQKDVLIYDLVHNFNNFGKRCYQAKETQFPDSALEFQRFYSSIYNKLNKNRENPILL
LIDECDKIAPNNSRMPGGLAELNDLHRHAKFNTSFVCVARRPATLNTNLKEIADYIIFFKLTGKNDKSFLNDLHKDLYNA
VESLNAEEHEFIIYDMPMSKIYKSKLDLNINFKK
>WP_013179460.1 hypothetical protein [Methanococcus voltae]
MTKTINGINFKAMGIVTISKVVGEQVLTPIIGNGTVKSGLPKILGAVLLAGTKNTYAKYAGTGLAVDGIEDILMGSGILS
KLGAVAGAKTTAGTGNTNNIDIM
>WP_013179461.1 hypothetical protein [Methanococcus voltae]
MAVVKPGNGDPSVLGLNDFEFNAKGDTIKAGRWTDIYKFTVPVQEQVAIGSDDNGNVGILYGIIKDNSETPAEVSGVIRI
SRRPARENVSDRQLEVRTEMIKDTMTDRLKAYFLPVKRNKRIGENSKLVIEFMPDTDFVLGDSVLQIPITRW

Example output

G,6.2152758802848975,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
A,5.358317495592757,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
L,9.04238847295845,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
M,2.5514448269790093,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
F,3.7227292323199204,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
W,0.5889817901403234,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
K,9.248430899887264,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
Q,2.067103444227877,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
E,7.731246192364286,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
S,6.133911385564481,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
P,3.1477129897808624,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
V,6.35290736389157,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
I,9.564983312182612,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
C,1.347169345420616,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
H,1.4429041862234933,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
N,7.343449247378424,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
R,2.929526608416165,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
D,5.753603212480747,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
T,5.148429483092579,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
Y,4.309484630813666,GCF_000006175.1_ASM617v2_protein.faa,Methanococcus voltae
Improve this question
\$\endgroup\$
2
  • \$\begingroup\$ Can your provide a sample faa file and the output? \$\endgroup\$
    – hjpotter92
    Jul 5, 2018 at 9:18
  • 1
    \$\begingroup\$ I have edited to add example input, and output and pasted the raw-code, with comments (are comments allowed here?) which definitely works! \$\endgroup\$
    – Biomage
    Jul 5, 2018 at 9:30

1 Answer 1

Reset to default
10
\$\begingroup\$

  1. Add functions

    Everything is currently in the global namespace, you can add a few functions to split up the code and make it more readable.

    You could use a count_occurence function and a write_output function.

  2. Avoid working with global variables, it can make the code less readable, secondly these are not really global variables.

  3. Make use of a dictionary, to count the amino acids.

  4. Remove unused imports.

  5. Remove unused variables:

    GAP = 0
SEQCOUNTERPERALIGNMENT = 0
AMINOACIDSPERSTRAIN = 0
NONE = 0

    You only use the AMINOACIDSPERTRAIN variable, the rest is not used in this script.

  6. Use a if __name__ == '__main__' guard, it is standard idiom in Python.

  7. Use str.format(x) or f"string{x}" instead of string concatenation.

  8. Count the percentages only after the entire file is parsed, you are only writing the output after the entire file is parsed.

  9. You never close the FILE, you can use with open(filename) as f: context manager to let Python close the file automatically after the file is used.

Revised script

from Bio import SeqIO
import glob
import time

def count_occurence(f, amino_acids_to_track='ACDEFGHIKLMNPQRSTVWY'):
    amino_acids_per_train = 0
    amino_acids = dict.fromkeys(amino_acids_to_track, 0)
    # or amino_acids = { k: 0 for k in amino_acids_to_track}

    for record in SeqIO.parse(f, "fasta"):
        for char in record.seq:
            if char in amino_acids:
                amino_acids_per_train += 1
                amino_acids[char] += 1

    percentages = { k : (v * 100.0) / amino_acids_per_train for k, v in amino_acids.items()}
    speciesname = record.description.split('[', 1)[1].split(']', 1)[0]
    return percentages, speciesname

def write_output(out_filename, from_file, percentages, species_name):
    with open(out_filename, 'a+') as out_file:
        for k, v in percentages.items():
            out_file.write('{} is {}% in {}, {}\n'.format(k, v, from_file, species_name))

if __name__ == '__main__':
    current_time = time.strftime("%d.%m.%y %H:%M", time.localtime())
    out_file= 'test#{}.txt'.format(current_time)
    files = glob.glob('*.faa')
    for f in files:
        percentages, species_name = count_occurence(f)
        write_output(out_file, f, percentages, species_name)
Improve this answer
\$\endgroup\$
13
  • 1
    \$\begingroup\$ Does the revised script work for you? I'm getting this error which I'm not sure I understand: File "propgen3.py", line 30, in <module> percentages, species_name = count_occurence(f) File "propgen3.py", line 16, in count_occurence percentages = { k : (v * 100.0) / amino_acids_per_train for k, v in amino_acids.iteritems()} AttributeError: 'dict' object has no attribute 'iteritems' \$\endgroup\$
    – Biomage
    Jul 5, 2018 at 10:14
  • 1
    \$\begingroup\$ @Biomage This is a Python 2 / Python 3 issue. This code assume Python 2, you seems to run Python 3. Replace iteritems with items. \$\endgroup\$
    – 301_Moved_Permanently
    Jul 5, 2018 at 10:15
  • 1
    \$\begingroup\$ Change to iteritems to items i assumed you were using python2.7 \$\endgroup\$
    – Ludisposed
    Jul 5, 2018 at 10:16
  • 1
    \$\begingroup\$ @Biomage You should have a look at argparse and/or fileinput \$\endgroup\$
    – 301_Moved_Permanently
    Jul 5, 2018 at 10:29
  • 1
    \$\begingroup\$ @Biomage Click is super nice \$\endgroup\$
    – Martin Thoma
    Jul 5, 2018 at 11:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.