3
\$\begingroup\$

I am a beginner in python programming, and I am trying to implement multiprocessing to the code I use daily for my data treatment. Essentially, the input file has around 10 - 30 millions of reads where each read contains 4 lines of information. There are two elements in the first line of every four lines that I am looking at. First I search if the first element is in the dictionary. If not, then I add to the dictionary the first element is the key while the second element is the value. But if the first element is already in the dictionary, I proceed to compare whether the second element is within the value list corresponding to the key. Similarly, if the second element is not already present in the list of values, it will be added, otherwise it is considered a duplicate. Now as you can see, the dictionary is going to grow as more lines are read. In the end, it will just be too slow for a single process to complete the task efficiently. For my data it sometimes can take a week to do finish. Now I have a simple understanding that if I can implement something such multiple processes can be utilized to simultaneously compare with my dictionary, it would speed things up significantly (I can reserve as many as 16 CPU threads dedicated to my task). But unfortunately, multiprocessing/multithreading seems to be a complicated thing. I would appreciate if someone can direct me a path on what should I be focusing on. Below is the code that I am currently using:

def sequence_mismatch (input_seq, input_dict, input_linker_length, j):
    local_smm_list=[]
    for stored in input_dict[input_seq[0:j:]]:
        local_smm=0
        for x in range(0, 7):
            if input_seq[input_linker_length+j+x] != stored[x]:
                local_smm += 1
        local_smm_list.append(local_smm)
    if min(local_smm_list) == 0:
        return 1
    else:
        return 0



# Main
from datetime import datetime
import itertools
import sys
try:
    input_fastq = sys.argv[1]
    linker_length_str = sys.argv[2]  
    output_file_name = sys.argv[3]
    randomer_length = sys.argv[4]
except:
    print("Usage: ./PCRDupRm.py your_fastq.fastq linker_length output_file_name randomer_length")
    sys.exit(1)
startTime = datetime.now()
linker_length=int(linker_length_str)
i = int(randomer_length)
output_file = open(output_file_name, 'w+')
scanned_dict = {}
with open (input_fastq, 'r') as seq_data:
    for line1, line2, line3, line4 in itertools.izip_longest(*[seq_data]*4):
        if line2[0:i:] not in scanned_dict:
            scanned_dict.update({line2[0:i:]: [line2[(linker_length + i):(linker_length + i + 7):]]})
            output_file.write(line1)
            output_file.write(line2)
            output_file.write(line3)
            output_file.write(line4)
        else:
            comparison = sequence_mismatch (line2, scanned_dict, linker_length, i)
            if comparison == 1:
                continue
            else:
                scanned_dict[line2[0:i:]].extend([line2[(linker_length + i):(linker_length + i + 7):]])
                output_file.write(line1)
                output_file.write(line2)
                output_file.write(line3)
                output_file.write(line4)

output_file.close()
print (datetime.now() - startTime)
\$\endgroup\$
  • \$\begingroup\$ What is i in if line2[0:i:] not in scanned_dict? \$\endgroup\$ – vnp Mar 1 '17 at 18:17
  • \$\begingroup\$ Hi vnp, sorry about that, I missed a line of code stating i= int(randomer_length). now it is reflected in the code. \$\endgroup\$ – Yinzhou Zhu Mar 1 '17 at 18:22
1
\$\begingroup\$

I just figured out that I could have simply used a dictionary nested inside a dictionary instead of a list nested inside a dictionary. Now the script finishes within minutes for each sample.

\$\endgroup\$

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.