text: user2141: 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??
  
  user92: 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.
  
  user1310: 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.
  
original_text: 2141: <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>
   92: <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>
   1310: <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>
  