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 le cours/TP “variant-calling et annotation” du niveau 1.

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

3 Lire le fichier vcf sous R

3.1 vcfR package: Manipulate and Visualize VCF Data

Le package R vcfR 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 / équivalent de module load
library(vcfR)

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()

Exécuter la fonction read.vcfR()

vcf <- vcfR::read.vcfR(file = vcf_file, verbose = FALSE)

L’objet vcf appartient à quelle classe ?

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

Trois sections (slots) : * 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

La liste des slots :

slotNames(vcf) 
## [1] "meta" "fix"  "gt"

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

Rappel * @meta => VCF Header * @fix => VCF INFO * @gt => VCF GENOTYPE

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

vcf.tidy.list <- 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 liste de 3 tibbles fix, gt, meta

is(vcf.tidy.list)
## [1] "list"   "vector"

Pour accéder aux noms de ces tibbles (les noms des élements d’une liste), on utilise la foction name.

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

Pour accéder à un tibble (un élement de la liste), on utilise $.

# vcf.tidy.list$meta
head(vcf.tidy.list$meta)

La fonction head permet d’afficher que les 6 premières lignes d’un tableau. Le tibble meta contient les informations présentent dans l’entête du fichier vcf.

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.2 Fusionner les deux tibbles $fix et $gt :

Pour cela on va utiliser la fonction full_join() du package dplyr pour futionner 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.3 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)

3.4 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)

3.5 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 afin de rajouter une colonne avec le type de variations.

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

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

3.6 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 definir les colonnes ”gt_ref.AD” et “gt_alt.AD” comme numériques, pour pouvoir calculer le ratio.

vcf.u.tidy <- dplyr::mutate(vcf.u.tidy, gt_ref.AD = as.integer(gt_ref.AD), gt_alt.AD = as.integer(gt_alt.AD) )
dplyr::select(vcf.u.tidy, ID, gt_AD, gt_ref.AD, gt_alt.AD, gt_DP, gt_GT, Indiv) %>% 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 = gt_alt.AD/gt_DP)
select(vcf.u.tidy, ID, gt_AD, 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))

4 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.

4.1 Score qualité (colonne QUAL)

tb <- dplyr::select(vcf.u.tidy, ID, QUAL)
head(tb)

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

library(ggplot2)

ggplot() est une fonction clé du package et permet de initialiser la figure.

p <- ggplot(data = tb, aes(x = QUAL))
p

La fonction geom_density() permet de représenté la densité de nos scores de qualité

p <- p + geom_density()
p

Le fonction xlim() permet de définir la feunettre à visualiser et les fonctions xlab()/ylab() permettent de changer les noms des axes X et Y.

p <- p +  xlim(0,500) +  xlab("score de qualité") 
p

Il est possible de superposer des graphes sur le même fond. La fonction geom_vline() nous permet de tracer un ligne verticale pour représenter le cut-off qui sera utilisé pour filtrer.

p <- p + geom_vline(xintercept = 40, color = "red", linetype = "longdash") 
p

4.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).

4.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.

  • 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)
head(vcf.q.tidy)

5 Statistiques discriptives

5.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).

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 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")

5.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)
# transformer la matrice de données en format "wide" -> ID en lignes et Indiv en colonne
vcf.mat.AR <- reshape2::dcast(vcf.q.tidy, ID~Indiv, value.var="gt_AR", fun.aggregate = sum)
# l'option (fun.aggregate = sum) permet de remplacer NA par 0.
head(vcf.mat.AR)

Pour obtenir une matrice binaire avec 0/1, on utilise la fonction if_else() vu précédemment et la fonction mutate_at() du package dplyr pour remplacer les AR supérieurs à 0 par 1.

vcf.mat.bin <- 
  dplyr::mutate_at(vcf.mat.AR, .vars = c("SRR1205973", "SRR1205992", "SRR1262731"), 
                   .funs = function(x){dplyr::if_else(x == 0, 0, 1)})
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")

5.3 Clustering hiérarchique des échantillons (BONUS)

On va reprendre labelau “vcf.mat.AR”.

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   13.76642           
## SRR1262731   15.61776   16.05623

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")

6 Recherche 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 contient les informations d’annotation des 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() 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 str_detect() 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

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

write.table(file = "QTLminus_mutation.xls", x = resultats, quote = FALSE, sep = "\t", row.names = FALSE, col.names = TRUE)