What's the fastest way to find all occurences of the substring in large string in python
<p>Edited to clarify input/output. I think this will be somewhat slow no matter what, but up until now I haven't really considered speed in my python scripts, and I'm trying to figure out ways to speed up operations like these.</p>

<p>My input is pickled dictionaries of the genome sequences. I'm currently working with two genomes, the budding yeast genome (11.5 MB on disk), and the human genome (2.8 GB on disk). These dictionaries have the form:</p>

<pre><code>seq_d = { 'chr1' : 'ATCGCTCGCTGCTCGCT', 'chr2' : 'CGATCAGTCATGCATGCAT', 
        'chr3' : 'ACTCATCATCATCATACTGGC' }
</code></pre>

<p>I want to find all single-base instances of a nucleotide(s) in both strands of the genome. Where the <code>+</code> strand refers to the sequence in the above dictionary, and the <code>-</code> strand is the reverse complement of the sequences. My output is a nested dictionary, where the top level <code>keys</code> are <code>+</code> or <code>-</code>, the nested <code>keys</code> are chromosome names, and the values are lists of 0-indexed positions:</p>

<pre><code>nts = 'T'
test_d = {'+': {'chr3': [2, 5, 8, 11, 14, 17], 'chr2': [3, 7, 10, 14, 18], 
           'chr1': [1, 5, 9, 12, 16]}, '-': {'chr3': [0, 4, 7, 10, 13, 15], 
           'chr2': [2, 5, 9, 13, 17], 'chr1': [0]}}
</code></pre>

<p><code>test_d</code> defines a set of positions to examine in a large Illumina sequencing dataset later in the script.</p>

<p>My first attempt uses <code>enumerate</code>, and iteration.</p>

<pre><code>import time
import numpy as np

rev_comps = { 'A' : 'T', 'T' : 'A', 'G' : 'C', 'C' : 'G', 'N' : 'N'}
test_d = { '+' : {}, '-' : {}}
nts = 'T'

s = time.time()
for chrom in seq_d:
    plus_pos, minus_pos = [], []
    chrom_seq = seq_d[chrom]
    for pos, nt in enumerate(chrom_seq):
        if nt in nts:
            plus_pos.append(pos)
        if rev_comps[nt] in nts:
            minus_pos.append(pos)

    test_d['+'][chrom] = plus_pos
    test_d['-'][chrom] = minus_pos
e = time.time()
print 'The serial version took {} minutes...'.format((e-s)/60)
</code></pre>

<p>The output for yeast:</p>

<pre><code>The serial version took 0.0455190300941 minutes...
</code></pre>

<p>The output for human:</p>

<pre><code>The serial version took 10.1694028815 minutes...
</code></pre>

<p>I tried using numpy arrays rather than <code>enumerate()</code> and iteration:</p>

<pre><code>s = time.time()
for chrom in seq_d:
    chrom_seq = np.array(list(seq_d[chrom]))
    nts = list(nts)
    rev_nts = [rev_comps[nt] for nt in nts]
    plus_pos = list(np.where(np.in1d(chrom_seq, nts) == True)[0])
    minus_pos = list(np.where(np.in1d(chrom_seq, rev_nts) == True)[0])

    test_d['+'][chrom] = plus_pos
    test_d['-'][chrom] = minus_pos
e = time.time()
print 'The numpy version took {} minutes...'.format((e-s)/60)
</code></pre>

<p>The output for yeast:</p>

<pre><code>The numpy version took 0.0309354345004 minutes...
</code></pre>

<p>The output for human:</p>

<pre><code>The numpy version took 9.86174853643 minutes...
</code></pre>

<p>Why does the numpy version lose its performance advantage for the larger human genome? Is there a faster way to do this? I tried implementing a version using <code>multiprocessing.Pool</code>, but that's slower than either version:</p>

<pre><code>def getTestHelper(chrom_seq, nts, rev_comp):

    rev_comps = { 'A' : 'T', 'T' : 'A', 'G' : 'C', 'C' : 'G', 'N' : 'N'}
    if rev_comp:
        nts = [rev_comps[nt] for nt in nts]
    else:
        nts = list(nts)
    chrom_seq = np.array(list(chrom_seq))
    mask = np.in1d(chrom_seq, nts)
    pos_l = list(np.where(mask == True)[0])
    return pos_l

s = time.time()
pool = Pool(4)
plus_pos = pool.map(functools.partial(getTestHelper, nts=nts, rev_comp=False), seq_d.values())
minus_pos = pool.map(functools.partial(getTestHelper, nts=nts, rev_comp=True), seq_d.values())
e = time.time()
print 'The parallel version took {} minutes...'.format((e-s)/60)
</code></pre>

<p>I haven't run this on the human genome, but the yeast version is slower:</p>

<pre><code>The parallel version took 0.652778700987 minutes...
</code></pre>

<p>To be honest, I think you've been doing the right things.</p>

<p>There are a few more tweaks you can make to your code, though.  In general, when performance is key, only do the bare minimum in your innermost loops.  Looking through your code, there are still some quick optimizations left on this front:</p>

<ol>
<li>Use if...elif instead of if...if.</li>
<li>Don't use lists where you don't have to - e.g. just a single string is sufficient for nts and the reverse.</li>
<li>Don't evaluate the same result multiple times - e.g. the reverse lookup.</li>
</ol>

<p>I'm guessing that your problem with multi-processing is down to the serialization of these very large strings, offsetting any performance gain you might have from running in parallel.  However, there may be another way to do this - see <a href="http://stackoverflow.com/questions/11442191/parallelizing-a-numpy-vector-operation">Parallelizing a Numpy vector operation</a>.  I can't verify as I am having difficulty installing <code>numexpr</code>.</p>

<p>Putting them together and trying out some of the other suggestions in this trail, I get the following results:</p>

<pre><code>$ python test.py
Original serial took 0.08330821593602498 minutes...
Using sets took 0.09072601397832235 minutes...
Using built-ins took 0.061421032746632895 minutes...
Using regex took 0.11649663050969442 minutes...
Optimized serial took 0.05909080108006795 minutes...
Original numpy took 0.04050511916478475 minutes...
Optimized numpy took 0.03438538312911987 minutes...
</code></pre>

<p>My code is as follows.</p>

<pre><code>import time
import numpy as np
from random import choice
import re

# Create single large chromosome for the test...
seq = ""
for i in range(10000000):
    seq += choice("ATGCN")
seq_d = {"Chromosome1": seq}


rev_comps = { 'A' : 'T', 'T' : 'A', 'G' : 'C', 'C' : 'G', 'N' : 'N'}
nts = 'T'

# Original serial implementation
def serial():
    test_d = { '+' : {}, '-' : {}}
    for chrom in seq_d:
        plus_pos, minus_pos = [], []
        chrom_seq = seq_d[chrom]
        for pos, nt in enumerate(chrom_seq):
            if nt in nts:
                plus_pos.append(pos)
            if rev_comps[nt] in nts:
                minus_pos.append(pos)

        test_d['+'][chrom] = plus_pos
        test_d['-'][chrom] = minus_pos

# Optimized for single character tests
def serial2():
    test_d = {'+': {}, '-': {}}
    rev_nts = rev_comps[nts]
    for chrom in seq_d:
        plus_pos, minus_pos = [], []
        chrom_seq = seq_d[chrom]
        for pos, nt in enumerate(chrom_seq):
            if nt == nts:
                plus_pos.append(pos)
            elif nt == rev_nts:
                minus_pos.append(pos)

        test_d['+'][chrom] = plus_pos
        test_d['-'][chrom] = minus_pos

# Use sets instead of lists
def set_style():
    test_d = { '+' : {}, '-' : {}}
    for chrom in seq_d:
        plus_pos, minus_pos = set(), set()
        chrom_seq = seq_d[chrom]
        for pos, nt in enumerate(chrom_seq):
            if nt in nts:
                plus_pos.add(pos)
            if rev_comps[nt] in nts:
                minus_pos.add(pos)

        test_d['+'][chrom] = plus_pos
        test_d['-'][chrom] = minus_pos

# Use regex to find either nucleotide...
def regex_it():
    test_d = { '+' : {}, '-' : {}}
    search = re.compile("(T|A)")
    for chrom in seq_d:
        pos = 0
        plus_pos, minus_pos = [], []
        chrom_seq = seq_d[chrom]
        for sub_seq in search.split(chrom_seq):
            if len(sub_seq) == 0:
                continue
            if sub_seq[0] == 'T':
                plus_pos.append(pos)
            elif sub_seq[0] == 'A':
                minus_pos.append(pos)
            pos += len(sub_seq)

        test_d['+'][chrom] = plus_pos
        test_d['-'][chrom] = minus_pos

# Use str.find instead of iteration
def use_builtins():
    test_d = { '+' : {}, '-' : {}}
    for chrom in seq_d:
        plus_pos, minus_pos = [], []
        chrom_seq = seq_d[chrom]
        start = 0
        while True:
            pos = chrom_seq.find("T", start)
            if pos == -1:
                break
            plus_pos.append(pos)
            start = pos + 1

        start = 0
        while True:
            pos = chrom_seq.find("A", start)
            if pos == -1:
                break
            minus_pos.append(pos)
            start = pos + 1

        test_d['+'][chrom] = plus_pos
        test_d['-'][chrom] = minus_pos

# Original numpy implementation
def numpy1():
    test_d = { '+' : {}, '-' : {}}
    for chrom in seq_d:
        chrom_seq = np.array(list(seq_d[chrom]))
        for_nts = list(nts)
        rev_nts = [rev_comps[nt] for nt in nts]
        plus_pos = list(np.where(np.in1d(chrom_seq, for_nts) == True)[0])
        minus_pos = list(np.where(np.in1d(chrom_seq, rev_nts) == True)[0])

        test_d['+'][chrom] = plus_pos
        test_d['-'][chrom] = minus_pos

# Optimized for single character look-ups
def numpy2():
    test_d = {'+': {}, '-': {}}
    rev_nts = rev_comps[nts]
    for chrom in seq_d:
        chrom_seq = np.array(list(seq_d[chrom]))
        plus_pos = np.where(chrom_seq == nts)
        minus_pos = np.where(chrom_seq == rev_nts)

        test_d['+'][chrom] = plus_pos
        test_d['-'][chrom] = minus_pos

for fn, name in [
        (serial, "Original serial"),
        (set_style, "Using sets"),
        (use_builtins, "Using built-ins"),
        (regex_it, "Using regex"),
        (serial2, "Optimized serial"),
        (numpy1, "Original numpy"),
        (numpy2, "Optimized numpy")]:
    s = time.time()
    fn()
    e = time.time()
    print('{} took {} minutes...'.format(name, (e-s)/60))
</code></pre>

python