Présentation des données de RNAseq
Vous allez travailler sur des données non publiées (elles sont sur moodle et non disponibles sur ce site du cours), merci de ne pas les partager et ne pas les déposer en ligne.
Il s’agit de RNAseq qui ont été effectués pour un sauvage Col-0 et le
mutant arp6-1 qui ont poussé avec soit une faible (0.5mM) soit
une forte (10mM) concentration en azote. Il y a donc au total:
- 3 réplicats du sauvage Col-0 avec une
faible (0.5 mM) concentration en azote
- 3 réplicats du sauvage Col-0 avec une
forte (10 mM) concentration en azote
- 3 réplicats du mutant arp6-1 avec une
faible (0.5 mM) concentration en azote
- 3 réplicats du mutant arp6-1 avec une
forte (10 mM) concentration en azote
La protéine ARP6 est impliquée dans l’intégration du variant d’histone H2A.Z dans la chromatine. Le mutant arp6-1 comporte donc une faible quantité de H2A.Z dans la chromatine, remplacé par l’histone H2A.
L’objectif de ce projet est d’étudier le rôle possible de H2A.Z dans la réponse transcriptionnelle à la quantité d’azote dans le milieu de culture des plantes
Une grosse partie du code est indiquée pour les premières étapes, mais plus le projet avance plus vous aurez à trouver le code par vous-même à l’aide des indications dans le matériel et de ce que vous trouvez dans les aides R et en ligne
Préparation de l’environnement et chargement des données
Tout d’abord téléchargez le support du cours et les données.
Ensuite ouvrez Rstudio et créez un nouveau document R Markdown (ou script R) que vous enregistrez dans votre dossier de travail (là où sont vos données et où vous voulez enregistrer les figures).
Attention, si vous utilisez un script R, il faudra indiquer le dossier de travail avec la fonctionsetwd(). Cette étape n’a pas besoin d’être effectuée si vous utilisez un R Markdown.Chargez les librairies dont vous aurez besoin en utilisant la fonction suivante (exemple pour une librairie).
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsVoici la liste des librairies à charger:
- tidyverse
- visdat
- scTenifoldNet
- factoextra
- Enfin, chargez les données et mettez les dans un objet à l’aide du code suivant (modifiez en fonction du chemin vers votre fichier):
RNAseq_data <- read_tsv("../../../data/Raw_featureCounts_WT_arp6.txt")
RNAseq_data## # A tibble: 36,795 × 13
## Geneid arp6_0.5mM_rep1 arp6_0.5mM_rep2 arp6_0.5mM_rep3 WT_10mM_rep1
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 AT1G01010 428 393 344 762
## 2 AT1G01020 689 662 747 739
## 3 AT1G01030 106 234 208 202
## 4 AT1G01040 2048. 2256. 2358. 2448.
## 5 AT1G01050 2743 2936. 3645 3990
## 6 AT1G01060 196 136 228 84
## 7 AT1G01070 422 500 362 390
## 8 AT1G01080 2204 2686 2178 4011
## 9 AT1G01090 7250 10233 13556 10228
## 10 AT1G01100 15991 20832 20542 22295
## # ℹ 36,785 more rows
## # ℹ 8 more variables: WT_10mM_rep2 <dbl>, WT_10mM_rep3 <dbl>,
## # arp6_10mM_rep1 <dbl>, arp6_10mM_rep2 <dbl>, arp6_10mM_rep3 <dbl>,
## # WT_0.5mM_rep1 <dbl>, WT_0.5mM_rep2 <dbl>, WT_0.5mM_rep3 <dbl>Vous avez une ligne par gène et une colonne par échantillon. Pour chaque gène dans chaque échantillon, nous avons une valeur d’expression. C’est le nombre de reads qui ont été séquencé pour ce gène dans cet échantillon.
Vérification de l’import des données
Avant de faire la normalisation des données et leur vérification, nous devons nous assurer que l’import c’est bien passé.
Pour cela, vérifiez que nous avons bien le nombre de lignes et colonnes attendu, que les colonnes sont bien au bon format et qu’il n’y a pas de données manquantes.
Exemple de fonctions que vous pouvez utiliser:
dim(),vis_dat
Nous voyons ici que le format des données est attendu:
- Il n’y a pas de données manquantes.
- La colonne avec le nom des gènes est bien du texte.
- Les colonnes avec les valeurs d’expression sont bien numériques.
Normalisation par la taille de la librarie
Afin de pouvoir comparer les données dans les différens échantillons, il faut les normaliser pour corriger par la taille de la librarie. La taille de la librairie est le nombre total de reads séquencés pour chaque échantillon. Plus ce nombre total de reads est élevé pour un échantillon, plus la valeur d’expression des gènes est élevé. Il faut donc corriger ces différences entre échantillons.
Pour cela nous utilisons la fonction
cpmNormalizationde la librairiescTenifoldNet.Cette fonction a besoin d’un tableau uniquement numérique. Nous modifions donc le tableau pour que la colonne des gènes devienne le nom des lignes avec la fonction suivante:
RNAseq_data <- column_to_rownames(RNAseq_data, var="Geneid")- Normalisation avec la fonction
cpmNormalizationpour avoir le nombre de reads de chaque gène par million de reads de la librarie:
CPMnorm_RNAseq_data <- cpmNormalization(RNAseq_data)- Comme le format de sortie de la fonction n’est pas un data.frame (qui est le format nécessaire pour la suite), nous devons changer le format avec la fontion suivante:
CPMnorm_RNAseq_data <- as.data.frame.matrix(CPMnorm_RNAseq_data)## arp6_0.5mM_rep1 arp6_0.5mM_rep2 arp6_0.5mM_rep3 WT_10mM_rep1
## AT1G01010 10.052817 8.073799 6.019492 12.978821
## AT1G01020 16.183156 13.600140 13.071398 12.587072
## AT1G01030 2.489716 4.807300 3.639693 3.440580
## AT1G01040 48.091454 46.357576 41.252772 41.687225
## AT1G01050 64.427281 60.327508 63.782122 67.959970
## AT1G01060 4.603626 2.793986 3.989664 1.430736
## WT_10mM_rep2 WT_10mM_rep3 arp6_10mM_rep1 arp6_10mM_rep2
## AT1G01010 11.275950 14.355056 7.092418 8.887736
## AT1G01020 15.680861 13.855750 12.332926 13.688946
## AT1G01030 2.450619 3.162273 2.442944 3.188590
## AT1G01040 42.017262 44.063780 34.319422 41.735709
## AT1G01050 60.187518 57.316202 61.585828 65.512694
## AT1G01060 2.140414 6.407764 2.876369 1.099514
## arp6_10mM_rep3 WT_0.5mM_rep1 WT_0.5mM_rep2 WT_0.5mM_rep3
## AT1G01010 8.835543 11.968238 11.3271061 14.596911
## AT1G01020 12.290042 15.863643 17.0578159 14.106847
## AT1G01030 2.391576 3.313739 3.2235243 3.360440
## AT1G01040 39.726729 45.757798 46.9761505 55.981078
## AT1G01050 58.748614 66.371721 67.3134547 61.161756
## AT1G01060 3.543075 1.938890 0.9849657 1.750229Nous avons maintenant un tableau avec les niveaux d’expression normalisés par la taille de la librarie pour chaque gène.
ACP pour s’assurer que les réplicats sont suffisamment similaires
Une des premières vérifications biologique que nous faisons sur nos données est de s’assurer que les 3 réplicats pour chaque génotype et condition sont suffisamment similaire. Pour cela nous utilisons une méthode de visualisation: l’Analyse en Composante Principale (ACP). Avant d’aller plus loin, si vous ne connaissez pas le principe de l’ACP, je vous recommande de regarder cette courte vidéo d’introduction à l’ACP.
- Tout d’abord nous enlevons les gènes non exprimés dans 5 échantillons ou plus (avec une valeur de 0 dans 5 échantillons ou plus):
CPMnorm_RNAseq_data_expressed <- filter( CPMnorm_RNAseq_data,
rowSums(CPMnorm_RNAseq_data[1:12]==0)<5)En suivant la même logique, enlevez aussi les gènes faiblement exprimés (avec une moyenne d’expression inférieure à 2). Dans ce cas, utilisez
rowMeansà la place derowSums.
- Nous pouvons maintenant effectuer l’ACP (notez que nous utilisons
les options
scale. = TRUEetcenter = TRUEafin d’avoir des données réduites (écart type de 1) et centrées autour de zéro). Pour cela nous utilisons la fonction suivante:
pca_res_CPMnorm_RNAseq_data_expressed <- prcomp(t(CPMnorm_RNAseq_data_expressed),
scale. = TRUE, center = TRUE) Nous faisons l’ACP sur les données transformées t() avec
les gènes en colonne et les échantillons en ligne, car nous souhaitons
analyser les échantillons à l’aide de l’information donnée par les
gènes. Nous aurions gardé les données telle quelle si nous avions voulu
représenter les gènes analysés en fonction de l’information des
échantillons (ce que je ne recommande pas car dans ce cas nous avons
trop peu de variables explicatives (12 échantillons) par rapport au
nombre de gènes (autour de 20000 gènes)).
- L’étape suivante consiste à regarder le pourcentage de la variance expliqué par chaque dimension de l’ACP:
fviz_eig(pca_res_CPMnorm_RNAseq_data_expressed)
Nous voyons que les deux première dimensions de l’ACP expliquent environ
50% de la variance. Dans le graphique suivant, quand nous présenterons
les données en fonction de ces deux premières dimensions (ce qui est
généralement fait), gardons en tête que 50% de la variance reste non
représentée.
- Faisons maintenant le graphique de l’ACP:
fviz_pca_ind(pca_res_CPMnorm_RNAseq_data_expressed, pointsize=3)
Il est difficile de comprendre cette ACP car le texte n’est pas très lisible.
- Afin de faciliter l’interprétation de l’ACP, ajoutons un code couleur pour les échantillons:
# chargement d'un tableau expliquant à quel génotype et quelle
# condition environnementale correspond chaque échantillon
data_info <- read_tsv("../../../data/Data_info.txt")
# Maintenant refaisons l'ACP avec le code couleur en utilisant
# les paramères suivants:
# col.ind=data_info$genotype_KNO3 pour ajouter le code couleur
# repel=TRUE pour faire un sorte que le texte soit lisible
# invisible="quali" pour enlever la moyenne de chaque groupe
fviz_pca_ind(pca_res_CPMnorm_RNAseq_data_expressed,
pointsize=3, col.ind=data_info$genotype_KNO3,
repel=TRUE, invisible="quali") 
Que concluez vous de cette ACP?
Est-ce que les réplicats sont suffisamment similaires?
A quoi correspond la dimension 1?
A quoi correspond la dimension 2?
Etude de l’expression de quelques gènes pour lesquels nous avons des résultats préliminaires
L’analyse par ACP semble encourageante. Cependant pour aller plus loin dans la vérification des données biologique, nous pouvons aussi vérifier que l’expression de certains gènes se comporte commme attendu dans les données de RNA-seq.
Les gènes pour lesquels nous savons quelle doit être l’expression
dans les échantillons sont:
- ARP6 (AT3G33520)
- NRT2.1 (transporteur de nitrate, AT1G08090)
Pour faire le graphique pour chaque gène, suivez les étapes ci-dessous:
- Passer le nom des lignes en une colonne “gene” à l’aide de la fonctionrownames_to_column
- Ne garder que la ligne pour notre gène d’intérêt en utilisant la fonctionfilter
- Transformer le format de manière à avoir une colonne avec le nom de l’échantillon et une colonne avec les valeurs d’expression plutôt qu’une colonne par échantillon en utilisantpivot_longer
- Séparez la colonne des échantillons de manière à avoir 3 nouvelle colonnes: une pour le génotype, une pour la condition en nitrate et une pour le réplicat. Pour cela utilisez la fonctionseparate
- Utilisezggplotpour faire un graphique avec l’expression du gène en fonction du génotype et les points colorés en fonction de la condition en nitrate
Graphiqe attendu pour ARP6
