Statistiques pour l’imagerie cellulaire
Proposition de résolution des exercices
Avril 2019
Introduction
Nous vous proposons dans ce document des éléments de résolution des exercices que nous avons vus ensemble. Pour les intéressé-e-s, cette solution a été écrite avec le langage R Mardown : le code source se trouve dans le fichier Correction.Rmd, que nous avons placé dans un dossier compressé contenant également les données nécessaires à son fonctionnement. Attention, pour pouvoir exécuter ce code, il vous faudra peut-être installer des librairies supplémentaires. Si vous n’êtes pas intéressé par le langage R Markdown, vous pouvez simplement lire la correction.
1. COMPARAISON DE LA LONGUEUR DU SEGMENT INITIAL DE L’AXONE
On veut comparer des longueurs d’AIS entre deux traitements : le traitement d’intérêt et une situation contrôle. On suppose que la longueur des AIS est influencée par ce traitement, mais nous voudrions le vérifier au moyen d’un test statistique. Les données se trouvent dans le fichier 01_AIS.tsv.
Fig. 1 : Mesure des AIS sur une coupe de cervelet
Les données peuvent être visualisées sous la forme d’histogrammes et de diagrammes de dispersion. On voit que la distribution empirique des longueurs d’AIS n’est pas symétrique, et le nombre de mesures par groupe est faible (<30). Nous appliquons donc un test non-paramétrique de Wilcoxon, dont le rapport est présenté ci-dessous.
##
## Wilcoxon rank sum test
##
## data: AIS.Length by Groupe
## W = 4, p-value = 0.02597
## alternative hypothesis: true location shift is not equal to 0La conclusion de ce test est qu’il y a bien une différence significative entre les deux groupes.
2. DENSITE DE SYNAPSES
Dans cette expérience, on compte le nombre de synapses le long de la membrane cellulaire et on en déduit une densité par unité de longueur de membrane. Les données se trouvent dans le fichier 02_SYNAPSE.tsv.
Fig. 2 : Comptage de synapses autour d’un corps cellulaire
Cette fois-ci encore, les distributions ne sont pas symétriques. Nous appliquons donc préférentiellement un test de Wilcoxon.
##
## Wilcoxon rank sum test
##
## data: Dens by Type
## W = 63, p-value = 0.0001748
## alternative hypothesis: true location shift is not equal to 0Le rapport du test de Wilcoxon nous indique que la différence entre les deux types est significative. Pour les besoins de l’exercice, nous allons comparer ces résultats avec ceux donnés par un test de Student :
##
## Welch Two Sample t-test
##
## data: Dens by Type
## t = 7.2166, df = 13.972, p-value = 4.501e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 2.312790 4.269418
## sample estimates:
## mean in group Type 1 mean in group Type 2
## 15.19013 11.89903On remarque que le test choisi par défaut par la fonction t.test est un test dit de Welch, qui est adapté à une situation d’hétéroscédasticité (i.e. de variances différentes entre les groupes). Le résultat de ce test est bien sûr significatif mais il est intéressant de constater que la p-valeur est d’un ordre de grandeur inférieur à la p-valeur du test de Wilcoxon, ce qui est souvent le cas lorsque les échantillons sont petits.
3. COMPTAGE DE NEURONES
Nous prenons un peu de recul avec ce nouvel exemple dans lequel ce sont des neurones qui sont comptés dans une zone d’intérêt sur une coupe de cerveau. Les données se trouvent dans le fichier 03_NEURONES.tsv.
Fig. 3 : Comptage de neurones dans une région d’intérêt
La distribution des données, comme l’illustre la figure ci-dessous, n’est pas non plus symétrique, mais la taille de l’échantillon (environ 30 comptages) nous permet sans souci d’appliquer directement un test de Student.
##
## Welch Two Sample t-test
##
## data: Compte by Type
## t = -39.224, df = 56.991, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -109.42617 -98.79605
## sample estimates:
## mean in group Type 1 mean in group Type 2
## 147.4444 251.5556Le résultat du test est, sans surprise par rapport au diagramme de dispersion, que la différence entre les groupes est (très !) significative.
4. POSITION DE VESICULES PROCHES DE SYNAPSES
Ici encore, on s’intéresse à comparer deux échantillons de populations différentes, mais dans cette expérience, nous allons faire plusieurs mesures par objet, nous avons donc des données répétées. Les données se trouvent dans le fichier 04_VESICULES.tsv.
Fig. 4 : Distance entre les vésicules au centre et la synapse, en haut à droite de l’image
Dans cet exemple, il est très important de se rendre compte qu’il y a plusieurs objets mesurés par “individu”, ici, par image. On peut voir cet effet “individu” sur le graphe ci-dessous.
Si on ne prend pas en compte cet effet, par exemple en appliquant directement un test de Student, on obtient une p-valeur artificiellement petite (cf. ci-dessous) !
##
## Welch Two Sample t-test
##
## data: DIST by Type
## t = 6.3004, df = 377.93, p-value = 8.266e-10
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.5041024 0.9614930
## sample estimates:
## mean in group Type 1 mean in group Type 2
## 11.43525 10.70246Une façon de faire, qui n’est pas optimale, pourrait être de calculer les moyennes par individu et de comparer les deux groupes de moyennes à l’aide d’un test de Student. Cette approche n’est pas correcte au sens où on diminue la reproductibilité des résultats et on risque de comparer des grandeurs qui ne sont pas vraiment comparables en mélangeant des moyennes qui sont calculées sur des nombres d’objets différents. Ceci étant, le résultat du test, présenté ci-dessous, donne déjà une idée de la différence entre les deux groupes d’images.
##
## Welch Two Sample t-test
##
## data: DIST by Type
## t = 3.0395, df = 16.096, p-value = 0.007766
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 0.2192938 1.2286920
## sample estimates:
## mean in group Type 1 mean in group Type 2
## 11.42592 10.70192La procédure correcte est d’utiliser un modèle linéaire à effets mixtes. Le résultat de cette analyse est présenté pour information ci-dessous. On remarque que le résultat de cette analyse est plus proche de celui de l’analyse sur les moyennes que de celui de la première analyse.
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DIST ~ Type + (1 | ID)
## Data: vesicule
##
## REML criterion at convergence: 1167.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.63444 -0.70174 -0.06745 0.72600 2.63047
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.2249 0.4742
## Residual 1.1046 1.0510
## Number of obs: 386, groups: ID, 20
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 11.4274 0.1681 17.9265 67.960 < 2e-16 ***
## TypeType 2 -0.7256 0.2378 17.9169 -3.052 0.00689 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## TypeType 2 -0.707Dans un cas comme celui-ci, nous vous conseillons de vous rapprocher d’un-e statisticien-ne pour vérifier la validité de vos analyses et choisir la méthode la plus adaptée.
5. RECOLONISATION
Après avoir séparé artificiellement une couche cellulaire à confluence (“scratch”), on mesure la vitesse de recolonisation sous plus de deux conditions différentes. Les données se trouvent dans le fichier 05_RECOLO.tsv.
Fig. 5 : Recolonisation d’une zone après “scratch”
Dans cet exercice, il s’agit de comparer les vitesses moyennes de recolonisation dans 3 conditions différentes. Une représentation adaptée des données est présentée ci-après.
On peut constater qu’une vitesse négative apparaît dans le Type 2, il s’agit vraisemblablement d’une valeur aberrante. En pratique, il faudrait enquêter pour comprendre d’où vient cette vitesse négative, pour pouvoir décider s’il faut la supprimer. Ici, une solution est d’appliquer un test de Kruskal-Wallis suivi d’un test de Mann-Whitney-Wilcoxon pairwise pour les comparaisons multiples. En effet, cette approche non paramétrique est plus robuste en présence de valeurs aberrantes que son équivalent paramétrique (analyse de variance) et permet de traiter l’ensemble du jeu de données sans avoir à recourir à la suppression de cette observation.
##
## Kruskal-Wallis rank sum test
##
## data: Vitesse by Type
## Kruskal-Wallis chi-squared = 6.284, df = 2, p-value = 0.0432Le test de Kruskal-Wallis, basé sur les rangs, est significatif et conclut qu’il y a une différence significative entre les 3 types. Ensuite, les comparaisons pairwise permettent de savoir quelles moyennes sont significativement différentes 2 à 2.
##
## Pairwise comparisons using Wilcoxon rank sum test
##
## data: recolo$Vitesse and recolo$Type
##
## Type 1 Type 2
## Type 2 0.586 -
## Type 3 0.090 0.085
##
## P value adjustment method: holmOn peut remarquer que les p-valeurs des comparaisons multiples sont ajustées sur la multiplicité des tests, car 3 comparaisons sont réalisées. Ici, on a utilisé la méthode de correction de holm, qui est choisie par défaut par la fonction pairwise.wilcox.test. Sur la sortie ci-dessus, on voit qu’aucune des comparaisons 2 à 2 n’est significative, alors que le test global l’est. Cette situation peut se rencontrer lorsque les différences entre les groupes sont peu marquées. Le test global de Kruskal-Wallis étant plus puissant que le test de Mann-Whitney-Wilcoxon, celui-ci peut détecter une différence globale entre plusieurs moyennes, alors que les comparaisons multiples ne permettront pas de conclure.
6. COLOCALISATION
L’objectif est ici de localiser les objets doublement marqués et de quantifier dans chaque objet les deux marquages :
- Marquage des noyaux au DAPI (bleu),
- Histone H3 phosphorylée (rouge).
Les données se trouvent dans le fichier 06_COLOC.tsv.
Fig. 6 : Mesure de fluorescence en vue de calculer la colocalisation
Ici, on cherche à savoir si les objets marqués en bleu le sont aussi en rouge. On cherche donc à étudier la corrélation entre les 2 marquages. Sur le nuage de points ci-après (figure de gauche), on constate que de nombreux objets ne sont pas marqués en rouge alors qu’ils le sont en bleu. On peut déjà constater visuellement que le lien entre les deux marquages n’est pas franc, voire supposer qu’il est absent sur ces données. Cependant, on peut calculer un coefficient de corrélation pour le nuage de points, après filtrage des points non marqués en rouge (figure de droite).
##
## Pearson's product-moment correlation
##
## data: Blue and Red
## t = 1.2598, df = 78, p-value = 0.2115
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.08101427 0.35007132
## sample estimates:
## cor
## 0.1412163Sur la sortie ci-dessus, on voit que la corrélation de Pearson entre les marquages bleu et rouge, après filtrage des objets non marqués, est de 0.14. Celle-ci n’est pas significative. On ne peut donc pas conclure sur l’existence d’une corrélation entre ces 2 marquages.
9. Translocation nucleo-cytoplasmique de NFkB. Effet dose du TNF une étude en imagerie HTS
L’objectif de l’étude est de déterminer si la réponse au TNF est différente en fonction du type cellulaire (MCF7 ou A549).
La quantité de NF\(\kappa\)B adressée au noyau dépendant de la concentration de TNF, on a mesuré sa fluorescence dans les deux types cellulaires et on va essayer de déterminer s’il y a une différence ou même un effet d’interaction.
Les données se trouvent dans le fichier 09_NFKB.tsv.
Fig. 7 : Effet dose du TNF.
Une représentation adaptée des données est la suivante (figure de gauche).
Ici, on remarque sur la représentation des données que l’effet de la dose semble varier en fonction du type cellulaire, ce qui est signe d’une potentielle interaction entre les 2 facteurs (dose et type cellulaire) dont il va falloir tester l’existence. Une méthode d’analyse adaptée est d’appliquer une analyse de variance à 2 facteurs avec interaction. En revanche, on remarque que la variance est beaucoup plus faible dans la dose 9 que dans la dose 2, ce qui implique que l’hypothèse d’homoscédasticité (variances égales entre les différents groupes) de l’ANOVA n’est pas respectée. Une solution est de réaliser l’ANOVA en utilisant la fonction lm appliquée aux données transformées en log (voir figure de droite, des données transformées). Cette transformation, très souvent utilisée en biologie, permet de stabiliser la variance entre les différents groupes. Pour réaliser les tests de Fisher des effets globaux du modèle, on peut utiliser la fonction Anova du package car.
## Anova Table (Type III tests)
##
## Response: log(1 + NCRatio)
## Sum Sq Df F value Pr(>F)
## (Intercept) 6.3724 1 5490.4139 < 2.2e-16 ***
## Dose 0.5547 1 477.9007 4.924e-11 ***
## G 0.0266 1 22.9575 0.0004399 ***
## Dose:G 0.0027 1 2.2995 0.1553039
## Residuals 0.0139 12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sur cette sortie, on voit que l’interaction n’est pas significative. Comme on s’intéresse à toutes les comparaisons entre groupes, on choisit de la laisser dans le modèle. On voit aussi que l’effet associé à la dose et l’effet associé au type cellulaire sont significatifs, ce qui signifie que ces 2 facteurs ont un impact sur le NCRatio. Ensuite, on utilise la fonction emmeans du package du même nom pour réaliser les comparaisons multiples entre les différents groupes.
## $`emmeans of Dose, G`
## Dose G emmean SE df lower.CL upper.CL
## 2 A549 1.262 0.017 12 1.225 1.299
## 9 A549 0.736 0.017 12 0.698 0.773
## 2 MCF7 1.147 0.017 12 1.110 1.184
## 9 MCF7 0.672 0.017 12 0.635 0.709
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $`pairwise differences of Dose, G`
## contrast estimate SE df t.ratio p.value
## 2,A549 - 9,A549 0.5266 0.0241 12 21.861 <.0001
## 2,A549 - 2,MCF7 0.1154 0.0241 12 4.791 0.0021
## 2,A549 - 9,MCF7 0.5904 0.0241 12 24.508 <.0001
## 9,A549 - 2,MCF7 -0.4112 0.0241 12 -17.070 <.0001
## 9,A549 - 9,MCF7 0.0638 0.0241 12 2.647 0.0866
## 2,MCF7 - 9,MCF7 0.4750 0.0241 12 19.716 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 4 estimatesOn voit que toutes les comparaisons sont significatives, à l’exception de la comparaison entre A549 et MCF7 pour la dose 9. On note que là encore, les p-valeurs pour ces comparaisons multiples ont été corrigées sur la multiplicité des tests, par la méthode de Tukey.
Bonus - fluorescence intracellulaire
On considère des cellules en culture portant un récepteur membranaire. On enregistre toutes les 5 minutes la fluorescence de chaque cellule pendant 80 minutes. Pour cet exercice, on considère les données pour 1 cellule par condition, elles sont disponibles dans le fichier 10_BONUS_FLUO.tsv. On considère 2 conditions expérimentales, dans la 1ère, un ligand (fluorescent) du récepteur est ajouté au milieu de culture. Dans la seconde condition, on ajoute au ligand un antagoniste qui bloque l’accès du ligand à son récepteur. L’hypothèse est la suivante : dans la 1ere condition (LF), le ligand serait internalisé donc la fluorescence augmenterait dans la cellule ; dans la seconde condition (LFA), l’antagoniste ajouté bloquerait l’internalisation du ligand et la fluorescence augmenterait moins.
Remarque cette expérience a été répliquée pour 1 seule cellule par condition. Les conclusions de cet exercice ne sont donc valables que pour les 2 cellules utilisées lors de l’expérience. Pour pouvoir étendre ces conclusions, il faudrait répliquer l’expérience sur un grand nombre de cellules. On aurait alors un design imbriqué et il faudrait utiliser des approches statistiques plus élaborées (modèle linéaire à effet aléatoire).
Fig. 8 : Fluorescence en fonction du temps, dans 2 conditions différentes
Une représentation adaptée des données est la suivante :
Pour répondre à la question biologique, une solution est d’estimer une droite de régression dans chaque condition expérimentale (voir figure ci-après). Si l’antogoniste a un effet sur la fluorescence de la cellule, alors les pentes dans les 2 conditions expérimentales seront significativement différentes. On peut déjà le suspecter à l’oeil nu grâce à la représentation des données : la pente de la droite associée à la condition ligand seul semble plus importante que celle associée à la condition ligand + antagoniste.
Une manière de procéder est d’utiliser la fonction lm() ajustée sur les effets de la condition (LF/LFA) et du temps puis y ajouter une interaction entre le temps et la condition. Cette interaction permet d’estimer une droite de régression pour chaque condition (voir figure suivante) dans laquelle la pente de la droite dépend de la condition expérimentale. Pour conclure qu’il y a un effet de l’antagoniste sur la fluorescence, on peut alors tester la nullité du coefficient de cette interaction.
Sur la sortie suivante, on voit que:
Le terme
(Intercept), colonneEstimatecorrespond à l’ordonnée à l’origine de la droite de régression dans la condition LF. Ce terme est significatif (p-valeur < 2e-16) ce qui signifie qu’au temps 0, la fluorescence n’est pas nulleLe terme
ConditionLFA, colonneEstimatecorrespond à l’écart d’ordonnée à l’origine entre la condition LF et LFA. Ce paramètre est négatif: on voit que l’ordonnée à l’origine dans la condition LFA est un peu inférieure à celle de la condition LF. Ce coefficient n’est pas significativement différent de 0, ce qui laisse penser que les ordonnées à l’origine entre les 2 conditions ne sont pas significativement différentes.Le terme
Lapse, colonneEstimatecorrespond à la pente de la droite dans la condition LF, ce coefficient est significatif (et positif). Cela signifie que la fluorescence augmente avec le temps dans la condition LFEnfin, le terme d’interaction
ConditionLFA:Lapsecorrespond à la différence de pente entre les 2 conditions. Ce coefficient est significatif. Cela implique que l’effet du temps est différent dans la condition LF et dans la condition LFA. Autrement dit, l’ajout de l’antagoniste modifie bien la variation de la fluorescence au cours du temps. On voit même que l’antagoniste ralentit la croissance de la fluorescence (coefficient de la colonneEstimate = -10749.1négatif).Ensuite, on voit tout en bas que le test de Fisher est significatif (p-valeur < 2.2e-16). Ceci signifie que le modèle explique une part significative de l’information contenue dans les données. Le R\(^2\) associé à ce modèle est de 0.98 (part de variance expliquée par le modèle).
##
## Call:
## lm(formula = Fluo ~ Condition + Lapse + Condition:Lapse, data = fluo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -108937 -22342 8237 29592 117069
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2314231.3 28413.3 81.449 < 2e-16 ***
## ConditionLFA -24461.6 40182.4 -0.609 0.548
## Lapse 14645.9 587.7 24.921 < 2e-16 ***
## ConditionLFA:Lapse -10749.1 831.1 -12.933 2.49e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 54180 on 28 degrees of freedom
## Multiple R-squared: 0.9789, Adjusted R-squared: 0.9766
## F-statistic: 432.1 on 3 and 28 DF, p-value: < 2.2e-16