text: user2141: I am trying to write python script for customized filtering for fastq file (size >3 GB). My proposed script is as follows:
  def filtering(read):
      time.sleep(0.1)
  
      if len(read) >= 15 and \
         len(read) <= 30 and \
         np.mean(read.qual) >= 30 and \
         sum(i < 20 for i in read.qual) <=2:
             flag1 = True
             return  flag1
  
  pool = multiprocessing.Pool(5)
  pool.ncpus = 4
  for file in os.listdir('/moRNA_data/trimmed_fastq'):
      if file.endswith('.fq'):
          fastq_file = HTSeq.FastqReader( os.path.join('/moRNA_data/trimmed_fastq',file))
          name = str(file.split('.')[0]) + '.fastq'
          file_out = open(os.path.join('/moRNA_data/filtered_fastq',name), "w" )
          for read in tqdm(fastq_file):    
              flag1 = pool.apply_async(filtering,read)
              if flag1:
                  read.write_to_fastq_file( file_out )
  pool.close()
  pool.join()
  
  The problem with this script is that it works well with small file like upto 100-200 MB but when it goes to directory having multiple sample each having size more than 3 
  GB, it works well for first file but after that it hangs so severely that I had to restart the server [ I can't even kill the process]. I need suggestions to improve this script to make fastq file processing faster. Platform details is as follows:
  OS : Ubuntu-18.02
  Architecture:        x86_64
  CPU op-mode(s):      32-bit, 64-bit
  Byte Order:          Little Endian
  CPU(s):              12
  On-line CPU(s) list: 0-11
  Thread(s) per core:  2
  Core(s) per socket:  6
  Socket(s):           1
  NUMA node(s):        1
  Vendor ID:           GenuineIntel
  CPU family:          6
  Model:               158
  Model name:          Intel(R) Core(TM) i7-8700 
  CPU @ 3.20GHz
  Stepping:            10
  CPU MHz:             800.560
  CPU max MHz:         4600.0000
  CPU min MHz:         800.0000
  BogoMIPS:            6384.00
  Virtualization:      VT-x
  L1d cache:           32K
  L1i cache:           32K
  L2 cache:            256K
  L3 cache:            12288K
  NUMA node0 CPU(s):   0-11
  
  Thanks.
  
  user37: Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose.
  C++, Rust, Go, Nim and Julia can be as fast.
  // Download https://raw.githubusercontent.com/lh3/minimap2/master/kseq.h
  // Compile with: gcc -O2 -o myprog this-file.c -lz
  #include <zlib.h>  
  #include <stdio.h>  
  #include "kseq.h"
  // STEP 1: declare the type of file handler and the read() function  
  KSEQ_INIT(gzFile, gzread)  
  
  int main(int argc, char *argv[])  
  {  
      gzFile fp;  
      kseq_t *seq;  
      if (argc == 1) {  
          fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);  
          return 1;  
      }  
      fp = gzopen(argv[1], "r"); // STEP 2: open the file handler  
      seq = kseq_init(fp); // STEP 3: initialize seq  
      while (kseq_read(seq) >= 0) { // STEP 4: read sequence  
          int i, sum, low;
          if (seq->qual.l == 0) continue; // no quality
          if (seq->seq.l < 15 || seq->seq.l > 30) continue;
          for (i = sum = low = 0; i < seq->seq.l; ++i) {
              int qual = seq->qual.s[i] - 33; // assuming Sanger quality
              sum += qual;
              if (qual < 20) ++low;
          }
          if (sum < seq->seq.l * 30.0 || low > 2) continue;
          printf("@%s\n", seq->name.s); // output can be made faster
          puts(seq->seq.s);
          printf("+\n");
          puts(seq->qual.s);
      }  
      kseq_destroy(seq); // STEP 5: destroy seq  
      gzclose(fp); // STEP 6: close the file handler  
      return 0;  
  }  
  
  To apply the program to all *.fq in a directory (again, not tested):
  mkdir -p filtered
  ls *.fq | xargs -i echo ./myprog {} \> filtered/{} | sh
  
  You can also replace the final sh with parallel, but the above should be fast enough.
  
original_text: 2141: <p>I am trying to write python script for customized filtering for fastq file (size >3 GB). My proposed script is as follows:</p>
  
  <pre><code>def filtering(read):
      time.sleep(0.1)
  
      if len(read) &gt;= 15 and \
         len(read) &lt;= 30 and \
         np.mean(read.qual) &gt;= 30 and \
         sum(i &lt; 20 for i in read.qual) &lt;=2:
             flag1 = True
             return  flag1
  
  pool = multiprocessing.Pool(5)
  pool.ncpus = 4
  for file in os.listdir('/moRNA_data/trimmed_fastq'):
      if file.endswith('.fq'):
          fastq_file = HTSeq.FastqReader( os.path.join('/moRNA_data/trimmed_fastq',file))
          name = str(file.split('.')[0]) + '.fastq'
          file_out = open(os.path.join('/moRNA_data/filtered_fastq',name), "w" )
          for read in tqdm(fastq_file):    
              flag1 = pool.apply_async(filtering,read)
              if flag1:
                  read.write_to_fastq_file( file_out )
  pool.close()
  pool.join()
  </code></pre>
  
  <p>The problem with this script is that it works well with small file like upto 100-200 MB but when it goes to directory having multiple sample each having size more than 3 
  GB, it works well for first file but after that it hangs so severely that I had to restart the server [ I can't even kill the process]. I need suggestions to improve this script to make fastq file processing faster. Platform details is as follows:</p>
  
  <pre><code>OS : Ubuntu-18.02
  Architecture:        x86_64
  CPU op-mode(s):      32-bit, 64-bit
  Byte Order:          Little Endian
  CPU(s):              12
  On-line CPU(s) list: 0-11
  Thread(s) per core:  2
  Core(s) per socket:  6
  Socket(s):           1
  NUMA node(s):        1
  Vendor ID:           GenuineIntel
  CPU family:          6
  Model:               158
  Model name:          Intel(R) Core(TM) i7-8700 
  CPU @ 3.20GHz
  Stepping:            10
  CPU MHz:             800.560
  CPU max MHz:         4600.0000
  CPU min MHz:         800.0000
  BogoMIPS:            6384.00
  Virtualization:      VT-x
  L1d cache:           32K
  L1i cache:           32K
  L2 cache:            256K
  L3 cache:            12288K
  NUMA node0 CPU(s):   0-11
  </code></pre>
  
  <p>Thanks.</p>
   37: <p>Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose.</p>
  
  <p>C++, Rust, Go, Nim and Julia can be as fast.</p>
  
  <pre class="lang-cpp prettyprint-override"><code>// Download https://raw.githubusercontent.com/lh3/minimap2/master/kseq.h
  // Compile with: gcc -O2 -o myprog this-file.c -lz
  #include &lt;zlib.h&gt;  
  #include &lt;stdio.h&gt;  
  #include "kseq.h"
  // STEP 1: declare the type of file handler and the read() function  
  KSEQ_INIT(gzFile, gzread)  
  
  int main(int argc, char *argv[])  
  {  
      gzFile fp;  
      kseq_t *seq;  
      if (argc == 1) {  
          fprintf(stderr, "Usage: %s &lt;in.seq&gt;\n", argv[0]);  
          return 1;  
      }  
      fp = gzopen(argv[1], "r"); // STEP 2: open the file handler  
      seq = kseq_init(fp); // STEP 3: initialize seq  
      while (kseq_read(seq) &gt;= 0) { // STEP 4: read sequence  
          int i, sum, low;
          if (seq-&gt;qual.l == 0) continue; // no quality
          if (seq-&gt;seq.l &lt; 15 || seq-&gt;seq.l &gt; 30) continue;
          for (i = sum = low = 0; i &lt; seq-&gt;seq.l; ++i) {
              int qual = seq-&gt;qual.s[i] - 33; // assuming Sanger quality
              sum += qual;
              if (qual &lt; 20) ++low;
          }
          if (sum &lt; seq-&gt;seq.l * 30.0 || low &gt; 2) continue;
          printf("@%s\n", seq-&gt;name.s); // output can be made faster
          puts(seq-&gt;seq.s);
          printf("+\n");
          puts(seq-&gt;qual.s);
      }  
      kseq_destroy(seq); // STEP 5: destroy seq  
      gzclose(fp); // STEP 6: close the file handler  
      return 0;  
  }  
  </code></pre>
  
  <p>To apply the program to all <code>*.fq</code> in a directory (again, <em>not</em> tested):</p>
  
  <pre class="lang-sh prettyprint-override"><code>mkdir -p filtered
  ls *.fq | xargs -i echo ./myprog {} \&gt; filtered/{} | sh
  </code></pre>
  
  <p>You can also replace the final <code>sh</code> with <code>parallel</code>, but the above should be fast enough.</p>
  