text: user16945: 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 :)
  
  user298: 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.
  
original_text: 16945: <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>
   298: <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>
  