Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-02-12T09:44:34.680Z
Answers.Text: 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
Answers.UserDisplayName: Damianos
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-02-11T18:28:41.077Z
Answers.LastActivityDate: 2023-02-11T18:28:41.077Z
Answers.UserDisplayName: terdon
Body: <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>
  
Tags: <python><trimming><primer><multi-fasta>
Title: Primer trimming-fasta files
ThreadText: # 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:
  
  > **CGAGAAGACCCTATGGAGCT**TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG**GTTGGGGCGACCGCGG**AACGT
  
  Seq2:
  
  > ATGGCCCATCCCTCACAGCTAGGATTTCAAGATGCAGCTTCCCCAGTTATAGAAGAACTTCTCCACTTTCACGACCATGCCCTAATAATCGTTTTTTTAATTAGTACACTAGTA
  
  Seq3:
  
  > CATAAGA**CGAAAAGACCCTATGGAGCT**TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG**GTTCGGGCGACCGCGG**
  
  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 (CGA**A**AA...) and in Reverse (GTT**C**GG...). I don't need the primers to have 100% identity. 85-90% is OK.
  
  Issues:
  
  1. primer sequences can be anywhere in each sequence, not just in the beginning.
  2. The sequence between the primers is not constant. it is approximately 80 bases but not exactly 80.
  3. I don't know how to put the primers' identity parameter.
  4. 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:
  
  ```py
  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](https://cutadapt.readthedocs.io/en/stable/installation.html) 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.  Σε ευχαριστώ θερμά! @username_1  
  
  ---
  
  