Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-22T05:03:07.670Z
Answers.LastActivityDate: 2019-07-22T05:03:07.670Z
Answers.UserDisplayName: finswimmer
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-22T04:31:01.207Z
Answers.Text: 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.
Answers.UserDisplayName: Lot_to_learn
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-22T04:33:24.797Z
Answers.Text: @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`
Answers.UserDisplayName: Scot
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-07-20T22:44:00.320Z
Answers.LastActivityDate: 2019-07-25T02:01:22.823Z
Answers.LastEditDate: 2019-07-25T02:01:22.823Z
Answers.UserDisplayName: Scot
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-07-20T13:46:19.980Z
Comments.Text: 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` ?
Comments.UserDisplayName: glenn jackman
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-07-20T14:37:55.873Z
Comments.Text: Either `zcat | grep` when you're trying to pick reads based on sequences, or `zgrep`. This code should not work as-is.
Comments.UserDisplayName: Ram RS
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-07-21T06:31:17.063Z
Comments.Text: @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.
Comments.UserDisplayName: Lot_to_learn
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-07-21T06:33:32.353Z
Comments.Text: @RamRS : Can you tell me why? and how to get this done using other ways either awk, perl,R or python.
Comments.UserDisplayName: Lot_to_learn
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-07-22T13:56:06.560Z
Comments.Text: 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.
Comments.UserDisplayName: Ram RS
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-07-23T07:06:26.967Z
Comments.Text: 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.
Comments.UserDisplayName: Lot_to_learn
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-07-25T11:39:32.413Z
Comments.Text: You say that you are "trying awk to make it faster". What else have you tried that was slower?
Comments.UserDisplayName: winni2k
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2019-07-25T11:43:17.517Z
Comments.Text: 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?
Comments.UserDisplayName: winni2k
Tags: <fastq><bash><awk>
Title: Filter reads belonging to unique sequences with threshold
ThreadText: # 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**: @username_1: 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**: @username_2 : 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
  
  ```
  
  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](http://mywiki.wooledge.org/ParsingLs), 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.
  
  > **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**: @username_0  `zcat | grep` may be a better option depending on your setup, as the `--no-group-separator` option may not work properly with `zgrep`  
  
  ---
  
  