<p>I'm working on a Python script that shuffles DNA sequences while preserving their protein-coding potential. My goal is to adjust the codon usage in such a way that the ACGT (adenine, cytosine, guanine, thymine) content of the sequence remains as close as possible to the original. I'm using Python's <code>multiprocessing</code> module to speed up the shuffling process by performing multiple attempts in parallel.</p>
<h3>What the Code Does:</h3>
<p>The code shuffles DNA sequences while trying to preserve the ACGT content and maintain the protein-coding potential of the sequence. It uses multiprocessing to perform multiple shuffles in parallel and selects the best result based on how close the ACGT content of the shuffled sequence is to the original sequence.</p>
<h3>What I Want to Improve:</h3>
<p>Speed: Although the code uses multiprocessing, I'm looking for ways to further optimize the performance, especially when dealing with large sequences.</p>
<p>Accuracy: The goal is to get the ACGT difference as close to zero as possible. Currently, some sequences don't achieve this, and I'd like to improve the likelihood of achieving an exact match.</p>
<h3>My Questions:</h3>
<ol>
<li>Are there any obvious optimizations I could apply to this code to make it faster or more efficient?</li>
<li>How can I improve the accuracy to ensure that more sequences have an exact ACGT match after shuffling?</li>
<li>Is there a better approach to achieving these goals, either algorithmically or through different Python libraries?</li>
</ol>
<p>Here's a simplified version of the core functions I'm using:</p>
<pre class="lang-py prettyprint-override"><code>from collections import Counter
from Bio.Data import CodonTable
from multiprocessing import Pool, cpu_count
import random
import logging

def adjust_acgt_content(shuffled_codons, original_counter, codon_table, remaining_nucleotides):
    # Adjust the shuffled sequence to match the ACGT content of the original sequence.
    for base in 'ACGT':
        attempt_count = 0
        shuffled_counter = Counter(''.join(shuffled_codons))
        while original_counter[base] != shuffled_counter[base]:
            attempt_count += 1
            if attempt_count &gt; 100:  # Prevent infinite loops
                break
            difference = original_counter[base] - shuffled_counter[base]
            if difference &gt; 0:
                for i in range(len(shuffled_codons)):
                    for j in range(3):
                        if shuffled_codons[i][j] != base:
                            current_codon = shuffled_codons[i]
                            new_codon = list(current_codon)
                            new_codon[j] = base
                            new_codon_str = ''.join(new_codon)
                            if codon_table.forward_table.get(new_codon_str) == codon_table.forward_table.get(current_codon):
                                shuffled_codons[i] = new_codon_str
                                shuffled_counter = Counter(''.join(shuffled_codons))
                                break
                    if original_counter[base] == shuffled_counter[base]:
                        break
            elif difference &lt; 0:
                for i in range(len(shuffled_codons)):
                    for j in range(3):
                        if shuffled_codons[i][j] == base:
                            current_codon = shuffled_codons[i]
                            new_codon = list(current_codon)
                            new_codon[j] = random.choice('ACGT'.replace(base, ''))
                            new_codon_str = ''.join(new_codon)
                            if codon_table.forward_table.get(new_codon_str) == codon_table.forward_table.get(current_codon):
                                shuffled_codons[i] = new_codon_str
                                shuffled_counter = Counter(''.join(shuffled_codons))
                                break
                    if original_counter[base] == shuffled_counter[base]:
                        break
    return ''.join(shuffled_codons) + remaining_nucleotides

def shuffle_once(codons, sequence, original_counter, codon_table, codons_by_aa, remaining_nucleotides):
    # Shuffle codons while preserving the ACGT content as much as possible.
    shuffled_codons = []
    shuffled_counter = Counter('')
    for codon in codons:
        aa = codon_table.forward_table.get(codon)
        if aa in codons_by_aa:
            possible_codons = [c for c in codons_by_aa[aa] if c != codon]
            if possible_codons:
                min_diff = min(sum(abs(shuffled_counter[base] + c.count(base) - original_counter[base]) for base in 'ACGT') for c in possible_codons)
                best_codons = [c for c in possible_codons if sum(abs(shuffled_counter[base] + c.count(base) - original_counter[base]) for base in 'ACGT') == min_diff]
                best_codon = random.choice(best_codons)
                shuffled_codons.append(best_codon)
                shuffled_counter.update(best_codon)
            else:
                shuffled_codons.append(codon)
                shuffled_counter.update(codon)
        else:
            shuffled_codons.append(codon)
            shuffled_counter.update(codon)
    return adjust_acgt_content(shuffled_codons, original_counter, codon_table, remaining_nucleotides)

def worker(codons, sequence, original_counter, codon_table, codons_by_aa, remaining_nucleotides):
    shuffled_sequence = shuffle_once(codons, sequence, original_counter, codon_table, codons_by_aa, remaining_nucleotides)
    shuffled_codons = [shuffled_sequence[i:i+3] for i in range(0, len(shuffled_sequence) - len(shuffled_sequence) % 3, 3)]
    shuffled_counter = Counter(shuffled_sequence[:len(shuffled_codons) * 3])
    acgt_diff = sum(abs(original_counter[base] - shuffled_counter[base]) for base in 'ACGT')
    num_different_codons = sum(1 for o, s in zip(codons, shuffled_codons) if o != s)
    return shuffled_sequence, acgt_diff, num_different_codons

def shuffle_coding(sequence, transcript_id, num_shuffles=10):
    # Load the codon table for standard DNA coding
    codon_table = CodonTable.unambiguous_dna_by_id[1]
    full_codons_length = len(sequence) - len(sequence) % 3
    codons = [sequence[i:i+3] for i in range(0, full_codons_length, 3)]
    remaining_nucleotides = sequence[full_codons_length:]
    codons_by_aa = {}
    for codon, aa in codon_table.forward_table.items():
        codons_by_aa.setdefault(aa, []).append(codon)
    original_counter = Counter(sequence[:full_codons_length])
    with Pool(processes=cpu_count()) as pool:
        results = pool.starmap(worker, [(codons, sequence, original_counter, codon_table, codons_by_aa, remaining_nucleotides) for _ in range(num_shuffles)])
    best_sequence, min_acgt_diff, max_diff_codons = min(results, key=lambda x: (x[1], -x[2]))
    procentage_acgt_diff = (min_acgt_diff / len(sequence)) * 100
    if procentage_acgt_diff &gt; 1.0:
        logging.warning(f&quot;ACGT difference after shuffling for {transcript_id}: {procentage_acgt_diff:.2f}%&quot;)
    return best_sequence
</code></pre>
<p>Thank you in advance for any suggestions or insights!</p>

CC BY-SA 4.0
2024-08-23T12:04:26.437Z
Have you identified which of these functions is the slowest? And how does the script call them? Is the main body of the script simply calling each function once for each input sequence or multiple times?
terdon
CC BY-SA 4.0
2024-08-23T12:26:53.757Z
What exactly do you mean by "maintain the protein-coding potential of the sequence"? Do you mean that the shuffled sequence still needs to code for protein (so no in frame stop codons), or that it should code for the _same_ protein? Or is it more about protein coding models and codon usage so it isn't enough to code but it needs to use the specific species' codon frequency? And how long are your sequences?
terdon
CC BY-SA 4.0
2024-08-23T17:23:53.823Z
Thanks, but please add all of this to the question, [edit] and add it (no need for a special "edit" section, just insert into the flow of the question) so that everyone can see it. Comments are hard to read, easy to miss and can be deleted without warning.
terdon
<python><biopython><codon><benchmarking>
Improving DNA Sequence Shuffling Code for Speed and Accuracy
# Improving DNA Sequence Shuffling Code for Speed and Accuracy

(Asked by: username_0 on 2024-08-23)

I'm working on a Python script that shuffles DNA sequences while preserving their protein-coding potential. My goal is to adjust the codon usage in such a way that the ACGT (adenine, cytosine, guanine, thymine) content of the sequence remains as close as possible to the original. I'm using Python's multiprocessing module to speed up the shuffling process by performing multiple attempts in parallel.
What the Code Does:
The code shuffles DNA sequences while trying to preserve the ACGT content and maintain the protein-coding potential of the sequence. It uses multiprocessing to perform multiple shuffles in parallel and selects the best result based on how close the ACGT content of the shuffled sequence is to the original sequence.
What I Want to Improve:
Speed: Although the code uses multiprocessing, I'm looking for ways to further optimize the performance, especially when dealing with large sequences.
Accuracy: The goal is to get the ACGT difference as close to zero as possible. Currently, some sequences don't achieve this, and I'd like to improve the likelihood of achieving an exact match.
My Questions:

Are there any obvious optimizations I could apply to this code to make it faster or more efficient?
How can I improve the accuracy to ensure that more sequences have an exact ACGT match after shuffling?
Is there a better approach to achieving these goals, either algorithmically or through different Python libraries?

Here's a simplified version of the core functions I'm using:
from collections import Counter
from Bio.Data import CodonTable
from multiprocessing import Pool, cpu_count
import random
import logging

def adjust_acgt_content(shuffled_codons, original_counter, codon_table, remaining_nucleotides):
    # Adjust the shuffled sequence to match the ACGT content of the original sequence.
    for base in 'ACGT':
        attempt_count = 0
        shuffled_counter = Counter(''.join(shuffled_codons))
        while original_counter[base] != shuffled_counter[base]:
            attempt_count += 1
            if attempt_count > 100:  # Prevent infinite loops
                break
            difference = original_counter[base] - shuffled_counter[base]
            if difference > 0:
                for i in range(len(shuffled_codons)):
                    for j in range(3):
                        if shuffled_codons[i][j] != base:
                            current_codon = shuffled_codons[i]
                            new_codon = list(current_codon)
                            new_codon[j] = base
                            new_codon_str = ''.join(new_codon)
                            if codon_table.forward_table.get(new_codon_str) == codon_table.forward_table.get(current_codon):
                                shuffled_codons[i] = new_codon_str
                                shuffled_counter = Counter(''.join(shuffled_codons))
                                break
                    if original_counter[base] == shuffled_counter[base]:
                        break
            elif difference < 0:
                for i in range(len(shuffled_codons)):
                    for j in range(3):
                        if shuffled_codons[i][j] == base:
                            current_codon = shuffled_codons[i]
                            new_codon = list(current_codon)
                            new_codon[j] = random.choice('ACGT'.replace(base, ''))
                            new_codon_str = ''.join(new_codon)
                            if codon_table.forward_table.get(new_codon_str) == codon_table.forward_table.get(current_codon):
                                shuffled_codons[i] = new_codon_str
                                shuffled_counter = Counter(''.join(shuffled_codons))
                                break
                    if original_counter[base] == shuffled_counter[base]:
                        break
    return ''.join(shuffled_codons) + remaining_nucleotides

def shuffle_once(codons, sequence, original_counter, codon_table, codons_by_aa, remaining_nucleotides):
    # Shuffle codons while preserving the ACGT content as much as possible.
    shuffled_codons = []
    shuffled_counter = Counter('')
    for codon in codons:
        aa = codon_table.forward_table.get(codon)
        if aa in codons_by_aa:
            possible_codons = [c for c in codons_by_aa[aa] if c != codon]
            if possible_codons:
                min_diff = min(sum(abs(shuffled_counter[base] + c.count(base) - original_counter[base]) for base in 'ACGT') for c in possible_codons)
                best_codons = [c for c in possible_codons if sum(abs(shuffled_counter[base] + c.count(base) - original_counter[base]) for base in 'ACGT') == min_diff]
                best_codon = random.choice(best_codons)
                shuffled_codons.append(best_codon)
                shuffled_counter.update(best_codon)
            else:
                shuffled_codons.append(codon)
                shuffled_counter.update(codon)
        else:
            shuffled_codons.append(codon)
            shuffled_counter.update(codon)
    return adjust_acgt_content(shuffled_codons, original_counter, codon_table, remaining_nucleotides)

def worker(codons, sequence, original_counter, codon_table, codons_by_aa, remaining_nucleotides):
    shuffled_sequence = shuffle_once(codons, sequence, original_counter, codon_table, codons_by_aa, remaining_nucleotides)
    shuffled_codons = [shuffled_sequence[i:i+3] for i in range(0, len(shuffled_sequence) - len(shuffled_sequence) % 3, 3)]
    shuffled_counter = Counter(shuffled_sequence[:len(shuffled_codons) * 3])
    acgt_diff = sum(abs(original_counter[base] - shuffled_counter[base]) for base in 'ACGT')
    num_different_codons = sum(1 for o, s in zip(codons, shuffled_codons) if o != s)
    return shuffled_sequence, acgt_diff, num_different_codons

def shuffle_coding(sequence, transcript_id, num_shuffles=10):
    # Load the codon table for standard DNA coding
    codon_table = CodonTable.unambiguous_dna_by_id[1]
    full_codons_length = len(sequence) - len(sequence) % 3
    codons = [sequence[i:i+3] for i in range(0, full_codons_length, 3)]
    remaining_nucleotides = sequence[full_codons_length:]
    codons_by_aa = {}
    for codon, aa in codon_table.forward_table.items():
        codons_by_aa.setdefault(aa, []).append(codon)
    original_counter = Counter(sequence[:full_codons_length])
    with Pool(processes=cpu_count()) as pool:
        results = pool.starmap(worker, [(codons, sequence, original_counter, codon_table, codons_by_aa, remaining_nucleotides) for _ in range(num_shuffles)])
    best_sequence, min_acgt_diff, max_diff_codons = min(results, key=lambda x: (x[1], -x[2]))
    procentage_acgt_diff = (min_acgt_diff / len(sequence)) * 100
    if procentage_acgt_diff > 1.0:
        logging.warning(f"ACGT difference after shuffling for {transcript_id}: {procentage_acgt_diff:.2f}%")
    return best_sequence

Thank you in advance for any suggestions or insights!

**Tags:** python, biopython, codon, benchmarking
> **Comments:**
> **username_1**: Have you identified which of these functions is the slowest? And how does the script call them? Is the main body of the script simply calling each function once for each input sequence or multiple times?  
> **username_1**: What exactly do you mean by "maintain the protein-coding potential of the sequence"? Do you mean that the shuffled sequence still needs to code for protein (so no in frame stop codons), or that it should code for the _same_ protein? Or is it more about protein coding models and codon usage so it isn't enough to code but it needs to use the specific species' codon frequency? And how long are your sequences?  
> **username_1**: Thanks, but please add all of this to the question, [edit] and add it (no need for a special "edit" section, just insert into the flow of the question) so that everyone can see it. Comments are hard to read, easy to miss and can be deleted without warning.  

---

## 0 answers

