I have written a very basic algorithm, that counts the amount of times sub-string appears in a string.
The algorithm is given:
- seq - a string sequence
- k - length of a sub-string
- L - a window of a string to search in
- t - times that sub-string, length k has to be present to be added to results
from collections import defaultdict
import concurrent.futures
def kmer_clumps(seq, k, L, t):
"""Returns a set of L-t kmers that form a clump. If a kmer of length [k] apperas [t] times in window [L], it is added to a set"""
print("Starting a job")
clump_dict = defaultdict(list)
clumps = set()
for pos in range(len(seq) - k + 1):
current_pos = seq[pos:pos + k]
clump_dict[current_pos].append(pos)
if len(clump_dict[current_pos]) >= t:
if ((clump_dict[current_pos][-1] + (k - 1)) - clump_dict[current_pos][-t]) < L:
clumps.add(current_pos)
if clumps:
return clumps
else:
return None
Running this algorithm in one thread/core is easy, just pass all of the parameters and it scans the string and returns a list.
Now, this is a perfect candidate for running on multiple cores as sometimes I need to provide a large amount of data to this algorithm (text file 20mb+ in size or even more).
Below, I have used Python's import multiprocessing
and/or import concurrent.futures
, but this does not really matter as the question is: what is the best way of splitting the data to be passed to multiple threads/cores to avoid duplicate data, large amounts of memory usage, etc... I am hoping that people can share their experience in preparing data in this manner and best practices for this.
seq = "CGGACTCGACAGATGTGAAGAACGACAATGTGAAGACTCGACACGACAGAGTGAAGAGAAGAGGAAACATTGTAA"
k = 5
L = 50
t = 4
cores = 2
def sub_seq_clumps_multicore_job(seq, k, L, t, cores=1):
seqLength = len(seq)
if cores == 1:
print("Using only 1 core as specified or data is to small.")
return kmer_clumps(seq, k, L, t)
else:
print("Preparing data...")
# Basic logic to split the data to provide to multiple jobs/threads
jobSegments = seqLength // cores
extraSegments = seqLength % cores
overlapSize = (L // cores)
# Creating multiple parts from a single string and adding same arguments to all
# Amount of parts is based on core/thread count
parts = []
for part in range(0, seqLength, jobSegments + extraSegments):
tmpList = [seq[part:part + jobSegments +
extraSegments + overlapSize + 1], k, L, t]
parts.append(tmpList)
print(f"Processing data on {cores} cores...")
resultSet = set()
# Starting jobs/threads by passing all parts to a thread pool
with concurrent.futures.ProcessPoolExecutor() as executer:
jobs = [
executer.submit(kmer_clumps, *sec) for sec in parts
]
# Just collecting results
for f in concurrent.futures.as_completed(jobs):
resultSet.update(f.result())
# Returning collected results, which are sub-strings of length k,
# that were found in all windows length L
return resultSet
print(sub_seq_clumps_multicore_job(seq, k, L, t, cores))
Some sample output:
Data, segmented and passed to kmer_clumps(seq, k, L, t)
to run on 1 core/thread:
"CGGACTCGACAGATGTGAAGAACGACAATGTGAAGACTCGACACGACAGAGTGAAGAGAAGAGGAAACATTGTAA", 5, 50, 4
Data, segmented and passed to kmer_clumps(seq, k, L, t)
to run on 2 core/thread:
Thread 1:
['CGGACTCGACAGATGTGAAGAACGACAATGTGAAGACTCGACACGACAGAGTGAAGAGAAGAGG', 5, 50, 4]
Thread 2:
['CGACACGACAGAGTGAAGAGAAGAGGAAACATTGTAA', 5, 50, 4]
Data, segmented and passed to kmer_clumps(seq, k, L, t)
to run on 4 core/thread:
Thread 1:
['CGGACTCGACAGATGTGAAGAACGACAATGTGAA', 5, 50, 4]
Thread 2:
['ACGACAATGTGAAGACTCGACACGACAGAGTGAA', 5, 50, 4]
Thread 3:
['ACGACAGAGTGAAGAGAAGAGGAAACATTGTAA', 5, 50, 4]
Thread 4:
['GAAACATTGTAA', 5, 50, 4]]
Results in running time, when reading a huge data file (4.5 MB):
File reading complete!
Using only 1 core as specified or data is to small.
Starting a job
real 0m4.536s
File reading complete!
Preparing data...
Processing data on 2 cores...
Starting a job
Starting a job
real 0m2.438s
File reading complete!
Preparing data...
Processing data on 4 cores...
Starting a job
Starting a job
Starting a job
Starting a job
real 0m1.413s
File reading complete!
Preparing data...
Processing data on 6 cores...
Starting a job
Starting a job
Starting a job
Starting a job
Starting a job
Starting a job
real 0m1.268s
Because of the nature of the algorithm, when running on multiple threads/cores there is more data generated as we need to make sure we cover the whole data string and take a window L into the account. So there are data overlaps happening:
Let's take two examples:
With this post I have two goals:
1) To have experienced people look at it and suggest any improvements to how to segment the data for threading like that, to avoid data duplicates, memory overload, etc... any suggestions would be very helpful.
2) To learn from amazing people and so others can learn from this example.
Thank you all very much!