Returning substrings (motifs) and original strings (sequences) from a file of strings (sequences)

Question

I would like to get some help/tips on writing better and more pythonic code, as well as variable naming.

Code info:

• .backtranscribe() is just a method to convert the U's in the sequence to T's (ATUG -> ACTG)
• a - 1-based index where substring starts
• b - size of substring
• c - text input of names of sequences to get substrings of (the data/sorted...fa file has hundreds of name/sequence pairs)

def comp(a, b, c):
    mirna_set = set(c.data.splitlines())
    motif_size = b
    results = []
    motif_set = set()
    handle = open("data/sorted.human.mature.fa", "rU")
    for record in SeqIO.parse(handle, "fasta"):
        if record.id in mirna_set: # check if db mirna is one in user input
            index_5p = a-1
            flag_motif_unique = False
            motif = ""
            while(flag_motif_unique == False):
                motif = record.seq.back_transcribe()[index_5p:index_5p+motif_size]
                if motif not in motif_set:
                    flag_motif_unique = True
                    motif_set.add(motif)
                index_5p += 1
            pos = "pos:" + str(index_5p) + "-" + str(index_5p+10)
            results.append(">" + record.id + " " + motif + " " + pos)
            results.append(record.seq.back_transcribe())
        else:
            continue
    handle.close()
    return(results)

Full app here (with input/output)

Answer 1

Disclaimer: I know very little about bioinformatics, so I’m reviewing this on the quality of the Python code. I’m not qualified to comment on the scientific accuracy.

So here are a few suggestions:

• You need better comments and documentation. Ideally a docstring on the function(s), and some more comments that explain why the code was written this way – how does it relate to the underlying scientific principles? As it stands, it’s fairly opaque, and so it would be hard to another scientist to verify that it’s a correct implementation.

  As a sidebar to this, you should pick much better variable names as the inputs to your function – a, b and c are utterly useless. Likewise the function name isn’t very helpful – lots of functions compute something, but what does this function compute?

  Perhaps something like:

  def gen_uniq_subsequences(index_5p, motif_size, mirna_list):
  

  More characters, but substantially easier to read.

• The output isn’t a particularly useful format. You’re returning a list of strings, which is fine for printing, but what if somebody wants to get the individual components of those results. Say, break out the position and the motif separately from the record ID?

  I’d have this function return some chunk of structured data (say, a namedtuple), and then the pretty printing goes off to a separate function. Makes it much easier to reuse the outputs of this function.

  Which leads me to my next comment…

• Your function is doing too much work. Ideally a function should have a single responsibility, an action that can be described in a single sentence. But your function does quite a lot:

  • Unpacks the third input from a string into a set.
  • Opens the sorted.human.mature.fa file, and starts parsing it.
  • Constructs the results list, which is some sort of scientific computation.
  • Assembles that into rigid strings for pretty-printing.

  Most of those actions don’t need to be in this function – pushing them out will allow you to simplify this part of the code, which makes it easier to review, read and maintain.

  For example, suppose I already have a set that contains my miRNAs. It would be much cleaner if I could pass that in directly, rather than packing it into a string so that you can unpack it again.

• You can simplify your looping constructs. Two places immediately jump out at me:

  • You don’t need the else: continue; if you omit this the you get the same behaviour.
  • If you slightly rearrange the while loop and add a break, you can get rid of the flag_motif_unique sentinel variable.

  Neither of these are ground-breaking, but they make your code a bit easier to write.