Atelier Variant

niveau 2

Nadia BESSOLTANE - INRAe

Vivien DESHAIES -APHP

École de bioinformatique niveau 2 AVIESAN-IFB-INSERM 2023

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)

du Fastq aux VCFs (ebaiin1)

Fastq Quality Control

-- FastQC --

Mapping

-- Bwa --

Reads (Fastq)

Reference genome (Fasta)

Processing Post Alignement

-- GATK/Picard --

small variants

-- GATK --

Structural variants

-- Delly --

Annotation

-- SnpEff --

Variant Calling

UNIX

VCF (variant call format)

du VCFs aux marqueurs (ebaiin2)

du Fastq aux VCFs (ebaiin1)

5

Fastq Quality Control

-- FastQC --

Mapping

-- Bwa --

Reads (Fastq)

Reference genome (Fasta)

Processing Post Alignement

-- GATK/Picard --

small variants

-- GATK --

Structural variants

-- Delly --

Annotation

-- SnpEff --

Variant Calling

QC & Filtres

stats descriptives

Recherche de marqueurs/QTL

Genome Wide Association Study (GWAS)

analyse phylogénétique

….

Dépend de la question biologique

UNIX

6

##fileformat=VCFv4.2

##FILTER=<ID=LowQual,Description="Low quality">

##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">

##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtere

##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VC

##GATKCommandLine=<ID=HaplotypeCaller,CommandLine="HaplotypeCaller --min-base-quality-score 18 --emit-ref-confidence NONE

##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">

##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">

##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities

...

##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read positio

##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">

##contig=<ID=6,length=119458736>

##source=HaplotypeCaller

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SRR1262731 SRR1262732

6 2 . T A 67.64 . AC=1;AF=0.500;... GT:AD:DP:GQ:PL 0/1:3,2:5:75:75,0,105 0/1:3,2:5:75:75,

6 4 . GT G 58.60 . AC=1;AF=0.500;... GT:AD:DP:GQ:PL 0/1:1,2:3:28:66,0,28 0/1:1,2:3:28:66,

6 9 . C CA 55.60 . AC=1;AF=0.500;... GT:AD:DP:GQ:PL 0/1:7,2:9:63:63,0,279 0/1:7,2:9:63:63,

​

​

VCF header

Body

Lire un VCF sous R

​

METADATA

INFO

GENOTYPE

Lire un VCF sous R

vcfR package : Manipulate and Visualize VCF Data

• vcfR documentation : https://knausb.github.io/vcfR_documentation/index.html

> ## installer le package

> install.packages(“vcfR”)

> ## charger le package

> library(vcfR)

> # lire le fichier vcf

> my.vcf <- read.vcfR(“pool.vcf”)

> # l’objet vcf appartient à quelle class

> is(my.vcf)

[1] "vcfR"

> # la liste des slots (sections)

> slotNames(my.vcf)

[1] "meta" "fix" "gt"

Objet de la classe vcfR

Trois sections :

• meta-information : entête du vcf
• Fixed information : information par variant mais commune à tous les échantillons (position, allèles, qualité…)
• Genotype information : information de génotypage par échantillon

8

##fileformat=VCFv4.2

##FILTER=<ID=LowQual,Description="Low quality">

##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">

##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtere

##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">

##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">

##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VC

##GATKCommandLine=<ID=HaplotypeCaller,CommandLine="HaplotypeCaller --min-base-quality-score 18 --emit-ref-confidence NONE

##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">

##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">

##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities

...

##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read positio

##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">

##contig=<ID=6,length=119458736>

##source=HaplotypeCaller

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SRR1262731 SRR1262732

6 2 . T A 67.64 . AC=1;AF=0.500;... GT:AD:DP:GQ:PL 0/1:3,2:5:75:75,0,105 0/1:3,2:5:75:75,

6 4 . GT G 58.60 . AC=1;AF=0.500;... GT:AD:DP:GQ:PL 0/1:1,2:3:28:66,0,28 0/1:1,2:3:28:66,

6 9 . C CA 55.60 . AC=1;AF=0.500;... GT:AD:DP:GQ:PL 0/1:7,2:9:63:63,0,279 0/1:7,2:9:63:63,

​

​

VCF header

Body

objet de la classe vcfR

​

@meta

@fix

@gt

Rappel : Filtre des variants

• De nombreux filtres peuvent être appliqués sur le VCF

→ type de variants à garder (SNVs seulement, Indels...)

→ région d’intérêt

→ filtres sur la qualité (seuils arbitraires : profondeur, génotype (0/1, 1/1), ratio allélique…)

​

​

• Filtres difficilement transposables entre analyse :

→ dépendent de la question biologique

→ dépendent des outils utilisés

​

​

​

9

→ d’où l'intérêt de faire ses propres scripts

Rappel : Faux positif

• De nombreux variants Faux Positifs peuvent survenir des étapes précédentes :
• Artéfacts issus des cycle PCR pendant la préparation des échantillons
• Artéfacts liés à la technologie de séquençage (PacBio, HiSeq, NextSeq, … )
• Difficultées d’alignement (régions d’ADN répétées)
• Erreurs de lecture lors du “BaseCalling”

​

​

Rappel :

Faux positifs vs Filtres qualité

True Positive (polymorphisms,mutations)

False

Positive

False Negative

True Negative

Filtres/seuils

résultat final

Variant calling results

Plus on est stringent plus on va éliminer les faux positifs mais avec le risque de perdre de vrais variants

​

max de vrais variants

min de faux positifs

​

• Analyse post-VCF des petits variants (SNV & InDel) → 4h
• Analyse post-VCF des variants structuraux (SV) → 2h

12

Objectifs�

Petites Variations Génomiques

SNV and InDels

Jeux de données #1:

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.

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.

Quelle mutation est responsable de ce QTL ?

Pour le TP nous disposons des résultats du variant calling de 3 échantillons (en Multi-VCF annoté).

TP 1 : recherche de mutation dans un QTL

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

1- Se connecter sur Rstudio via JupyterLab

​

2- Copier le matériel de TP dans le dir TP_variants

> file.copy(from = “/shared/projects/2315_ebaii_2/DNAseq/TP_small_variants”,

to = “~/”, recursive=TRUE)

​

2- Positionner l’espace du travail

> setwd(“~/TP_small_variants”)

​

5- Changer le workspace sur l’interface RStudio.

Préparation des données

RMD

Variations Structurales

17

18

Principe de détection des SVs

Short reads (Illumina) : selon l’outil et la qualité des données

→ faible sensibilité : 10 à 70% des SVs détectés

→ faible précision : jusqu’à 90% de Faux Positifs

→ Difficulté à caractériser des SVs complexes (alignement imprécis dans les régions répétées et faible résolution)

/!\ Un calling consensus avec plusieurs outils de détection peut être utile avec des données short reads /!\

​

​

​

​

​

19

Short reads ou long reads?

Long reads (PacBio / Nanopore) :

→ Meilleure caractérisation des altérations des régions répétées

→ Une faible profondeur de couverture suffit (15-30x)

Sedlazeck, F. J. et al.Nat. Methods 15, 461–468 (2018).

• La détection des SVs manque de précision et engendre des faux négatifs

→ Nécessité de croiser différents outils / technologies

→ Nécessité de bien utiliser les métriques des outils

→ Nécessité d’une bonne profondeur (variant hétérozygote)

​

• Vérifier visuellement les résultats sur IGV permet d’augmenter la confiance dans les SVs détectés

20

Conclusion

​

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 Nanopore

​

21

Your turn !

Retrouvez les délétions de grande taille

Jeux de données #2 : SVs

Data : souche de Zymoseptoria tritici séquencées à la fois en Illumina et en Nanopore.

→ chaque set de reads a été aligné sur le génome de référence avec les outils dédiés

→ les données ont été réduites aux premiers 500kb du chr10

Tools :

• Delly (Bioinformatics, Volume 28, Issue 18, 15 September 2012, Pages i333–i339, https://doi.org/10.1093/bioinformatics/bts378)
• Sniffles (Nature Methods volume 15, pages 461–468 (2018) , https://www.nature.com/articles/s41592-018-0001-7) with NGMLR mapping

​

​

​

�

​

22

Partie TP

1- Se connecter sur Rstudio via JupyterLab

​

2- Copier le matériel de TP dans le dir TP_variants

> file.copy(from = “/shared/projects/2315_ebaii_2/DNAseq/TP_structural_variants”,

to = “~/”, recursive=TRUE)

​

2- Positionner l’espace du travail

> setwd(“~/TP_R”)

​

5- Changer le workspace sur l’interface RStudio.

Préparation des données

RMD