how to parse reads containing a particular sequence in any orientation

[PoorSeq]

Good point gringer! however, because I have already flash-stitched the reads, I expect that the sequence I'm looking for will be in the same orientation in all reads. Still, I would see if I have any potholes. Also, I chose 14nt because my samples are from a bacteria with a genome size of 4.2 million, so I expect anything equal or more than 12nt will be unique.

[gringer]

bioawk -c fastx '/SEQUENCE/ {print "@"$name; print $seq; print "+"; print $qual }' inut.fq > output.fq

Hmm, bioawk seems pretty neat.

Note that what you've got there won't work for reverse complement orientation, so you'll need to have both forward and reverse included. Also, picking any 14+nt substring of SEQUENCE (or its reverse-complement) will be a bit trickier to implement.

[PoorSeq]

Thanks all for the answers, specially dariober for the code. I have also developed a bioawk code later which I was able to use. The code is:

bioawk -c fastx '/SEQUENCE/ {print "@"$name; print $seq; print "+"; print $qual }' inut.fq > output.fq

Yes, my aim is to collect all the reads that contains a sequence (or subsequence, at least 14nt) and make a new file.

[dariober]

I'm pretty new to bioinformatics and sorry if it's not worthy of asking here.

I have a FLASH stitched fastq file from my paired end data, from which I want to sort the reads containing a particular sequence or part of that sequence in any orientation, and make a new fastq file with them. Is there any easy tool/code to do that?

Hi- As gringer pointed out, you should clarify what you are trying to do. Anyway... This sequence of unix commands reads a fastq file and outputs the reads in fastq format matching a regular expression. See if it helps:

## Get reads containing substring AAA or its revcomp TTT
gunzip -c fastq.fq.gz \
| paste - - - - \
| grep -P '^@.*?\t(.*?AAA.*?)|(.*?TTT.*?)\t\+' \
| tr '\t' '\n' \
| gzip > sub.fq.gz

## Example input fastq:
@seq1
ACTGAAACTG
+comment
IIIIIIIIII
@seq2
ACTGNNNCTGTTT
+comment
BBBBBBBBBBBBB
@seq3
CCCCCCCCCCCCC
+comment
BBBBBBBBBBTTT
@seq4
AAACCCCCCCCCC
+comment
BBBBBBBBBBTTT

## Output sub.fq.gz
@seq1
ACTGAAACTG
+comment
IIIIIIIIII
@seq2
ACTGNNNCTGTTT
+comment
BBBBBBBBBBBBB
@seq4
AAACCCCCCCCCC
+comment
BBBBBBBBBBTTT

If your input is unzipped use "paste - - - - < fastq.fq" instead of "gunzip -c fastq.fq.gz \
| paste - - - -"

[yueluo]

If you are trying to get reads that align to a particular sequence, try bowtie2. Not sure if that's your purpose though?

[gringer]

I have a FLASH stitched fastq file from my paired end data, from which I want to sort the reads containing a particular sequence

Okay, that's doable with grep, and very quick.

or part of that sequence

Wait, what? You want any subsequence? Will one base pair do? What are your limits?

[PoorSeq]

I'm pretty new to bioinformatics and sorry if it's not worthy of asking here.

I have a FLASH stitched fastq file from my paired end data, from which I want to sort the reads containing a particular sequence or part of that sequence in any orientation, and make a new fastq file with them. Is there any easy tool/code to do that?