Python program to check substring match locations with a lot of permutations of the substring

Question

This is my first code here and also my first program of this size. I don't really know what is expected from a programmer who writes "good, readable" code. This is my first program which will be used in a real-world application. Also I'm extremely new to Python. So while reviewing please be kind enough to give constructive criticism about how this code or my code in general can be better with respect to both Python and programming in general. I'll try to explain the problem in the following paragraph in the best way possible. If any clarifications are needed about my code/logic/the problem, feel free to ask in the comments, I'll try my best to clear the doubts.

The problem -

• Consider two files.
• Each containing a list of strings.
• One list has strings of some combination of 'a','t','g', and 'c'
• One list has strings of some combination of 'A',U','G', and 'C'
• I have to convert the strings from the capitalized list as a's to t's, c's to g's, u's to a's and g's to c's [ a-t, c-g, g-c, u-a ]. And another special condition is that there can be at most of two instances where u's convert to g's and/or g's convert to t's [ u-g, g-t ]
• The conversions only need to be done for four regions of the strings, indices 2-7(6 characters), 2-8(7 characters), 1-7(7 characters) and 1-8(8 characters), provided the starting index is 1
• After generating all possible conversions, I have to check each of them against all the strings in the other list and find out the locations where they match.

If you are looking for an output of sorts, I won't be able to provide it yet, since the comparisons I need to do are about (38869 * 2588 * all possibble combinatons of each of the 2588) + time taken to generate all the permutations. So my machine is extremely inadequate of doing something like that.

My Program -

## Date : 2017-08-10
## Author : dadyodevil
## Contact : [email protected]
##
## A python program to detect all indices of complimentary Micro-RNA(miRNA) target sites on Messenger-RNAs(mRNA)
##
## As an input, this program needs two lists - 
##  1. A list of mRNAs where each entry is represented in a two line format:
##      >hg19_refGene NM_032291 range=chr1:67208779-67210768...
##      Sequence of mRNA
##   2. A list of miRNAs where each entry is represented in a two line format:
##      >hsa-miR-576-3p MIMAT000...
##      Sequence of miRNA
##
##  Pre-requisites for the reader -  
##  1. Understanding of programming concepts
##  2. A moderate understanding of the Python programming language version 2.7
##  3. Knowledge of terms regarding miRNA-mRNA target detection


import re


def extractSeed(miRNA):

    ## There are 4 seed regions with indices from 2-7, 2-8, 1-7 and 1-8

    miRNAfor6mer.append(miRNA[1:7][::-1])
    miRNAfor7mer.append(miRNA[1:8][::-1])
    miRNAfor7a1.append(miRNA[:7][::-1])
    miRNAfor8mer.append(miRNA[0:9][::-1])


def createCompliment(allCompliments, miRNA, wobbleCount, compliment):

    ## For the compliment, the convertions include a:t, u:a, g:c, c:g and for Wobble-Pairs, u:g and g:u

    if wobbleCount == 2:        
        for letter in miRNA:
            if letter == 'a':
                compliment += 't'
            elif letter == 'c':
                compliment += 'g'
            elif letter == 'g':
                compliment += 'c'
            else:
                compliment += 'a'
        allCompliments.append(compliment)

    else:
        for index, letter in enumerate(miRNA):
            if letter == 'a':
                compliment += 't'
            elif letter == 'c':
                compliment += 'g'
            elif letter == 'g':
                createCompliment(allCompliments, miRNA[index+1:], wobbleCount + 1, compliment + "t")
                createCompliment(allCompliments, miRNA[index+1:], wobbleCount + 1, compliment + "c")
                compliment += 'c'
            elif letter == 'u':
                createCompliment(allCompliments, miRNA[index+1:], wobbleCount + 1, compliment + "g")
                createCompliment(allCompliments, miRNA[index+1:], wobbleCount + 1, compliment + "a")
                compliment += 'a'

    ## Now that all possibilities are generated, the duplicates need to be removed
    allCompliments = sorted(list(set(allCompliments)))


def checkForMatch(miRNACompliments, seedRegion, miRNAname):

    ## Each miRNA that is recived by this function will be compared against the whole list of mRNAs and the matching indices will be saved

    ## Since the mRNA sequences are in alternate lines the sequences will be extracted as such and the when matches are found, the name of the mRNA will be extracted from teh index just before the current one

    for index in range(1, len(mRNA_List), 2):
        for entry in miRNACompliments:
            mRNA =  mRNA_List[index]
            matchesStart = [m.start() for m in re.finditer(entry, mRNA)]

            if (len(matchesStart) > 0):             
                mRNAname = mRNA_List[index-1][14:mRNA_List[index-1].find(" ",15)]
                matchesEnd = []             
                for index2 in range(0, len(matchesStart)):
                    matchesEnd.append(matchesStart[index2] + len(entry))
                allindices = zip(matchesStart, matchesEnd)
                complimentarySiteList.append([miRNAname, mRNAname, seedRegion, allindices])


def prepareForMatch(miRNA, miRNAname):

    global miRNAfor6mer, miRNAfor7mer, miRNAfor7a1, miRNAfor8mer
    miRNAfor6mer, miRNAfor7mer, miRNAfor7a1, miRNAfor8mer = [], [], [], []

    ## First the seed sites will be extracted and reversed
    extractSeed(miRNA)

    ## Empty lists will be generated to store all the compliments
    miRNAfor6mer.append([])
    miRNAfor7mer.append([])
    miRNAfor7a1.append([])
    miRNAfor8mer.append([])

    ## Then the compliments will be generated from the seed regions along with atmost of two Wobble-Pairs
    miRNAfor6mer.append(createCompliment(miRNAfor6mer[1], miRNAfor6mer[0], 0, ""))
    miRNAfor7mer.append(createCompliment(miRNAfor7mer[1], miRNAfor7mer[0], 0, ""))
    miRNAfor7a1.append(createCompliment(miRNAfor7a1[1], miRNAfor7a1[0], 0, ""))
    miRNAfor8mer.append(createCompliment(miRNAfor8mer[1], miRNAfor8mer[0], 0, ""))

    ## After generating all possible compliments, they will be checked for matching sites
    checkForMatch(miRNAfor6mer[1], "6mer", miRNAname)
    checkForMatch(miRNAfor7mer[1], "7mer", miRNAname)
    checkForMatch(miRNAfor7a1[1], "7A1", miRNAname)
    checkForMatch(miRNAfor8mer[1], "8mer", miRNAname)


def Main():

    global mRNA_List, miRNA_List, complimentarySiteList
    miRNA_List = open('miRNA_list.txt').read().splitlines()
    mRNA_List = open('mRNA_list.txt').read().splitlines()
    complimentarySiteList = []

    ## Since the sequences are in every alteRNAte lines, the 'index' needs to be incremeted by 2 to access only the sequences

    ## The miRNA lengths are also checked whether they are atleast 8 neucleotides long, if they are not, they will not be checked

    for index in range(1,len(miRNA_List),2):

        miRNAname = miRNA_List[index-1][5:miRNA_List[index-1].find(' ')]

        if (len(miRNA_List[index]) < 8):
            print "%s at %d has insufficient length." %(miRNAname, index)

        else:
            prepareForMatch(miRNA_List[index].lower(), miRNAname)

    for entry in complimentarySiteList:
        print entry


if __name__ == '__main__':
    Main()

Answer 1

• Please stop abusing lists. Mutations are hard to reason with, and so when a function can mutate the data it makes the function a lot harder to understand.
• The function createCompliment is doesn't return anything, instead it mutates allCompliments. This mutation doesn't work when you assign to it, such as on the last line. allCompliments = sorted(list(set(allCompliments))). This does nothing, as you don't use it afterwards.

• Rather than using two slow for loops in createCompliment, you could instead do all the standard conversions, and then handle the special conversions, by looping through all the combinations of the special indexes.

  By this, if you perform the loop when if wobbleCount == 2:, first, so that you get the basic conversion, then you don't have to care about them, when doing the special conversions. By this, if you have the input ccagaa, then you convert it to ggtctt, without caring about g. The simplest way to do this would be to use str.translate.

  After that, you want to convert the special conversions, which are g -> t and u -> g. However, since we performed the above, they are c -> t and a -> g. To convert these, you want to get the indexes of c and a. Which you can do with a list comprehension. [i for i, c in enumerate(rna) if c in 'ac'].

  After this you want to convert all combinations of one to two occurrences of these characters. This means we can use itertools.combinations to loop through all the combinations we want to change.

  Finally, we have to convert the values, and so making a copy of the list using [:], which is a slice over the entire list. Then looping through the indexes, we can convert the value in the list, and yield the string version of the list.

  An example of this is:

  We start with rna = 'cagu', we convert it to the basic conversion, gtca. After this, we get all the indexes of the special characters, which are [2, 3]. We then go through all combinations of these, which are [(2,), (3,), (2, 3)], and yield the converted word, being gtta, gtcg, and gttg.

  The simplest way to think of yield is as array.append. So the following to functions, are kinda the same:

  def fn_1():
      yield 1
      yield 2
  
  def fn_2():
      array = []
      array.append(1)
      array.append(2)
      return array
  
  def fn_3():
      return [1, 2]
  
  list(fn_1()) == fn_2() == fn_3() # True
  

• To simplify check_for_match, as you're looping through a flattened list of 2d tuples, you want to un-flatten the list. So [0, 1, 2, 3] would become [(0, 1), (2, 3)], allowing for the simpler looping method of for a, b in ....

  To do this, you could use the grouper recipe:

  def grouper(iterable, n, fillvalue=None):
      "Collect data into fixed-length chunks or blocks"
      # grouper('ABCDEFG', 3, 'x') --> ABC DEF Gxx
      args = [iter(iterable)] * n
      return izip_longest(fillvalue=fillvalue, *args)
  

  The way this works is as [item] * n doesn't perform a copy on item, and by exploiting how iterators work. The former means that [item] * 2 is the same as [item, item], rather than say [item, copy(item)]. This is important, as this ensures that both items are the same iterator.

  Using a single iterator multiple times is important as zip basically uses [(next(it), next(it)), (next(it), next(it)), ...], it's a little more complex, as it works with any size it, and also knows when it stops. However it's pretty much how it works.

• You should follow PEP8.

And so I'd change your code to:

import re
import string
import itertools

TRANS = string.maketrans('acgu', 'tgca')
CONVS = {'a': 'g', 'c': 't'}

SEEDS = [
    "6mer",
    "7mer",
    "7A1",
    "8mer"
]


def create_compliments(rna):
    rna = rna.translate(TRANS)
    yield rna
    all_indexes = [i for i, c in enumerate(rna) if c in CONVS]
    rna = list(rna)
    for n in (1, 2):
        for indexes in itertools.combinations(all_indexes, n):
            t = rna[:]
            for index in indexes:
                t[index] = CONVS[t[index]]
            yield ''.join(t)


def grouper(iterable, n, fillvalue=None):
    "Collect data into fixed-length chunks or blocks"
    # grouper('ABCDEFG', 3, 'x') --> ABC DEF Gxx"
    args = [iter(iterable)] * n
    return itertools.izip_longest(*args, fillvalue=fillvalue)


def check_for_match(mi_RNAs, seed, mi_RNA_name, m_RNA_list):
    mi_RNAs = list(mi_RNAs)
    for m_RNA_name, m_RNA in grouper(m_RNA_list, 2):
        m_RNA_name = m_RNA_name[14:m_RNA_name.find(" ", 15)]
        for entry in mi_RNAs:
            matches = [m.start() for m in re.finditer(entry, m_RNA)]
            if matches:
                all_indices = tuple(
                    (match, match + len(entry))
                    for match in matches
                )
                yield mi_RNA_name, m_RNA_name, seed, all_indices


def prepare_for_match(mi_RNA, mi_RNA_name, m_RNA_list):
    mi_RNAs = [
        mi_RNA[1:7][::-1],
        mi_RNA[1:8][::-1],
        mi_RNA[:7][::-1],
        mi_RNA[0:9][::-1]
    ]

    for mi_RNA, seed in zip(mi_RNAs, SEEDS):
        for entry in check_for_match(create_compliments(mi_RNA), seed, mi_RNA_name, m_RNA_list):
            yield entry


def main():
    mi_RNA_list = open('miRNA_list.txt').read().splitlines()
    m_RNA_list = open('mRNA_list.txt').read().splitlines()

    for index in range(1, len(mi_RNA_list), 2):
        mi_RNA_name = mi_RNA_list[index-1][5:mi_RNA_list[index-1].find(' ')]
        if (len(mi_RNA_list[index]) < 8):
            print "{} at {} has insufficient length.".format(mi_RNA_name, index)
        else:
            for entry in prepare_for_match(mi_RNA_list[index].lower(), mi_RNA_name, m_RNA_list):
                print tuple(entry)


if __name__ == '__main__':
    main()

Answer 2

I will add some stuff to what @Peilonrayz has already mentioned.

Reading from a file here is done OK. But, after that you iterate only every other line. You can read only the odd lines from the file. This will make your code more readable and easier to understand.

So this:

miRNA_List = open('miRNA_list.txt').read().splitlines()

for index in range(1,len(miRNA_List),2):
    miRNAname = miRNA_List[index-1][5:miRNA_List[index-1].find(' ')]

becomes this:

with open('miRNA_list.txt') as f:
    mi_rna_list = [line.strip() for i, line in enumerate(f) if i % 2 == 1]

And regular iteration over mi_rna_list will give you what you need.

Please notice that if print %s at %d has insufficient length. is not absolutely essential then all of the length < 8 micro-rna can be removed ahead of time and save you some time iterating over dead-end cases. This doesn't seem a lot but when it comes to runtime the if statement won't be checked which will shave some time off you runtime.

A little note about this ^: This looks like I'm just pushing the if around the code but the effects of this is obvious when you remove the print. A simple example on my system cut the runtime by half.

In this case the code will look like:

with open('miRNA_list.txt') as f:
    mi_rna_list = [line.strip() for i, line in enumerate(f) if i % 2 == 1 and len(line) >= 8]

So now there is no need for the if condition checking since all of the micro-rna s are OK for check.

So now we are at:

for mi_rna in mi_rna_list:
    mi_rna_name = mi_rna[5:mi_rna.index(' ')]
    prepare_for_match(mi_rna.lower(), mi_rna_name)

but I really prefer returning a result rather then mutating a global variable. So let's make this happen:

for mi_rna in mi_rna_list:
    mi_rna_name = mi_rna[5:mi_rna.index(' ')]
    resutls.append(prepare_for_match(mi_rna.lower(), mi_rna_name))

One thing I think should be different from @Peilonrayz's answer is that I think you should return the result (per specific micro-rna) as a dictionary. That will allow you to "save" the "metadata" regarding that result. What I mean here, assuming that will be helpful to you, is that the results would look something like:

{
"miRNA name": "bla bla",
"compliments": [
    {
    "region": "2-7",
    "possible permutations" : [
        "first permutation",
        "second permutation",
        ...
        ]
    },
    {
    "region": "1-7",
    "possible permutations" : [
        "first permutation",
        "second permutation",
        ...
        ]
    },
    ...
    ]
}

Now the results are human readable and can also be used as an input for another function \ script \ whatever.

So our modified main function would look something like:

def main():

    results = []
    with open('miRNA_list.txt') as f:
        mi_rna_list = [line.strip() for i, line in enumerate(f) if i % 2 == 1 and len(line) >= 8]

    for mi_rna in mi_rna_list:
        mi_rna_name = mi_rna[5:mi_rna.index(' ')]
        results.append(prepare_for_match(mi_rna.lower(), mi_rna_name))

    for entry in results:
        print json.dumps(entry, indent=4)

The json.dumps part will arrange results variable in a human readable JSON form. Before you do that, don't forget to import json at the top of the file.

One extra little thing

Since this program has a lot of computing to do I suggest you try to make this program multi-threaded \ multi-processed. A good way to start that implementation would be with multiprocessing pool. I didn't ran tests of it on my system but I think this would also improve your runtime.