Filter reads belonging to unique sequences with threshold

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??

• I don't know fastq, but at the start of the pipeline you zcat $fq and in the middle you grep -B 1 -A 2 -f tmp.txt$fq -- do you need to zcat $fq | grep -B1 -A1 -f tmp.txt ? – glenn jackman Jul 20 '19 at 13:46 • Either zcat | grep when you're trying to pick reads based on sequences, or zgrep. This code should not work as-is. – Ram RS Jul 20 '19 at 14:37 • @glennjackman: In tmp.txt, I am collecting all unique sequences and finally using this file to count the number of occurrences in the original file. I tried to make it continuous without creating any tmp.txt but that didn't work. Rest other arguments of grep are used to take a complete chunk of 4 lines in fastq files it satisties the criteria of having the number of occurrences more than 10. For more information, you can check the fastq file format. – Lot_to_learn Jul 21 '19 at 6:31 • @RamRS : Can you tell me why? and how to get this done using other ways either awk, perl,R or python. – Lot_to_learn Jul 21 '19 at 6:33 • Your file $fq is either gzipped or not gzipped. If it is gzipped, grep ... $fq should not work. If it's not gzipped, zcat$fq should not work. I don't think both of those statements can work with the same $fq file. – Ram RS Jul 22 '19 at 13:56 2 Answers 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. 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


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.

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

2. 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

3. 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.

• Thanks for your response. I tried zgrep, but I am not getting any output. Although the script is working till generating tmp.txt. But after that, I didn't see any output file in fastq2. So I tried zcat|grep and it was working. I'll check whether it is giving a correct response or not after completing the pipeline. – Lot_to_learn Jul 22 '19 at 4:31
• @Lot_to_learn zcat | grep may be a better option depending on your setup, as the --no-group-separator option may not work properly with zgrep – Scot Jul 22 '19 at 4:33