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??
zcat $fq
and in the middle yougrep -B 1 -A 2 -f tmp.txt $fq
-- do you need tozcat $fq | grep -B1 -A1 -f tmp.txt
? $\endgroup$zcat | grep
when you're trying to pick reads based on sequences, orzgrep
. This code should not work as-is. $\endgroup$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 anytmp.txt
but that didn't work. Rest other arguments ofgrep
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. $\endgroup$$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. $\endgroup$