Workflow - Alignement et post-processing

École de bioinformatique INRAe-IFB-INSERM 2025

• Vivien DESHAIES - AP-HP

I - Alignement

BWA est un logiciel d'alignement dont la particularité est d'utiliser la transformée de Burrows-Wheeler (Burrows-Wheeler Aligner) ce qui lui permet d'utiliser peu de mémoire.

On commence par charger le module bwa. Comment trouver le module bwa ? (allez sur le terminal)

# Cette commande ne fonctionne pas ici il faut aller dans un terminal pour la lancer
#module avail -l bwa 

module load bwa/0.7.17

# Afficher l'aide de l'algorithme mem et se rendre dans le dossier tp_variant
bwa | cat

1. Indexation du génome

BWA a besoin de genome de ref indexé. pour cela on va utiliser la fonction index de bwa.

# Créer l'index du génome pour bwa
# bwa index # pour afficher l'aide
# usage: bwa index <fasta_file>
# remplacer <fasta_file> par le fasta du génome de référence.
bwa index genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa

2. Alignement des reads sur le génome de référence

BWA propose plusieurs algorithme d'alignement ... Celui qui est recommandé pour les données short-reads est mem (Maximum exact match)

Afficher d'abord l'aide de bwa mem

bwa mem | head

# usage: bwa mem -t 4 -R <readGroup> <fasta_file> <fastq1> <fastq2> > <sam_output>
# -R "@RG\tID:1\tPL:Illumina\tPU:PU\tLB:LB\tSM:SRR1262731"
bwa mem -t 4  \
    genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa  \
    fastq/SRR1262731_extract_R1.fq.gz  \
    fastq/SRR1262731_extract_R2.fq.gz  >  SRR1262731_extract.sam

Le caractère "\" permet de revenir à la ligne sans interrompre la commande. Il faut veiller à ne rien écrire après, y compris des espaces.

Quelle est la taille du fichier SRR1262731_extract.sam

ls -lh SRR1262731_extract.sam

2.1 Alignement long-reads

Pour les long-read l'outil d'alignement le plus courament utilisé est minimap2. Il utilise un sytème d'index et un algorithme mieux adapté pour gérer les homopolymères les taux d'erreurs plus élevés liés aux long-read.

L'outil propose des "preset" de paramètres en fonction de la technologie (PacBio, Nanopore ...) : https://lh3.github.io/minimap2/minimap2.html#8

# exemple
minimap2 -ax map-ont ref.fasta reads.fastq.gz > aln.sam

3. Trie et indexation de l'alignement (outil : samtools)

L'outil samtools permet de lire et de manipuler les fichiers d'alignement. Ses sous-commandes permettent de convertir, de calculer des métriques de couverture ...

# Chargement du module samtools en version 1.13
module load samtools/1.13

# Afficher l'aide
samtools | cat

a) Convertir le fichier sam (volumineux) en fichier bam (compressé).

samtools view -Sh -bo SRR1262731_extract.bam SRR1262731_extract.sam

Quelle est la taille du fichier SRR1262731_extract.bam

ls -lh SRR1262731_extract.bam SRR1262731_extract.sam

b) Ajout de la provenance des échantillons (readGroup)

• ReadGroups (RG) : permet d'associer des informations sur la provenance des reads
  • Identité : numéro de run ou d'échantillon
  • Séquençage, type de librairie, ...

Ces informations sont nécessaires à la recherche de variants

Plus d'informations sur les ReadGroups

Comment vérifier la présence de ReadGroups dans un fichier BAM ?

# Afficher l'entête du bam (-H)
samtools view -H SRR1262731_extract.bam

Comment pouvons-nous l'ajouter ?

# usage : samtools addreplacerg -r <readGroup> -o <output_bam> <input_bam>
samtools addreplacerg -@ 4 -r "@RG\tID:1\tPL:Illumina\tPU:PU\tLB:LB\tSM:SRR1262731" \
    -o SRR1262731_extract.rg.bam \
    SRR1262731_extract.bam

Vérifions en filtrant avec grep

samtools view -H SRR1262731_extract.rg.bam | grep '^@RG'

b) Trie des reads alignés par coordonnées génomique

# usage: samtools sort -@ nb_cpu -o <sorted_bam>  <input_bam>
samtools sort -@ 4 -o SRR1262731_extract.sort.bam SRR1262731_extract.rg.bam

c) Créaction de l'index

# usage: samtools index <sorted_bam>
samtools index SRR1262731_extract.sort.bam

Qu'est-ce qui s'est passé ?

ls -l 

d) Visualiser le contenu du BAM

samtools view -h SRR1262731_extract.sort.bam | head -n6

e) Contrôle qualité des données alignées

Quelles informations regarder une fois le mapping effectué ?

• Pourcentage total de reads alignés
• Pourcentage de reads pairés “proprement”

Quels outils ?

• Samtools flagstat/stat
• Qualimap [optionnel]
• multiqc [optionnel]

Samtools :

module load samtools/1.13
samtools flagstat | cat # affiche l'aide
samtools stat | cat # affiche l'aide

# usage: samtools flagstat <bam_file>
samtools flagstat SRR1262731_extract.sort.bam

Stocker le résultat dans un fichier txt:

samtools flagstat SRR1262731_extract.sort.bam > SRR1262731_extract.sort.flagstat.txt

samtools stats SRR1262731_extract.sort.bam > SRR1262731_extract.sort.stats.txt

Qualimap :

module load qualimap/2.2.2b
qualimap bamqc | cat # affiche l’aide

qualimap bamqc -nt 4 -outdir SRR1262731_extract_qualimap_report --java-mem-size=4G -bam SRR1262731_extract.sort.bam

MultiQC:

module load multiqc/1.11
multiqc -h | cat # affiche l'aide

mkdir -p log
multiqc -f . 2>log/multiqc_stdout.txt

# Création d'un archive zip
zip -r SRR1262731_multiqc.zip multiqc_report.html multiqc_data

Télécharger le zip pour visualiser le rapport sûr votre machine

f) Optimisation de l'espace

Il est recommandé de supprimer tous les fichiers intermediares qui sont volumineux.

Qu'est-ce qu'il est possible de supprimer ?

ls -lh SRR1262731_extract*

rm SRR1262731_extract.sam
rm SRR1262731_extract.bam
rm SRR1262731_extract.rg.bam

II - Processing post alignement

Il est nécessaire de préparer les données avant de pouvoir détecter les variants

Les étapes que nous allons voir durant le TP :

Rappel du workflow :

1. Marquage des duplicats de PCR

Identifier les reads provenant d'une même molécule issus de :

• PCR duplicates : amplification PCR durant la préparation de la librairie
• Optical duplicates : Cluster Illumina identifié comme deux clusters

• Garder les duplicats : probabilité importante de confondre les duplicats avec des fragments biologiques issus du même locus
• Marquer les duplicats mais les conserver dans le fichier BAM
• Supprimer les duplicats du fichier BAM : Certains outils les supprimeront par défaut (samtools, GATK, ...)

Avec l'outil MarkDuplicates de la suite PicardTools intégrée à la suite GATK4

# Charger le module gatk4 en 4.2.3.0
module load gatk4/4.2.3.0

# Afficher l'aide de gatk MarkDuplicates
gatk MarkDuplicates --help | head 

# Lancement de l'outils
gatk MarkDuplicates --java-options '-Xmx8G' \
    --INPUT SRR1262731_extract.sort.bam \
    --OUTPUT SRR1262731_extract.sort.md.bam \
    --METRICS_FILE SRR1262731_extract_metrics_md.txt \
    --VALIDATION_STRINGENCY SILENT  2>log/MarkDuplicates.log

# Afficher les 20 premières lignes du fichier
head SRR1262731_extract_metrics_md.txt

# Rechercher le % de duplicat de PCR
grep -A1 "LIBRARY" SRR1262731_extract_metrics_md.txt

2. Filtres sur les alignements

Restreindre le fichier BAM en fonction de métriques d'alignements :

• Qualité de mapping (MAPQ) suffisante
• Retrait des reads non mappés

# Suppression des reads non mappés et filtre sur les reads avec MAPQ < 30
samtools view -bh -F 4 -q 30 SRR1262731_extract.sort.md.bam \
    > SRR1262731_extract.sort.md.filt.bam

Pour utiliser le paramètre -F : plus d'information sur les SAM flags

# Résumé de l'ensemble des flag présents dans notre BAM
samtools flagstat SRR1262731_extract.sort.md.filt.bam \
    > SRR1262731_extract.filt.flagstat.txt

# Visualiser le fichier
cat SRR1262731_extract.filt.flagstat.txt

Restreindre le fichier BAM en fonction de métriques d'alignements (suite) :

• Alignements intersectant les régions d'intérêt
• En fonction du nombre de mismatchs, de la taille d'insert, de paires mappées sur des chromosomes différents, ...

# Conservation des alignements dans les régions ciblées
module load bedtools/2.29.2

# Afficher l'aide
bedtools intersect --help | head

# Lancement de bedtools intersect
bedtools intersect \
    -a SRR1262731_extract.sort.md.filt.bam \
    -b additionnal_data/QTL_BT6.bed \
    > SRR1262731_extract.sort.md.filt.onTarget.bam

# Indexation du bam
samtools index SRR1262731_extract.sort.md.filt.onTarget.bam

3. Analyse de la couverture

Contrôle qualité de l'enrichissement de ma capture :

• Est ce que ma région est couverte par suffisamment de reads ?
• Cette couverture est-elle homogène sur toute la région ?

# Afficher l'aide de la fonction Samtools
samtools depth --help | head

# Lancement de la commande
samtools depth -b additionnal_data/QTL_BT6.bed \
    SRR1262731_extract.sort.md.filt.onTarget.bam \
    > SRR1262731_extract.onTarget.depth.txt

head SRR1262731_extract.onTarget.depth.txt

# Compter les positions avec une profondeur inférieure à 3
awk '{if($3<3)print}' SRR1262731_extract.onTarget.depth.txt | wc -l

4. Contrôle qualité des données alignées après filtres

A vous de jouer