3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Make sure that filtering() actually returns something, though that's likely not the root of the problem. $\endgroup$ – Devon Ryan Jun 12 at 8:56
  • $\begingroup$ @DevonRyan just to clarify, tqdm only reports time in such cases it does not do a list conversion as you indicated- you can easily test with with an infinite generator pastebin.com/NJvwf5DT $\endgroup$ – Chris_Rands Jun 12 at 8:58
  • $\begingroup$ Take out the time.sleep() call for starters. and I think you're missing an else: flag1 = False clause. If you still have issues, remove the multiprocessing code to assist your debugging $\endgroup$ – Chris_Rands Jun 12 at 9:07
  • $\begingroup$ Did you try it on just one CPU or two just to test if the problem is with the threading and other applications on the server? Also how much RAM and memory is available? $\endgroup$ – llrs Jun 12 at 9:49
  • 1
    $\begingroup$ If you can't kill the process and your only option is to restart your server, it seems that you exceed your available RAM. So you could monitor the RAM usage at first and check how to clear the memory. $\endgroup$ – Mr_Z Jun 12 at 13:41
6
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ :Thanks for your response. Till now I have done all my analysis in either python or R. I have never touched C/C++. Although I studies those language but I left using this before 8 or more years. But still I think I should try your way [although it is very much time taking for me to start heavy text file processing in C/C++]. For time being, Do you have any python alternative solution for this problem? $\endgroup$ – Lot_to_learn Jun 14 at 1:52
  • $\begingroup$ @Lot_to_learn The inner loop mimics your script. It won't take you much time to modify it for other purposes even if you don't know C well. For tasks like this, C will be times faster than python whatever you do. $\endgroup$ – user172818 Jun 14 at 2:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.