Atelier Variant

Introduction

Nadia Bessoltane - INRAE

Elodie Girard - Institut Curie

(Maria Bernard - INRAE)

(Olivier Rué - INRAE)

​

École de bioinformatique AVIESAN-IFB-INSERM 2021

Olivier Quenez - INSERM

Mathieu Charles - INRAE

Odile Rogier - INRAE

​

• Introduction, qualité de lectures et alignement (Elodie)
• Utilisation du visualiseur IGV (Olivier Q.)
• Pré- processing des alignements (Odile)
• SNVs et indels de petite taille, à l’aide de 3 outils : GATK, Mpileup/Varscan et discoSNP (Nadia, Mathieu)
• Introduction à R (Nadia)
• Variations Structurales (SV) (Olivier Q.)
• Utilisation de R pour visualiser des métriques obtenus (Elodie)
• Workflow/Conclusion (Odile)

2

Objectifs�

JupyterHub @ IFB

https://jupyterhub.cluster.france-bioinformatique.fr

​

​

​

3

doc cluster

Cluster de l’IFB

L’Institut Français de Bioinformatique met à disposition de la communauté un cluster de calculs

Your turn! Se connecter au cluster

​

​

​

4

# Sous Windows avec MobaXterm

​

Session : ssh

​

Host : core.cluster.france-bioinformatique.fr

​

Specify username : coché et complété

# Sous Mac avec Cyberduck

​

Open connexion : SFTP

Server : core.cluster.france-bioinformatique.fr

Username/Password : à compléter

doc cluster

Cluster de l’IFB

• /!\ Connexion initiale : tout le monde sur le noeud maître sur lequel il ne faut pas travailler /!\
• Lancement de “jobs” ou d’une session interactive sur le cluster
• [Vidéo] : The 5 minutes IFB Core Cluster tutorial

​

5

# Chargement de l’environnement dédié à chaque outil (exemple pour varscan)

$ module avail -l | grep varscan

$ module load varscan/2.4.3 # ou module load varscan

​

# Nous aurons besoin au cours du TP de ressources CPU et mémoire (RAM)

$ sbatch --cpus=4 --mem=16G -J toolName_<user_name> --wrap="tool command line"

​

# Pour suivre vos “jobs” soumis sur le cluster, 2 solutions

​

​

Cluster de l’IFB

Remember : Tous les jobs doivent être lancés sur un noeud du cluster !

$ squeue -u <user_name>

$ scontrol show job <job_id>

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.

​

​

​

​

6

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

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 ?

​

​

​

​

7

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

Jeux de données #1 : SNVs/Indels

Zymoseptoria tritici : Champignon ascomycète, pathogène du blé tendre, responsable d’une maladie foliaire (septoriose).

• Principale maladie du blé (jusqu’à 50% de perte de rendement).
• Haploïde, génome de 40 Mb séquencé en 2011 : 13 chromosomes essentiels + 8 chromosomes accessoires
• Souche séquencée avec deux technologies : Illumina et MinIon

​

8

Your turn !

Retrouvez les délétions de grande taille

Jeux de données #2 : SVs

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

→ /shared/projects/form_2021_26/data/atelier_variant/variants

• Jeux de données #2 : SVs

→ /shared/projects/form_2021_26/data/atelier_variant/sv

​

Cheatsheet :

→ Version html : /shared/projects/form_2021_26/data/atelier_variant/EBAII2021_variants.html

9

Emplacement des données brutes

Workflow

10

Fastq Quality Control

-- FastQC --

Mapping

-- Bowtie2/Bwa --

Reads (Fastq)

Reference genome (Fasta)

Processing Post Alignement

Variant Calling

Filtrage & Annotation

Variants Structuraux

Copy Number

Copie du jeu de données #1

11

# Listing des fichiers FASTQ, Genome et BAM

$ ls -lh /shared/projects/form_2021_26/data/atelier_variant/variants/fastq

$ ls -lh /shared/projects/form_2021_26/data/atelier_variant/variants/genome

​

​

# Copie des fichiers dans notre home

$ mkdir -p ~/tp_variant

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

​

# Se déplacer dans le dossier optional

$ mkdir -p ~/tp_variant/optional

$ cd ~/tp_variant/optional

​

​

​

​

​

Les lectures (raw reads) au format fastq

​

12

# Vous pouvez utiliser la commande less pour visualiser le contenu du fichier

# q pour quitter

​

$ less -S ~/tp_variant/fastq/SRR1262731_extract_R1.fq.gz

​

​

​

​

​

​

Header

Sequence

+ (optional header)

Quality

@QSEQ32.249996 HWUSI-EAS1691:3:1:17036:13000#0/1 PF=0 length=36

GGGGGTCATCATCATTTGATCTGGGAAAGGCTACTG

+

=.+5:<<<<>AA?0A>;A*A################

Le score de qualité Sanger

​

• Une valeur de score Sanger est attribuée à chaque base séquencée
• Basée sur p, la probabilité d’erreur (i.e. que la base soit fausse)
• Q_Sanger = -10*log₁₀(p)
• p = 0.1 ⇔ Q_Sanger10
• p = 0.01 ⇔ Q_Sanger20
• p = 0.001 ⇔ Q_Sanger30
• …
• Les scores sont encodés en ASCII 33
• Objectif : compresser les données en diminuant le nombre de caractères utilisés pour encoder la qualité
• Le score de qualité Sanger varie entre 0 et 40

​

​

13

Le score de qualité Sanger

​

• ! correspond à 0
• “ correspond à 1
• # correspond à 2
• $ correspond à 3
• …
• I correspond à 40

​

14

Contrôle qualité des données brutes

15

$ module load fastqc/0.11.9

$ fastqc --version # affiche la version (v0.11.9)

$ fastqc --help # affiche l’aide

​

$ mkdir -p Fastqc/logs

$ cd ~/tp_variant/optional/Fastqc

​

$ sbatch -J FastQC_SRR1262731_R1 -o logs/FastQC_SRR1262731_R1.out -e logs/FastQC_SRR1262731_R1.err --cpus-per-task=2 --wrap=" \

fastqc --threads 2 --outdir . ~/tp_variant/fastq/SRR1262731_extract_R1.fq.gz"

​

$ sbatch -J FastQC_SRR1262731_R2 -o logs/FastQC_SRR1262731_R2.out -e logs/FastQC_SRR1262731_R2.err --cpus-per-task=2 --wrap=" \

fastqc --threads 2 --outdir . ~/tp_variant/fastq/SRR1262731_extract_R2.fq.gz"

​

# ouvrir les fichiers html via jupyter

​

​

​

​

​

​

Trimmer les lectures

• Une étape de pré-processing
• Les reads en entrée sont rognés afin d’éliminer des extrémités de mauvaises qualités
• En fonction de la capacité de l’outil à faire des alignements locaux ou globaux et de la qualité intrinsèque des données, cette étape peut être cruciale
• Risque: peu de reads alignés
• Quelques logiciels existants
• Sickle-trim (sliding window-based trimming)
• FASTX-Toolkit (cut a defined number of nucleotides)
• Trimmomatic
• Cutadapt

16

Principe de cutadapt

• Objectif :
• Supprimer les extrémités de mauvaise qualité
• Solution:
• Parcourir le read avec un fenêtre coulissante de droite à gauche. Calculer la qualité moyenne dans chaque fenêtre
• Si la valeur de qualité chute en dessous d’une valeur seuil q, déléter l’extrémité 3’.
• Si la taille restante du read est inférieure à une longueur seuil l, déléter le read.

​

​

17

ACTCGCTCGCTGGTTAATCGATGATCGTGCAGTCGTACTCGTAGCTAGCTAGTCGTAACATAGCTAGTC

Retrait des séquences de mauvaises qualité

18

$ module load cutadapt/2.10

$ cutadapt --version # affiche la version (v0.2.10)

$ cutadapt --help # affiche l’aide

​

$ mkdir -p ~/tp_variant/optional/Cutadapt/logs

$ cd ~/tp_variant/optional/Cutadapt

​

$ sbatch -J Cutadapt_SRR1262731 -o logs/Cutadapt_SRR1262731.out -e logs/Cutadapt_SRR1262731.err --cpus-per-task=2 --wrap=" \

cutadapt --cores 2 --trim-n --max-n 0.3 --error-rate 0.1 -q 30,30 \

--minimum-length 50 --pair-filter both \

--paired-output SRR1262731_extract_R2.trimmed.fq \

--output SRR1262731_extract_R1.trimmed.fq \

~/tp_variant/fastq/SRR1262731_extract_R1.fq.gz \

~/tp_variant/fastq/SRR1262731_extract_R2.fq.gz \

> SRR1262731_extract_trimming_stats.txt"

​

​

​

​

​

​

​

​

Workflow

19

Fastq Quality Control

-- FastQC --

Mapping

-- Bowtie2/Bwa --

Reads (Fastq)

Reference genome (Fasta)

Processing Post Alignement

Variant Calling

Filtrage & Annotation

Variants Structuraux

Copy Number

Aligner les reads

• Objectif
• Trouver la région du génome qui a produit les read
• Trouver dans le génome le mot correspondant au read

​

​

20

Ref. Genome

Reads

• Une position dans le génome
• Plusieurs positions (éléments répétés, région de basse complexité)

L’approche seed & extend

• Une extrémité du read est interrogée (la graine = the seed)
• On cherche ses régions correspondantes sur le génome (à l’aide d’un index créé initialement) avec ou sans mismatch

/!\ Ne faire qu’une seule fois par génome d’intérêt et version majeure /!\

• On teste si le reste du read s’aligne avec la séquence
• Les données générées sont au format SAM ou BAM

​

​

21

✅

❌

❌

❌

L’ajout de Read Group

• Associe une identification de provenance à chaque read
• Utile dans les analyses multi-échantillons ou de reséquençage
• Obligatoire pour utiliser certains outils (comme GATK)
• Un Read Group (RG) est défini par :
• ID : Read group IDentifier (barcode)
• PU : Platform Unit
• SM : Sample Biological NaMe
• PL : PLatform/Technology utilisée (e.g.: Illumina)
• LB : préparation de la LiBrary

​

22

Indexation du génome pour BWA

23

$ module load bwa/0.7.17

$ module load samtools/1.13

$ module load gatk4/4.2.3.0

​

​

​

$ cd ~/tp_variant/genome/

$ mkdir -p logs

​

$ sbatch -J BWA_index -o logs/BWA_index.out -e logs/BWA_index.err --wrap="bwa index Bos_taurus.UMD3.1.dna.toplevel.6.fa"

​

$ sbatch -J samtools_index -o logs/samtools_index.out -e logs/samtools_index.err --wrap="samtools faidx Bos_taurus.UMD3.1.dna.toplevel.6.fa"

​

$ sbatch -J GATK_index -o logs/GATK_index.out -e logs/GATK_index.err --wrap=" \

gatk CreateSequenceDictionary --REFERENCE Bos_taurus.UMD3.1.dna.toplevel.6.fa \

--OUTPUT Bos_taurus.UMD3.1.dna.toplevel.6.dict"

​

​

​

​

​

​

Alignement des données

24

$ bwa # affiche la version et l’aide (v0.7.17-r1188)

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

​

$ cd ~/tp_variant/optional/

$ mkdir -p alignment_bwa/logs

$ cd alignment_bwa

​

$ sbatch -J SRR1262731_mapping -o logs/SRR1262731_mapping.out -e logs/SRR1262731_mapping.err --cpus-per-task=4 --mem=16G --wrap=" \

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

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

~/tp_variant/fastq/SRR1262731_extract_R1.fq.gz \

~/tp_variant/fastq/SRR1262731_extract_R2.fq.gz \

| samtools view -Sh - -bo SRR1262731_extract.bam"

​

# Visualiser le contenu du BAM

​

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

​

​

​

​

​

​

​

​

Visualiser le contenu du BAM

25

# Visualiser le contenu du BAM

​

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

​

Tri et indexage du BAM

26

# On trie le fichier BAM par coordonnées et on crée un index (.bai)

$ sbatch -J SRR1262731_mappingSort -o logs/SRR1262731_mappingSort.out -e logs/SRR1262731_mappingSort.err --cpus-per-task=4 --mem=16G --wrap=" \

samtools sort -@ 4 --write-index \

-o SRR1262731_extract.sort.bam##idx##SRR1262731_extract.sort.bam.bai \

SRR1262731_extract.bam"

​

samtools index SRR1262731_extract.sort.bam

​

# On produit les statistiques d'alignement

$ sbatch -J SRR1262731_flagstat -o logs/SRR1262731_flagstat.out -e logs/SRR1262731_flagstat.err --wrap=" \

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

​

$ cat SRR1262731.flagstat.txt

​

​

​

​

​

​

​

​

​