HAZAN_Jeremy
Jeremy Hazan
03 mai 2020
library(tidyverse)
library(ggmosaic)
library(genetics)
library(LDheatmap)
library(qqman)1ere partie : Calculs sur pedigree
Exercice 1
On se pose le problème d’identifier le modèle sous-jacent à une maladie donnée par analyse de ségrégation. Pour cela on considère le pedigree de la figure en bas dans lequel les individus en gris sont atteints par la maladie. On note Gi le génotype de l’individu i et Pi son phénotype (1 pour malade et 0 pour sain).
On fait l’hypothèse que la maladie est une maladie génétique due à mutation dans un seul locus d’allèles S,s. On suppose que la fréquence de l’allèle de susceptibilité S est q et que il y a équilibre de Hardy-Weinberg.
On commence par considérer un modèle de maladie dominant à pénétrance complète.
Question 1
Ecrire les fonctions de pénétrances \(ℙ(P=1|G=g)\) pour tous les trois génotypes possibles.
Nous avons 3 types de génotypes possibles : SS,Ss,ss, on a donc trois fonctions de pénétrances qui sont les suivantes :
fSS=ℙ(P=1|G1=SS)=1 car l’individu sera malade vu qu’il possède l’allèle de susceptibilité S. fSs=ℙ(P=1|G1=Ss)=1 fss=ℙ(P=1|G1=ss)=0 car l’individu ne sera pas malade
Question 2
Trouver les génotypes possibles pour chaque individu.
On rappelle que S est l’allèle de susceptibilité, c’est donc l’allèle associé à la maladie. De plus, nous somme dans le cas de maladie dominant à pénétrance complète. D’après le pedigree, nous avons une famille avec les individus 1 et 2 qui sont les parents et l’individu 3 qui est leur enfant. Les individus 2 et 3 sont malades tandis que l’individu 1 ne l’est pas. Grâce à ces informations :
Individu 1 : Il n’est pas atteint par la maladie, il ne peut donc pas posséder l’allèle S. Il ne peut donc avoir que le génotype ss.
Individu 2 : Il est atteint par la maladie, il possède forcément au moins un allèle S. Il peut donc avoir deux génotype Ss ou SS.
Individu 3 : Il est atteint par la maladie, il possède forcément au moins un allèle S. Or, il est le fruit de la fécondation de l’individu 1 et de l’individu 2. En effet, il récupère de l’individu 1 un allèle s et comme il est malade il récupère forcément un allèle S de l’individu 2. Il ne peut donc avoir qu’un seul génotype qui est Ss.
En conclusion, nous avons : g1∈{ss}, g2∈{SS,Ss} et g3∈{Ss}.
Question 3
Pour tout génotype possible g1 pour l’individu 1, on calcule ℙ(G=g). Même question pour l’individu 2
On rappelle que la fréquence de S=q, donc la fréquence de s=1−q. D’après l’équilibre de HW, on a :
Pour l’individu 1 : ℙ(G1=g1)=ℙ(G1=ss)=(1−q)2
Pour l’individu 2 : ℙ(G2=g2)=ℙ(G2∈Ss,SS)=ℙ(G2=Ss)+ℙ(G2=SS)=2q(1−q)+q2
Question 4
Pour chaque combinaison de génotypes possibles (g1,g2,g3), on calcule ℙ(G3=g3|G1=g1,G2=g2).
Les génotypes des individus 1 et 3 sont unique, réciproquement : ss et Ss. Il faut donc calculer ℙ(G3=Ss|G1=ss,G2=SS(ou)Ss).
Génotype G2=Ss : ℙ(G3=Ss|G1=ss,G2=Ss) = 12. En effet, la probabilité que l’individu 3 reçoive un allèle s de l’individu 1 vaut 1 car il ne peut donner qu’un allèle s. De plus, la probabilité que l’individu 3 reçoive un allèle S de l’individu 2 vaut 12. En multipliant ces probabilité, on obtient le résultat ci-dessus.
Génotype G2=SS : ℙ(G3=Ss|G1=ss,G2=SS)=1. En effet, de la même manière que précédemment, l’individu 3 reçoit forcément un allèle s de l’individu 1 et cette probabilité vaut 1 . De plus, l’individu 3 reçoit forcément un allèle S de l’individu 2 et cette probabilité vaut 1. En multipliant ces probabilité, on obtient le résultat ci-dessus.
Question 5
Ecrire la vraisemblance ℙ(P1=0,P2=1,P3=1) du trio sous ce modèle de maladie, à l’aide des questions précédentes.
Notons V la vraisemblance du modèle.
V=ℙ(P1=0,P2=1,P3=1)=ℙ(P1=0)ℙ(P2=1|P1=0)ℙ(P3=1|P1=0,P2=1)
Or, d’après les questions précédentes, les individus 1 et 2 sont indépendants vu que ce sont les parents de l’individu 3.
D’où V=ℙ(P1=0)ℙ(P2=1)ℙ(P3=1|P1=0,P2=1)
Or : ℙ(P1=0)=ℙ(P1=0|G1=ss)ℙ(G1=ss)=(1−q)2
Cas où G2=SS : ℙ(P2=1)=ℙ(P2=1|G2=SS)ℙ(G2=SS)=ℙ(G2=SS)=q2 ℙ(P3=1|P1=0,P2=1)=ℙ(P3=1|G3=Ss)ℙ(G3=Ss|G1=ss,G2=SS)=1
Cas où G2=Ss : ℙ(P2=1)=ℙ(P2=1|G2=Ss)ℙ(G2=Ss)=ℙ(G2=Ss)=2q(1−q) ℙ(P3=1|P1=0,P2=1)=ℙ(P3=1|G3=Ss)ℙ(G3=Ss|G1=ss,G2=Ss)=12
Donc on obtient : V=(1−q)2[1×q2+12×2q(1−q)]=(1−q)2(q2+q(1−q))
Question 6
Calculer la vraisemblance du trio, si q= \(\frac{1}{2}\).
q = 1/2
Vrais = (1-q)^2 * (q^2 + q*(1-q))
Vrais #1/8[1] 0.125Dans un deuxième temps, on considère un modèle de maladie récessif à pénétrance complète.
Question 7
Quelle est la vraisemblance du trio ?
Afin de répondre correctement à cette question, nous devons faire les mêmes étapes que précédemment. On se situe maintenant dans un modèle de maladie récessif à pénétrance complète : il faut obligatoirement la présence de deux allèles S pour que la maladie se manifeste. Nous avons donc les fonctions de pénétrances suivantes: fSS=ℙ(P=1|G1=SS)=1 fSs=ℙ(P=1|G1=Ss)=0 fss=ℙ(P=1|G1=ss)=0
Les génotypes ne sont plus les mêmes que précédemment. En effet, l’individu 1 peut posséder un allèle S car S est considérer comme récessif. De plus, pour que les individus soient malades, il faut qu’ils possèdent deux allèles S. Nous avons donc : g1∈{Ss,ss}, g2∈{SS} et g3∈{SS}.
Passons maintenant au calcul de la vraisemblance : Notons V2 la vraisemblance du modèle récessif et que les individus 1 et 2 sont indépendants. De plus, rappelons que les probabilités d’avoir le génotype SS,Ss,ss sont données par l’équilibre de HW. V2=ℙ(P1=0,P2=1,P3=1)=ℙ(P1=0)ℙ(P2=1)ℙ(P3=1|P1=0,P2=1)
Or, ℙ(P2=1)=ℙ(P2=1|G2=SS)ℙ(G2=SS)=ℙ(G2=SS)=q2
Cas où G1=ss : ℙ(P1=0)=ℙ(P1=0|G1=ss)ℙ(G1=ss)=ℙ(G1=ss)=(1−q)2 ℙ(P3=1|P1=0,P2=1)=ℙ(P3=1|G3=SS)ℙ(G3=SS|G1=ss,G2=SS)=ℙ(G3=SS|G1=ss,G2=SS)=0 (L’individu 3 ne va jamais recevoir un allèle S de l’individu 1, vu qu’il a pour génotype G1=ss).
Cas où G1=Ss : ℙ(P1=0) = ℙ(P1=0|G1=Ss)ℙ(G1=Ss) = ℙ(G1=Ss)=2q(1−q) ℙ(P3=1|P1=0,P2=1)=ℙ(P3=1|G3=SS)ℙ(G3=SS|G1=Ss,G2=SS)=ℙ(G3=SS|G1=Ss,G2=SS)=12
Donc on obtient : V2=q2(0×(1−q)2)(12×2q(1−q))=q2(q(1−q))=q3(1−q)
Pour q= \(\frac{1}{2}\), on obtient V2:
q = 1/2
V2 = q^3*(1-q)
V2 #1/16[1] 0.0625Enfin, on suppose que la maladie n’a pas de composante génétique. On note F la fréquence de la maladie dans la population.
Question 8
Ecrire la vraisemblance du trio observé sous cette nouvelle hypothèse.
Afin de répondre correctement à cette question, il faut expliquer notre modèle. Nous sommes dans le cas où la fréquence de tomber malade est de F. On peut donc modéliser ce modèle par une loi de Bernoulli(F) : La probabilité d’être malade : ℙ(P=1)=F et la probabilité de ne pas être malade ℙ(P=0)=1−F. Ici, nous ne parlons pas d’allèles, on peut donc dire que nos individus sont indépendants.
On calcule alors la vraisemblance du modèle : V3=ℙ(P1=0,P2=1,P3=1)=ℙ(P1=0)ℙ(P2=1)ℙ(P3=1)=(1−F)×F×F=F2(1−F)
Question 9
Calculer la vraisemblance du trio, si F= \(\frac{1}{20}\)
F = 1/20
V3 = F^2*(1-F)
V3 #19/8000[1] 0.0023752ème partie : Etude Human Genetic Diversity Project (HGDP)
Le projet HGDP de l’Université de Stanford a pour but la caractérisation de la diversité génétique dans les populations humaines. Pour cela, les génotypes de 650 000 marqueurs (SNPs) ont été identifiés pour plus de 1000 individus dans de nombreux pays. Nous allons travailler avec les données relative aux SNPs se trouvant dans le gène AKT1.
Exercice 2 : Explorations préliminaires - covariables
Question 1 : Chargement des données
data_hgdp = read.table("/Users/jeremyhazan/Desktop/HGDP_AKT1.txt" , header = TRUE, sep = '\t', fill = TRUE)
head(data_hgdp) Well ID Gender Population Geographic.origin Geographic.area
1 B12 HGDP00980 F Biaka Pygmies Central African Republic Central Africa
2 A12 HGDP01406 M Bantu Kenya Central Africa
3 E5 HGDP01266 M Mozabite Algeria (Mzab) Northern Africa
4 B9 HGDP01006 F Karitiana Brazil South America
5 E1 HGDP01220 M Daur China China
6 H2 HGDP01288 M Han China China
AKT1.C0756A AKT1.C6024T AKT1.G2347T AKT1.G2375A
1 CA CT TT AA
2 CA CT TT AA
3 AA TT TT AA
4 AA TT TT AA
5 AA TT TT AA
6 AA TT TT AAQuestion 2
Trouver le nombre d’observations (individus) et le nombre de SNPs disponibles (les SNPs dans le gène AKT1 sont dénotés avec le préfixe AKT1).
print(dim(data_hgdp))[1] 1064 10names(data_hgdp) [1] "Well" "ID" "Gender"
[4] "Population" "Geographic.origin" "Geographic.area"
[7] "AKT1.C0756A" "AKT1.C6024T" "AKT1.G2347T"
[10] "AKT1.G2375A" Notre jeu de données est composée de 1064 observations et 10 variables. De plus, on remarque que nos 4 dernières variables possèdent le préfixe AKT1 : ce sont nos 4 SNPs.
Question 3
Combien de femmes et d’hommes y a-t-il dans l’étude? Combien de populations? Pour combien de zones géographiques regroupant plusieurs pays a-t-on des données?
summary(data_hgdp$Gender) F M
380 684 summary(data_hgdp$Population) Adygei Balochi Bantu Bedouin Biaka Pygmies
17 25 20 49 11
Biaka Pygmies Brahui Burusho Cambodian Colombian
25 25 25 11 13
Dai Daur Druze French French Basque
10 10 48 29 24
Han Hazara Hezhen Japanese Kalash
45 25 10 31 25
Karitiana Lahu Makrani Mandenka Maya
24 10 25 24 25
Mbuti Pygmies Miaozu Mongola Mozabite NAN Melanesian
15 10 10 30 22
Naxi North Italian Orcadian Oroqen Palestinian
10 14 16 10 51
Papuan Pathan Pima Russian San
17 25 25 25 7
Sardinian She Sindhi Surui Tu
28 10 25 21 10
Tujia Tuscan Uygur Xibo Yakut
10 8 10 9 25
Yizu Yoruba
10 25 length(summary(data_hgdp$Population))[1] 52summary(data_hgdp$Geographic.area) Central Africa Central America China Israel Japan
119 50 184 148 31
New Guinea Northern Africa Northern Europe Pakistan Russia
17 30 16 200 67
South Africa South America Southeast Asia Southern Europe
8 58 11 125 length(summary(data_hgdp$Geographic.area))[1] 14On peut voir que dans notre étude qu’il y a 380 femmes et 684 hommes. De plus, il y a 52 populations différentes et 14 zones géographiques.
Question 4
Quelles sont les populations et les zones géographiques le plus représentées? Et le moins? Représenter graphiquement ces deux variables à l’aide d’un diagramme en bâtons.
ggplot(data_hgdp) +
geom_bar(mapping = aes(x = Geographic.area, fill = Population)) +
ggtitle('Diagramme en bâtons des zones géographiques avec les populations')+
xlab('Zone géographique')+
coord_flip()
ggplot(data_hgdp)+
geom_bar(mapping = aes(x = Population, fill = Geographic.area)) +
ggtitle('Diagramme en bâtons des populations avec les zones géographiques')+
xlab('Population')+
coord_flip()
On peut voir à l’aide de la question 3 et des diagrammes en bâtons ci-dessus que les zones géographiques les plus représentés sont le Pakistan, la Chine et Israël. En outre, les populations les plus représentées sont les Palestiniens, les Bédouins ainsi que les Druzes.
De même, on peut voir que les zones géographiques les moins représentées sont le sud de l’Afrique, le sud de l’Asie et le nord de l’Europe; quant aux populations, celles qui sont les moins représentées sont les Sans, les Toscans et les Xibos.
Exercice 3 : Exploration préliminaires - génotypes
Question 1
Estimer les fréquences alléliques et génotypiques du SNP AKT1.C6024T. Pour cela, on fait l’hypothèse que ces fréquences peuvent être estimées sans prendre en compte le fait que certaines données sont manquantes. Expliquer pourquoi il s’agit d’une hypothèse assez forte en général et la raison pour la quelle elle est justifiée dans ce cas. Vous pouvez utiliser le package genetics ou faire directement le calcul des fréquences.
Génotypes du SNP AKT1.C6024T :
summary(data_hgdp$AKT1.C6024T) CC CT TT NA's
719 294 50 1 Calcul des fréquences alléliques et génotypiques à l’aide de la fonction genotype du package genetics :
genotype_hgdp = genotype(data_hgdp$AKT1.C6024T, sep='')
summary(genotype_hgdp)
Number of samples typed: 1063 (99.9%)
Allele Frequency: (2 alleles)
Count Proportion
C 1732 0.81
T 394 0.19
NA 2 NA
Genotype Frequency:
Count Proportion
C/C 719 0.68
C/T 294 0.28
T/T 50 0.05
NA 1 NA
Heterozygosity (Hu) = 0.3021008
Poly. Inf. Content = 0.2563692Question 2
Afficher les proportions des génotypes de AKT1.C6024T pour chaque zone géographique à l’aide d’un ‘mosaic plot’.
ggplot(data_hgdp)+
geom_mosaic(mapping = aes(x = product(AKT1.C6024T), fill = Geographic.area), na.rm = TRUE)+
labs(x = 'Génotype', y = 'Zone géographique' ) +
ggtitle('Mosaic plot')
ggplot(data_hgdp)+
geom_mosaic(mapping = aes(x = product(AKT1.C6024T), fill = Geographic.area), na.rm = TRUE)+
labs(x = 'Génotype', y = 'Zone géographique' ) +
ggtitle('Mosaic plot')
ggplot(data_hgdp)+
geom_mosaic(mapping = aes(x = product(Geographic.area), fill = AKT1.C6024T), na.rm = TRUE)+
labs(x = 'Zone géographique', y = 'Génotype') +
ggtitle('Mosaic plot')
Question 3
Au vu du mosaic plot, qu’a-t-on envie de conclure concernant l’homogénéité des fréquences genotypiques de AKT1.C6024T? Que faudrait-t-il faire pour vérifier la réponse précédente d’un point de vue statistique?
Au vu du mosaic plot, on voit que le génotype CC est beaucoup plus présent que les deux autres génotypes. On pourrait donc penser que les fréquences génotypiques ne sont pas homogènes.
Pour confirmer cette hypothèse, d’un point de vue statistique, nous pouvons faire un test du \(\chi^2\) : en effet ce test statistique est possible car nous avons deux variables qualitatives.
On souhaite donc tester :
H0 : Les génotypes sont homogènes H1 : Les génotypes ne sont pas homogènes (donc ils sont hétérogènes)
chisq.test(table(data_hgdp$AKT1.C6024T,data_hgdp$Geographic.area))
Pearson's Chi-squared test
data: table(data_hgdp$AKT1.C6024T, data_hgdp$Geographic.area)
X-squared = 88.036, df = 26, p-value = 1.148e-08On obtient comme p-value 1,14.10−8 qui est significativement plus petit que \(α=0.05\) : on rejette donc H0 avec un seuil de confiance à 5% et on décide H1. En conclusion, d’après le test du \(\chi^2\), les génotypes sont hétérogènes.
Question 4
En considérant l’ensemble de la population, calculer la mesure D′ de déséquilibre de liaison (LD) pour toute paire de SNPs dans le gène AKT1. Commenter le résultat.
On commence par calculer pour chacun de nos SNPs les fréquences alléliques et génotypiques à l’aide de la fonction genotype, on construit ensuite un dataframe constitué des résultats précédents.
genotype1 = genotype(data_hgdp$AKT1.C0756A, sep='')
genotype2 = genotype(data_hgdp$AKT1.C6024T, sep='')
genotype3 = genotype(data_hgdp$AKT1.G2347T, sep='')
genotype4 = genotype(data_hgdp$AKT1.G2375A, sep='')genotype_akt = data.frame(genotype1,genotype2,genotype3,genotype4)Ensuite, à l’aide de la fonction LD, on récupère la mesure de déséquilibre de liaison D’ :
LD(genotype_akt)$"D'" genotype1 genotype2 genotype3 genotype4
genotype1 NA 0.9934369 0.9481195 0.9843420
genotype2 NA NA 0.9429031 0.9842787
genotype3 NA NA NA 0.9995238
genotype4 NA NA NA NAOn observe que nos déséquilibre de liaison sont proche de 1 ce qui indique un fort déséquilibre de liaison.
LDheatmap(genotype_akt, LDmeasure="D'")
3ème partie : Etude Functional Single Nucleotide Polymorphisms Associated with Human Muscle Size and Strength (FAMuSS)
Exercice 4 : HWE
Question 1 : Chargement des données
Exclure des données le seul individu d’origine amérindienne (cette observation crée des problèmes quand on essaie automatiser l’analyse sur les strates).
fms_data = read.table("http://w3.mi.parisdescartes.fr/~vperduca/data/FMS_data.txt", header = TRUE, sep = '\t', fill = TRUE)On souhaite enlever l’individu d’origine amérindienne :
table(fms_data$Race)
African Am Am Indian Asian Caucasian Hispanic Other
44 1 97 791 52 49 amerindien = fms_data[which(fms_data$Race =="Am Indian"),]
head(amerindien) id acdc_rs1501299 ace_id actn3_r577x actn3_rs540874 actn3_rs1815739
1107 UC-356 <NA> <NA> <NA> <NA> <NA>
actn3_1671064 ardb1_1801253 adrb2_1042713 adrb2_1042714 adrb2_rs1042718
1107 <NA> <NA> <NA> <NA> <NA>
adrb3_4994 agrp_5030980 akt1_t22932c akt1_g15129a akt1_g14803t
1107 <NA> <NA> TT AG <NA>
akt1_c10744t_c12886t akt1_t10726c_t12868c akt1_t10598a_t12740a
1107 CC CC TA
akt1_c9756a_c11898t akt1_t8407g akt1_a7699g akt1_c6148t_c8290t
1107 CT TT GG <NA>
akt1_c6024t_c8166t akt1_c5854t_c7996t akt1_c832g_c3359g akt1_g288c
1107 CT <NA> <NA> <NA>
akt1_g1780a_g363a akt1_g2347t_g205t akt1_g2375a_g233a akt1_g4362c
1107 <NA> GT GG GG
akt1_c15676t akt1_a15756t akt1_g20703a akt1_g22187a akt1_a22889g
1107 <NA> <NA> <NA> GA <NA>
akt1_g23477a akt2_7254617 akt2_rs892118 akt2_2304186 akt2_969531
1107 <NA> <NA> <NA> <NA> <NA>
ampd1_c34t ankrd6_a550t ankrd6_q122e ankrd6_m485l ankrd6_p636l
1107 <NA> <NA> <NA> <NA> <NA>
ankrd6_t233m ankrd6_i128l ankrd6_g197805a ankrd6_a545t ankrd6_k710x
1107 <NA> <NA> <NA> <NA> <NA>
apoe_c472t b2b bcl6_4686467 bcl6_3774298 bcl6_17797517 bcl6_1056932
1107 <NA> <NA> <NA> <NA> <NA> <NA>
bmp2_rs15705 bmp2_235768 c8orf68_rs6983944 carp_exon3_snp1 carp_exon3_snp2
1107 <NA> <NA> <NA> <NA> <NA>
carp.exon3_snp3 carp_exon5_snp1 carp_exon5_snp2 carp_exon6_ag
1107 <NA> <NA> <NA> <NA>
carp_exon8_ct carp_3utr carp_a8470g carp_c105t cast_rs754615
1107 <NA> <NA> <NA> <NA> <NA>
cast_rs7724759 cav2_q130e cntf_g6a cox5b_11904110 cox5b_17022045
1107 <NA> <NA> <NA> <NA> <NA>
ctsf_572846 ddit_rs1053227 esr1_rs1801132 esr1_rs1042717 esr1_rs2228480
1107 <NA> <NA> <NA> <NA> <NA>
esr1_rs2077647 esr1_rs9340799 esr1_rs2234693 esr2_1256112 esr2_3020450
1107 <NA> <NA> <NA> <NA> <NA>
fatso fbox32_rs6690663 fbox32_rs3739287 fbox32_rs4871385 fbox32_rs2891770
1107 <NA> <NA> <NA> <NA> <NA>
fbox32_rs6983944 fst_722910 fst_1423560 fstl1_rs9631455 galntl1_rs8005903
1107 <NA> <NA> <NA> <NA> <NA>
gapd_7307229 gapd_7971637 gdf8_t55a gdf8_1225t gdf8_p198a gdf8_e164k
1107 <NA> <NA> <NA> <NA> <NA> <NA>
gdf8_k153r gdf8_rs1805086 gdf8_rs1805085 gnb3_rs5443 gs_s287nga
1107 <NA> <NA> <NA> <NA> <NA>
hsd11b1_rs4844880 hsd11b1_rs846910 hsd11b1_rs846907 hsd11b1_rs860185
1107 <NA> <NA> <NA> <NA>
hsd11b1_rs846911 hsd11b1_rs3753519 hsd11b1_rs846909 igf1_pro igf1_2288377
1107 <NA> <NA> <NA> <NA> <NA>
igf1_t1245c igf2_rs680 igf2_2230949 igf2_3213221 igf2_3213220 igf2_3213216
1107 <NA> <NA> <NA> <NA> <NA> <NA>
igfbp3_6670 igfbp3_2132570 il6_g572c il6_c174g il15ra_3136618
1107 <NA> <NA> <NA> <NA> <NA>
il15ra_2228059 il15_1057972 il15_1589241 il15_rs2296135 insig2_rs7566605
1107 <NA> <NA> <NA> <NA> <NA>
irs1_g972r kchj11_rs5219 lep_2167270 lepr_1137100 lepr_1137101
1107 <NA> <NA> <NA> <NA> <NA>
lepr_8179183 lepr_1805096 lrp5_snp2 lrp5_snp3 lrp5_snp4 mgst3_4147542
1107 <NA> <NA> <NA> <NA> <NA> <NA>
mgst3_7549530 mlck_c49t mlck_c37885a mlck_g91689t mulk_3757470
1107 <NA> <NA> <NA> <NA> <NA>
mylk_c37885a mylk_g91689t myod1_rs2249104 nnmt_g5082t nos3_rs1799983
1107 <NA> <NA> <NA> <NA> <NA>
nr3c1_rs10482616 nr3c1_rs10482614 nr3c1_rs4634384 nr3c1_rs4582314
1107 <NA> <NA> <NA> <NA>
nr3c1_rs6195 opg_2073618 opg_3102735 p2rx2_6560891 p2rx2_7966242
1107 <NA> <NA> <NA> <NA> <NA>
p2ry2_2511241 p2ry2_rs1783596 pai1_4g5g pcr13_snp1 pcr15_snp1 pcr15_snp2
1107 <NA> <NA> <NA> <NA> <NA> <NA>
pcr15_snp3 pcr15_snp4 pcr5_snp1 pcr5_snp2 pcr8_snp1 pgc1_g76039a
1107 <NA> <NA> <NA> <NA> <NA> <NA>
pgm1_2284998 pik3cg_4727666 pik3_rs3173908 ppara_1800206 ppara_135539
1107 <NA> <NA> <NA> <NA> <NA>
ppara_4252778 ppar_gp12a pparg_c1a_rs8192678 prdx6_4382766 pygm_620006
1107 <NA> <NA> <NA> <NA> <NA>
rankl_17458177 rankl_17536280 rankl_17639305 rankl_4531631 resistin_c30t
1107 <NA> <NA> <NA> <NA> <NA>
resistin_c398t resistin_g540a resistin_c980g resistin_c180g resistin_a537c
1107 <NA> <NA> <NA> <NA> <NA>
rnf28_c16242t rs302964 rs11630261 rs849409 runx_12526462
1107 <NA> <NA> <NA> <NA> <NA>
slc35f1_rs12110370 slc35f1_rs10484290 syngr2_c886t tcfl72_12255372
1107 <NA> <NA> <NA> <NA>
tcfl72_7903146 tcfl72_rs12255372 tcfl72_rs7903146 tnfa_g308a tnni1_2799672
1107 <NA> <NA> <NA> <NA> <NA>
tnni1_3738287 tpd52l1_3778458 tpd52l1_514096 tpd52l1_4896782
1107 <NA> <NA> <NA> <NA>
tdp52l1_9321028 tpd52l1_3799736 ucp2_ct vdr_taq1 vdr_fok1 vdr_bsm1
1107 <NA> <NA> <NA> <NA> <NA> <NA>
vdr_rs731236 vim_6602186 visfatin_2041681 visfatin_6947766
1107 <NA> <NA> <NA> <NA>
visfatin_3801272 visfatin_10487820 visfatin_929604 visfatin_10953502
1107 <NA> <NA> <NA> <NA>
Status Center Term Gender Age Race Pre_NDRM_Max Post_NDRM_Max
1107 0 UC 02-3 Male 38 Am Indian 22.5 25
NDRM_DIFF NDRM.CH Pre_DRM_Max Post_DRM_Max DRM_DIFF DRM.CH Pre.weight
1107 2.5 11.1 22.5 27.5 5 22.2 191
Pre.height pre.BMI SBP DBP Pre.HR Pre_RBi_Avg Pre_RTri_Avg Pre_LBi_Avg
1107 67.5 29.47 160 100 86 9.33 16.33 10
Pre_LTri_Avg Dom.Arm Post.weight Post.Height Calc.post.BMI Post_RBi_Avg
1107 16.33 R 200 67.5 30.859 29.67
Post_RTri_Avg Post_LBi_Avg Post_LTri_Avg V1_ND1 V1_ND2 V1_ND3 V1_ND_AVG
1107 42 29.67 45 172.7 178.3 189 180
V2_ND1 V2_ND2 V2_ND3 V2_ND_AVG V3_ND1 V3_ND2 V3_ND3 V3_ND_AVG V23_ND_AVG
1107 168.7 201.1 197 188.93 129.5 156.8 134.4 140.23 164.58
V123_ND_AVG Post1_ND1 Post1_ND2 Post1_ND3 Post1_ND_AVG Post2_ND1 Post2_ND2
1107 169.72 145.9 134.1 146.6 142.2 NA NA
Post2_ND3 Post2_ND_AVG Post_ND_avg ND23_DIFF ND23.CH V1_D1 V1_D2 V1_D3
1107 NA NA 142.2 -22.38 -13.6 196.6 183.8 210.2
V1_D_AVG V2_D1 V2_D2 V2_D3 V2_D_AVG V3_D1 V3_D2 V3_D3 V3_D_AVG V23_D_AVG
1107 196.87 195.1 190.1 207 197.4 185.8 186.3 155.8 175.97 186.68
V123_D_AVG Post1_D1 Post1_D2 Post1_D3 Post1_D_AVG Post2_D1 Post2_D2
1107 190.08 126.3 133.7 131.7 130.57 NA NA
Post2_D3 Post2_D_AVG Post_D_avg D23_DIFF D23.CH FLGU TG CHOL HDL_C
1107 NA NA 130.57 -56.12 -30.06 94 76 226 76
CHOL_HDL_C VLDL_TG LDL_C FINS CRP Mean_BP HOMA Met_syn CK_day_0 CK_day_2
1107 2.97 15 135 5 0.4 120 1.15 0 NA NA
CK_day_4 ND_Pre_WA_VOL ND_Pre_WM_VOL ND_Pre_F_VOL ND_Pre_BM_VOL
1107 NA 688330.6 665629.5 292496.2 22701.16
ND_Pre_M_VOL ND_Pre_B_VOL ND_Post_WA_VOL ND_Post_WM_VOL ND_Post_F_VOL
1107 5490.972 17210.19 867765.7 844335.9 436569.2
ND_Post_BM_VOL ND_Post_M_VOL ND_Post_B_VOL D_Pre_WA_VOL D_Pre_WM_VOL
1107 23429.77 5073.678 18356.09 713159.2 690650.5
D_Pre_F_VOL D_Pre_BM_VOL D_Pre_M_VOL D_Pre_B_VOL D_Post_WA_VOL
1107 289098.7 22508.68 5909.655 16599.02 837716.4
D_Post_WM_VOL D_Post_F_VOL D_Post_BM_VOL D_Post_M_VOL D_Post_B_VOL
1107 814765.6 424385.5 22950.82 3446.034 19504.79On remarque que les informations concernant l’individu d’origine amérindienne sont à la ligne 1107. On décide donc d’enlever la ligne 1107 de notre dataframe :
fms_data = fms_data[-1107,]table(fms_data$Race)
African Am Am Indian Asian Caucasian Hispanic Other
44 0 97 791 52 49 Question 2
Tester si le SNP akt1_t10726c_t12868c est en HWE dans l’ensemble de la population.
hwe_geno = genotype(fms_data$akt1_t10726c_t12868c, sep = '')
summary(hwe_geno)
Number of samples typed: 1244 (89.1%)
Allele Frequency: (2 alleles)
Count Proportion
C 2063 0.83
T 425 0.17
NA 304 NA
Genotype Frequency:
Count Proportion
C/C 881 0.71
C/T 301 0.24
T/T 62 0.05
NA 152 NA
Heterozygosity (Hu) = 0.2833949
Poly. Inf. Content = 0.2431569On teste donc à l’aide du test du \(\chi^2\) :
H0 : Il y a équilibre de HW H1 : Il n’y a pas équilibre de HW
HWE.chisq(hwe_geno)
Pearson's Chi-squared test with simulated p-value (based on 10000
replicates)
data: tab
X-squared = 26.467, df = NA, p-value = 9.999e-05On obtient comme p-value 9,99.10−5 qui est significativement plus petit que \(α=0.05\) : on rejette donc H0 avec un seuil de confiance à 5% et on décide H1. En conclusion, d’après le test du \(\chi^2\), il n’y a pas équilibre de HW.
Question 3
Tester si le SNP akt1_t10726c_t12868c est en HWE dans chaque strate de la variable Race. Comment expliquez-vous le résultat du point 2 au vu de ces résultats?
ggplot(fms_data)+
geom_mosaic(mapping = aes(x = product(akt1_t10726c_t12868c), fill = Race), na.rm = TRUE)+
labs(x = 'Génotype', y = 'Origine éthnique' )
On utilise la fonction genotype pour afficher les fréquences alléliques et génotypiques pour notre SNP en fonction de chacune de nos origines ethniques :
geno_african = genotype(fms_data$akt1_t10726c_t12868c[fms_data$Race == 'African Am'], sep = '')
geno_asian = genotype(fms_data$akt1_t10726c_t12868c[fms_data$Race == 'Asian'], sep = '')
geno_caucasian = genotype(fms_data$akt1_t10726c_t12868c[fms_data$Race == 'Caucasian'], sep = '')
geno_hispanic = genotype(fms_data$akt1_t10726c_t12868c[fms_data$Race == 'Hispanic'], sep = '')
geno_other = genotype(fms_data$akt1_t10726c_t12868c[fms_data$Race == 'Other'], sep = '')On teste maintenant à l’aide du test du χ2 pour chaque origine ethnique si il y a équilibre de HW ou non : C’est à dire qu’on teste :
H0 : Il y a équilibre de HW H1 : Il n’y a pas équilibre de HW
HWE.chisq(geno_african)
Pearson's Chi-squared test with simulated p-value (based on 10000
replicates)
data: tab
X-squared = 2.0144, df = NA, p-value = 0.2329HWE.chisq(geno_asian)
Pearson's Chi-squared test with simulated p-value (based on 10000
replicates)
data: tab
X-squared = 2.3611, df = NA, p-value = 0.1668HWE.chisq(geno_caucasian)
Pearson's Chi-squared test with simulated p-value (based on 10000
replicates)
data: tab
X-squared = 0.0030898, df = NA, p-value = 1HWE.chisq(geno_hispanic)
Pearson's Chi-squared test with simulated p-value (based on 10000
replicates)
data: tab
X-squared = 0.70512, df = NA, p-value = 0.5784HWE.chisq(geno_other)
Pearson's Chi-squared test with simulated p-value (based on 10000
replicates)
data: tab
X-squared = 0.60706, df = NA, p-value = 0.6528On remarque que pour chaque test du \(\chi^2\) que nous avons effectué, toutes les p-values (0.224,0.1673,1,0.5779,0.6599) sont strictement supérieure à α=0.05 : on décide donc H0 avec seuil de confiance de 5%. En conclusion, pour chaque origine ethnique il y a équilibre de HW. On a donc pas équilibre de HW pour la population globale de notre jeu de données, mais une fois qu’elle est stratifiée (via les origines ethniques) il y a équilibre de HW sur l’ensemble de la population.
Exercice 5 : Association, tests multiples
On s’intéresse à la variable NDRM.CH, le changement en pourcentage de la force du bras non dominant avant et après le programme d’entraînement physique prévu dans l’étude. On se demande si NDRM.CH est associée à un ou plusieurs SNPs.
Question 1
A partir de NDRM.CH construire la variable aléatoire binaire Y qui vaut 1 si NDRM.CH>60 et 0 autrement.
Y = fms_data$NDRM.CH
Y[fms_data$NDRM.CH > 60] = 1
Y[fms_data$NDRM.CH <= 60] = 0Question 2
Tester l’association entre tous les SNPs et Y. On effectue donc un test du \(\chi^2\). On teste :
H0 : Le snp est indépendant de NDRM.CH H1 : Le snp est lié à NDRM.CH
#On commence par enlever les colonnes qui ne sont pas des snps
data_snp = fms_data[, -c(1,3, 53, 165, 227:347)]
snps = names(data_snp)
pval = apply(data_snp, 2, function(x)chisq.test(table(x,Y))$p.value)
#On décide de faire le test à un seuil de 5% :
snp_sign = pval[pval< 0.05]
snp_sign ardb1_1801253 akt1_g4362c <NA> carp_exon3_snp1
2.982194e-02 2.480870e-02 NA 3.556811e-04
carp_exon3_snp2 carp.exon3_snp3 carp_exon5_snp1 cox5b_17022045
4.652582e-04 4.652582e-04 6.051723e-04 4.786887e-02
fbox32_rs4871385 gdf8_1225t gdf8_p198a gdf8_e164k
1.763484e-02 5.082999e-04 5.082999e-04 5.082999e-04
hsd11b1_rs4844880 <NA> mylk_g91689t myod1_rs2249104
4.231747e-02 NA 2.233967e-02 4.932344e-02
nr3c1_rs4634384 pygm_620006 resistin_c30t resistin_c398t
3.822049e-02 4.799825e-03 4.703417e-02 1.122223e-02
resistin_g540a resistin_c180g rs849409 vdr_taq1
3.059336e-03 2.094796e-02 9.735323e-14 1.693855e-02
visfatin_3801272 visfatin_10487820
1.936865e-02 3.351943e-02 On remarque que pour chaque test du χ2 que nous avons effectué, qu’il y a 24 SNPs dont leurs pvaleurs sont strictement inférieure à α=0.05 : on rejette H0 et on décide donc H1 avec seuil de confiance de 5%. En conclusion, parmi nos 222 SNPs de notre jeu de données, seulement 24 sont associés au changement en pourcentage de la force du bras non dominant avant et après le programme d’entraînement physique prévu dans l’étude.
Question 3, 4 et 5
datafms= data.frame(SNP = snps, CHR = 1:length(snps), BP = 1:length(snps), P = pval)
datafms2 = datafms[-c(40,130),] #il faut supprimer les valeurs manquantes
manhattan(datafms2, highlight = names(snp_sign), annotatePval = 0.0001)
Sur le graphique ci-dessous, nous avons utilisé la library tidyverse et la fonction ggplot pour afficher un Manhattan Plot. On affichera en vert les 24 SNPs qui sont liés à la maladie.
datafms= data.frame(SNP = snps, P = pval)
df_sign = filter(datafms, P < (0.05/222)) #0.05/222 correspond à la correction de tests multiples
df_sign2 = filter(datafms, P < 0.05)
ggplot(datafms, aes(x = SNP, y = -log10(P))) +
geom_point()+
geom_point(data = df_sign2, col = 'green')+
geom_hline(aes(yintercept = -log10((0.05/1397)), linetype = 'Correction de Bonferroni'),col = 'red')+
geom_hline(aes(yintercept = -log10(0.05/222), linetype = 'Correction tests multiples'), col = 'blue')+
scale_linetype_manual(name = 'Légende', values = c(2, 2),
guide = guide_legend(override.aes = list(color = c('red', 'blue'))))+
geom_label(data = df_sign, aes(label = SNP), size = 3,vjust = 1.5, hjust = 0.5)+
ggtitle('Manhattan Plot')+
theme_minimal()
Après correction pour tests multiples (ligne pointillée bleue), on remarque que le seul SNP qui semble avoir une association significative avec Y est le SNP rs849409. On pourrai conclure à une association entre Y et ce SNP, mais afin d’être plus précis il faudrait effectuer une régression linéaire et un calcul de l’OR pour confirmer que le SNP influe sur Y.
4ème partie : autres données
Exercice 6 : Modélisation, intéraction
Question 1
Charger les données. Comme dans l’exercice précédent, tester l’association entre toutes les SNP et M.
modelisation = read.table('http://w3.mi.parisdescartes.fr/~vperduca/data/exo5.txt')Dans notre jeu de données, on trouve :
M: phénotype, 1 pour les cas et 0 pour les témoins E: covariable d’exposition environnementale, 1 pour les exposés, 0 pour les non exposés. Exemple: 1 pour les fumeurs et 0 pour les non fumeurs SNP1 - SNP100: génotype 0,1,2 pour 100 SNP bialléliques. La valeur du génotype indique le nombre d’allèles rares Testons l’allocation entre toutes les SNPs et M, pour cela on effectue donc un test du \(\chi^2\). On teste :
H0 : Le snp a un lien avec le phénotype H1 : Le snp n’a aucun lien avec le phénotype
data_snp = modelisation[,-(1:2)]
snp = names(data_snp)
pval = apply(data_snp, 2, function(x)chisq.test(table(x,modelisation$M))$p.value)#On décide de faire le test à un seuil de 0.05% :
snp_sign = pval[pval < 0.05]
snp_sign SNP29 SNP31 SNP40 SNP41 SNP42 SNP56
4.426461e-02 6.633656e-03 5.718667e-06 1.100742e-08 5.032654e-09 5.035117e-03
SNP59 SNP71 SNP81
2.938665e-02 1.491203e-02 2.325385e-02 On remarque que pour chaque test du \(\chi^2\) que nous avons effectué, qu’il y a 9 SNPs dont leurs pvaleurs sont strictement inférieure à \(α=0.05\) : on rejette H0 et on décide donc H1 avec seuil de confiance de 5%. En conclusion, parmi nos 100 SNPs de notre jeu de données, seulement 9 sont associés au phénotype : ce sont les SNPS 29, 31, 40, 41, 42, 56, 59, 71 et 81.
Question 2
Afficher le Manhattan plot avec la barre horizontale correspondante à la correction de Bonferroni. Quels sont les SNPs associé à la maladie?
datamodelisation = data.frame(SNP = snp, CHR = 1:length(snp), BP = 1:length(snp), P = pval )
manhattan(datamodelisation, highlight = names(snp_sign), annotatePval = 0.01)
De la même façon qu’à l’exercice 5, nous allons construire le Manhattan Plot avec la fonction ggplot. On affichera en vert les SNPs associés à la maladie.
datamodelisation= data.frame(SNP = snp, P = pval)
df_sign = filter(datamodelisation, P < (0.05/1191))
df_sign2 = filter(datamodelisation, P < 0.05)
ggplot(datamodelisation, aes(x = SNP, y = -log10(P))) +
geom_point()+
geom_point(data = df_sign2, col = 'green')+
geom_hline(aes(yintercept = -log10((0.05/1191)), linetype = 'Correction de Bonferroni'),col = 'red')+
geom_hline(aes(yintercept = -log10(0.05), linetype = 'Correction à 5%'),col = 'blue') +
scale_linetype_manual(name = 'Légende', values = c(2, 2),
guide = guide_legend(override.aes = list(color = c('blue', 'red'))))+
geom_label(data = df_sign, aes(label = SNP), size = 2.5,vjust = 0, hjust = 1.25)+
ggtitle('Manhattan Plot')+
theme_minimal()
A l’aide du Manhattan plot, on se rend compte qu’après correction de Bonferonni, les SNPs associés au phénotype sont les SNPs 40,41 et 42.
Question 3
On considère le SNP X pour le quel le signal d’association est plus forte (c’est à dire pour le quel le test d’association a donné la plus petite p-valeur). Ecrire l’équation du modèle de régression logistique de M sur X (utiliser le codage additif pour X). Estimer les paramètres du modèle, tester leur significativité et les interpréter en terme de Odds Ratios.
On rappelle que nos 9 SNPs associés à la maladie :
names(snp_sign)[1] "SNP29" "SNP31" "SNP40" "SNP41" "SNP42" "SNP56" "SNP59" "SNP71" "SNP81"names(which.min(snp_sign))[1] "SNP42"On remarque que celui qui a la plus petite pvaleur est le SNP42.
On souhaite maintenant savoir si le SNP42 influe sur M, c’est à dire s’il influe sur le phénotype. M est une variable binaire (elle prend les valeurs 0 ou 1) qui peut être modélisé par une loi de Bernouilli de paramètre pi∈[0,1]. On étudie alors le modèle suivant : \(logit(pi)=β0+β1X1,i\) Avec : - β0 qui correspond à l’ordonnée à l’origine du modèle - β1 qui est associé à la variable explicative - X1,i qui correspond à la variable explicative qui est donc notre SNP pour les i individus
X = modelisation$SNP42
reglog = glm(M~X, data = modelisation, family = binomial(logit))
summary(reglog)
Call:
glm(formula = M ~ X, family = binomial(logit), data = modelisation)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.506 -1.045 -1.045 1.089 1.316
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.32044 0.07918 -4.047 5.19e-05 ***
X 0.53265 0.08981 5.931 3.01e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1651.1 on 1190 degrees of freedom
Residual deviance: 1614.6 on 1189 degrees of freedom
AIC: 1618.6
Number of Fisher Scoring iterations: 4D’après les résultats de notre régression logistique : - On peut écrire le modèle sous la forme suivante : \(Y=−0.42+0.53X\) - On a dans la colonne estimate, les estimations de nos coefficients de notre régression logistique - On a dans la colonne \(Pr(>|z|)\), les pvaleurs du test de Student où l’on test : \(H0:"βk=0"\) \(H1:"βk≠0"\) Ce test permet de savoir si nos coefficients sont significatifs ou non, c’est à dire si la variable explicative associée joue un rôle dans la détermination du phénotype d’un individu. On voit que nos pvalues sont très faibles : nous allons donc rejeter H0 et ainsi accepter H1. En conclusion, nos coefficients sont significatifs : le SNP42 influe sur la détermination du phénotype d’un individu.
Regardons maintenant les odds ratio obtenu avec cette régression :
coeff = reglog$coefficients
exp(coeff)(Intercept) X
0.7258286 1.7034435 On remarque que notre OR entre X et M est supérieur à 1 : X influe bien sur M.
Question 4
On suspecte que l’effet de X sur M dépend de l’exposition E. Proposer un modèle pour vérifier cette hypothèse, estimer ses paramètres, et conclure.
Comme à la question précédente nous allons effectuer un modèle logistique, la seule différence avec celui ci-dessus est que nous rajoutons une variable qualitative au modèle afin de comprendre si l’effet de X sur M dépend de l’exposition E. On écrit donc le modèle 2 de la façon suivante : \(logit(pi)=β0+β1X1,i+β2X2,i\) Avec : - \(β0\) qui correspond à l’ordonnée à l’origine du modèle - \(β1\) et β2 qui sont associés aux variables explicatives - X1,i et X2,i qui correspondent aux variables explicatives qui sont donc notre SNP et E la covariable d’exposition environnementale pour les i individus
reglog2 = glm(M~X+E, data = modelisation, family = binomial(logit))
summary(reglog2)
Call:
glm(formula = M ~ X + E, family = binomial(logit), data = modelisation)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.678 -1.190 -1.000 1.128 1.366
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.43202 0.08695 -4.969 6.74e-07 ***
X 0.54875 0.09042 6.069 1.29e-09 ***
E 0.46218 0.14345 3.222 0.00127 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1651.1 on 1190 degrees of freedom
Residual deviance: 1604.1 on 1188 degrees of freedom
AIC: 1610.1
Number of Fisher Scoring iterations: 4D’après les résultats de notre régression logistique : - On peut écrire le modèle sous la forme suivante : \(Y=−0.43+0.55X1+0.46X2\) - On a dans la colonne estimate, les estimations de nos coefficients de notre régression logistique - On a dans la colonne \(Pr(>|z|)\), les pvaleurs du test de Student où l’on test : \(H0:"βk=0"\) \(H1:"βk≠0"\) Ce test permet de savoir si nos coefficients sont significatifs ou non, c’est à dire si les variables explicatives associées joue un rôle dans la détermination du phénotype d’un individu. On voit que nos pvalues sont très faibles : nous allons donc rejeter H0 et ainsi accepter H1. En conclusion, nos coefficients sont significatifs : le SNP42 ainsi que l’exposition environnementale influent sur la détermination du phénotype d’un individu.
Cette conclusion nous amène à penser que E est un facteur de confusion si E et X ne sont pas indépendants. Commençons par faire un teste du \(\chi^2\) : On teste : H0 : X est indépendant de E H1 : X est lié à E
chisq.test(X, modelisation$E)
Pearson's Chi-squared test
data: X and modelisation$E
X-squared = 1.9467, df = 2, p-value = 0.3778On obtient comme p-value 0.3778 qui est significativement plus grande que α=0.05 : on accepte donc H0 avec un seuil de confiance à 5%. En conclusion, X et E sont indépendants. Nous allons maintenant faire une régression logistique et calculer l’OR correspondant pour confirmer nos propos :
reglog3 = glm(E~X, data = modelisation, family = binomial(logit))
summary(reglog3)
Call:
glm(formula = E ~ X, family = binomial(logit), data = modelisation)
Deviance Residuals:
Min 1Q Median 3Q Max
-0.7350 -0.7350 -0.6892 -0.6456 1.8280
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.17073 0.09225 -12.69 <2e-16 ***
X -0.14585 0.10645 -1.37 0.171
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1260.1 on 1190 degrees of freedom
Residual deviance: 1258.2 on 1189 degrees of freedom
AIC: 1262.2
Number of Fisher Scoring iterations: 4coeff3 = reglog3$coefficients
exp(coeff3)(Intercept) X
0.3101410 0.8642857 On voit que le coefficient lié à la variable explicative X n’est pas significatif. De plus, on remarque que l’OR entre X et E est proche de 1, E n’a donc aucun effet sur X et on peut dire que X et E sont bien indépendants. En conclusion, E la covariable d’exposition environnementale n’influe par sur X, notre SNP42 : E n’est donc pas un facteur de confusion.
On a montré à la question précédente que X influe sur M, le phénotype de l’individu i. Il reste donc à savoir si M et E sont indépendants ou non, pour cela nous ferons un test du χ2. On teste : H0 : E est indépendant de M H1 : E est lié à M
chisq.test(modelisation$E,modelisation$M)
Pearson's Chi-squared test with Yates' continuity correction
data: modelisation$E and modelisation$M
X-squared = 8.2695, df = 1, p-value = 0.004032On obtient comme p-value 0.004 qui est significativement plus petite que α=0.05 : on rejette H0 et donc on accepte H1 avec un seuil de confiance à 5%. En conclusion, E est lié à M. Nous allons maintenant faire une régression logistique et calculer l’OR correspondant pour confirmer nos propos :
reglog4 = glm(M~E, data = modelisation, family = binomial(logit))
summary(reglog4)
Call:
glm(formula = M ~ E, family = binomial(logit), data = modelisation)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.316 -1.138 -1.138 1.217 1.217
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.09284 0.06576 -1.412 0.15801
E 0.41375 0.14095 2.935 0.00333 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1651.1 on 1190 degrees of freedom
Residual deviance: 1642.4 on 1189 degrees of freedom
AIC: 1646.4
Number of Fisher Scoring iterations: 3coeff4 = reglog4$coefficients
exp(coeff4)(Intercept) E
0.9113402 1.5124740 On voit que le coefficient lié à la variable explicative E est significatif. De plus, on remarque que l’OR entre M et E est supérieur à 1 : E a donc un effet sur M et on peut dire que M et E sont bien dépendants. En conclusion, E la covariable d’exposition environnementale influe par M, le phénotype d’un individu i.
Pour conclure sur la question de savoir si l’effet de X sur M dépend de l’exposition E : - On peut dire que E n’influe pas X : ils ne sont pas liés - On peut dire que E influe sur M