Searching for and extracting sequences within a fastQ file

[rtleenay]

Hello,

So I'm trying to extract information from a nucleotide library I have sequenced on an Illumina MiSeq.

In short, I amplified an ~150bp region, which contains a 4 nucleotide library about 70 nucleotides into the read. This means that about 90% of each read is unnecessary, and all I need is the small 4 bp region.

Previously, I have grabbed out my region of interest using the following code:

grep 'TTCATTAAAAATTGAATTGACATTAACCTATAAAAATAGGCGTCGAGGCCCTTTCGTCTTC[TCAG][TCAG][TCAG][TCAG]GTCGAGTGCA' sample.fastq | cut -c 62-65 | sort | uniq -c | sort -nr | less > sample.txt

This code cuts out the library region of interest based on the character number (see cut -c) and then sorts and counts how many times that sequence occurred.

My problem is that now I have done this again, except with a slightly modified sequence. This time, there are a few homopolmeric regions on either side of my library, so occasionally they are miscalled as one nucleotide longer or shorter.

The sequence looks like this: GGGNNNNGGG

Instead of cutting by character, I would just like to cut out this region from each read, and then put it through the rest of my counting pipeline. However, I cannot seem to find a good way to cut out only this sequence from within the fastq file.

Does anyone have any suggestions to fix this?

Thank you!

-rtleenay

[mastal]

changing to something like:

grep -P 'GGG{min, max}[ACGT]{4}GGG{min,max} sample.fastq |

might work.

where min and max could be set to one less (for min), or one more (max) than the number of nucleotides in your homopolymer stretch.

[rtleenay]

But then once I do use grep to grab that line, what command (or tool) can I use to grab that stretch by its sequence instead of by the character number?

[mastal]

In perl, you could capture the sequence that matches the [ACGT]{4} using regular expressions, and print it to a file for downstream processing.