Filter reads belonging to unique sequences with threshold

Question

I have CLL samples as fastq files and I want to remove those reads which have a unique sequence with less than 10 read counts each. I tried this by following way using awk to make it faster:

for fq in $(ls fastq1/*.fq.gz); do
    basename=$(basename "$fq" | cut -f1 -d '.')
    zcat $fq | awk "/^@/{getline; print}" |sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt; grep -B 1 -A 2 -f tmp.txt $fq | grep -v ^--'| pigz -p 10 > fastq2/$basename.fastq.gz

done

rm tmp.txt

But this is not working well as many reads with more than 10 read counts are also missing. I think it's OK to stick with bash/awk for fast fastq file processing. What am I missing here, or is this command wrong. Can anyone suggest me the right and fast way to do it??

Answer 1

Checking for line start to find the read header is not the best option. Even if unlikely it is possible that the qualitys starts with a @. It's better to compare the current line number. Also you can use awk directly to count and print the result.

zcat $fq|awk 'NR%4==2 {count[$0]++} END { for(sequence in count) { if(count[sequence]>=10) print sequence }}'

As the other statet out, the second command should be fine if you use zcat instead of cat.

Answer 2

Your algorithm generally works, you just had one minor issue, and a few stylistic changes that could be made.

I'd adjust your script to the following:

for fq in fastq1/*.fq.gz; do
    basename=$(basename "$fq" .fq.gz)
    zcat "$fq" | awk "/^@/{getline; print}" | sort | uniq --count | awk '{if($1>= 10) print $2}' > tmp.txt
	zcat | grep --no-group-separator -B 1 -A 2 -f tmp.txt "$fq" | pigz -p 10 > "fastq2/$basename.fastq.gz"
done

rm tmp.txt

1. Instead of

for fq in $(ls fastq1/*.fq.gz)

use this

for fq in fastq1/*.fastq.gz

Iterating over ls is generally not ideal, and can give unexpected results. Better to just rely on bash globbing.

2. Purely stylistic, but I'd simplify and use basename to remove the suffix, instead of cut.

3. Rather than removing the separators that grep will introduce by filtering them out, use the following option to not create them: grep --no-group-separator

4. And, minor bug, as others have mentioned, you should use zgrep or zcat | grep since your files are compressed.

Oh - and quote your variables.

Finally, while this does work, it may not be the most efficient method. The file is being decompressed twice, and may be better to decompress once to a temporary file that can be used and then deleted.