Commit f55201d8 authored by Theo Serralta's avatar Theo Serralta

final code to find CNV with multithreading

parent 70f8dbd4
...@@ -2,10 +2,13 @@ import pysam ...@@ -2,10 +2,13 @@ import pysam
import concurrent.futures import concurrent.futures
import sys import sys
import getopt import getopt
import logging
#Options #Options
output_file = ""
try: try:
opts, args = getopt.getopt(sys.argv[1:], 'b:w:s:t:o:e') opts, args = getopt.getopt(sys.argv[1:], 'b:w:s:t:o:e:')
for opt, arg in opts: for opt, arg in opts:
if opt in ("-b"): if opt in ("-b"):
bamfile = arg bamfile = arg
...@@ -16,15 +19,18 @@ try: ...@@ -16,15 +19,18 @@ try:
if opt in ("-t"): if opt in ("-t"):
num_threads = int(arg) num_threads = int(arg)
if opt in ("-o"): if opt in ("-o"):
output_file = open(arg, "w") output_file = open(arg, "a")
if opt in ("-e") : if opt in ("-e"):
logfile = open(arg, "w") logfile = arg
except getopt.GetoptError: except getopt.GetoptError:
print('Invalid argument') print('Invalid argument')
sys.exit(1) sys.exit(1)
def calcul_depth_seq(bamfile, seq, window_size, step_size): logging.basicConfig(filename = '%s' % (logfile), filemode = 'a', level = logging.INFO, format = '%(asctime)s %(levelname)s - %(message)s')
logging.info('start')
def calcul_depth_seq(bamfile, seq, window_size, step_size, numero_de_thread):
sys.stderr.write("\t entering calcul_depth_seq : thread numero %s\n" % numero_de_thread)
samfile = pysam.AlignmentFile(bamfile, "rb") samfile = pysam.AlignmentFile(bamfile, "rb")
seq_length = samfile.lengths[samfile.references.index(seq)] #Permet d'obtenir la longueur de la sequence seq_length = samfile.lengths[samfile.references.index(seq)] #Permet d'obtenir la longueur de la sequence
...@@ -36,30 +42,31 @@ def calcul_depth_seq(bamfile, seq, window_size, step_size): ...@@ -36,30 +42,31 @@ def calcul_depth_seq(bamfile, seq, window_size, step_size):
window_depths.append(pileupcolumn.n)#ajoute le nb de lecture couvrant cette position à la liste window_depths 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 if window_depths:#calcul la profondeur moyenne pour chaque sequence
sys.stderr.write("\t\t entering wd thread numero %s\n" % numero_de_thread)
avg_depth = sum(window_depths) / len(window_depths) avg_depth = sum(window_depths) / len(window_depths)
print(f"{seq}:{pos_start}-{pos_end}: {avg_depth}") sys.stdout.write(f"{seq}:{pos_start}-{pos_end}: {avg_depth}\n")
output_file.write(f"{seq}:{pos_start}-{pos_end}: {avg_depth}\n")
samfile.close() sys.stderr.write("\t leaving calcul_depth_seq thread numero %s\n" % numero_de_thread)
def calcul_depth_all(bamfile, window_size, step_size, num_threads): def calcul_depth_all(bamfile, window_size, step_size, num_threads):
sys.stderr.write("entering calcul_depth_all \n")
samfile=pysam.AlignmentFile(bamfile, "rb") samfile=pysam.AlignmentFile(bamfile, "rb")
with concurrent.futures.ThreadPoolExecutor(max_workers = num_threads) as executor: with concurrent.futures.ThreadPoolExecutor(max_workers = num_threads) as executor:
futures = [] #liste pour stocker les futures futures = [] #liste pour stocker les futures
n=0
for seq in samfile.references: for seq in samfile.references:
future = executor.submit(calcul_depth_seq, bamfile, seq, window_size, step_size) n+=1
future = executor.submit(calcul_depth_seq, bamfile, seq, window_size, step_size, n)
futures.append(future) futures.append(future)
concurrent.futures.wait(futures) concurrent.futures.wait(futures)
samfile.close() samfile.close()
if output_file != sys.stdout: #if output_file != sys.stdout:
output_file.close() # output_file.close()
if logfile:
logfile.close()
#Programme principal #Programme principal
if len(sys.argv) <2 or bamfile is None: if len(sys.argv) <2 or bamfile is None:
...@@ -68,3 +75,5 @@ if len(sys.argv) <2 or bamfile is None: ...@@ -68,3 +75,5 @@ if len(sys.argv) <2 or bamfile is None:
bamfile_path = bamfile bamfile_path = bamfile
calcul_depth_all(bamfile_path, window_size, step_size, num_threads) calcul_depth_all(bamfile_path, window_size, step_size, num_threads)
logging.info('end')
samfile.close()
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