Rapport d’analyses bioinformatiques
Numéro de devis: 00000
2023-01-12
Rapport du projet : Analyse transcriptomique de 6
échantillons de Genre espèce correspondant à 2 conditions et 1
comparaison
A l’attention de : NOM PRENOM , ENTREPRISE - VILLE, PAYS
Analyses/Rédaction : NOM PRENOM, Bioinformatics Engineer
Corrections/Validation : NOM PRENOM, Bioinformatics deputy managerRésultats clés
L’analyse d’expression différentielle a permis de rechercher les
gènes différentiellement exprimés entre vos différentes conditions. Dans
votre cas, nous avons trouvés des différences d’expression
significatives pour votre comparaison UHR vs HBR avec 289 gènes
différentiellement exprimés lorsque l’on choisit un seuil de faux
positif de 5% et 129 gènes différentiellement exprimés lorsque l’on se
place au seuil de 0.1%. Nous vous conseillons donc dans un premier temps
d’examiner les résultats obtenus avec le filtre de p-value ajustée de
0.001 (seuil de faux positif de 0.1%) qui est le plus stringeant afin de
sélectionner les gènes dont le niveau de modulation vous intéresse ou
correspond à vos attentes ou hypothèses scientifiques.
Introduction
Rapport_analyse_differentielle.Rmd Les analyses reportées dans ce report font suite au rapport sur le processing des données de RNA-seq. L’objectif de ce projet est de trouver des features différentiellement exprimés pour les comparaisons des conditions biologiques dans le Tableau 1 ci-dessous. Les analyses statistiques incluent l’exploration des données, la normalisation, les tests statistiques entre les conditions à comparer ainsi que la liste des features différentiellement exprimées entre les conditions.
Ci-dessous le workflow pour l’analyse d’expression différentielle à partir de la table de comptages brute :
Figure 1 - Pipeline
Description des données
Les tables de comptages et les conditions biologiques associées à chaque échantillon sont listées dans le Tableau 2 suivant :
Comptages nuls
La Figure 2 montre le pourcentage de features avec des comptages nuls. Nous nous attendons à ce que ce pourcentage soit similaire dans certaines conditions. Les features avec des comptages nuls dans tous les échantillons sont conservés dans les données mais ne sont pas pris en compte pour l’analyse différentielle avec DESeq2. Ici, 1031 features (36.34%) sont dans cette situation (ligne pointillée). Pour ces features, les valeurs de fold-change et des p-valeurs sont définis comme NA dans les fichiers de résultats.
Figure 2 - Densite de distribution des comptages
Densité
La Figure 3 montre la distribution des nombres de reads pour chaque échantillon (sur une échelle logarithmique pour améliorer la lisibilité). Encore une fois, nous nous attendons à ce que les réplicats aient des distributions similaires. De plus, cette figure montre si le nombre de lectures est de préférence faible, moyen ou élevé. Cela dépend des organismes ainsi que des conditions biologiques considérées.
Figure 3 - Densité de distribution des comptages
Exploration de la variabilité de l’expérience
La principale variabilité au sein de l’expérience devrait provenir des différences biologiques entre les échantillons. Cela peut être vérifié de différentes manières.
La première est de réaliser un calcul de matrice de distance et de clusteriser les résultats. Pour cette mesure de corrélation, la distance euclidienne est utilisée comme mesure de dissimilarité et la méthode de “ward.D2” pour le clustering hiérarchique. Le résultat obtenu est une heatmap avec un clustering hiérarchique : Figure 3 dans l’onglet Dissimilarité des échantillons. Notons ici que cette analyse est effectuée après une transformation des données de comptage avec une transformation stabilisatrice de variance (VST).
La seconde consiste à regarder les premières composantes d’une ACP (Analyse en Composantes Principales) comme sur la Figure 4 dans l’onglet ACP. Cette analyse est également effectuée après transformation de VST.
La troisième est une méthode complémentaire à la première et qui consiste à partir des données brutes non transformées à calculer la statistique de SERE (Simple Error Ratio Estimate). En effet, les méthodes classiques de calculs de similarités ne sont pas de parfaits indicateurs pour estimer la similarité entre des réplicats. La statistique SERE s’avère être un bon indicateur de la similarité entre les échantillons de données RNA-seq. Elle peut ainsi déterminer si deux librairies de RNA-seq sont des réplicats fidèles ou si elles sont globalement différentes. Voir la Figure 5 dans l’onglet Statistique SERE
Définition : VST
VST est une transformation des données qui les rend homoscédastiques, c’est-à-dire que la variance est alors indépendante de la moyenne. Elle est réalisée en deux étapes : (i) une relation moyenne-variance est estimée à partir des données avec la même fonction qui est utilisée pour normaliser les données de comptage et (ii) à partir de cette relation, une transformation des données est effectuée pour rendre la variance indépendante de la moyenne. L’homoscédasticité est une condition préalable à l’utilisation de certaines méthodes d’analyse de données, telles que le clustering hiérarchique ou l’Analyse en Composantes Principales (ACP).Dissimilarité des échantillons
La Figure 4 représente ainsi la heatmap avec le clustering hiérarchique de la dissimilarité des échantillons. Plus la valeur entre deux échantillons est proche de 0 plus les échantillons sont similaires. Réciproquement, au plus la valeur est grande, au moins les échantillons sont similaires.
Figure 4 - Heatmap de la distance entre echantillons
ACP
Sur la Figure 5, la première composante principale (PC1) est censée séparer les échantillons des différentes conditions biologiques, ce qui signifie que la variabilité biologique est la principale source de variance dans les données.
Figure 5 - ACP des échantillons de l’experience
Statistique SERE
Les valeurs de la statistique SERE sont :
- 0 : lorsque les échantillons sont identiques (aucune variabilité) : peut arriver en cas de duplication de l’échantillon
- 1 : pour les réplicats techniques (variabilité technique suivant une distribution de Poisson)
- >1 : pour les réplicats biologiques et les échantillons provenant de différentes conditions biologiques (la variabilité biologique est plus importante que la variabilité technique). Plus le coefficient de la valeur SERE est élevé, moins les échantillons sont similaires. Il est également attendu que cette valeur soit aux alntours de 1 pour des réplicats et soit beaucoup plus grande pour des échantillons provenant de différentes conditions biologiques.
Figure 6 - Heatmap des coefficients SERE
Analyse différentielle
L’analyse d’expression différentielle avec DESeq2 comporte plusieurs étapes. En bref, DESeq2 modélise les comptages bruts, en utilisant des facteurs de normalisation (size factor) pour tenir compte des différences de profondeur de librairie. Puis, il détecte les outliers (valeurs aberrantes) du modèle pour lesquels au moins un échantillon ne semble pas lié au plan d’expérience. Ensuite, il estime les dispersions par gène et réduit ces estimations pour générer des estimations plus précises. L’avant dernière étape consiste à effectuer un filtrage indépendant afin d’augmenter le puissance de détection des features exprimés différentiellement en se basant sur la moyenne des comptes normalisés, indépendamment de la condition biologique. Enfin, DESeq2 ajuste le modèle binomial négatif et effectue des tests d’hypothèse en utilisant le test de Wald.
Normalisation
La normalisation a pour but de corriger les biais systémiques dans les données et pour rendre comparable les comptages de reads à travers les échantillons. La normalisation réalisée avec DESeq2 repose sur l’hypothèse que la plupart des features ne sont pas différentiellement exprimés. Elle calcule ainsi un facteur d’ajustement (size factor) pour chaque échantillon. Les comptages de reads normalisés sont obtenus en divisant les comptes bruts par le facteur d’ajustement associé à chaque échantillon.
- Un size factor autour de 1 signifie qu’aucune normalisation n’est requise
- Un size factor < 1 produira des comptes normalisées plus élevés que les comptes bruts
- Un size factor > 1 produira des comptes normalisées plus bas que les comptes bruts
| Size_factor | |
|---|---|
| UHRMix1_REP1 | 1.53 |
| UHRMix1_REP2 | 0.91 |
| UHRMix1_REP3 | 1.25 |
| HBRMix2_REP1 | 0.82 |
| HBRMix2_REP2 | 0.94 |
| HBRMix2_REP3 | 0.93 |
Des boxplots sont souvent utilisés comme une mesure qualitative de la qualité de la normalisation, cela montre comment les distributions sont globalement affectés après le processus de normalisation. Il est attendu que ce processus de normalisation stabilise les distributions sur l’ensemble des échantillons. La Figure 7 montre les boxplots avant (à gauche) et après normalisation (à droite) des données.
Figure 7 - Boxplots des comptages bruts (gauche) et normalises (droite)
Comparaison : UHR vs HBR
Tests statistiques pour l’expression différentielle
Une fois l’estimation de la dispersion et l’ajustement du modèle sont effectués, DESeq2 peut faire les tests statistiques. La Figure 8 montre les distributions des p-valeurs ajustées, calculées par le test statistique pour la comparaison effectuée. Cette distribution devrait être un mélange d’une distribution uniforme sur [0,1] et d’un pic autour de 0 correspondant aux features significatifs exprimées de manière différentielle.
L’ajustement des p-valeurs est effectué pour prendre en compte les tests multiples et contrôler le taux de faux positifs à un niveau α choisi. Pour cette analyse, un ajustement de la p-valeur a été effectué avec la méthode Banjamini-Hochberg (BH) et le niveau du taux de faux positifs α a été fixé à 0,05.
Le Tableau 4 montre le nombre de features différentiellement exprimés en fonction des seuils α correspondant. Aucun filtre sur le log2FoldChange n’a été appliqué, tous les features avec des log2FoldChange < 0 ou > 0 sont comptés dans ce tableau.| 0.05 | 0.01 | 0.005 | 0.001 | |
|---|---|---|---|---|
| LFC_Up | 139 | 100 | 90 | 62 |
| LFC_Down | 150 | 107 | 95 | 67 |
Visualisation des résultats
Pour chaque comparaison, nous avons généré un tableau dynamique lié à un volcano plot et un MA-plot
A gauche, le volcano plot pour la comparaison réalisée et les features différentiellements exprimés sont mis en surbrillance (bleu pour les down-regulated et rouge pour les up-regulated) are still highlighted in red. Un volcano plot représente le log de la p-valeur ajustée en fonction du log ratio de la différence d’expression.
A droite, le MA-plot pour la comparaison réalisée avec les mêmes critères de surbrillance que pour les volcano plots. Un MA-plot représente le log ratio de la différence d’expresion en fonction de la moyenne d’expression de chaque feature.
Le tableau contient les informations importantes pour la comparaison réalisée. Voir les fichiers complets correspondant
Note: comment lire un volcano-plot et un MA-plot ?
Un volcano plot est un type de diagramme de dispersion qui montre la signification statistique (p-valeur) par rapport à l’ampleur du changement (fold change). Il permet d’identifier rapidement et visuellement les gènes dont les fold change sont importants et statistiquement significatifs. Il peut s’agir des gènes les plus importants sur le plan biologique. De plus, les gènes les plus régulés vers le haut (up_regulated) sont à droite, les gènes les plus régulés vers le bas (down-regulated) sont à gauche et les gènes les plus significatifs sur le plan statistique sont en haut du graphique.
Le MA-plot est aussi un type de graphe fréquemment utilisé pour comparer deux échantillons, une fois les comptes obtenus. La lettre M représente les log-ratios, et la lettre A est pour “Average” (la moyenne). Comme son nom l’indique, le MA-plot est un plot de la moyenne (en x) contre le ratio (en y) de deux valeurs: l’expression d’un gène dans deux conditions, chaque gène étant un point du graphe. Il s’interprète comme suit: plus un point est à droite du graphe, plus le gène est globalement fortement exprimé, et inversement. Plus un point est en haut du graphe, plus il est surexprimé dans la première condition par rapport à la seconde, et inversement. Les gènes “intéressants” se trouvent donc plutôt vers la droite (ratios plus convaincants si les comptes son élevés) et le plus possible en haut ou en bas du graphe, visiblement séparés de la masse des autres points.
Note: interprétation des noms de colonnes du tableau
- Id : identifiant unique du feature
- baseMean : moyenne de base sur tous les échantillons
- Test et Ref : moyennes (arrondies) des comptes normalisés des conditions biologiques
- log2FoldChange : log2(FC) tel qu’estimé par le modèle GLM. Il
reflète l’expression différentielle entre Test et Ref et peut être
interprété comme log2(Test/Ref)
- valeur est autour de 0 : l’expression du feature est similaire dans les deux conditions
- positive : la caractéristique est régulée vers le haut (Test>Ref)
- négative : la caractéristique est régulée à la baisse (Test<Ref)
- pvalue : p-valeur brute du test statistique
- padj : p-valeur ajustée sur laquelle le cut-off α est appliqué
- sig : statut du feature : significatif (up ou down) ou non significatif (not sig)
Ici, seule la première comparaison est affiché dans le rapport. Les résultats des autres comparaisons sont disponibles au format HTML, chaque nom de fichier contient le nom de la comparaison correspondante. Les comparaisons réalisées et les fichiers HTML correspondant, contenant les volcano-plots, MA-plots et les tableaux d’expression différentielle sont les suivants : UHR_vs_HBR.html.
Conclusion
L’analyse d’expression différentielle a permis de rechercher les gènes différentiellement exprimés entre vos différentes conditions. Dans votre cas, nous avons trouvés des différences d’expression significatives votre comparaison UHR vs HBR avec 289 gènes différentiellement exprimés lorsque l’on choisit un seuil de faux positif de 5% et 129 gènes différentiellement exprimés lorsque l’on se place au seuil de 0.1%. Nous vous conseillons donc dans un premier temps d’examiner les résultats obtenus avec le filtre de p-value ajustée de 0.001 (seuil de faux positif de 0.1%) qui est le plus stringeant afin de sélectionner les gènes dont le niveau de modulation vous intéresse ou correspond à vos attentes ou hypothèses scientifiques.
Méthodologie détaillée
Contrôle qualité
La principale variabilité au sein de l’expérience devrait provenir des différences biologiques entre les échantillons. Ceci peut être vérifié de deux manières. La première consiste à effectuer un regroupement hiérarchique de l’ensemble des échantillons. Ceci est effectué après une transformation des données de comptage avec une transformation stabilisatrice de la variance (VST).
Une VST est une transformation des données qui les rend homoscédastiques, ce qui signifie que la variance est alors indépendante de la moyenne. Elle s’effectue en deux étapes : (i) une relation moyenne-variance est estimée à partir des données avec la même fonction qui est utilisée pour normaliser les données de comptage et (ii) à partir de cette relation, une transformation des données est effectuée afin d’obtenir un ensemble de données dans lequel la variance est indépendante de la moyenne. L’homoscédasticité est une condition préalable à l’utilisation de certaines méthodes d’analyse des données, telles que le regroupement hiérarchique ou l’analyse en composantes principales (ACP).
Analyse différentielle
Modélisation
DESeq2 vise à ajuster un modèle linéaire par feature. Pour ce projet, le modèle utilisé est celui des comptages ~ conditions et le but est d’estimer les coefficients des modèles qui peuvent être interprétés comme le log2(FC). A partir de ces coefficients, des tests statistiques sont ensuite réalisés pour obtenir des p-valeurs.
Détection des outliers
Les valeurs aberrantes (outliers) du modèle sont des features pour lesquels au moins un échantillon ne semble pas lié au plan d’expérience. Pour chaque feature et pour chaque échantillon, la distance de Cook reflète la façon dont l’échantillon correspond au modèle. Une valeur élevée de la distance de Cook indique un nombre aberrant et les p-valeurs ne sont pas calculées pour le feature correspondant.
Estimation de la dispersion
Le modèle DESeq2 suppose que les données de comptage suivent une distribution binomiale négative qui est une alternative robuste à la loi de Poisson lorsque les données sont trop dispersées (la variance est supérieure à la moyenne). La première étape de la procédure statistique consiste à estimer la dispersion des données. Son but est de déterminer la forme de la relation moyenne-variance. Par défaut, on applique une méthode basée sur le GLM (Generalized Linear Model). Ensuite, DESeq2 impose une maximisation de la vraisemblance de profil ajustée par la méthode de Cox Reid [McCarthy, 2012] et utilise le maximum a posteriori (MAP) de la dispersion [Wu, 2013].
Filtre indépendant
DESeq2 peut effectuer un filtrage indépendant afin d’augmenter le pouvoir de détection des features exprimées de manière différentielle avec la même erreur de type I à l’échelle de l’expérience. Étant donné que les features présentant des comptes très faibles ne sont pas susceptibles de présenter des différences significatives en raison d’une dispersion élevée, il définit un seuil sur la moyenne des comptes normalisés, indépendamment de la condition biologique. Les p-valeurs ajustées des features écartés sont ensuite fixées à NA.
1 rue du Pr. Calmette,
59000 LILLE, France
+33 (0) 3 62 26 37 77
bioinfo@genoscreen.fr
© GenoScreen 2023