Workflow - Alignement et post-processing¶
École de bioinformatique IFB-INSERM 2024
- 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 -l SRR1262731_extract.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 -l SRR1262731_extract.bam
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
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 # affiche l'aide
# samtools stat # 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 # 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 : Il faut mettre à jour ce shema!!!!!!!
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 :)