pattern finding algorithm

[kdt]

Hi,

I'm checking to see if you guys may be able to help me with an algorithm for finding patterns. I have around 2000 short sequences (of length 9) that are aligned. I want to be able to extract all common patterns on the same positions and report the number of occurrences.

For example in the following:

ACGCATTCA
ACTGGATAC
TCAGCCATC

I would like the following output (where a full stop represents any character)

(AC....T..) 2 occurrences (pattern between sequence 1 and 2)
(.C.G...C) 2 occurrences (pattern between sequence 2 and 3)
(.C.......) 2 occurrences (pattern between sequence 1 and 3)

As you can see, the way that I am planning on doing this now requires sum(n-1...1) comparisons. Is there a more efficient way of doing this with less comparisons?

Thanks

[bvdet]

Hi,

A quick question on the above algorithm, it seems to performing more comparisons than necessary, though I can't seem to find where.

I decided to experiment to get the number of comparisons required to compare each line with every other line - and not against itself (line 1 to 1) or lines already compared ( as comparing line 1 to 2 is the same as comparing line 2 to 1).

line_num = 0
counter = 0
 
for line in myList:
    elmt_num = 0
 
    for elmt in myList:
        i = 0
        if line_num < elmt_num # As I don't want to compare twice and against the same seq
            while i<9:
                counter += 1
                i += 1
 
        elmt_num += 1
 
    line_num += 1
 
print "%d comparisons were made between %d lines" % (counter, line_num)
 

I get the following output for my dataset

7587459 comparisons were made between 1299 lines

Which is exactly what I would expect from comparing 9 times (n**2 -n)/2 where n is the number of sequences.

The output I get from bvdet's algorithm is however 8430510 comparisons (10 time (n**2-n)/2) I fail to see where the extra 1 comparison each time is coming from.

If someone can let me know, I would be thankful.

Cheers

kdt - I ran the following function and counted the comparisons:

def patt_match(sList):
    global count
    count = 0
    sList = sList[:]
    patt = re.compile('[ACGT]')
    dd = {}
    indx = 0
    sList = strList[:]
    while len(sList) > 0:
        s1 = sList[0]
        for j, item in enumerate(sList[1:]):
            res = ''
            for i, s in enumerate(s1):
                count += 1
                if s == item[i]:
                    res += s
                else:
                    res += '.'
            if patt.search(res):
                if dd.has_key(res):
                    dd[res].append([indx, j+1+indx])
                else:
                    dd[res] = [[indx, j+1+indx], ]
        indx += 1      
        sList.pop(0)
    return dd

Output:

>>> dd = patt_match(strL ist)
>>> count
7587459
>>>

[kdt]

hi bvdet,

thanks for getting back to me on this. Unfortunately when I run your script I still get the same results. Would it be possible for you to check the number of lines by adding another counter after the line

while len(sList) > 0:

Thanks

[bartonc]

hi bvdet,

thanks for getting back to me on this. Unfortunately when I run your script I still get the same results. Would it be possible for you to check the number of lines by adding another counter after the line

while len(sList) > 0:

Thanks

# line 2 & 3 would read
    global count linecount
    count = linecount = 0
# insert at the level of (current line 10) - this looks like the lines to me -
        linecount += 1

[bvdet]

hi bvdet,

thanks for getting back to me on this. Unfortunately when I run your script I still get the same results. Would it be possible for you to check the number of lines by adding another counter after the line

while len(sList) > 0:

Thanks

Here are the results:

>>> len(strList)
1299
>>> dd = patt_match(strL ist)
>>> len(dd)
53939
>>> count
7587459
>>> linecount
843051
>>> itemcount
1299
>>>

Here is the function:

def patt_match(sList):
    global count, linecount, itemcount
    count = linecount = itemcount = 0
    sList = sList[:]
    patt = re.compile('[ACGT]')
    dd = {}
    indx = 0
    while len(sList) > 0:
        s1 = sList[0]
        for j, item in enumerate(sList[1:]):
            res = ''
            for i, s in enumerate(s1):
                count += 1
                if s == item[i]:
                    res += s
                else:
                    res += '.'
            if patt.search(res):
                if dd.has_key(res):
                    dd[res].append([indx, j+1+indx])
                else:
                    dd[res] = [[indx, j+1+indx], ]
            linecount += 1
        indx += 1
        itemcount += 1
        sList.pop(0)
    return dd