python counting number of presence and absence of substrings in list of sequences

Excellent question. This is a classic computer science problem. Yes, there is indeed a better algorithm. Yours processes each long string 16384 times. A better way is to process each long string only once.

Rather than searching for each motif within each long string, instead you should just record which motifs appear in each long string. For example, if you were searching for length 2 motifs in the following string:

s = 'ACGTAC'

then you could run a loop over the length 2 substrings and record which ones are present in a dict:

motifAppearances = {}for i in range(len(s)-1):    motif = s[i:i+2]                   # grab a length=2 substring    if motif not in motifAppearances:        motifAppearances[motif] = 0    # initialize the count    motifAppearances[motif] += 1       # increment the count

Now you've processed the entire string exactly once and found all the motifs present in it. In this case the resulting dict would look like:

motifAppearances = {'AC':2, 'CG':1, 'GT':1, 'TA':1}

Doing something similar for your case should cut down your run time by exactly 16384-fold.

A clean and very fast way (~15s with OP's data) would be to use the CountVectorizer of scikits-learn as it uses numpy under the hood, for example:

from sklearn.feature_extraction.text import CountVectorizerdef make_chunks(s):    width = 2    return [s[i:i+width] for i in range(len(s)-width+1)]l = ['ATTGCGGCTCACGAA', 'ACCTAGATACGACGG', 'CCCCTGTCCATGGTA']vectorizer = CountVectorizer(tokenizer=make_chunks)X = vectorizer.fit_transform(l)

Now X is a sparse matrix having all possible chunks as columns and sequence as rows, where every value is the number of occurrences of the given chunk in each sequence:

>>> X.toarray()# aa ac ag at ca cc cg ...[[1  1  0  1  1  0  2  1 1 2 1 0 0 1 1 1]      # ATTGCGGCTCACGAA [0  3  1  1  0  1  2  1 2 0 1 0 2 0 0 0]      # ACCTAGATACGACGG [0  0  0  1  1  4  0  1 0 0 1 2 1 1 2 0]]     # CCCCTGTCCATGGTA>>> (X.toarray()>0).astype(int)  # the same but counting only once per sequence[[1 1 0 1 1 0 1 1 1 1 1 0 0 1 1 1] [0 1 1 1 0 1 1 1 1 0 1 0 1 0 0 0] [0 0 0 1 1 1 0 1 0 0 1 1 1 1 1 0]]>>> vectorizer.get_feature_names()     # the columns(chunks)[u'aa', u'ac', u'ag', u'at', u'ca', u'cc', u'cg', u'ct', u'ga', u'gc', u'gg', u'gt', u'ta', u'tc', u'tg', u'tt']

Now you can sum along the columns, mask non-values or whatever manipulation you need to do, for example:

>>> X.sum(axis=0)[[1 4 1 3 2 5 4 3 3 2 3 2 3 2 3 1]]

Finally to find how many occurrences a given motif has, you must find the index of the corresponding motif/chunk and then evaluate in the previous sum:

>>> index = vectorizer.vocabulary_.get('ag')    # 'ag' is your motif2   # this means third column

In your case you would have to divide your list in two parts(positive and negative values) in order to include the down condition. I made a quick test with the list from DSM's answer and it takes around 3-4 seconds in my computer to run. If I use 12 000 length 4000 sequences, then it takes a little less than a minute.

EDIT: The whole code using OP's data can be found here.

There are several odd things about your code.

1. What you're calling "permutations" looks more like the Cartesian product, which can be computed using itertools.product.

2. Because Python is zero-indexed, the first element of a string is at index 0, and so the comparison like i[2].find(sMotif) < 1 will return True if the string is there right at the start, which seems a little odd.

3. Your OddsRatio, PValue and Enrichment calculations are inside the loop, but neither the zeroing of the counts nor the print is, which means you're computing them cumulatively for each new row but not doing anything with that information.

4. You recompute i[2].find(sMotif) multiple times in the typical case. That result isn't cached.

Assuming that I understand the numbers you're trying to compute -- and I could well be wrong, because there are several things you're doing that I don't understand -- I'd flip the logic. Instead of looping over each motif and trying to count it in each row, loop over each row and see what's there. That will be roughly 7*the number of rows instead of the number of motifs * the number of rows.

For example:

import randomfrom itertools import productfrom collections import defaultdict, CounterN = 12000datalength = 400listoflists = [[str(i), random.uniform(-1, 1),                 ''.join([random.choice('AGCT') for c in range(datalength)])]               for i in range(N)]def chunk(seq, width):    for i in range(len(seq)-width+1):        yield seq[i:i+width]def count_motifs(datatriples, width=7):    motif_counts_by_down = defaultdict(Counter)    nonmotif_counts_by_down = defaultdict(Counter)    all_motifs = set(''.join(p) for p in product('AGCT',repeat=width))    for symbol, value, sdata in datatriples:        down = value < -0.5        # what did we see?        motifs_seen = set(chunk(sdata, width))        # what didn't we see?        motifs_not_seen = all_motifs - motifs_seen        # accumulate these        motif_counts_by_down[down].update(motifs_seen)        nonmotif_counts_by_down[down].update(motifs_not_seen)    return motif_counts_by_down, nonmotif_counts_by_down

(I lowered the linelength just to make the output faster; if the line is 10 times longer, the code takes 10 times as long.)

This produces on my slow laptop (after inserting some linebreaks):

>>> %time mot, nomot = count_motifs(listoflists, 7)CPU times: user 1min 50s, sys: 60 ms, total: 1min 50sWall time: 1min 50s

So I'd figure about 20 minutes for the full problem, which isn't bad for so little code. (We could speed up the motifs_not_seen part by doing the arithmetic instead, but that'd only get us a factor of two anyway.)

On a much smaller case, where it's easier to see the output:

>>> mot, nomot = count_motifs(listoflists, 2)>>> motdefaultdict(<class 'collections.Counter'>, {False: Counter({'CG': 61, 'TC': 58, 'AT': 55, 'GT': 54, 'CA': 53, 'GA': 53, 'AC': 52, 'CT': 51, 'CC': 50, 'AG': 49, 'TA': 48, 'GC': 47, 'GG': 45, 'TG': 45, 'AA': 43, 'TT': 40}), True: Counter({'CT': 27, 'GT': 26, 'TC': 24, 'GC': 23, 'TA': 23, 'AC': 22, 'AG': 21, 'TG': 21, 'CC': 19, 'CG': 19, 'CA': 19, 'GG': 18, 'TT': 17, 'GA': 17, 'AA': 16, 'AT': 16})})>>> nomotdefaultdict(<class 'collections.Counter'>, {False: Counter({'TT': 31, 'AA': 28, 'GG': 26, 'TG': 26, 'GC': 24, 'TA': 23, 'AG': 22, 'CC': 21, 'CT': 20, 'AC': 19, 'GA': 18, 'CA': 18, 'GT': 17, 'AT': 16, 'TC': 13, 'CG': 10}), True: Counter({'AA': 13, 'AT': 13, 'GA': 12, 'TT': 12, 'GG': 11, 'CC': 10, 'CA': 10, 'CG': 10, 'AG': 8, 'TG': 8, 'AC': 7, 'GC': 6, 'TA': 6, 'TC': 5, 'GT': 3, 'CT': 2})})