Filtrage des reads selon le phred quality d'une position génomique

Cet outil filtre les reads provenant d'un fichier BAM en fonction des scores de qualité Phred à des positions génomiques spécifiques. Il est conçu pour s'assurer que les reads conservés dans le fichier de sortie respectent un seuil de qualité spécifié à deux positions distinctes.

Cet outil est créé par mes soins

Utilisation

Pour utiliser ce script, vous devez avoir Python et la bibliothèque pysam installés. Vous pouvez exécuter le script depuis la ligne de commande en utilisant la syntaxe suivante :

python bam_quality_filter.py <chemin_bam_entree> <chemin_bam_sortie> <chromosome> <position1> <position2> <qualite_min> <chemin_fichier_log>

Paramètres

• chemin_bam_entree : Chemin vers le fichier BAM d'entrée.

• chemin_bam_sortie : Chemin vers le fichier BAM de sortie où les reads filtrés seront sauvegardés.

• chromosome : Chromosome à partir duquel les reads sont analysés (ex. : chr1).

• position1 : Première position génomique (base 1) pour vérifier la qualité du read.

• position2 : Seconde position génomique (base 1) pour vérifier la qualité du read.

• qualite_min : Qualité Phred minimale requise aux deux positions pour conserver un read.

• chemin_fichier_log : Chemin vers le fichier log où les résultats du filtrage seront écrits.

Sortie

Le script créera un fichier BAM de sortie avec les reads ayant un score de qualité Phred supérieur au minimum spécifié aux deux positions indiquées. Il génère également un fichier log détaillant le statut de traitement de chaque read.

Dépendances

• Python 3.x

• pysam

Installez pysam via pip :

pip install pysam

Exemple

Exécutez le script avec la commande suivante :

python bam_quality_filter.py exemple.bam filtre.bam chr1 100000 100500 30 log_qualite.log

Cette commande filtre les reads en fonction de leur qualité aux positions 100000 et 100500 sur le chromosome 1, en exigeant une qualité Phred minimale de 30.

Le programme python

import pysam
import sys

input_bam_path = sys.argv[1]
output_bam_path = sys.argv[2]
chromosome = sys.argv[3]
position1 = int(sys.argv[4])  # Première position génomique, base 1
position2 = int(sys.argv[5])  # Seconde position génomique, base 1
min_quality = int(sys.argv[6])  # Qualité Phred minimale pour conserver un read
log_file_path = sys.argv[7]

input_bam = pysam.AlignmentFile(input_bam_path, "rb")
output_bam = pysam.AlignmentFile(output_bam_path, "wb", template=input_bam)

with open(log_file_path, 'w') as log_file:
    log_file.write("ReadName\tBaseAtPosition1\tQuality1\tBaseAtPosition2\tQuality2\tStatus\n")

    # Parcourir tous les reads qui couvrent au moins une des positions
    start = min(position1, position2) - 1
    end = max(position1, position2)
    for read in input_bam.fetch(chromosome, start, end):
        if not read.is_unmapped:
            qualities = []
            bases = ['N', 'N']
            positions = [position1, position2]
            status = "Kept"

            for i, pos in enumerate(positions):
                read_pos = None
                for qpos, rpos in read.get_aligned_pairs():
                    if rpos == pos - 1:
                        read_pos = qpos
                        break

                if read_pos is not None:
                    if read.query_qualities is not None and 0 <= read_pos < len(read.query_qualities):
                        qualities.append(read.query_qualities[read_pos])
                    else:
                        qualities.append(None)

                    if read.query_sequence is not None and 0 <= read_pos < len(read.query_sequence):
                        bases[i] = read.query_sequence[read_pos]
                else:
                    qualities.append(None)

            if any(q is not None and q < min_quality for q in qualities):
                status = "Removed"
            else:
                output_bam.write(read)

            log_file.write(f"{read.query_name}\t{bases[0]}\t{qualities[0]}\t{bases[1]}\t{qualities[1]}\t{status}\n")

input_bam.close()
output_bam.close()