3
\$\begingroup\$

The script takes two input files of protein sequences and runs an external program (installed in linux/MacOS). The result provides a text output file example output.. Identity percentage is extracted from the matched keyword "Identity" in the line. Also, the script returns the result text file.

Can you help me to improve this script?

import re
import subprocess

def cmdline(command):
    process = subprocess.Popen(
        args=command,
        stdout=subprocess.PIPE,
        shell=True
    )
    return process.communicate()[0]

def blast_two_sequences(file1, file2):

    try:
        # fetch the identity percentage
        #https://www.biostars.org/p/200256/

        cmd = 'needle -datafile EBLOSUM62 -auto Y' + ' -asequence ' + file1 +' -bsequence ' + file2 +  ' -sprotein1 Y -sprotein2 Y ' + ' -auto -stdout'
        # cmd = 'needle -asequence ' + file1 +' -bsequence ' + file2 + ' -auto -stdout'

        results = cmdline(cmd)

        results = results.decode("utf-8")


        identity = re.search(r"\d{1,3}\.\d*\%", results)
        identity = identity.group()
        identity = identity.replace('%', '')

        return identity, results

    except:
        return None


file1 = "/Users/catuf/Desktop/file1.txt"
file2 = "/Users/catuf/Desktop/file2.txt"

if blast_two_sequences(file1, file2) is not None:
    identity, result = blast_two_sequences(file1, file2)
\$\endgroup\$
5
\$\begingroup\$

Not looking bad for as far as I can see. If the example file is accurate for the lengths of the input files, then I don't forsee any real problems, though others may of course disagree.

Naming:

cmdline is quite a... short name for a function. I'd think cmd_line for snake_case convention. However, what it does is create what's basically a text file in a bytes buffer by means of running a program. While you could call that a command line execution if you're working from a terminal, I'd name it after the invoked program, or something more generic like get_input_data. Since generic is always survivable at best, I'd pick something like get_needle_results or something like that.

For a module like subprocess which you can guess does creative things dependent on what platform you're on, it's always advisable to grab the docs:

args should be a sequence of program arguments or else a single string. By default, the program to execute is the first item in args if args is a sequence. If args is a string, the interpretation is platform-dependent and described below. See the shell and executable arguments for additional differences from the default behavior. Unless otherwise stated, it is recommended to pass args as a sequence.

And your platform matters - a lot. If you never switch platform or needle implementation, then you can do what you're doing, but if not you're risking a collossal headache with potentially giving it wrong arguments. Therefore, you might want to restructure it to a list of arguments.

And if you don't, use f-strings, the best invention since the list comprehension:

cmd = 'needle -datafile EBLOSUM62 -auto Y' + ' -asequence ' + file1 +' -bsequence ' + file2 +  ' -sprotein1 Y -sprotein2 Y ' + ' -auto -stdout'
# Becomes:
cmd = f'needle -datafile EBLOSUM62 -auto Y -asequence {file1} -bsequence {file2} -sprotein1 Y -sprotein2 Y  -auto -stdout'

Which is shorter and more readable.

        results = cmdline(cmd)

        results = results.decode("utf-8")

These can be combined:

        results = cmdline(cmd).decode("utf-8")

So can these:

        identity = re.search(r"\d{1,3}\.\d*\%", results)
        identity = identity.group()
        identity = identity.replace('%', '')
# New:
        identity = re.search(r"\d{1,3}\.\d*\%", results).group().replace('%', '')
if blast_two_sequences(file1, file2) is not None:
    identity, result = blast_two_sequences(file1, file2)

And ouch... you're actually running your entire program twice. If your needle executions are counted, you might have noticed this. You should refactor this to:

output = blast_two_sequences(file1, file2)
if output is not None:
    identity, result = output 
# Or if you work with Python 3.8+
if output := blast_two_sequences(file1, file2) is not None:
    identity, result = output 

You should also put a guard around the actual execution, so that a later scripter can import your functions and use them without executing your script on import:

if __name__ == "__main__":
    file1 = "/Users/catuf/Desktop/file1.txt"
    file2 = "/Users/catuf/Desktop/file2.txt"
    if output := blast_two_sequences(file1, file2) is not None:
        identity, result = output 

Exceptions:

There's a small issue you have - you have a bare except in your code. That means that if weird stuff happens, your program will just ignore it and return none. In your case, I would recommend removing the try: statement completely, and you'll be informed about any problems that arise. If you're watching out for specific Exceptions, narrow your except: to those, and preferably restrict your try/except to directly around the origin lines of these Exceptions.

This will obsolete your None check on the outside. Exceptions showing up on your screen is a good thing - it means you'll be notified when the unexpected happens and you have some unknown error - because now you can go fix it.

\$\endgroup\$
  • \$\begingroup\$ Thank you very much for doing this. It help me a lot. \$\endgroup\$ – catuf Sep 12 at 12:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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