Tool to remove a PCR contamination in NGS data

Question

BACKGROUND

I work on NGS data (illumina paired ends reads) coming from a full extract of RNA (metagenomic). We are interested in the viral fraction of this extract.

I observed a contamination with a PCR amplicon. This contamination has been identified as coming from a PCR result that we used to make near the library preparation. So I know exactly the primers used in the PCR.

Here you can see a tablet visualization illustrating the contamination (2 contaminations). I obtained that image with a mapping on a known reference of my virus.

PROBLEM

The problem I have is that these PCR products are exactly the type of sequence I want from my extract (identifying part of the virus). So I have some "false positive" and some misidentification of viruses present in the extract.

I have thought about removing all reads containing primer sequences but it leads to a gap in my data and I can't identify the viruses in my extract anymore. I am thinking that all of the read pairs that begin before or end after my primers are "true-positives" and I would really like to keep them (even having the primer sequence in them).

Do you know any tool that can help me in that situation? i.e. removing reads containing a sequence but only if it begins with it.

If the tool can also handle the degeneration of my primer it would be perfect.
An example of primer I used is

CNTGGGAGGGCGATCGCAA

Complementary information (asked in comments)
I have around 30 million read pairs (60 millions reads in total). Contamination reads are between 200.000 to a 1 millions reads (depending on sample). "true" reads (not from the contamination) represent around 0.5% of the contamination.

Answer 1

This bash script will get you all of the reads that do not start with your primer sequence. It also handles degeneracy. It requires the tool bioawk.

It works by looking at the first bases in the read and only prints out the headers that do not have a match with the degenerate primer. Then, comm makes sure that the read was retained in both the R1 and R2 filtering. This allows reads to be filtered whether or not the primers start at the beginning of the R1 or R2 read.

seqtk subseq is where you recover the reads that don't start with the primer.

Be warned though, any off-by-one or sequencing errors will allow some of the sequences with primers at the beginning. Also, this script requires that the headers for the R1 and R2 reads are the same and aren't appended with anything like "\1", "\2", et cetera.

#!/bin/bash
get_primer_free_starts () {
    bioawk -cfastx 'BEGIN{seqs[0] = "CATGGGAGGGCGATCGCAA"; seqs[1] = "CCTGGGAGGGCGATCGCAA"; \
        seqs[2] = "CGTGGGAGGGCGATCGCAA"; seqs[3] = "CTTGGGAGGGCGATCGCAA"} \
        {thisseq = substr($seq, 1, length(seqs[0])); printme=1\
        for (i=0; i<=3; i++) { \
            if (thisseq != seqs[i]){ \
                printme=0 \
            }
        }; if (printme==1){ \
            print($seq)
        }}' "${INPUT} | sort > "${OUTPUT}
}

INPUT="R1.fastq.gz"
OUTPUT="R1_OK.txt"
get_primer_free_starts

INPUT="R2.fastq.gz"
OUTPUT="R2_OK.txt"
get_primer_free_starts

comm -12 R1_OK.txt R2_OK.txt | sort | uniq > passing_reads.txt

seqtk subseq R1.fastq.gz passing_reads.txt > R1.filtered.fq
seqtk subseq R2.fastq.gz passing_reads.txt > R2.filtered.fq