<p>This is beyond a single answer. 'Refactoring' (code efficiency, high performance code) is a huge topic. Like its own project.</p>
<p>The first thing I would personally absolutely target first is <code>re</code> here:</p>
<pre><code>    if re.search(id_,s, re.I):
</code></pre>
<p>Python is not Perl and the <code>re</code> is EXTREMELY slow and in particular using <code>search</code> is especially slow, most notably when uncompiled (see notes and point 3 below).</p>
<h3>How to avoid <code>re.search</code></h3>
<p>Could you revise this to a direct <code>C</code> style operation? This is the best solution if you get it working.
For example,</p>
<ol>
<li>Using <code>==</code>:</li>
</ol>
<p>e.g.</p>
<pre><code>if str(id_) == str(s):
</code></pre>
<p>if the case is important (like they're different) then...</p>
<pre><code>if str(id_).lower() == str(s).lower():
</code></pre>
<p>You could dump the <code>str()</code> it ain't doing much its just bug insurance.</p>
<p>The other approach if you want to stick with <code>re</code> is to switch to <code>match</code>, that will be quicker.</p>
<ol start="2">
<li>Using <code>match</code></li>
</ol>
<p>e.g.</p>
<pre><code>if re.match(id_, s, re.I):
</code></pre>
<p><code>match</code> anchors  the search at the first character, i.e. id_[0], <code>search</code> goes everywhere. <code>match</code> is faster.</p>
<ol start="3">
<li>Sticking with <code>search</code> but combining with compile:</li>
</ol>
<p>e.g.</p>
<pre><code>id_ = id_.strip()
x = re.compile (id_) 
if x.search(s, re.I): # don't use re. x replaces it
</code></pre>
<p>or</p>
<pre><code>id_ = id_.strip()
x = re.compile (id_)
if x.match(s, re.I):
</code></pre>
<p>I have not put these in order of preference. My preference is for 1. but it will risk throwing logical errors (missing stuff).</p>
<p><em>If you are using <code>re</code> you must compile the search term</em>. Its a pain but whilst compiled <code>re</code> is slow, uncompiled <code>re</code> is really, really slow. Disclaimer, I'm from a Perl background and from this perspective <code>re</code> is the pits because previously relied on the efficient regex (Perl).</p>
<p>The <code>==</code> and <code>match</code> could fall over if <code>&gt;</code> is in one variable and not in the other, but the <code>search</code> will return the correct 'match' (i.e. not the output of 'match') at the cost of algorithm efficiency.</p>
<p>This isn't an exhaustive analysis of <code>re</code> but thats where I would start.</p>
<p>Please note this is a starting point from prior experience to pick a 'big issues'. It is useful for future, particularly for questions, to place time checks across the code to identify the bottlenecks.</p>
<hr />
<p><strong>Notes</strong></p>
<ul>
<li>The thing we don't know is what the fasta headers look like.</li>
<li>For every query in <code>search</code> there are <span class="math-container">$1.10^6$</span> fasta headers. So if you've <span class="math-container">$1.10^5$</span> queries thats <span class="math-container">$1.10^{11}$</span> searches for <code>search</code> thats why it can be taking up a lot of time.</li>
<li>Rewriting the code using Biopython objects throughout is good practice. <code>SeqIO</code> in particular are cool. <code>AlignIO</code> are algorithmically expensive however.</li>
</ul>

CC BY-SA 4.0
2022-08-06T23:42:58.050Z
2022-08-07T03:47:08.363Z
2022-08-07T03:47:08.363Z
M__
<h1></h1>
<p>I wrote simple Python script which has an excessive run time, notably when using large data sets. For example a data set of 1 000 000 sequences.</p>
<p>I am seeking assistance refactoring the code and would be very grateful for any assistance.</p>
<pre><code>from Bio import SeqIO
import re

def move_item(old, new):
    k, v = old.popitem()
    new[k] = v

def map_ids_to_aa_sequence_list(id_):
    id_ = id_.strip()
    for s in aa:
        if re.search(id_,s, re.I):
            return s

file = open(r'list_of_ids', 'r').readlines()
names = [for name.strip() in file if len(name) &gt; 1]
        
handle = open('sequences.fasta','r')
sequences = list(SeqIO.parse(handle, 'fasta'))

seq_dict = {}

for s in sequences:
    seq_dict[str(s.description)] = str(s.seq)
    
aa = []

for k, v in seq_dict.items():
    sequence = '&gt;' + k +'\n' + v + '\n'
    aa.append(sequence)
    
separated_aminoacid_sequences_set = set(map(map_ids_to_aa_sequence_list, names))

separated_aminoacid_sequences_list = list(separated_aminoacid_sequences_set)
    
with open('separated_sequences.fasta','w') as fs:
    for s in separated_aminoacid_sequences_list:
        fs.write(str(s))
</code></pre>

CC BY-SA 4.0
2022-08-06T19:33:32.370Z
On a brief glance the problem could be `re` but I'd need to look. Could you kindly look at this code here https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example 
Are there any parallels?
M__
CC BY-SA 4.0
2022-08-07T07:39:49.897Z
Mapping millions of sequences to reference sequences is a common problem in bioinformatics, and has already been solved efficiently many times over. What is the problem that you are trying to solve? Are you creating your own code for educational purposes?
gringer
CC BY-SA 4.0
2022-08-07T17:46:14.453Z
Hi @gringer. Thanks, the 'mapping' in the question I think (could be wrong) is the `map()` Python function. I edited the title - so lock of clarity is perhaps my error. The OP is using `map()` to try and speed up the function. (I personally didn't know whether `yield` would work better here). 'Recodifying' via `map()` is known practice. Lets see what the OP says. I'm happy to reedit accordingly.
M__
CC BY-SA 4.0
2022-08-07T23:05:53.730Z
The code is mapping amino acid sequences to reference sequences, but the ultimate goal of the problem is not clear from the question. Either way, more clarification is needed. If the goal is to make a hash lookup more efficient, then the question would be better answered in another StackExchange community. If the goal is to find some way to efficiently map amino acids to sequences, then plenty of other better existing solutions exist. Educational / learning is one reason to get answers about code optimisation.
gringer
<python><algorithms><benchmarking>
Refactoring a script mapping ID to sequence
# Refactoring a script mapping ID to sequence

(Asked by: username_0 on 2022-08-06)

I wrote simple Python script which has an excessive run time, notably when using large data sets. For example a data set of 1 000 000 sequences.
I am seeking assistance refactoring the code and would be very grateful for any assistance.
from Bio import SeqIO
import re

def move_item(old, new):
    k, v = old.popitem()
    new[k] = v

def map_ids_to_aa_sequence_list(id_):
    id_ = id_.strip()
    for s in aa:
        if re.search(id_,s, re.I):
            return s

file = open(r'list_of_ids', 'r').readlines()
names = [for name.strip() in file if len(name) > 1]
        
handle = open('sequences.fasta','r')
sequences = list(SeqIO.parse(handle, 'fasta'))

seq_dict = {}

for s in sequences:
    seq_dict[str(s.description)] = str(s.seq)
    
aa = []

for k, v in seq_dict.items():
    sequence = '>' + k +'\n' + v + '\n'
    aa.append(sequence)
    
separated_aminoacid_sequences_set = set(map(map_ids_to_aa_sequence_list, names))

separated_aminoacid_sequences_list = list(separated_aminoacid_sequences_set)
    
with open('separated_sequences.fasta','w') as fs:
    for s in separated_aminoacid_sequences_list:
        fs.write(str(s))

**Tags:** python, algorithms, benchmarking
> **Comments:**
> **username_1**: On a brief glance the problem could be `re` but I'd need to look. Could you kindly look at this code here https://bioinformatics.stackexchange.com/questions/19506/improving-list-speed-glycosylation-example  Are there any parallels?  
> **username_2**: Mapping millions of sequences to reference sequences is a common problem in bioinformatics, and has already been solved efficiently many times over. What is the problem that you are trying to solve? Are you creating your own code for educational purposes?  
> **username_1**: Hi @gringer. Thanks, the 'mapping' in the question I think (could be wrong) is the `map()` Python function. I edited the title - so lock of clarity is perhaps my error. The OP is using `map()` to try and speed up the function. (I personally didn't know whether `yield` would work better here). 'Recodifying' via `map()` is known practice. Lets see what the OP says. I'm happy to reedit accordingly.  
> **username_2**: The code is mapping amino acid sequences to reference sequences, but the ultimate goal of the problem is not clear from the question. Either way, more clarification is needed. If the goal is to make a hash lookup more efficient, then the question would be better answered in another StackExchange community. If the goal is to find some way to efficiently map amino acids to sequences, then plenty of other better existing solutions exist. Educational / learning is one reason to get answers about code optimisation.  

---

## 1 answers

### Answer by username_1 on 2022-08-06. Score: 2 | Accepted answer

This is beyond a single answer. 'Refactoring' (code efficiency, high performance code) is a huge topic. Like its own project.
The first thing I would personally absolutely target first is re here:
    if re.search(id_,s, re.I):

Python is not Perl and the re is EXTREMELY slow and in particular using search is especially slow, most notably when uncompiled (see notes and point 3 below).
How to avoid re.search
Could you revise this to a direct C style operation? This is the best solution if you get it working.
For example,

Using ==:

e.g.
if str(id_) == str(s):

if the case is important (like they're different) then...
if str(id_).lower() == str(s).lower():

You could dump the str() it ain't doing much its just bug insurance.
The other approach if you want to stick with re is to switch to match, that will be quicker.

Using match

e.g.
if re.match(id_, s, re.I):

match anchors  the search at the first character, i.e. id_[0], search goes everywhere. match is faster.

Sticking with search but combining with compile:

e.g.
id_ = id_.strip()
x = re.compile (id_) 
if x.search(s, re.I): # don't use re. x replaces it

or
id_ = id_.strip()
x = re.compile (id_)
if x.match(s, re.I):

I have not put these in order of preference. My preference is for 1. but it will risk throwing logical errors (missing stuff).
If you are using re you must compile the search term. Its a pain but whilst compiled re is slow, uncompiled re is really, really slow. Disclaimer, I'm from a Perl background and from this perspective re is the pits because previously relied on the efficient regex (Perl).
The == and match could fall over if > is in one variable and not in the other, but the search will return the correct 'match' (i.e. not the output of 'match') at the cost of algorithm efficiency.
This isn't an exhaustive analysis of re but thats where I would start.
Please note this is a starting point from prior experience to pick a 'big issues'. It is useful for future, particularly for questions, to place time checks across the code to identify the bottlenecks.

Notes

The thing we don't know is what the fasta headers look like.
For every query in search there are $1.10^6$ fasta headers. So if you've $1.10^5$ queries thats $1.10^{11}$ searches for search thats why it can be taking up a lot of time.
Rewriting the code using Biopython objects throughout is good practice. SeqIO in particular are cool. AlignIO are algorithmically expensive however.


---

