2
$\begingroup$

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

$\endgroup$
8
  • 2
    $\begingroup$ 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 ? $\endgroup$ Jul 20, 2019 at 13:46
  • $\begingroup$ Either zcat | grep when you're trying to pick reads based on sequences, or zgrep. This code should not work as-is. $\endgroup$
    – Ram RS
    Jul 20, 2019 at 14:37
  • $\begingroup$ @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. $\endgroup$ Jul 21, 2019 at 6:31
  • $\begingroup$ @RamRS : Can you tell me why? and how to get this done using other ways either awk, perl,R or python. $\endgroup$ Jul 21, 2019 at 6:33
  • $\begingroup$ 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. $\endgroup$
    – Ram RS
    Jul 22, 2019 at 13:56

2 Answers 2

1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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.

  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.

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ Jul 22, 2019 at 4:31
  • $\begingroup$ @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 $\endgroup$
    – Scot
    Jul 22, 2019 at 4:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.