3
$\begingroup$

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. twin towers and world trade center

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

enter image description here

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.

$\endgroup$
5
  • $\begingroup$ Welcome to the site Untitpoi. Could you please clarify how many read do you remove when you "remove all reads containing primer sequences"? Also how many reads you have in total? Could you quantify the contamination (i.e 50% of your reads are contaminated)? Also how did you try to remove these sequence did you use a bash script or something like that? $\endgroup$
    – llrs
    Oct 10, 2018 at 15:34
  • $\begingroup$ Have you searched this site yet? There are other similar questions such as bioinformatics.stackexchange.com/questions/5163/…. Also what tools have you tried? $\endgroup$
    – Bioathlete
    Oct 10, 2018 at 16:58
  • $\begingroup$ Can you post an example of what you reasonably know to be a PCR contaminant and a reasonably similar read that you believe not to be but gets removed using the methods you've tried? That would help us fully grasp exactly what you're dealing with. $\endgroup$
    – Devon Ryan
    Oct 10, 2018 at 17:15
  • $\begingroup$ @DevonRyan I added some graph to explain my problem, I hope this is clearer now. I didn't try anything yet but I will try a cut-adapt method today and the bioawk method proposed by conchoecia. $\endgroup$
    – Untitpoi
    Oct 11, 2018 at 7:49
  • $\begingroup$ I just saw your post with the exact sequence of the forward and reverse primers. Is it true that you expect the beginning of one read to start with the forward primer, and the beginning of its paired read to start with the reverse primer? If so, I should change my answer to reflect that and to better filter your data. $\endgroup$
    – conchoecia
    Oct 12, 2018 at 21:18

1 Answer 1

2
$\begingroup$

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.