1 Objectif

L’objectif de ce TP est de vous familiariser avec les packages R utilisés pour la manipulation et la visualisation de variations génomiques à partir de fichiers VCF.

2 Jeux de données

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.

Nom d’échantillon Phénotype
SRR1262731 QTL-
SRR1205992 QTL+
SRR1205973 QTL+

Quelle mutation est responsable du phénotype QTL+?

Le point d’entrée de ce TP sera le fichier multi-VCF généré par le workflow décrit dans les cours/TP “Atelier Variants, niveau 1”.

3 Extraire les variations génomiques du VCF (vcf -> tabelau)

L’objectif est d’extraire les variations génomiques à partir d’un fichier VCF, puis de parser et de normaliser les informations du champ INFO ainsi que les annotations qui y sont associées, afin de les représenter sous forme de table structurée : chaque ligne correspondant à un variant, et chaque colonne à une annotation ou une métrique spécifique liée à ce variant.

Définir le nom (+ chemin) du fichier VCF (one-sample, multi-samples)

# le chemain vers le vcf
vcf_file <- "data/pool_GATK_annot.vcf"

On va utiliser le package R vcfR, qui est un ensemble d’outils conçu pour lire, écrire, manipuler et analyser les données au format VCF. vcfR documentation

Installer le package

#install.packages("vcfR")

Charger le package

library(vcfR)

3.1 Lire le fichier VCF

La fonction read.vcfR() du package vcfR permet de lire un fichier VCF/multi-VCF et de retourner un objet de la classe vcfR.

# pour chercher l'aide : help(read.vcfR) / ?read.vcfR()
?read.vcfR()
# Lire le fichier VCF
vcf <- vcfR::read.vcfR(file = vcf_file, verbose = FALSE)

L’objet vcf appartient à quelle classe :

is(vcf) # ou class(vcf)
## [1] "vcfR"

Cet objet est composé de :

# Les slots correspondent aux éléments composant un objet de classe S4.
slotNames(vcf) 
## [1] "meta" "fix"  "gt"
  • meta (meta-information) : entête du vcf
  • fix (fixed information) : information par variant mais commune à tous les échantillons (position, allèles, qualité…)
  • gt (genotype information) : information de génotypage par échantillon

Ces slots corresponds aux différents parties d’un fichier vcf.

3.2 Convertir l’objet vcfR en objet tibble

La fonction vcfR2tidy() du package vcfR permet de convertir un objet vcfR en un data frame de type tibble, un format compatible avec les fonctions du tidyverse et plus facile à manipuler.

vcf.tidy.list <- vcfR::vcfR2tidy(x = vcf)

On peut spécifier les champs que l’on souhaite conserver (voir l’exemple ci-dessous).

#vcf.tidy.list <- vcfR2tidy(x = vcf.snv, 
#                           format_fields = c("GT", "AD", "DP"), 
#                           info_fields = c("AC", "AN", "MQ"))

Cette fonction retourne une :

class(vcf.tidy.list)
## [1] "list"

Pour accéder aux noms des éléments de cette liste :

names(vcf.tidy.list)
## [1] "fix"  "gt"   "meta"

Pour accéder à un des éléments (tibbles) de cette liste :

# vcf.tidy.list[["meta"]]
head(vcf.tidy.list[["meta"]])
head(vcf.tidy.list$meta)

La fonction head() permet d’afficher uniquement les 6 premières lignes d’un tableau.

Le tibble meta contient les informations présentes dans l’en-tête du fichier VCF.

!!! Attention : Les outils de variant calling n’utilisent pas toujours les mêmes noms de champs dans la colonne ID.

Le tibble fix contient les informations liées à chaque variation génomique.

head(vcf.tidy.list[["fix"]])

Le nombre de ligne de ce tableau correspond au nombre de variations génomiques.

dim(vcf.tidy.list[["fix"]]) # nrow() permet d'afficher que le nombre de lignes
## [1] 2826   28

Le tibble gt contient les inforamtions liées au génotype de chaque échantillon et pour chaque variation génomique.

head(vcf.tidy.list[["gt"]])

Le nombre de ligne de ce tableau correspond au nombre de variations génomiques multiplier par le nombre d’échantillons.

nrow(vcf.tidy.list[["gt"]]) # dim() permet d'afficher le nombre de lignes et le nombre de colonnes.
## [1] 8478

3.3 Fusionner les deux tibbles $fix et $gt :

Pour cela on va utiliser la fonction full_join() du package dplyr pour fusionner les deux tableaux. Il faut commencer par definir les colonnes communes entre les deux tibbles. Dans notre cas c’est la colonne ChromKey et la colonne POS.

library(dplyr) # manipulation des data
# appliquer le merge
vcf.tidy <- dplyr::full_join(x = vcf.tidy.list[["fix"]], 
                             y = vcf.tidy.list[["gt"]],  
                             by = c("ChromKey", "POS"))

Afficher les 6 premières lignes du tibble fusionné.

head(vcf.tidy)

La taille de ce tibble

dim(vcf.tidy)
## [1] 8478   41

3.4 Cas de plusieurs fichiers VCF

Dans le cas où l’on dispose de plusieurs fichiers VCF (un fichier par échantillon).

On utilise une boucle for ou lapply.

Extraire les noms (+chemin) des fichiers vcf

vcf.list <- list.files(path = "data/", pattern = ".vcf", full.names = TRUE)

Extraire une liste de tibbles par échantillons (avec la boocle lapply)

tidy.list.tmp <- lapply(vcf.list, function(vcf_file){
  vcf <- read.vcfR(file = vcf_file, verbose = FALSE)
  vcf.tidy.list <- vcfR2tidy(x = vcf)
  full_join(x = vcf.tidy.list[["fix"]], y = vcf.tidy.list[["gt"]],
            by = c("ChromKey", "POS"))
})

Extraire une liste de tibbles par échantillons (avec la boocle for)

tidy.list.tmp <- list()
for(i in 1:length(vcf.list)){
  vcf <- read.vcfR(file = vcf.list[[i]], verbose = FALSE)
  vcf.tidy.list <- vcfR2tidy(x = vcf)
  tidy.list.tmp[[i]] <- full_join(x = vcf.tidy.list[["fix"]], y = vcf.tidy.list[["gt"]],
                                  by = c("ChromKey", "POS"))

}

Fusionner les tableaux avec la fonction de base rbind(), qui concatène deux matrices ou data frames, par les lignes. !!! Attention : assurez-vous que les noms des échantillons sont distincts afin d’éviter toute confusion lors de l’agrégation.

# tidy.list <- rbind(tidy.list.tmp[[1]], tidy.list.tmp[[2]])
# tidy.list <- rbind(tidy.list,          tidy.list.tmp[[3]])
# ...
# la fonction de base Reduce permet d'appliquer la fonction rbind successivement
# 2 par 2 à partir d'une liste de matrice
tidy.list <- Reduce(rbind, tidy.list.tmp)

4 Processer le tableau des variants

Cette étape permet de rendre le tableau des variants plus simple à manipuler.

4.1 Retirer les variations multi-alleliques : (Qu’est ce que c’est ?)

Dans ce TP on va s’intéresser uniquement aux variations avec un allele alternatif. Dans le cas des variants multi-alléliques les valeurs de la colonne ALT contiennent une “,” pour séparer les différents allèles alternatifs.

Pour detecter ces cas on va utiliser la fonction str_detect() du package stringr qui permet de detecter un pattern dans un vecteur de chaines de caractères.

library(stringr)

Example :

fruits <- c("banane", "pomme", "mangue", "raisin", "kiwi")
index <- stringr::str_detect(string = fruits , pattern = "a")
index
## [1]  TRUE FALSE  TRUE  TRUE FALSE

La fonction str_detect() renvoie un vecteur de bouléens.

Pour récupérer la liste des fruits dont le nom contient la lettre ‘a’ :

fruits[index]
## [1] "banane" "mangue" "raisin"

Pour récuperer les noms de fruits qui ne contiennent pas la lettre ‘a’, on utilise le signe ! pour inverser la condition.

fruits[!index]
## [1] "pomme" "kiwi"

On va appliquer le même principe pour récupérer les variations génomiques avec un allèle unique. On va d’abord repérer les variations multi-alleliques grâce à la présence de “,” dans la colonne ALT.

index <- stringr::str_detect(string = vcf.tidy[["ALT"]] , pattern = ",")

On va ensuite les retirer du tableau des variants. La fonction filter() du package dplyr permet de filtrer les lignes selon une ou plusieurs conditions.

vcf.u.tidy <- dplyr::filter(vcf.tidy, !index)
head(vcf.u.tidy)

4.2 Définir un identifiant unique par variation :

Le champs ID renseigne un identifiant unique pour chaque variant. S’il n’est pas renseigné il est recommandé d’en créer un pour faciliter la manipulation des tableaux.

Une variation est définie par sa position et par son allèle.

On va utiliser la fonction unite() du package tidyr, pour créer un identifiant qui est le résultat de la concaténation des colonnes : CHROM, POS, REF et ALT.

# on commence par charger le package
library(tidyr) 
?unite()
# unite(data, col, ..., sep = "_", remove = TRUE, na.rm = FALSE)
vcf.u.tidy <- 
  tidyr::unite(data = vcf.u.tidy, col="ID", CHROM, POS, REF, ALT, sep="_", remove = FALSE)
head(vcf.u.tidy)

4.3 Ajouter une colonne avec le type des variation (SNV/InDel)

Certains outils de variant-calling donnent l’info du type de variants dans le VCF. Ce n’est pas le cas de tous les outils!

Pour différencier les SNV des INDEL on teste le nombre de bases dans les colonnes REF et ALT. Si le nombre de bases est supéreur à 1 cela signifie que le variant est de type INDEL, sinon c’est un SNV. Pour effectuer ce test on va utiliser la fonction if_else() du package dplyr combiner à la fonction mutate() de dplyr afin de rajouter une colonne avec le type de variations.

Pour calculer le nombre de bases (caractères) de l’allèle alternatif ou de l’allèle de référence, on utilise la fonction str_length() du package stringr.

Exemple d’utilisation de if_else() et de str_length()

?if_else()
#if_else(condition, true, false, ...)

alt <- c("AA", "C", "AT", "G")
if_else(str_length(alt) == 1, "?", "INDEL")
## [1] "INDEL" "?"     "INDEL" "?"

On va appliquer ce principe sur nos données. On va rajouter une colonne nommée var_type avec le type de variations.

vcf.u.tidy <- 
  dplyr::mutate(
    vcf.u.tidy, 
    var_type = if_else(str_length(REF) == 1 & str_length(ALT) == 1, "SNV", "INDEL"))

La fonction select() du package dplyr permet de selectionner les colonnes que on souhaite garder ou visualiser dans ce cas.

dplyr::select(vcf.u.tidy, ID, var_type, Indiv) |> head()

4.4 Ajouter une colonne avec les allèles ratio (AR) des variants

L’AR à une position donnée est le ratio entre le nombre de reads (depth) qui couvrent l’allèle variant/alternatif et le nombre de reads total qui couvrent cette position.

L’information liée au nombre de reads par allèle est stockée dans la colonne gt_AD (Allele depth) sous forme de (ref AD,alt AD). Pour accèder à “alt AD”, on va utiliser la fonction separate() du package tidyr pour separer la colonne gt_AD en 2 colonnes, gt_ref.AD et gt_alt.AD.

vcf.u.tidy <- tidyr::separate(data = vcf.u.tidy, 
                              col = gt_AD, into = c("gt_ref.AD", "gt_alt.AD"), 
                              sep = ",", remove=FALSE)
dplyr::select(vcf.u.tidy, ID, Indiv, gt_AD, gt_ref.AD, gt_alt.AD) |> head()

On va utiliser la fonction mutate() aussi pour ajouter une colonne avec l’allèle ratio à chaque ligne

vcf.u.tidy <- dplyr::mutate(vcf.u.tidy, gt_AR = as.numeric(gt_alt.AD)/gt_DP)
select(vcf.u.tidy, ID, gt_ref.AD, gt_alt.AD, gt_DP, gt_AR, gt_GT, Indiv) |> head()

BONUS : gt_DP peut contenir des 0 ou des NA. Dans ces cas là on ne pourra pas calculer le ratio. La solution serait :

# dplyr::mutate(vcf.u.tidy, gt_AR = if_else(gt_DP %in% c(0, NA), 0, gt_alt.AD/gt_DP))

5 Contrôle qualité

Le champs INFO contient plusieurs métriques qui vont nous permettre d’évaluer la qualité de nos variants, comme par exemple QUAL, MQ, DP…etc

Le meilleur moyen d’accèder à ces metriques est d’afficher leur distribution sous forme de graphe. Pour cela on va utiliser les fonctions du package ggplot2. La visualisation des différents métriques nous permet de voir la qualité de nos données mais aussi de définir les filtres et les seuils à appliquer.

5.1 Score qualité (colonne QUAL)

dplyr::select(vcf.u.tidy, ID, QUAL, Indiv)

Pour visualiser la distribution de la qualité on va utiliser les fonctions du package ggplot2.

library(ggplot2)

La fonction geom_density() permet de représenté la densité de nos scores de qualité. Il est possible de superposer des graphes sur le même fond. La fonction geom_vline() nous permet de tracer une ligne verticale pour représenter le cut-off qui sera utilisé pour filtrer.

p <- ggplot(data = vcf.u.tidy, aes(x = QUAL))
p <- p + geom_density()
p <- p + xlim(0,500) +  xlab("score de qualité") 
p <- p + geom_vline(xintercept = 40, color = "red", linetype = "longdash") 
p

5.2 Couverture

  • Distribution de la profondeur de couverture de toutes les positions variantes.

On utilise les même fonction que ci-dessus.

p <- ggplot(data=vcf.u.tidy, aes(x = gt_DP, color=Indiv)) + 
  geom_density() + 
  xlab("Variant position coverage") + xlim(c(0,100)) + 
  geom_vline(xintercept = 4, color = "red", linetype = "longdash")
p

  • Distribution de l’allèle ratio des variants
p <- 
  ggplot(data=vcf.u.tidy, aes(gt_AR, color=Indiv)) + 
  geom_density() + 
  labs(y="Alternative allele ratio")
p

Ches les organismes diploide hétérozygote on observe un pic à 0.5 (qui correspond au variants hétérozygotes) et 2 pics à 1 et 0 (qui corresponent au variants homozygotes).

5.3 Filtrer les variations génomiques de moins bonnes qualité

Il y a plusieurs filtre qui pourrraient être appliqué pour réduire la présence de faux positifs et de variations de mauvaises qualités. Ci-dessous on vous donne quelques exemples.

Nombre de variants avant le filtre

length(unique(vcf.u.tidy[["ID"]]))
## [1] 2766
  • Filtrer les variations avec un score qualité inférieur à 40

On va utiliser la fonction filtre() pour filtrer que les variations sur la base de la colonne QUAL.

# filtrer les variations avec un score qualité inférieur à 30
dplyr::filter(vcf.u.tidy, QUAL >= 40) |> head()
  • Filtrer les variation avec une profondeur de couverture supérieur ou égale à 4

L’information de la couverture est stockée dans la colonne gt_DP.

dplyr::filter(vcf.u.tidy, gt_DP >= 4) |> head()
  • Filtrer les variation avec une profondeur de couverture de l’allèle alternatif supérieur ou égale à 2

Cette information est stockée dans la colonne gt_alt.AD.

dplyr::filter(vcf.u.tidy, gt_alt.AD >= 2) |> head()

Il est possible de regrouper plusieurs conditions avec la fonction filter().

vcf.q.tidy <- dplyr::filter(vcf.u.tidy, 
                            QUAL      >= 30, 
                            gt_DP     >= 4, 
                            gt_alt.AD >= 2, 
                            gt_AR     >= 0.2)

Nombre de variants après le filtre

length(unique(vcf.q.tidy[["ID"]]))
## [1] 1321

6 Statistiques discriptives

6.1 Nombre de variants

Compter le nombre de variations génomiques par échantillon.

La fonction group_by() du package dplyr permet de grouper un tibble afin d’appliquer une opération sur chacun des groupe. La fonction count() du même package permet de compter le nombre de ligne par groupe (par individu dans notre cas).

# |> permet d'enchainer les commandes
tb.nbr.var <- dplyr::group_by(vcf.q.tidy, Indiv) |> dplyr::count(name = "nb_var")
tb.nbr.var

Pour représenter ces comptages sous forme de barres on va utiliser la fonction geom_bar() du package ggplot.

ggplot(tb.nbr.var, aes(x = Indiv, y = nb_var)) + 
  geom_bar(stat = "identity") + 
  labs(y="nombre de variants par individu", x="")

Dans l’échantillon ‘SRR1262731’, il y a moins de variations que dans les autres. Dans notre cas, cela est dû à la profondeur de couverture qui est plus faible pour cet échantillon. Dans d’autres cas, cela peut être dû à la distance entre les variants de l’espèce et le génome de référence.

BONUS

Il est possible de compter le nombre de variations par type de variants et par échantillon. Cette fois on groupe le tableau sur deux niveaux (par undividu et par type de variant).

tb.nbr.var <- 
  dplyr::group_by(vcf.q.tidy, Indiv, var_type) |> 
  dplyr::count(name = "nb_var")

tb.nbr.var

Le paramêtre fill de ggplot permet de distinguer les deux types de variants

ggplot(tb.nbr.var, aes(x = Indiv, y = nb_var, fill = var_type)) + 
  geom_bar(stat="identity") + 
  labs(y="nombre de variants par individu", x="", fill = "Type de variants")

BONUS Une autre façon de faire :)

# position=position_dodge() permet de mettre les barres côte à côte
# http://www.sthda.com/french/wiki/ggplot2-barplots-guide-de-demarrage-rapide-logiciel-r-et-visualisation-de-donnees
ggplot(tb.nbr.var, aes(x = Indiv, y = nb_var, fill = var_type)) + 
  geom_bar(stat = "identity", position=position_dodge()) + 
  labs(y="nombre de variants par individu", x="", fill = "Type de variants")

6.2 Intersection des variations génomiques entre les 3 échantillons

On peut comparer les variants des 3 échantiilons avec un diagramme de venn.

Pour porduire un diagramme de venn, on a besoin d’extraire une matrice binaire de presence absence des variants dans chaque échantillon. Avec les variants en lignes et les échantillons en colonnes. La fonction dcast() du package reshape2 permet de construire ce type de matrice.

library(reshape2)
?dcast()

Transformer la matrice de données en format “wide” -> ID en lignes et Indiv en colonne

vcf.mat.bin <-
  select(vcf.q.tidy, ID, Indiv) |> 
  mutate(gt_bin = 1) |>
  reshape2::dcast(ID~Indiv, value.var="gt_bin", fun.aggregate = sum) 
# (fun.aggregate = sum) permet de remplacer NA par 0.
head(vcf.mat.bin)

La fonction venn() du package venn permet de représenter les intersections entre les échantillons.

library(venn)
# mat[,-1] : retire la première colonne.
venn::venn(vcf.mat.bin[,-1], zcolor = "style", ilabels = "counts")

6.3 Clustering hiérarchique des échantillons (BONUS)

Pour produire un dendrogramme, on a besoin d’extraire une matrice de allèle ratio des variants dans chaque échantillon. Avec les variants en lignes et les échantillons en colonnes. La fonction dcast() du package reshape2 permet de construire ce type de matrice.

vcf.mat.AR <-
  select(vcf.q.tidy, ID, Indiv, gt_AR) |> 
  reshape2::dcast(ID~Indiv, value.var="gt_AR", fun.aggregate = sum) 
# (fun.aggregate = sum) permet de remplacer NA par 0.

head(vcf.mat.AR)

Calculer la distance euclidienne entre les lignes de la matrice

dist  <- dist(t(vcf.mat.AR[,-1]), method = "euclidean")
dist
##            SRR1205973 SRR1205992
## SRR1205992   12.76952           
## SRR1262731   13.31545   14.32771

Effectuer une classification hiérarchique des échantillons à partir des distances.

# clustering hiérarchique
hierarchical_result    <- hclust(dist)

Convertir la sortie de hclust() en objets de classe “dendrogram”.

dendrogram_result  <- as.dendrogram(hierarchical_result)

Afficher le dendrogramme avec la fonction plot()

plot(dendrogram_result, main = "Sample clustering using all variants")

7 Identification de la mutation associée au QLT-

La mutation recherchée :

  • se situe dans un gène codant pour une protéine,
  • c’est une variation missence (changement de l’acide aminé),
  • présente seulement sur le locus QTL- (SRR1262731)

On va commencer par appliquer les filtres nécessaires pour cerner la mutation d’interêt.

La colonne ANN, rajouter par l’outil d’annotation “snpEff”. contient les informations d’annotation associées aux variants.

vcf.q.tidy[1,"ANN"] # afficher la 1ere ligne de la colonne "ANN"

Pour filtrer les variants qui se trouve sur un gène codant pour une protéine, on va utiliser la fonction str_detect() ou (la fonction de base grepl) pour rechercher le motif protein_coding.

dplyr::filter(vcf.q.tidy, stringr::str_detect(string = ANN, pattern = "protein_coding")) |> head()

Pour filtrer les variants qui ont un effet “missense”, on va utiliser la même fonction grepl() pour rechercher le motif missense_variant dans la colonne ANN.

dplyr::filter(vcf.q.tidy, grepl("missense_variant", ANN)) |> head()

Il est possible de combiner plusieurs condition dans la même commande.

vcf.filt.tidy <- 
  dplyr::filter(vcf.q.tidy, grepl("protein_coding" ,ANN),  grepl("missense_variant" ,ANN))

Pour sélectionner que les variants qui portent l’allèle référent sur les échantillons QTL+ et l’allèle alternatif sur l’echantillons QTL-. On va utiliser de nouveau la fonction dcast() pour obtenir un tabeau au format large (matrice ID~Indiv).

vcf.filt.AR <- reshape2::dcast(vcf.filt.tidy, ID~Indiv, value.var="gt_AR", fun.aggregate = sum) 
vcf.filt.AR

Grace à ce format large et à l’aide de la fonction filter(), on va pouvoir identifier les mutations absentes (qui ont un AR égale à 0) chez les échantillons QTL+ et présentes (AR différent de 0) chez l’échantillons QTL-.

mutations <- dplyr::filter(vcf.filt.AR, SRR1262731 !=  0 & SRR1205973 == 0 & SRR1205992 == 0)
mutations

L’identifiant de la mutation recherchée :

mutations$ID
## [1] "6_38027010_A_C"

Il est possible d’extraire l’information complète de cette mutation du tableau d’origine :

resultats <- dplyr::filter(vcf.u.tidy, ID == mutations$ID)
resultats

Extraire le nom du gène impacté et l’effet de la mutation.

resultats <- 
  tidyr::separate(resultats, col = ANN, into = c("allele", "effect", "impact", "gene"), sep = "\\|") 
dplyr::select(resultats, ID, Indiv, gt_AR, gene)

Ecrire les résultats dans un fichier txt tabulé avec une extention “xls” pour le visualiser facilement sous excel.

dplyr::select(resultats, ID, CHROM, POS, REF, ALT, var_type, QUAL, effect, impact, gt_DP, gt_alt.AD, gt_AR ,gt_GT) |>
  write.table(file = "QTLminus_mutation.xls", quote = FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)