Commit 773d7dc8 authored by Theo Serralta's avatar Theo Serralta

Optimize function calcul_mappability

parent 6ed38039
......@@ -8,7 +8,6 @@ import pycuda.autoinit
from pycuda.compiler import SourceModule
import pycuda.gpuarray as gpuarray
from pycuda.autoinit import context
from bisect import bisect_right
# Options
try:
......@@ -157,16 +156,24 @@ def dico_mappabilite(mappability_file):
return mappability_dico #Dictionnaire avec les bornes de mappabilité en fonction des positions pour chaque chromosome.
def calcul_mappability(seq_length, mappability, chr):
map_data = np.zeros(seq_length, dtype = np.float32)
map_data = np.zeros(seq_length, dtype=np.float32)
sys.stderr.write("\t map_data =\n")
sorted_keys = sorted(mappability[chr].keys())
sys.stderr.write("\t sorted_keys =\n")
for i in range(seq_length):
bound_index = bisect_right(sorted_keys, i) - 1
nearest_bound = sorted_keys[bound_index] if bound_index >= 0 else 0
map_data[i] = mappability[chr][nearest_bound]
#print(f" i : {i}\t sorted_keys[i] : {sorted_keys[i]}\tbound_index : {bound_index}\t nearest_bound : {nearest_bound}\t map_data[i] : {map_data[i]}\t sorted_keys[bound_index] : {sorted_keys[bound_index]}")
prev_bound = 0
prev_mappability = 0
for bound in sorted_keys:
for i in range(prev_bound, min(seq_length, bound)):
map_data[i] = prev_mappability
prev_bound = bound
prev_mappability = mappability[chr][bound]
# Fill remaining positions if sequence length exceeds last bound
for i in range(prev_bound, seq_length):
map_data[i] = prev_mappability
return map_data
#############################################
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment