library(tidyverse)
library(kableExtra)
library(ggvenn)
kabex = function(o) o %>% kable(format="html", escape=F) %>% kable_styling(bootstrap_options = c("striped"))
var.n = function(x) 1/length(x) * sum( (x - mean(x))**2 )
t.stat = function(m1, m2, n1, n2, v1, v2) (m1-m2)/sqrt( (n1*v1 + n2*v2)  * (1/n1 + 1/n2) ) * sqrt(n1+n2-2)
rounding=2
correction=T
correction.start=function() ifelse(correction, '<font style="color: #F00000;">', "<!--")
correction.end=function() ifelse(correction, '</font>', "-->")

Contexte

La présence du marqueur (gène FLS2) induit une cascade de réactions résultant en réponses immunitaires de la plante Arabidopsis thaliana. Ici, vous voulez comparer le transcriptome entre plantes avec marqueur et plantes sans marqueur. Pour cela, vous marquez le transcriptome des plantes avec marqueur avec le fluorochrome vert (cy3) et le transcriptome des plantes sans marqueur avec le fluorochrome rouge (cy5). Vous disposez de 10 réplicats.

Test de l’expression différentielle d’un gène

Vous voulez tester l’expression différentielle d’un gène de réponse immunitaire. Pour cela, vous disposez des valeurs de sonde suivantes :

replicat=1:10
cy3=c(10, 30, 42, 12, 30, 47, 25, 27, 26, 41)
cy5=c(25, 47, 60, 47, 77, 56, 81, 70, 56, 68)
tibble(Réplicat=replicat, `Cy3 (G)`=cy3, `Cy5 (R)`=cy5) %>% t %>% kabex
Réplicat 1 2 3 4 5 6 7 8 9 10
Cy3 (G) 10 30 42 12 30 47 25 27 26 41
Cy5 (R) 25 47 60 47 77 56 81 70 56 68

Vous devez donc comparer (statistiquement) la moyenne de l’expression du gène des plantes avec et sans marqueur afin de déterminer si les moyennes sont différentes et donc de déterminer si le marqueur a un effet (significatif).

Elements de statistiques

On suppose que la population peut se résumer à (ou modéliser avec) une moyenne et un écart-type.

Exemple avec une population de moyenne \(µ = 30\) et d’écart-type \(\sigma = 10\) tous les 2 connus.

xi=seq(0,60,by=.01)
plot(xi, dnorm(xi, mean=30, sd=10), type='l', xlab='x', ylab='p(x)', main='Loi normale N(30,10)')

Quand on prend des individus au hasard dans la population, on a plus de chance de tomber sur quelqu’un proche de la moyenne.

En pratique, il est rarement possible de mesurer tous les individus d’une population. L’idée est d’estimer la moyenne et l’écart-type de la population à partir d’un échantillon, c’est-à-dire un certain nombre d’individus (→ taille de l’échantillon) pris au hasard. Il sera essentiel de pouvoir statistiquement quantifier la confiance que l’on peut accorder à ces estimation.

Moyenne de l’échantillon : \[\overline{X} = \frac{1}{n}\sum_{i=1}^{n} x_i\] avec \(n\) la taille de l’échantillon.

Ecart-type de l’échantillon noté \(S\) et variance \(S^2\) : \[S^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \overline{X} )^2\]

3 échantillons de taille 10 pour la loi N(moyenne=30, écart-type=10)

par(mfrow=c(3,1))
sample.size = 10
smp = round(rnorm(sample.size,30,10))
stripchart(smp, method="stack", xlim=c(5,55), frame.plot=F, axes=F)
axis(1, col=NA, col.ticks=1)
abline(v=mean(smp), lty="dashed")
text(mean(smp), 0, paste0("Moyenne de l'échantillon : ", round(mean(smp),rounding)))
smp = round(rnorm(sample.size,30,10))
stripchart(smp, method="stack", xlim=c(5,55), frame.plot=F, axes=F)
axis(1, col=NA, col.ticks=1)
abline(v=mean(smp), lty="dashed")
text(mean(smp), 0, paste0("Moyenne de l'échantillon : ", round(mean(smp),rounding)))
smp = round(rnorm(sample.size,30,10))
stripchart(smp, method="stack", xlim=c(5,55), frame.plot=F, axes=F)
axis(1, col=NA, col.ticks=1)
abline(v=mean(smp), lty="dashed")
text(mean(smp), 0, paste0("Moyenne de l'échantillon : ", round(mean(smp),rounding)))

2 échantillons de taille 100 pour la loi N(moyenne=30, écart-type=10)

par(mfrow=c(2,1))
sample.size = 100
smp = round(rnorm(sample.size,30,10))
stripchart(smp, method="stack", xlim=c(5,55), frame.plot=F, axes=F)
axis(1, col=NA, col.ticks=1)
abline(v=mean(smp), lty="dashed")
text(mean(smp), 0, paste0("Moyenne de l'échantillon : ", round(mean(smp),rounding)))
smp = round(rnorm(sample.size,30,10))
stripchart(smp, method="stack", xlim=c(5,55), frame.plot=F, axes=F)
axis(1, col=NA, col.ticks=1)
abline(v=mean(smp), lty="dashed")
text(mean(smp), 0, paste0("Moyenne de l'échantillon : ", round(mean(smp),rounding)))

2 échantillons de taille 1000 pour la loi N(moyenne=30, écart-type=10)

stripchart(round(rnorm(1000,30,10)), method="stack", xlim=c(5,55), ylim=c(0,15), frame.plot = F)

stripchart(round(rnorm(1000,30,10)), method="stack", xlim=c(5,55), ylim=c(0,15), frame.plot = F)

On remarque que plus la taille de l’échantillon est grande et plus on se rapproche de l’allure de la population.

Il en est de même pour la moyenne de l’échantillon par rapport à la moyenne de la population → plus l’échantillon est gros et plus sa moyenne est proche de la moyenne de la population (dans la plupart des cas, mais on peut tomber sur un échantillon peu représentatif de la population).

Moyennes de 10 échantillons de taille 10 pour N(30,10) :

r = 1:10
sample.size=10
samples = lapply(r, function(i) rnorm(sample.size,30,10))
sample.means = sapply(samples, mean ) %>% round(rounding)
sample.sds =  sapply(samples, sd ) %>% round(rounding)
tibble(moyenne=sample.means, 
       `écart-type`=sample.sds, 
       `erreur sur la moyenne`= sample.means-30, 
       ESM=round(sample.sds/sqrt(sample.size),rounding)
       ) %>% t %>% kabex
moyenne 29.56 32.27 27.09 34.20 29.02 31.54 32.10 21.36 30.62 30.47
écart-type 10.15 6.26 8.32 10.03 9.37 6.85 8.95 10.97 9.31 8.80
erreur sur la moyenne -0.44 2.27 -2.91 4.20 -0.98 1.54 2.10 -8.64 0.62 0.47
ESM 3.21 1.98 2.63 3.17 2.96 2.17 2.83 3.47 2.94 2.78

On observe bien que l’on a pas obtenu la moyenne de la population de 30, ni son écart-type de 10. La moyenne des erreurs obtenues est de -0.18.

Moyennes de 10 échantillons de taille 100 pour N(30,10) :

sample.size=100
samples = lapply(r, function(i) rnorm(sample.size,30,10))
sample.means = sapply(samples, mean ) %>% round(rounding)
sample.sds =  sapply(samples, sd ) %>% round(rounding)
tibble(moyenne=sample.means, 
       `écart-type`=sample.sds, 
       `erreur sur la moyenne`= sample.means-30, 
       ESM=round(sample.sds/sqrt(sample.size),rounding)
       ) %>% t %>% kabex
moyenne 30.09 28.07 30.72 28.42 30.10 29.91 29.66 30.81 29.27 30.61
écart-type 10.11 10.07 10.47 9.84 9.72 10.57 10.45 9.90 10.13 9.95
erreur sur la moyenne 0.09 -1.93 0.72 -1.58 0.10 -0.09 -0.34 0.81 -0.73 0.61
ESM 1.01 1.01 1.05 0.98 0.97 1.06 1.04 0.99 1.01 0.99

La moyenne des erreurs obtenues est de -0.23.

On observe donc une différence entre la moyenne de l’échantillon et la moyenne de la population. Cette différence se modélise (avec une loi normale) et dépend de la taille de l’échantillon → erreur standard à la moyenne : \[ESM = \frac{\sigma}{\sqrt{n}}\] Malheureusement, on ne connait pas toujours \(\sigma\) (de la population) et il est nécessaire d’utiliser son estimation \(S\) à partir de l’échantillon.

D’après cette modélisation, on a 2/3 ou 66% de chances que l’erreur de moyenne obtenue sur l’échantillon soit inférieure à \(ESM\) (et donc 1/3 de chances que l’erreur soit plus importante que \(ESM\)).

Moyennes de 10 échantillons de taille 1000 pour N(30,10) :

sample.size=1000
samples = lapply(r, function(i) rnorm(sample.size,30,10))
sample.means = sapply(samples, mean ) %>% round(rounding)
sample.sds =  sapply(samples, sd ) %>% round(rounding)
tibble(moyenne=sample.means, 
       `écart-type`=sample.sds, 
       `erreur sur la moyenne`= sample.means-30, 
       ESM=round(sample.sds/sqrt(sample.size),rounding)
       ) %>% t %>% kabex
moyenne 29.62 29.87 30.10 29.65 29.44 29.84 30.14 30.32 29.90 29.51
écart-type 9.87 9.82 10.02 10.16 10.29 10.27 9.98 10.28 10.13 9.85
erreur sur la moyenne -0.38 -0.13 0.10 -0.35 -0.56 -0.16 0.14 0.32 -0.10 -0.49
ESM 0.31 0.31 0.32 0.32 0.33 0.32 0.32 0.33 0.32 0.31

La moyenne des erreurs obtenues est de -0.16.

\(ESM\) nous donne donc une indication sur la différence entre la moyenne de l’échantillon et la moyenne de la population. En prenant 2 échantillons, on va pouvoir de la même manière modéliser la différence observée entre les moyennes des 2 échantillons (s’ils sont tirés d’une même population). Cette différence de moyenne permet de calculer la probabilité que les 2 échantillons soient tirés de la même population. La statistique calculée par le test de Student est :

\[T = \frac{ \text{différence entre les moyennes des 2 échantillons} } {\text{erreur standard de la différence de moyenne}} = \frac{DM}{ESDM}\] Cette statistique \(T\) suit une loi de Student dépendant de la taille des 2 échantillons.

Par exemple, pour 2 échantillons de taille 100, ayant pour moyennes 31.9 et 29.1 et pour écarts-types 9.5 et 10.3, \(T = 1.99\)

Avec 100 valeurs dans chaque échantillon, le nombre de degrés de liberté (noté ddl par la suite) est 198 (99 valeurs que l’on peut choisir comme on veut mais la dernière est fixée pour “retomber” sur la moyenne observée, pour chacun des échantillons → 198).

La fonction de densité de la statistique \(T\) pour 198 ddl est :

xi = seq(-4,4,by=0.1)
plot(xi, dt(xi, df=198), type='l', xlab='T', ylab='p(T)')
thr = qt(.975, df=198)
foo=sapply(seq(-thr, thr, by=.1), function(xi) lines(c(xi,xi), c(0.005,dt(xi,df=198))) )
points(1.99, dt(1.99, df=198), pch=3, col='darkgreen', cex=2)
text(0,0, "intervalle de confiance à 95%")

On observe donc que cette différence de moyennes s’écarte de 0 mais n’est pas extrême.

D’après la fonction de densité, 95% de la population a une valeur de \(T\) comprise entre -1.97 et 1.97.

Ainsi avec un seuil de 5% (seuil \(\alpha\) ou risque de 1ère espèce), on peut décider que la différence entre les moyennes des échantillons de 31.9 et 29.1 est significative, c’est-a-dire que les 2 échantillons sont issus de populations n’ayant pas la même moyenne.

Le risque de 1ère espèce ou erreur de type I est la probabilité de décider \(H_1\) alors que c’est \(H_0\) qui est vraie.

Ici, pour \(T = 1.99\), la p-valeur du test 0.048, c’est-à-dire la probabilité de \(p(T\lt-1.99 \cup T\gt1.99)\).

Ici, au seuil de 5%, nous décidons \(H_1\) (les moyennes sont différentes), et ceci, à tort puisque nous avons tiré les 2 échantillons de la même population (de moyenne 30 et d’écart-type 10).

Application du test de Student

tibble(Réplicat=replicat, `Cy3 (G)`=cy3, `Cy5 (R)`=cy5) %>% t %>% kabex
Réplicat 1 2 3 4 5 6 7 8 9 10
Cy3 (G) 10 30 42 12 30 47 25 27 26 41
Cy5 (R) 25 47 60 47 77 56 81 70 56 68

Hypothèses du test :

\(H_0\): \(\overline{X}_G = \overline{X}_R\) (ou bien \(\overline{X}_G - \overline{X}_R = 0\))

\(H_1\): \(\overline{X}_G \ne \overline{X}_R\)

Pour ce faire, vous devez calculer :

\[T = \frac{ \text{différence entre les moyennes des 2 échantillons} } {\text{erreur standard de la différence de moyenne}} = \frac{\overline{X}_G - \overline{X}_R}{ESDM}\]

Pour la moyenne d’un échantillon de taille \(n\) :

\[\overline{X} = \frac{1}{n}\sum_{i=1}^{n} x_i\]

pour l’erreur standard de la différence de moyenne \(ESDM\) :

\[ESDM = \sqrt{ \frac{S_p^2}{n_G} + \frac{S_p^2}{n_R}} \] que l’on peut réécrire : \[ESDM = \sqrt{ S_p^2 (\frac{1}{n_G} + \frac{1}{n_R} )} = S_p \sqrt{ \frac{1}{n_G} + \frac{1}{n_R}}\]

et où \(S_p\) représente la variance groupée (pooled en anglais) des 2 échantillons : \[S_p = \sqrt{ \frac{ n_G S^2_G + n_R S^2_R }{n_G + n_R - 2} } = \frac{ \sqrt{ n_G S^2_G + n_R S^2_R } }{ \sqrt{n_G + n_R - 2} } \]

Variance d’un échantillon :

\[S^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \overline{X} )^2\]

Calculez \(\overline{X_G}\) et \(\overline{X_R}\).

\(\overline{X_G}\) = moyenne(10, 30, 42, 12, 30, 47, 25, 27, 26, 41) = 29

\(\overline{X_R}\) = moyenne(25, 47, 60, 47, 77, 56, 81, 70, 56, 68) = 58.7

Calculez \(S^2_G\) et \(S^2_R\).

\(S_G^2\) = variance(10, 30, 42, 12, 30, 47, 25, 27, 26, 41) = 131.8

\(S_R^2\) = variance(25, 47, 60, 47, 77, 56, 81, 70, 56, 68) = 247.21

Calculez la statistique \(T\) du test. La formule réécrite dans le détail est :

\[T = \frac{\overline{X_G}-\overline{X_R}} { \sqrt{(n_G S^2_G + n_R S^2_R) (\frac{1}{n_G} + \frac{1}{n_R}) } } \sqrt{n_G + n_R - 2}\]

\(T\) = -4.5766992

A l’aide de la table ci-dessous, déterminer les bornes de l’intervalle des valeurs de \(T\) pour 95%, autrement dit les seuils pour un risque \(\alpha\) de 5%.

ddl = 1:30
tibble(ddl=ddl, 
       `α = 0.01` = -round(qt(0.01 /2, ddl),3), 
       `α = 0.02` = -round(qt(0.02 /2, ddl),3),
       `α = 0.025`= -round(qt(0.025/2, ddl),3),
       `α = 0.05` = -round(qt(0.05 /2, ddl),3),
       `α = 0.1`  = -round(qt(0.1  /2, ddl),3)
) %>% kabex
ddl α = 0.01 α = 0.02 α = 0.025 α = 0.05 α = 0.1
1 63.657 31.821 25.452 12.706 6.314
2 9.925 6.965 6.205 4.303 2.920
3 5.841 4.541 4.177 3.182 2.353
4 4.604 3.747 3.495 2.776 2.132
5 4.032 3.365 3.163 2.571 2.015
6 3.707 3.143 2.969 2.447 1.943
7 3.499 2.998 2.841 2.365 1.895
8 3.355 2.896 2.752 2.306 1.860
9 3.250 2.821 2.685 2.262 1.833
10 3.169 2.764 2.634 2.228 1.812
11 3.106 2.718 2.593 2.201 1.796
12 3.055 2.681 2.560 2.179 1.782
13 3.012 2.650 2.533 2.160 1.771
14 2.977 2.624 2.510 2.145 1.761
15 2.947 2.602 2.490 2.131 1.753
16 2.921 2.583 2.473 2.120 1.746
17 2.898 2.567 2.458 2.110 1.740
18 2.878 2.552 2.445 2.101 1.734
19 2.861 2.539 2.433 2.093 1.729
20 2.845 2.528 2.423 2.086 1.725
21 2.831 2.518 2.414 2.080 1.721
22 2.819 2.508 2.405 2.074 1.717
23 2.807 2.500 2.398 2.069 1.714
24 2.797 2.492 2.391 2.064 1.711
25 2.787 2.485 2.385 2.060 1.708
26 2.779 2.479 2.379 2.056 1.706
27 2.771 2.473 2.373 2.052 1.703
28 2.763 2.467 2.368 2.048 1.701
29 2.756 2.462 2.364 2.045 1.699
30 2.750 2.457 2.360 2.042 1.697

Avec un risque \(\alpha\) de 5%, acceptez-vous \(H_0\) ? Que concluez-vous ?

Pour 18 ddl et \(\alpha = 0.05\), on obtient un seuil de 2.101.

Ici, \(T=-4.5767\) et est très inférieure au seuil (de -2.101) donc on décide \(H_1\).

En R :

Normalité des échantillons ?

Test de Shapiro-Wilk avec \(H_0\) : adéquation de l’échantillon à la loi normale.

shapiro.test(cy3)
## 
##  Shapiro-Wilk normality test
## 
## data:  cy3
## W = 0.93676, p-value = 0.5175
shapiro.test(cy5)
## 
##  Shapiro-Wilk normality test
## 
## data:  cy5
## W = 0.95391, p-value = 0.7149

Les p-valeurs sont supérieures à 0.05, on décide \(H_0\) pour les 2 échantillons.

Homogénéité des variances ?

Test d’homoscédasticité, \(H_0\): les variances des 2 échantillons sont homogènes

var.test(cy3,cy5)
## 
##  F test to compare two variances
## 
## data:  cy3 and cy5
## F = 0.53315, num df = 9, denom df = 9, p-value = 0.3626
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1324269 2.1464586
## sample estimates:
## ratio of variances 
##            0.53315

p-valeur > 0.05 → on décide \(H_0\), les variances sont homogènes.

Test de Student

t.test(cy3, cy5, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  cy3 and cy5
## t = -4.5767, df = 18, p-value = 0.0002339
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -43.33371 -16.06629
## sample estimates:
## mean of x mean of y 
##      29.0      58.7

p-valeur < \(\alpha\), on décide \(H_1\), les moyennes sont différentes → le gène est différentiellement exprimé.

Ajustement dans le cas de tests multiples

Elements de statistique

Lorsque l’on réalise un seul test, comme précédemment, on rejette \(H_0\) si p-valeur \(< \alpha\).

Comme vu dans l’introduction, il arrive que l’on rejette \(H_0\) alors qu’elle est vraie (erreur de type I). Avec un seuil \(\alpha\) de 5%, on a donc une probabilité de 0.05 de rejeter \(H_0\) alors qu’elle est vraie.

Le génome d’arabette comporte ~25000 gènes. Pour chaque gène, le test de Student est réalisé pour décider s’il est différentiellement exprimé ou non. Supposons que pour 100 gènes, la p-valeur obtenue est inférieure à \(\alpha=0.05\). Sur ces 100 gènes, on peut donc s’attendre à tomber, pour certains, dans le cas où \(H_0\) est rejetée alors qu’elle est vraie. En moyenne ou en général, on peut donc s’attendre sur 100 p-valeurs significatives (à 5%) à ce que l’on rejette \(H_0\) à tort (donc ici 5 erreurs de type I parmi les 100 p-valeurs).

Il a été élaboré différentes méthodes pour ajustement le seuil \(\alpha\) dans le cadre de tests multiples. Nous évoquerons ici la méthode de Bonferroni et la False Discovery Rate (FDR) de Benjamini et Hochberg (1995).

La méthode de Bonferroni consiste à diviser le seuil par le nombre de tests effectués. Pour l’exemple d’arabette, le seuil ajusté devient donc 0.05 / 25000 = 0.000002. Cette méthode est connue pour être trop conservatrice ou trop stringente, c’est-à-dire que le seuil ajusté permet effectivement de limiter le nombre de faux positifs, mais malheureusement, il ne permet plus de retrouver les vrais positifs.

La méthode FDR est une méthode moins conservatrice qui va permettre d’adapter le seuil en tenant compte des p-valeurs obtenues. Elle consiste à trier les \(m\) p-valeurs obtenues de manière croissante (\(P_1 < P_2 < \dots < P_{m-1} < P_m\)) et d’ajuster le seuil selon le rang \(k\) de la p-valeur : \(P_k < \alpha \frac{k}{m}\).

Application de la FDR

On réalise le test de Student précédent pour 10 gènes de réponse immunitaire (G1, G2, …, G10). Pour chaque test, on calcule une p-valeur.

genes = paste0('G',1:10)
pvalues = c(0.07, .001, .54, .12, .85, .021, .0055, .01, .35, .98)
tibble(`gène`=genes, `p-valeur`=as.character(pvalues)) %>% t %>% kabex
gène G1 G2 G3 G4 G5 G6 G7 G8 G9 G10
p-valeur 0.07 0.001 0.54 0.12 0.85 0.021 0.0055 0.01 0.35 0.98

Quels sont les gènes différentiellement exprimés (sans ajustement) ?

tibble(`p-valeur`=pvalues, significative=pvalues<.05) %>% kabex
p-valeur significative
0.0700 FALSE
0.0010 TRUE
0.5400 FALSE
0.1200 FALSE
0.8500 FALSE
0.0210 TRUE
0.0055 TRUE
0.0100 TRUE
0.3500 FALSE
0.9800 FALSE

Quels sont les gènes différentiellement exprimés avec l’ajustement de la FDR ?

ps = sort(pvalues)
p.adj = round(ps * 10 / 1:10, 3)
as = 1:10*.05/10
tibble(k=1:10, `p-valeur`=ps, `seuil ajusté : 0.05 * k/m`=as, significative=ps<as, `p-valeur ajustée : Pk * m/k`=p.adj, `p.adj < 0.05`=p.adj<0.05) %>% kabex
k p-valeur seuil ajusté : 0.05 * k/m significative p-valeur ajustée : Pk * m/k p.adj < 0.05
1 0.0010 0.005 TRUE 0.010 TRUE
2 0.0055 0.010 TRUE 0.027 TRUE
3 0.0100 0.015 TRUE 0.033 TRUE
4 0.0210 0.020 FALSE 0.053 FALSE
5 0.0700 0.025 FALSE 0.140 FALSE
6 0.1200 0.030 FALSE 0.200 FALSE
7 0.3500 0.035 FALSE 0.500 FALSE
8 0.5400 0.040 FALSE 0.675 FALSE
9 0.8500 0.045 FALSE 0.944 FALSE
10 0.9800 0.050 FALSE 0.980 FALSE

Quel test utiliserait-on si on avait plus que 2 conditions à comparer ? et quelles sont les hypothèses de ce test ?

Test ANOVA (ANalysis Of VAriance)

\(H_0\): les moyennes des différents groupes sont égales

\(H_1\): au moins une des moyennes est différente

Le principe est de comparer la variance au sein des groupes à la variance entre les groupes. Si l’écart est trop grand, il doit y avoir au moins un des groupes pour lequel la variance diverge trop (et donc la moyenne).

Caractérisation d’un ensemble de gènes

Le test de Student est réalisé pour les 25000 gènes du génome d’A. thaliana. Seulement 50 gènes sont déclarés différentiellement exprimés après correction FDR.

On cherche à déterminer rapidement si les 50 gènes obtenus correspondent à des fonctions biologiques particulières.

Le principe va être de comparer l’ensemble des 50 gènes différentiellement exprimés à chacun des ensembles de gènes correspondantà une fonction biologique particulière.

Prenons un exemple fictif pour la fonction biologique respiration pour laquelle on aurait 250 gènes annotés avec ce terme de la Gene Ontology dans le génome d’arabette (soit 1% du génome). Serait-il surprenant d’en avoir 10 annotés respiration parmi les 50 qui sont différentiellement exprimés ?

sets = list(
  #`génome`= paste0('gene', 1:25000),
  `diff. exprimés` = paste0('gene', 1:50),
  respiration = paste0('gene', 41:290)
)
sets %>% 
  ggvenn(show_percentage=F) +
  geom_rect(aes(xmin=-2.5,ymin=-1.5,xmax=2.5,ymax=1.5), color='darkgreen', fill=NA) + 
  annotate("text", x=0, y=-1.25, label="génome 25000", color='darkgreen')

Intuitivement, on peut comparer la fréquence de ce GOTerm dans le génome à la fréquence de ce GOTerm dans l’ensemble de gènes différentiellement exprimés :

\(f_{\text{génome}}(\text{respiration}) = \frac{250}{25000} = 0.01\) vs. \(f_{\text{diff.expr}}(\text{respiration}) = \frac{10}{50} = 0.2\)

Il apparaît donc que ce GOTerm est 20x plus fréquent dans les gènes différentiellement exprimés que dans le génome. On dit que l’ensemble de gènes différentiellement exprimés est enrichi en gènes annotés avec le GOTerm respiration.

Il s’agit maintenant de réaliser un test statistique permettant de déterminer si l’enrichissement est significatif. Les plus couramment utilisés sont basés sur la loi du \(\chi^2\), la loi binomiale et la loi hypergéométrique.

Le test du \(\chi^2\) d’indépendance permet de déterminer s’il existe un lien entre 2 facteurs (variables qualitatives) ou s’ils sont indépendants.

Il s’agit de construire une table de contingence comme suit :

cont = matrix(
  c(10,     40,
    240,  24710), 
  ncol=2, byrow = T, dimnames=list(c("diff. exprimés","pas diff.exprimés"), c("respiration","pas respiration")))
cont %>% kabex
respiration pas respiration
diff. exprimés 10 40
pas diff.exprimés 240 24710

Si les 2 facteurs étaient indépendants, la probabilité d’être les 2 serait :

\[ p(\text{diff. exprimés} \cap \text{respiration}) = \frac{50}{25000} \times \frac{250}{25000} = 0.00002\]

et donc sur une population de 25000 individus, on en aurait \(25000 \times 0.00002 = 0.5\) plutôt que les 10 observés dans la case en haut à gauche. Il ne semble donc pas y avoir indépendance entre ces 2 facteurs.

Les tests basés sur la loi binomiale et la loi hypergéométrique consistent à considérer les 2 ensembles (l’un différentiellement exprimés, et l’autre annotés respiration) comme 2 échantillons tirés aléatoirement au sein d’une même population (génome). Ensuite, il s’agit de calculer la probabilité d’obtenir une intersection (les 10 gènes à la fois différentiellement exprimés et annotés respiration) aussi grande par hasard.

Pour la loi binomiale, il s’agit d’un tirage aléatoire avec remise (on peut tirer plusieurs fois le même gène) alors que pour la loi hypergéométrique il s’agit d’un tirage aléatoire sans remise.

Ainsi, si on tire 50 fois aléatoirement 1 gène parmi les 25000 possibles, et à chaque tirage (indépendant) la probabilité de 250/25000 = 0.01 d’en prendre un annoté respiration, le test binomial nous donne la probabilité d’en obtenir 10 par hasard qui soit annotés respiration :

dbinom(10, 50, prob=0.01) # 10 succès sur 50 tentatives avec une probabilité de succès de 0.01
## [1] 6.871864e-11

La probabilité d’en avoir exactement 10 parmi 50 est de \(6,9.10^{-11}\) donc presque nulle (\(P(X=x)\)). En pratique, il faut tester la probabilité d’en avoir 10 ou davantage sinon on pourrait avoir \(P(X=0)\) très faible (et donc significative alors que 0 gènes en commun n’est pas intéressant) : \(P(X \ge x)\) soit \(7,1.10^{-11}\).

pbinom(9, 50, prob=0.01, lower.tail = F) # strictement plus que 9 succès, sur 50 tentatives avec une proba de succès de 0.01
## [1] 7.132812e-11
x=0:50
plot(x, pbinom(x-1, 50, prob=0.01, lower.tail = F), type='l', ylim=c(0,1), ylab='p(X ≥ x)')

Pour les fonctions biologiques qui nous intéressent particulièrement, les résultats suivant ont été obtenus (p-valeur ajustée avec la FDR) :

rich = tibble(GOTerm=    c("Leaf growth", "Response to bacteria", "Stress response", "Response to PAMP", "Amino acid metabolism"),
       `p-valeur`=c( 0.8,        0.004,                0.001,           0.00005,          0.7 )) 
rich %>% t %>% kabex
GOTerm Leaf growth Response to bacteria Stress response Response to PAMP Amino acid metabolism
p-valeur 8e-01 4e-03 1e-03 5e-05 7e-01

Qu’en concluez-vous ?

rich %>%
  mutate(significative = `p-valeur`<0.05) %>%
  kabex
GOTerm p-valeur significative
Leaf growth 8e-01 FALSE
Response to bacteria 4e-03 TRUE
Stress response 1e-03 TRUE
Response to PAMP 5e-05 TRUE
Amino acid metabolism 7e-01 FALSE