Getting sequences from fastq file using Grep command

Question

[mohsamir2016]

I have been trying to get a sequence (e.g. GCGAGCCCCACATCGCCCCCCCGATTGTAATAAATAA) from a fastq file (file.fastq) I have and output a fq file. I have tried the command:

grep -A 2 -B 1 'GCGAGCCCCACATCGCCCCCCCGATTGTAATAAATAA' file.fastq | sed '/--/d' > output.fq

I got an output as 6- lines of reads that are different, while my target sequence is basically more in number. A part from that they might be just not there, Are the - A and - B is getting only the upstream and trailing sequences surroundng my target or do they get these along with the target sequences. I red the -- help page, but could not understand it. Thanks for any help.

[Pierre Lindenbaum]

use --no-group-separator instead of sed

[Joe]

Are you confident this sequence will only appear once?

[mohsamir2016]

Actually it is much longer I only gives example...

[Pierre Lindenbaum]

let's try another way. Show us the output of:

cat file.fastq | paste - - - - | grep -F GCGAGCCCCACATCGCCCCCCCGATTGTAATAAATAA 

[mohsamir2016]

Could you let me know what paste is doing ?

Thanks

[Joe]

Converting to a single line

Answer 1

[GenoMax]

Use bbduk.sh from BBMap suite in filter mode. grep has its place but a proper tool is foolproof.

bbduk.sh -Xmx2g in=your.fq outm=filtered.fq literal=GCGAGCCCCACATCGCCCCCCCGATTGTAATAAATAA

Add rcomp=f if you don't want to find reverse complemented sequence.

Answer 2

[shenwei356]

Try seqkit grep, which searches on both strands, you might need to add --only-positive-strand.

seqkit grep -s -p GCGAGCCCCACATCGCCCCCCCGATTGTAATAAATAA file.fastq -o out.fq.gz

Answer 3

[swbarnes2]

If things aren't working the way you think they should, go simpler to troubleshoot. Do

grep 'GCGAGCCCCACATCGCCCCCCCGATTGTAATAAATAA' file.fastq  | wc -l

Do you get more than 6 lines?

(Is your fastq really unzipped?)

[mohsamir2016]

It is actually unzipped and it gave me 6-lines when running your line

[swbarnes2]

Then that exact sequence is only in there 6 times.