This code aims to calculate arabidopsis thaliana methylation statistics it does not use chunk

Statistics_Arabidopsis_DeepSignalPlant.py → Statistics_Methylation_Arabidopsis_DeepSignalPlant.py

-###-------- CALCULATION OF METHYLATION STATISTICS FROM TSV FILE GENERATED BY DeepSignal-Plant --------####
+####-------- CALCULATION OF METHYLATION STATISTICS FROM TSV FILE GENERATED BY DeepSignal-Plant --------####
 # Autor: Fadwa EL KHADDAR
 # Lab : DIADE - IRD
 # University : Montpellier - France
@@ -6,6 +6,8 @@
 import os
 import subprocess
 import numpy as np
+import pandas as pd
+import  glob
 import sys
 import matplotlib
 import matplotlib.pyplot as plt
@@ -14,8 +16,6 @@ import argparse
 import re
 from tkinter import *
 import tkinter as tk
-import glob
-import pandas as pd
 from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
 matplotlib.use('TkAgg')
 # Arguments
@@ -60,8 +60,8 @@ def CpG_context():
                 rows.append(row)
     rows_CG_df = pd.DataFrame(rows)
     rows_CG = rows_CG_df.sort_values("chromosome")
-    #output_file = os.path.join(output_dir, filebasename + "_CG.csv")
-    #rows_CG.to_csv(output_file, index=False, sep='\t')
+    output_file = os.path.join(output_dir, filebasename + "_CG.csv")
+    rows_CG.to_csv(output_file, index=False, sep='\t')
     #if not rows_CG.empty:
     #    for chromosome in rows_CG['chromosome'].unique():
     #        positions_par_chromosome[chromosome] = list(rows_CG.loc[rows_CG['chromosome'] == chromosome, ['strand', 'position']])
@@ -186,13 +186,14 @@ def plotting():
     CHG = CHG_context()
     CHH = CHH_context()
 
-
-    fig, ax = plt.subplots(figsize=(10,6), subplot_kw=dict(aspect="equal"))
+    fig, ax = plt.subplots(figsize=(8,6), subplot_kw=dict(aspect="equal"))
     nb_unmethylated = nb_C - CpG - CHG - CHH
     counts = [CpG, CHG, CHH, nb_unmethylated]
 
     labels = ["CpG", "CHG", "CHH", "unmethylated"]
 
+    fig = plt.figure()
+    ax = fig.add_subplot(111)
     wedges, texts, autotexts = ax.pie(counts, autopct='%1.1f%%')
     bbox_props= dict(boxstyle="square, pad=0.3", fc="w", ec="k", lw=0.72)
     kw = dict(arrowprops=dict(arrowstyle="-"), bbox=bbox_props, zorder=0, va="center")
@@ -207,21 +208,24 @@ def plotting():
         ax.annotate(f"{labels[i]}: {counts[i]}", xy=(x, y), xytext=(1.35 * np.sign(x), 1.4 * y),
                     horizontalalignment=horizontalalignment, **kw)
 
-    plt.title(f"Methylation context amoung all cytosines in the sample : {filebasename}")
+    #ax.pie(counts, autopct='%1.1f%%', startangle=90, center = (-2,0))
+    #legend = ax.legend(labels, loc = 2)
+    #ax.axis('equal')
 
+    plt.title(f"Methylation context amoung all cytosines in the sample : {filebasename}")
     root = tk.Tk()
-    root.title("Methylation Statistics")
+    root.title("Graphique")
 
-    # Create a canvas to display the figure
+    # Créer une zone pour afficher le graphique
     canvas = FigureCanvasTkAgg(fig, master=root)
     canvas.draw()
     canvas.get_tk_widget().pack()
 
-    # Create a button to close the window
+    # Bouton pour fermer la fenêtre
     button = tk.Button(master=root, text="Fermer", command=root.quit)
     button.pack()
 
-    # Start the tkinter event loop
+    # Démarrer la boucle principale Tkinter
     tk.mainloop()
 
 print("####-------- CALCULATION OF METHYLATION STATISTICS FROM TSV FILE GENERATED BY DeepSignal-Plant  --------####")
@@ -240,12 +244,19 @@ output_dir = args.o
 
 # Running seqtk
 result = subprocess.run(f"seqtk comp {fasta} | awk '{{C+=$4}}END{{print C}}'", shell=True, capture_output=True, text=True)
-nb_C = int(result.stdout.strip())
-print(f"Le nombre total de cytosine : {nb_C}")
+if result.stdout.strip():
+    nb_C = int(result.stdout.strip())
+    print(f"Le nombre total de cytosine : {nb_C}")
+else:
+    print("Erreur : la commande n'a pas retourné de résultat.")
 
 result = subprocess.run(f"seqtk comp {fasta} | awk '{{CpG+=$10}}END{{print CpG}}'", shell=True, capture_output=True, text=True)
-nb_CpG = int(result.stdout.strip())
-print(f"Le nombre total de CpG : {nb_CpG}")
+if result.stdout.strip():
+    nb_CpG = int(result.stdout.strip())
+    print(f"Le nombre total de CpG : {nb_CpG}")
+else:
+    print("Erreur : la commande n'a pas retourné de résultat.")
+