Find CNV with Read Depth method on CPU with multithreading

CNV/Read_Depth_Multithreading.py

+import pysam
+import concurrent.futures
+import sys
+
+
+def calcul_depth_seq(bamfile, seq, window_size=100, step_size=10):
+
+    samfile = pysam.AlignmentFile(bamfile, "rb")
+    seq_length = samfile.lengths[samfile.references.index(seq)] #Permet d'obtenir la longueur de la sequence
+
+    for pos_start in range(0, seq_length - window_size +1, step_size): #Itère sur la sequence par "fenêtre" en fonction de la taille
+        pos_end = pos_start + window_size #pos_start + window_size donne la position de début de la fenêtre, puis grâce à min() et seq_length je peux définir position de fin. le min me permet d'être certain que pos_end < seq_lenght. C-a-d que si la fenêtre dépasse la taille de la sequence, j'aurais quand même pos_end
+        window_depths = [] #liste pour stocker la profondeur de chaque fenêtre
+
+        for pileupcolumn in samfile.pileup(seq, start=pos_start, stop=pos_end): #itère sur chaque colone de pileup
+            window_depths.append(pileupcolumn.n)#ajoute le nb de lecture couvrant cette position à la liste window_depths
+
+        if window_depths:#calcul la profondeur moyenne pour chaque sequence
+            avg_depth = sum(window_depths) / len(window_depths)
+            print(f"{seq}:{pos_start}-{pos_end}: {avg_depth}")
+           
+    samfile.close()
+    
+def calcul_depth_all(bamfile, window_size=100, step_size=10, num_threads=4):
+    samfile=pysam.AlignmentFile(bamfile, "rb")
+    
+    with concurrent.futures.ThreadPoolExecutor(max_workers = num_threads) as executor:
+        futures = [] #liste pour stocker les futures
+    
+        for seq in samfile.references:
+            future = executor.submit(calcul_depth_seq, bamfile, seq, window_size, step_size)
+            futures.append(future)
+            
+        concurrent.futures.wait(futures)
+        
+        samfile.close()
+
+#Programme principal
+if len(sys.argv) <2:
+    print("Usage : File alignement .bam")
+    sys.exit()
+bamfile_path = sys.argv[1]
+
+calcul_depth_all(bamfile_path, window_size=100, num_threads=4)