Statistics_methylation.py

import pandas as pd
import  glob

# Utilisation de glob pour pouvoir tourner le script sur tous les fichiers tsv.

path = r"/home/fadwa/Téléchargements/test"
filenames = glob.glob(path + '/*.tsv')

# Définition des fonctions

def CpG_contexte(methylated):

    """
    Cette fonction a pour but d'identifier la présence de la méthylation type CpG
    et de calculer le pourcentage

    """
    global CpG
    motif_counts = {}  # Création d'un dictionnaire pour stocker le comptage de chaque motif
    for motif in CpG:
        motif_counts[motif] = 0

    # Boucler à travers chaque k_mer dans df dataset

    for kmer in methylated["k_mer"][2:]:
        sub_kmer = kmer[2:4]
        sub_kmer_str = "".join(sub_kmer)

        # Comptage des occurences de chaque motif dans sub_kmer

        for motif in CpG:
            motif_count = sub_kmer_str.count(motif)
            motif_counts[motif] += motif_count

    # Calculer le total des occurences de chaque sub_kmer
    total_len = len(methylated["k_mer"])
    total_count = sum(motif_counts.values())
    pourcentage = (total_count / total_len) * 100

    print(f"Le pourcentage de la méthylation type CpG est  : {pourcentage: .2f}% ")



def CHG_context(methylated):
    """
    Cette fonction a pour but d'identifier la présence de la méthylation type CHG
    et de calculer le pourcentage

    """
    motif_counts = {}  # Création d'un dictionnaire pour stocker le comptage de chaque motif
    global CHG
    for motif in CHG:
        motif_counts[motif] = 0

    # Boucler à travers chaque k_mer dans methylated dataset

    for kmer in methylated["k_mer"][2:]:
        sub_kmer = kmer[2:]
        sub_kmer_str = "".join(sub_kmer)

        # Comptage des occurences de chaque motif dans sub_kmer

        for motif in CHG:
            motif_count = sub_kmer_str.count(motif)
            motif_counts[motif] += motif_count

    # Calculer le total des occurences de chaque sub_kmer
    total_len = len(methylated["k_mer"])
    total_count = sum(motif_counts.values())
    pourcentage = (total_count  / total_len) * 100

    print(f"Le pourcentage de la méthylation type CHG est: {pourcentage: .2f}% ")


def CHH_context(methylated):
    """
    Cette fonction a pour but d'identifier la présence de la méthylation type CpG
    et de calculer le pourcentage

    """
    global CHH
    motif_counts = {}  # Création d'un dictionnaire pour stocker le comptage de chaque motif
    for motif in CHH:
        motif_counts[motif] = 0

    # Boucler à travers chaque k_mer dans methylated dataset

    for kmer in methylated["k_mer"][2:]:
        sub_kmer = kmer[2:]
        sub_kmer_str = "".join(sub_kmer)

        # Comptage des occurences de chaque motif dans sub_kmer

        for motif in CHH:
            motif_count = sub_kmer_str.count(motif)
            motif_counts[motif] += motif_count

    # Calculer le total des occurences de chaque sub_kmer
    total_len = len(methylated["k_mer"])
    total_count = sum(motif_counts.values())
    pourcentage = (total_count  / total_len) * 100

    print(f"Le pourcentage de la méthylation type CHH est : {pourcentage: .2f}% ")


# Déclaration des trois contexte de méthylation
CpG = ["CG"]
CHG = ["CCG", "CAG", "CTG"]
CHH = ["CCC", "CAA", "CTT", "CCA", "CCT", "CAT", "CAC", "CTA", "CTC"]
seuil = 5



# Lecture de fichier
for file in filenames:
    print(f"Statistiques pour le fichier {file}")
    names = ["Chromosome", "position","strand", "pos_in_strand", "prob_0_sum", "prob_1_sum", "count_modified", "count_unmodified", "coverage", "modification_frequency", "k_mer"]
    call_mods = pd.read_csv(file, sep='\t', names=names)

# Création du DataFrame
df = pd.DataFrame(call_mods)

# Extraction des reads méthylés

methyl= df[df["count_modified"] >= seuil]
methylated=pd.DataFrame(methyl)

#methylated.to_csv("/home/fadwa/Téléchargements/FAS08151_methylated.tsv", sep='\t')

CpG_contexte(methylated)
CHG_context(methylated)
CHH_context(methylated)