Simple shell command for extracting indexed reads from fastq?

Question

[mbk0asis]

Hi.

All my reads in fastq files are supposed to have custom index sequences at 5th nucleotide from 5' end.

(for example, 5' NNNN-IndexSeq-restReadSequences)

I tried to demultiplexed them with several tools on the web, but I couldn't get consistent results from those tools.

For just "sanity check", I want to actually count or extract indexed reads with simple shell commands.

I'm not sure if I'm saying right, but most tools search indexes in the whole reads, but I want to look for indexes at specific position (i.e. 5th - 9th nucleotides).

I hope you guys could understand my poor English.

Thank you.

Answer 1

[Pierre Lindenbaum]

use awk. For example to filter the reads having bases 1-5 = "TCTCT"

gunzip -c your.fastq.gz | paste - - - - |\
awk -F '\t'  '(substr($2,1,5)=="TCTCT")' |\
tr "\t" "\n"

[dariober]

Shouldn't substr($2,1,5) be substr($2,5,5)? (Get substring from position 5 for length 5)

[Pierre Lindenbaum]

I said "bases 1-5"

Answer 2

[dariober]

Expanding Pierre's answer... You might want to check which are the most represented indexes in, say, the first 100000 reads:

zcat reads.fq.gz \
| head -n 100000 \
| paste  - - - - \
| awk '{idx=substr($2, 5, 5); print idx}' \
| sort \
| uniq -c \
| sort -k1,1nr

Remove the last three lines if you just want to see the indexes of each read.

Answer 3

[smilefreak]

Using grep before and after switches.

# extracts reads.
zcat reads.fq.gz | grep -A 2 -B 1 "^NNNNTCTCTC" > reads_out.txt
# count number of matches
zcat reads.fq.gz | grep "^NNNNTCTCTC" | wc -l