Extract specific reads from FASTQ files based on subsequence

Question

[Paul]

Dear all,

I have FASTQ files and on start of my read I have 7 nucleotides tag - I would like to extract reads with this specific tag.

I would like to search in first 15 nucleotides of my reads, if match - extract this read to new fastq files.

Thank you for any ideas or help.

[springpython]

any other tools for sam question?

[harold.smith.tarheel]

Why? Are all of these tools inadequate?

[Ram]

There's a wonderful quote in Jurassic Park: "...were so preoccupied with whether they could, they didn't stop to think if they should."

[Brian Bushnell]

Nice :) Jurassic Park has a lot of apropos quotables for genetic research...

Personally, I like "Life will find a way"

[Ram]

"Life ... uhh ... finds a way" :)

[Brian Bushnell]

My meme was more eloquent ;)

[Ram]

We need an /r/biostars

[springpython]

I am actually stuck in python code for same question. hence asked for suggestions.

[Ram]

Don't add answers unless you're answering the original question. Use Add Reply or Add Comment instead. Please read https://www.biostars.org/t/how-to/

I'm moving this to a comment on the top-level post now.

Answer 1

[5heikki]

Probably not the fastest option:

grep -B1 -A2 "^AGATCGG" file.fq | grep -v "^--$" > out.fq

Find lines that begin with "AGATCGG", grab one line before each hit and two lines after each hit, remove lines that are "--".

Note that it's possible (but extremely unlikely) that some quality value line begins like "AGATCGG", in which case the above command would mess up the output file. The likelihood of this is probably very close to zero but if you're processing a file with googolplex lines, maybe it could happen.

Answer 2

[Ram]

You can use Heng Li's bioawk:

To check if the tag is part of the first 15 bases

bioawk -c fastx 'substr($seq,0,15) ~ /$TAG/ { print }' reads.fq.gz

To match the first 7 bases to your tag,

bioawk -c fastx 'substr($seq,0,7) == $TAG { print }' reads.fq.gz

Answer 3

[michael.ante]

In order to stay old-school, you can use the FastX toolkit's barcode splitter.

You need a text file with your tag (let it be myTag.txt) and then you can run (with N mismatches):

cat sequence.fastq | /usr/local/bin/fastx_barcode_splitter.pl -Q33 --bcfile myTag.txt --bol --mismatches $N --prefix out --suffix .txt

In order to remove the tag-sequence, you can use the fastx_trimmer from the same toolkit

Using the Fastx-tools, do not forget to add the -Q33 option.

Answer 4

[Brian Bushnell]

You can also do this with BBDuk:

bbduk.sh -Xmx1g in=reads.fq out=filtered.fq k=7 mm=f rcomp=f restrictleft=15 literal=ACGTACG

If you want multiple tags, you can list them separated by commas; and if you want to allow 1bp mismatch, set hdist=1.