Primer trimming-fasta files

Question

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:

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:

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

Answer 1

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.