<p>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.</p>

<pre><code>zcat <span class="math-container">$fq|awk 'NR%4==2 {count[$</span>0]++} END { for(sequence in count) { if(count[sequence]&gt;=10) print sequence }}'
</code></pre>

<p>As the other statet out, the second command should be fine if you use <code>zcat</code> instead of <code>cat</code>.</p>

CC BY-SA 4.0
2019-07-22T05:03:07.670Z
2019-07-22T05:03:07.670Z
finswimmer
<p>Your algorithm generally works, you just had one minor issue, and a few stylistic changes that could be made.</p>

<p>I'd adjust your script to the following:</p>

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

rm tmp.txt
</code></pre>

<ol>
<li>Instead of</li>
</ol>

<p><code>for fq in $(ls fastq1/*.fq.gz)</code></p>

<p>use this</p>

<p><code>for fq in fastq1/*.fastq.gz</code></p>

<p>Iterating over <code>ls</code> is <a href="http://mywiki.wooledge.org/ParsingLs" rel="nofollow noreferrer">generally not ideal</a>, and can give unexpected results. Better to just rely on bash globbing.</p>

<ol start="2">
<li><p>Purely stylistic, but I'd simplify and use <code>basename</code> to remove the suffix, instead of <code>cut</code>.</p></li>
<li><p>Rather than removing the separators that <code>grep</code> will introduce by filtering them out, use the following option to not create them:
<code>grep --no-group-separator</code></p></li>
<li><p>And, minor bug, as others have mentioned, you should use <code>zgrep</code> or <code>zcat | grep</code> since your files are compressed.</p></li>
</ol>

<p>Oh - and quote your variables.</p>

<p>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.</p>

CC BY-SA 4.0
2019-07-22T04:31:01.207Z
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
CC BY-SA 4.0
2019-07-22T04:33:24.797Z
@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
CC BY-SA 4.0
2019-07-20T22:44:00.320Z
2019-07-25T02:01:22.823Z
2019-07-25T02:01:22.823Z
Scot
<p>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:</p>

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

done

rm tmp.txt
</code></pre>

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

CC BY-SA 4.0
2019-07-20T13:46:19.980Z
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
CC BY-SA 4.0
2019-07-20T14:37:55.873Z
Either `zcat | grep` when you're trying to pick reads based on sequences, or `zgrep`. This code should not work as-is.
Ram RS
CC BY-SA 4.0
2019-07-21T06:31:17.063Z
@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
CC BY-SA 4.0
2019-07-21T06:33:32.353Z
@RamRS : Can you tell me why? and how to get this done using other ways either awk, perl,R or python.
Lot_to_learn
CC BY-SA 4.0
2019-07-22T13:56:06.560Z
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
CC BY-SA 4.0
2019-07-23T07:06:26.967Z
My file `$fq` is surely gzipped. For this reason, I started with `zcat`. As per the answer suggested below, `zcat|grep` can work well with extracted unique sequences in the second step.
Lot_to_learn
CC BY-SA 4.0
2019-07-25T11:39:32.413Z
You say that you are "trying awk to make it faster". What else have you tried that was slower?
winni2k
CC BY-SA 4.0
2019-07-25T11:43:17.517Z
I have not encountered this problem before, and I am wondering if this is an [XY problem](http://xyproblem.info/). Why are you trying to do this?
winni2k
<fastq><bash><awk>
Filter reads belonging to unique sequences with threshold
# Filter reads belonging to unique sequences with threshold

(Asked by: username_0 on 2019-07-20)

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

**Tags:** fastq, bash, awk
> **Comments:**
> **username_1**: 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` ?  
> **username_2**: Either `zcat | grep` when you're trying to pick reads based on sequences, or `zgrep`. This code should not work as-is.  
> **username_0**: @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.  
> **username_0**: @RamRS : Can you tell me why? and how to get this done using other ways either awk, perl,R or python.  
> **username_2**: 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.  
> **username_0**: My file `$fq` is surely gzipped. For this reason, I started with `zcat`. As per the answer suggested below, `zcat|grep` can work well with extracted unique sequences in the second step.  
> **username_3**: You say that you are "trying awk to make it faster". What else have you tried that was slower?  
> **username_3**: I have not encountered this problem before, and I am wondering if this is an [XY problem](http://xyproblem.info/). Why are you trying to do this?  

---

## 2 answers

### Answer by username_4 on 2019-07-22. Score: 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 by username_5 on 2019-07-20. Score: 1

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


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.

Purely stylistic, but I'd simplify and use basename to remove the suffix, instead of cut.
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
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.

> **Comments:**
> **username_0**: 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.  
> **username_5**: @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`  

---

