# 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.

{% hint style="info" %}
Cet outil est créé par mes soins 
{% endhint %}

### 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 :

```bash
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 :

```sh
pip install pysam
```

### Exemple

Exécutez le script avec la commande suivante :

```sh
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

```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()

```


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://aleksandre80.gitbook.io/stage/filtrage-des-reads-selon-le-phred-quality-dune-position-genomique.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.