1 of 31

Atelier Variant

​

Nadia Bessoltane - INRAE

Vivien Deshaies - AP-HP

​

​

École de bioinformatique AVIESAN-IFB-INSERM 2022

2 of 31

Programme de l’atelier Variants

  • Détection des petites variations génomiques
  • Détection des variations structurales
  • Manipulation des variants avec R
  • Ecrire un script automatique

2

3 of 31

Introduction

4 of 31

https://www.genome.gov/Health/Genomics-and-Medicine/Polygenic-risk-scores#one

Définition

Il existe différents types de variations :

  • SNV : Single Nucleotide Variant
  • INDEL : INsertion ou DELetion
  • SV (Structural Variant)
  • CNV (Copy Number Variation)

5 of 31

Variant : variation génomique dans une séquence nucléotidique, en comparaison avec une séquence de référence

  • SNV : Single Nucleotide Variant

​

  • INDEL : INsertion ou DELetion d’une ou plusieurs bases

​

  • MNV (Multi-Nucleotide Variant) : plusieurs SNVs et/ou INDELS dans un bloc

​

  • SV (Structural Variant) : réarrangement génomique affectant > 50bp

5

Définition

6 of 31

  • SNV (Single Nucleotide Variant)

→ toute altération nucléotidique sans implication de fréquence populationnelle

​

  • SNP (Single Nucleotide Polymorphism)

→ implique qu’un variant est partagée dans la population (> 1%)

​

/!\ l’amalgame SNPs est souvent fait pour qualifier les SNVs /!\

6

SNV ≠ SNP

7 of 31

Workflow

7

Fastq Quality Control

-- FastQC --

Mapping

-- BWA --

Reads (Fastq)

Reference genome (Fasta)

Processing Post Alignement

-- GATK/samtools --

Small variant calling

Annotation

Structural variant calling

Variant analysis

Shell

R/Rstudio

8 of 31

Détection des petites variations génomiques

​

Vivien Deshaies - AP-HP

9 of 31

Depuis que l’homme fait de l’élevage, il essaie de faire en sorte de toujours améliorer sa production, que ce soit en quantité ou en qualité.

​

​

Les technologies de génotypage permettent maintenant de sélectionner les mâles reproducteurs en fonction du fond génétique qu’ils vont pouvoir transmettre à leur descendance.

​

​

​

​

9

Chez le bovin, il existe un locus de caractères quantitatifs (QTL) lié à la production de lait, situé sur le chromosome 6, et plus exactement sur une région de 700 kb, composée de 7 gènes.

Jeux de données #1 : SNVs/Indels

10 of 31

Les échantillons QTL+ sont caractérisés par une diminution de la production en lait et une augmentation des concentrations en protéine et lipide.

Vous aurez à votre disposition :

  • Un extrait des données de séquences d’un échantillon du projet 1000 génomes bovins, phénotypé comme QTL- : SRR1262731
  • Les résultats du variant calling pour deux échantillons phénotypés QTL+ : SRR1205992 et SRR1205973

​

Your turn !

Quelle mutation est responsable de ce QTL ?

​

​

​

​

10

Miri Cohen-Zinder et al, Genome Research 2005, DOI:10.1101/gr.3806705

Jeux de données #1 : SNVs/Indels

11 of 31

Copie du jeu de données #1

11

# Copie des fichiers dans notre home

$ mkdir -p ~/tp_variant

$ cp -r /shared/projects/form_2022_32/atelier_variant/variants/* ~/tp_variant/

$ ls -l

​

​

12 of 31

Détection des petites variations génomiques

Alignement et post-processing

13 of 31

Workflow - Alignement & Processing Post Alignement

  • Nécessité de préparer les données avant la détection des variants

13

Fastq Quality Control

-- FastQC --

Mapping

-- BWA --

Reads (Fastq)

Reference genome (Fasta)

Marquage des duplicats

-- MarkDuplicates --

Filtres sur l’alignement

-- samtools / bedtools --

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

-- Samtools / qualimap / R --

14 of 31

Indexation du génome

14

$ module load bwa/0.7.17

$ module load samtools/1.13

$ module load gatk4/4.2.3.0

​

​

$ # se déplacer dans le dossier genome

$ cd ~/tp_variant/genome/

$

$ bwa index Bos_taurus.UMD3.1.dna.toplevel.6.fa

$

$ samtools faidx Bos_taurus.UMD3.1.dna.toplevel.6.fa

$

$ gatk CreateSequenceDictionary --REFERENCE Bos_taurus.UMD3.1.dna.toplevel.6.fa --OUTPUT Bos_taurus.UMD3.1.dna.toplevel.6.dict

$

$ ls -l

15 of 31

Alignement des données avec l’outil BWA-mem

Mapping

15

$ bwa mem # affiche l’aide de l’algorithme mem

$ cd ~/tp_variant/

​

$ # Exécuter l’alignement (bwa mem -t 4 -R <readGroup> genome fastq1 fastq2 > sam)

$ bwa mem -t 4 -R "@RG\tID:1\tPL:Illumina\tPU:PU\tLB:LB\tSM:SRR1262731" \

genome/Bos_taurus.UMD3.1.dna.toplevel.6.fa \

fastq/SRR1262731_extract_R1.fq.gz \

fastq/SRR1262731_extract_R2.fq.gz > SRR1262731_extract.sam

​

​

​

​

16 of 31

Alignement des données avec l’outil BWA-mem

Trier et indexer l’alignement

16

$ module load samtools/1.13

$

$ # convertir le sam en bam

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

$

$ # On trie le fichier BAM par coordonnées

$ samtools sort -@ 4 -o SRR1262731_extract.sort.bam SRR1262731_extract.bam

$ # et on crée un index (.bai)

$ samtools index SRR1262731_extract.sort.bam

$ # Visualiser le contenu du BAM

$ samtools view -h SRR1262731_extract.bam | less -S

$ # supprimer le sam pour libérer de l’espace

$ rm SRR1262731_extract.sam

​

​

​

​

​

​

​

​

17 of 31

Ajout de la provenance des échantillons

  • ReadGRoups (RG) : associe des informations sur la provenance des reads

→ Identité : run/échantillon

→ Séquençage, librairie…

  • Nécessaire à la recherche de variants
  • Plus d’informations

​

17

$ samtools view -H SRR1262731_extract.bam | grep "^@RG"

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

$ module load gatk4/4.2.3.0

$ gatk AddOrReplaceReadGroups --help

  • Avec l’outil AddOrReplaceReadGroups de la suite PicardTools intégrée à GATK4

18 of 31

Workflow - Processing Post Alignement

  • Nécessité de préparer les données avant la détection des variants

18

19 of 31

Workflow - Alignement & Processing Post Alignement

  • Nécessité de préparer les données avant la détection des variants

19

Fastq Quality Control

-- FastQC --

Mapping

-- BWA --

Reads (Fastq)

Reference genome (Fasta)

Marquage des duplicats

-- MarkDuplicates --

Filtres sur l’alignement

-- samtools / bedtools --

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

-- Samtools / qualimap / R --

20 of 31

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

20

21 of 31

Marquage des duplicats de PCRa

→ 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

21

$ module load gatk4

$ gatk MarkDuplicates --help # affiche l’aide

​

$ gatk MarkDuplicates --java-options '-Xmx8G' \

-I SRR1262731_extract.sort.bam --VALIDATION_STRINGENCY SILENT \

-O SRR1262731_extract.sort.md.bam -M SRR1262731_extract_metrics_md.txt

22 of 31

Marquage des duplicats de PCR

→ 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…)

22

​

$ more SRR1262731_extract_metrics_md.txt

$ # % de pcrDup

$ grep -A1 "LIBRARY" SRR1262731_extract_metrics_md.txt

​

23 of 31

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

Pour utiliser le paramètre -F : plus d’information sur les SAM Flags

23

# 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

$ # https://broadinstitute.github.io/picard/explain-flags.html

​

$ samtools flagstat SRR1262731_extract.sort.md.filt.bam \

> SRR1262731_extract.filt.flagstat.txt

​

$ cat SRR1262731_extract.filt.flagstat.txt

24 of 31

Filtres sur les alignements

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

  • 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…

24

# Conservation des alignements dans les régions ciblées

$ module load bedtools/2.29.2

$ bedtools --version # affiche la version (v2.29.2)

$ bedtools intersect --help # affiche l’aide

​

$ bedtools intersect -a SRR1262731_extract.sort.md.filt.bam \

-b ~/tp_variant/additionnal_data/QTL_BT6.bed \

> SRR1262731_extract.sort.md.filt.onTarget.bam

$ samtools index SRR1262731_extract.sort.md.filt.onTarget.bam

25 of 31

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

  • Quelles informations regarder une fois l’alignement effectué ?

​

​

25

  • Quels outils ?
    • Samtools flagstat
    • Qualimap [optionnel]
    • MultiQC

→ Pourcentage total de reads alignés

→ Pourcentage de reads appariés “correctement”

26 of 31

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

26

# Lancement de samtools

$ samtools flagstat # affiche l'aide

​

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

​

$ cat SRR1262731.flagstat.txt # visualisation du résultat

​

$ samtools stats # affiche l'aide

​

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

​

$ cat SRR1262731.stats.txt # visualisation du résultat

27 of 31

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

​

  • Visualisation des contrôles qualité

​

27

# Lancement de Multiqc

$ module load multiqc/1.11

$ multiqc -h # affiche l’aide

​

$ multiqc -f .

​

# Téléchargement du fichier html à partir de jupyterhub

28 of 31

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

28

# Lancement de Qualimap

$ module load qualimap/2.2.2b

$ qualimap -h # affiche les outils disponibles (+ version)

$ qualimap bamqc # affiche l’aide

​

$ qualimap bamqc -nt 4 -outdir SRR1262731_extract_qualimap_report \

--java-mem-size=4G -bam SRR1262731_extract.sort.bam

​

# Création d’une archive zip

$ zip -r SRR1262731_extract_qualimap_report.zip \

SRR1262731_extract_qualimap_report

​

# Téléchargement du fichier zip à partir de jupyterhub

# Generation d’un rapport multiqc

$ multiqc -f .

29 of 31

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 ?

29

30 of 31

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 ?

30

# Calcul de la couverture avec samtools

$ samtools depth --help # affiche l’aide

​

$ samtools depth -b ~/tp_variant/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 position avec une profondeur inférieure à 3

$ awk '{if($3<3)print}' SRR1262731_extract.onTarget.depth.txt | wc -l

31 of 31

  • Jeux de données #1 : SNVs/Indels

→ /shared/projects/form_2022_32/atelier_variant/variants

​

​

Cheatsheet :

→ Version html : /shared/projects/form_2022_32/atelier_variant/EBAII2021_variants.html

31

Emplacement des données brutes