Pull out unknown primers from fastq file?

[Brian Bushnell]

Hi Jen,

My tools are all still unpublished, so please just cite my name and website. Thanks!

-Brian

[JenBarb]

Hi Brian,
How should I cite your tool in a manuscript in prep that I am doing? Do you have a reference or should I use your website?
Thanks,
Jen

[Brian Bushnell]

If you have not solved this problem yet, there's another option, using BBTools:

reformat.sh in=reads.fq out=trimmed.fq ftr=19
This will trim all but the first 20 bases (all bases after position 19, zero-based).

kmercountexact.sh in=trimmed.fq out=counts.txt fastadump=f mincount=10 k=20 rcomp=f
This will generate a file containing the counts of all 20-mers that occurred at least 10 times, in a 2-column format that is easy to sort in Excel. For example:

ACCGTTACCGTTACCGTTAC	100
AAATTTTTTTCCCCCCCCCC	85

...etc. If the primers are 20bp long, they should be pretty obvious.

[maubp]

e.g. FASTQC reports overrepresented sequences which ought to spot your primers:

http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/9%20Overrepresented%20Sequences.html

[JenBarb]

maubp,
I would love to find a tool that will spot and report possible primers. Can you be more specific?

I tried to sort and tally up the sequences and i am not finding them this way. Which tool are you referring to?

Thanks a bunch.
Jen

[maubp]

IIRC, some FASTQ quality control pipelines will spot and report possible primer sequences.

[cmbetts]

Something like that should be pretty easy to do with any scripting language with fastq parsing libraries (or heck maybe manually inspecting the fastq since there's only 6 primers if you're not keen on programming).

Example R psuedocode (only because that's what I'm comfortable with)
library("ShortRead"); #load the library for fastq manipulation
fq_data <- readFastq("reads.fastq.gz"); #read in fastq data
base_info <- sread(fq_data); #get just the base calls
first20 <- substring(base_info, 1, 20); #get the first 20bp of each read

then you could do something like
table(first20) to see the frequency of different 20bp sequences
alphabetFrequency(first20) to get consensus sequences

[JenBarb]

Hello,
I have fastq files from 16S sequencing data. The reads in these files have 6 different primers in them and I am wondering if anyone knows of a method where I can pull out the beginning of the reads, maybe 20-30 bps and then look for consensus sequences within these? This will then give me the primer sequence information for the 6 primers which is proprietary information from the company where the kits are made.

Can anyone think of a way to do this?
Thanks!