<p>The IUPAC ambiguity codes can be thought of as regular expression character classes. <code>R</code> matches any purine, so <code>[AG]</code>, <code>Y</code> matches any pyrimidine, so <code>[CT]</code>, and <code>N</code> matches anything at all, so <code>.</code> or <code>[ACTG]</code>. So those are easy to handle: instead of <code> CGAGAAGACCCTRTGRAGCT</code> you can use the regular expression <code>CGAGAAGACCCT[AG]TG[AG]AGCT</code>, and instead of <code>GTTGGGGYGACCNYGG</code> you can use <code>GTTGGGG[CT]GACC.[CT]GG</code>.</p>
<p>Allowing for random mismatches is trickier though and since you're open to using existing tools, I strongly urge you to use something like <code>cutadapt</code> instead of trying to reinvent the wheel. So, install <code>cutadapt</code> on your system (see <a href="https://cutadapt.readthedocs.io/en/stable/installation.html" rel="nofollow noreferrer">here</a> for instructions) and then run:</p>
<pre><code>cutadapt -a &quot;CGAGAAGACCCTRTGRAGCT...GTTGGGGYGACCNYGG;max_error_rate=0.15;&quot; --discard-untrimmed -o output.fa file.fa  
</code></pre>
<p>That assumes your input sequences are in fasta format in the file <code>file.fa</code> and  will create the output file <code>output.fa</code>:</p>
<pre><code>$ cat output.fa 
&gt;seq1
TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG
&gt;seq3
TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG
</code></pre>
<p>The options used are <code>-a $primer1..$primer2</code> which tells <code>cutudapt</code> these primers are linked adapters (<a href="https://cutadapt.readthedocs.io/en/stable/guide.html#linked-adapters" rel="nofollow noreferrer">https://cutadapt.readthedocs.io/en/stable/guide.html#linked-adapters</a>) and then we also add <code>max_error_rate=0.15;</code> to allow up to 85% identity matches. Cutadapt already knows how to handle IUPAC ambiguity codes so nothing special is needed for those.</p>
<p>I am not sure how this will handle cases where only one of the two sequences is present, you'll need to test for that. I expect that they will be discarded since we're specifying linked primers but I would double check to be sure. In any case, <code>cutadapt</code> is a very powerful tool, so I'm sure there will be a way to handle that in the way you want if you check the documentation.</p>

CC BY-SA 4.0
2023-02-12T09:44:34.680Z
Thank you a lot! Sequences with only 1 primer (and the 2nd is absent or more mismatches than the given threshold) aren't discarded.  cutadapt cuts the one primer which is present and let the sequence in the dataset which i dont want to.  But generally it works.  Σε ευχαριστώ θερμά! @terdon
Damianos
CC BY-SA 4.0
2023-02-11T18:28:41.077Z
2023-02-11T18:28:41.077Z
terdon
<p>My goal:
I want to trim off the primers (Forward : CGAGAAGACCCTRTGRAGCT, Reverse : GTTGGGGYGACCNYGG) from a fasta file with a lot of dna sequences allowing for some (e.g. 3) mismatches (identity).</p>
<p>Input sequences
Seq1:</p>
<blockquote>
<p><strong>CGAGAAGACCCTATGGAGCT</strong>TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG<strong>GTTGGGGCGACCGCGG</strong>AACGT</p>
</blockquote>
<p>Seq2:</p>
<blockquote>
<p>ATGGCCCATCCCTCACAGCTAGGATTTCAAGATGCAGCTTCCCCAGTTATAGAAGAACTTCTCCACTTTCACGACCATGCCCTAATAATCGTTTTTTTAATTAGTACACTAGTA</p>
</blockquote>
<p>Seq3:</p>
<blockquote>
<p>CATAAGA<strong>CGAAAAGACCCTATGGAGCT</strong>TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG<strong>GTTCGGGCGACCGCGG</strong></p>
</blockquote>
<p>Output sequences</p>
<p>Seq1:</p>
<blockquote>
<p>TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG</p>
</blockquote>
<p>Seq3:</p>
<blockquote>
<p>TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG</p>
</blockquote>
<p>Seq2 is absent because it has no match for the primers so it should be discarded.</p>
<p>In seq 3 there is a mismatch in Forward (CGA<strong>A</strong>AA...) and in Reverse (GTT<strong>C</strong>GG...).  I don't need the primers to have 100% identity.  85-90% is OK.</p>
<p>Issues:</p>
<ol>
<li><p>primer sequences can be anywhere in each sequence, not just in the beginning.</p>
</li>
<li><p>The sequence between the primers is not constant.  it is approximately 80 bases but not exactly 80.</p>
</li>
<li><p>I don't know how to put the primers' identity parameter.</p>
</li>
<li><p>How can I search for R or N nucleotides (IUPAC)? Because I don't know how, I wrote the primers with specific nucleotides and not R or N.</p>
</li>
</ol>
<p>This is the code I have so far, but it has the problems mentioned above:</p>
<pre class="lang-py prettyprint-override"><code>from Bio import SeqIO

def trim_adaptors(records, adaptorF,adaptorR):

    len_adaptorF = len(adaptorF)  
    len_adaptorR = len(adaptorR)

    for record in records:
        indexF = record.seq.find(adaptorF)
        indexR = record.seq.find(adaptorR)
        if indexR == -1  or indexF == -1:
            # adaptor not found, so won't trim
            print(record.id, &quot;no primer/s match&quot;)
            continue
        else:
            # trim off the adaptor
            yield record[indexF + len_adaptorF:indexR + 1]

original_reads = SeqIO.parse(&quot;trimming_testfas.fas&quot;, &quot;fasta&quot;)

trimmed_reads = trim_adaptors(original_reads,&quot;CGAGAAGACCCTATGGAGCT&quot; ,&quot;GTTGGGGCGACCGCGG&quot;)

count = SeqIO.write(trimmed_reads, &quot;trimmed.fasta&quot;, &quot;fasta&quot;)

print(&quot;Saved %i reads&quot; % count)

</code></pre>
<p>This doesn't have to be a script in Python, I am also open to using existing tools.</p>
<p>Thank you a lot in advance, hope you got me :)</p>

<python><trimming><primer><multi-fasta>
Primer trimming-fasta files
# Primer trimming-fasta files

(Asked by: username_0 on 2023-02-11)

My goal:
I want to trim off the primers (Forward : CGAGAAGACCCTRTGRAGCT, Reverse : GTTGGGGYGACCNYGG) from a fasta file with a lot of dna sequences allowing for some (e.g. 3) mismatches (identity).
Input sequences
Seq1:

CGAGAAGACCCTATGGAGCTTAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTGGTTGGGGCGACCGCGGAACGT

Seq2:

ATGGCCCATCCCTCACAGCTAGGATTTCAAGATGCAGCTTCCCCAGTTATAGAAGAACTTCTCCACTTTCACGACCATGCCCTAATAATCGTTTTTTTAATTAGTACACTAGTA

Seq3:

CATAAGACGAAAAGACCCTATGGAGCTTTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTGGTTCGGGCGACCGCGG

Output sequences
Seq1:

TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG

Seq3:

TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG

Seq2 is absent because it has no match for the primers so it should be discarded.
In seq 3 there is a mismatch in Forward (CGAAAA...) and in Reverse (GTTCGG...).  I don't need the primers to have 100% identity.  85-90% is OK.
Issues:

primer sequences can be anywhere in each sequence, not just in the beginning.

The sequence between the primers is not constant.  it is approximately 80 bases but not exactly 80.

I don't know how to put the primers' identity parameter.

How can I search for R or N nucleotides (IUPAC)? Because I don't know how, I wrote the primers with specific nucleotides and not R or N.


This is the code I have so far, but it has the problems mentioned above:
from Bio import SeqIO

def trim_adaptors(records, adaptorF,adaptorR):

    len_adaptorF = len(adaptorF)  
    len_adaptorR = len(adaptorR)

    for record in records:
        indexF = record.seq.find(adaptorF)
        indexR = record.seq.find(adaptorR)
        if indexR == -1  or indexF == -1:
            # adaptor not found, so won't trim
            print(record.id, "no primer/s match")
            continue
        else:
            # trim off the adaptor
            yield record[indexF + len_adaptorF:indexR + 1]

original_reads = SeqIO.parse("trimming_testfas.fas", "fasta")

trimmed_reads = trim_adaptors(original_reads,"CGAGAAGACCCTATGGAGCT" ,"GTTGGGGCGACCGCGG")

count = SeqIO.write(trimmed_reads, "trimmed.fasta", "fasta")

print("Saved %i reads" % count)


This doesn't have to be a script in Python, I am also open to using existing tools.
Thank you a lot in advance, hope you got me :)

**Tags:** python, trimming, primer, multi-fasta

---

## 1 answers

### Answer by username_1 on 2023-02-11. Score: 2 | Accepted answer

The IUPAC ambiguity codes can be thought of as regular expression character classes. R matches any purine, so [AG], Y matches any pyrimidine, so [CT], and N matches anything at all, so . or [ACTG]. So those are easy to handle: instead of  CGAGAAGACCCTRTGRAGCT you can use the regular expression CGAGAAGACCCT[AG]TG[AG]AGCT, and instead of GTTGGGGYGACCNYGG you can use GTTGGGG[CT]GACC.[CT]GG.
Allowing for random mismatches is trickier though and since you're open to using existing tools, I strongly urge you to use something like cutadapt instead of trying to reinvent the wheel. So, install cutadapt on your system (see here for instructions) and then run:
cutadapt -a "CGAGAAGACCCTRTGRAGCT...GTTGGGGYGACCNYGG;max_error_rate=0.15;" --discard-untrimmed -o output.fa file.fa  

That assumes your input sequences are in fasta format in the file file.fa and  will create the output file output.fa:
$ cat output.fa 
>seq1
TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG
>seq3
TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG

The options used are -a $primer1..$primer2 which tells cutudapt these primers are linked adapters (https://cutadapt.readthedocs.io/en/stable/guide.html#linked-adapters) and then we also add max_error_rate=0.15; to allow up to 85% identity matches. Cutadapt already knows how to handle IUPAC ambiguity codes so nothing special is needed for those.
I am not sure how this will handle cases where only one of the two sequences is present, you'll need to test for that. I expect that they will be discarded since we're specifying linked primers but I would double check to be sure. In any case, cutadapt is a very powerful tool, so I'm sure there will be a way to handle that in the way you want if you check the documentation.

> **Comments:**
> **username_0**: Thank you a lot! Sequences with only 1 primer (and the 2nd is absent or more mismatches than the given threshold) aren't discarded.  cutadapt cuts the one primer which is present and let the sequence in the dataset which i dont want to.  But generally it works.  Σε ευχαριστώ θερμά! @terdon  

---

