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:
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.
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 :)