tests de fisher mutiples et ajout des p_value sur graph

Version imprimable

Bonjour,

J'ai une data.frame qui est structurée comme ci dessous et j'aimerai pouvoir effectuer des tests de fisher entre chaque lignes et les placer sur mon barplot.
J"aimerai aussi pouvoir obtenir un tableau qui indique les p_value et les données comparées (A gene A vs B gene A) pour chaque test
Code:

1 2 3 4 5 6 7 8 9 10 11 12 data=data.frame(statut=c("A", "B", "C", "D", "A", "B", "C", "D"), nbr_in_list=c(5,15,7,8,4,3,5,8), nbr_not_in_list=c(9,7,5,4,8,6,4,5), in_list=c("geneA","geneA","geneA","geneA", "geneB", "geneB", "geneB", "geneB")) statut nbr_in_list nbr_not_in_list in_list 1 A 5 9 geneA 2 B 15 7 geneA 3 C 7 5 geneA 4 D 8 4 geneA 5 A 4 8 geneB 6 B 3 6 geneB 7 C 5 4 geneB 8 D 8 5 geneB
Pour faire un test de fisher manuellement, pas de pb mais je n'arrive pas à le faire sur chaque paire de ligne automatiquement, et je ne vois pas non plus comment obtenir un tableau résumé avec les p value. J'ai déjà essayé un apply mais je me suis vite retrouvée bloquée.

J'aimerai donc avoir un peu d'aide la dessus.

Ensuite quel seraient les outils à utiliser pour pouvoir ajouter une barre de comparaison (avec p value) uniquement pour les p value significatives qui seraient indiquée au dessus de chaque groupes significativement différents ?
Code:

1 2 3 4 5 6 library(ggplot2) bp_tab_pct<- ggplot(data, aes(x=in_list, y=nbr_in_list, fill=statut))+ geom_bar(stat = "identity",position="dodge")
Merci d'avance pour vos réponses,

Aline

04/06/2021, 13h02
olivier.decourt

Bonjour Aline.
Le package {ggpubr} propose des fonctionnalités en ce sens, y compris pour calculer les p-values. Voir la vignette https://rpkgs.datanovia.com/ggpubr/.
Les fonctions de ce package sont compatibles avec {ggplot2}, donc on peut imaginer enchaîner après votre geom_bar() des fonctions comme stat_compare_means avec un +, comme d'habitude.

2 pièce(s) jointe(s)

Bonjour,

Je vous remercie pour votre réponse, mais je n'ai rien trouvé concernant les test de fischer.
J'ai trouvé un autre moyen pour calculer les p_value en décomposant par gene mais maintenant je suis toujours bloquée pour rajouter les pvalues/odd ratio à mon barplot

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
 
data=data.frame(statut=c("A", "B", "C", "D", "A", "B", "C", "D"), nbr_in_list=c(5,15,7,8,4,3,5,8), nbr_not_in_list=c(9,7,5,4,8,6,4,5), in_list=c("geneA","geneA","geneA","geneA", "geneB", "geneB", "geneB", "geneB"))
data1=data %>% filter(in_list=="geneA")
 
 
idx = t(combn(seq_along(data1$statut),2))
 
res_A = lapply(1:nrow(idx),function(i){
  test = fisher.test(data1[idx[i,],c("nbr_in_list","nbr_not_in_list")])
  data.frame(
    group1 = data1$statut[idx[i,1]],
    group2 = data1$statut[idx[i,2]],
    odds_ratio = as.numeric(test$estimate),
    p = as.numeric(test$p.value)
  )
})
 
res_A = do.call(rbind,res_A)
res_A
 
 
data2=data %>% filter(in_list=="geneB")
 
idx = t(combn(seq_along(data1$statut),2))
 
res_B = lapply(1:nrow(idx),function(i){
  test = fisher.test(data1[idx[i,],c("nbr_in_list","nbr_not_in_list")])
  data.frame(
    group1 = data1$statut[idx[i,1]],
    group2 = data1$statut[idx[i,2]],
    odds_ratio = as.numeric(test$estimate),
    p = as.numeric(test$p.value)
  )
})
 
res_B = do.call(rbind,res_B)
 
res_A = cbind(res_A, in_list=c("geneA", "geneA", "geneA", "geneA", "geneA", "geneA"))
res_B = cbind(res_B, in_list=c("geneB", "geneB", "geneB", "geneB", "geneB", "geneB"))
 
res=rbind(res_A, res_B)
 
   group1 group2 odds_ratio          p in_list
1       A      B  0.2700299 0.08744158   geneA
2       A      C  0.4117416 0.43123974   geneA
3       A      D  0.2927966 0.23774398   geneA
4       B      C  1.5110621 0.71134315   geneA
5       B      D  1.0692477 1.00000000   geneA
6       C      D  0.7105367 1.00000000   geneA
7       A      B  0.2700299 0.08744158   geneB
8       A      C  0.4117416 0.43123974   geneB
9       A      D  0.2927966 0.23774398   geneB
10      B      C  1.5110621 0.71134315   geneB
11      B      D  1.0692477 1.00000000   geneB
12      C      D  0.7105367 1.00000000   geneB

code de mon barplot + image

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
data = data %>% mutate(tot=nbr_in_list+nbr_not_in_list)
 
res2=melt(as.data.table(data), id.vars = c("statut", "in_list", "tot")) 
res3= res2 %>% mutate(pct=value/tot*100) 
 
 
bp_tab_pct<- ggplot(res3, aes(x=statut, y=pct, fill=variable))+
  geom_bar(stat = "identity",position="dodge") + 
  theme(plot.title = element_text(color="black", size=15, face="bold.italic")) +
  labs(fill = "") +
  ylim(c(0,100))+
  ggtitle(label = "")+
  geom_text(aes(y = pct+2, label = paste(round(pct,2),"%", "\n","(",value,"/",tot,")")), color = "black",size=5, position = position_dodge(width = 0.9))+
  theme(plot.title = element_text(color="black", size=15, face="bold.italic"))+
  theme(axis.text.x = element_text(size=16), axis.text.y = element_text(size=14))+
  theme(axis.title.y = element_text(size=16), axis.title.x = element_blank(),) +
  ylab(paste("percentage of carriers (number)")) +
  labs(title = "% of carrier of exonic variations in HC_NDD genes")
 
bp_tab_pct + facet_grid(in_list ~ .)

Pièce jointe 599993