Fast processing of fastq data

Question

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.

Answer 1

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.