Evaluation Pratique 2
Voici l’evaluation pratique : 3 etapes:
Lire cette documentation et preparer le code avant l’evaluation finale.
Execution du code a partir des donnees fournies par Greg: lundi 20 avril (90 minutes).
Revision du travail (si necessaire): REMISE avant le 24 avril a minuit!
OBJECTIFS de cette evaluation:
- tester un jeu de donnees reel.
- se confronter a une nouvelle approche: RDA
- apprendre a resoudre les problemes.
Librairies et fonctions necessaires pour ce travail:
(“AnnotationHub”,“knitr”,“GenomicRanges”, “vegan”, “limma”, “biomaRt”)
Charger l’objet R dand votre session disponible pour l’evaluation.
1. Description des donnees: numeriques et label.
Quel est le nombre de replicat par groupe? groupe DG, groupe n_DG, DIETE, INSULINE, et combinaisons DG_DIETE et DG_INSULIN
dim(exprData)## [1] 46889 44head(colnames(exprData))## [1] "nDG_NA_1" "DG_INSULIN_2" "DG_DIETE_3" "DG_INSULIN_4"
## [5] "nDG_NA_5" "DG_INSULIN_6"head(rownames(exprData))## Ku8QhfS0n_hIOABXuE fqPEquJRRlSVSfL.8A ckiehnugOno9d7vf1Q
## "346389" "1" "1"
## x57Vw5B5Fbt5JUnQkI ritxUH.kuHlYqjozpE QpE5UiUgmJOJEkPXpc
## "29974" "29974" "29974"Utiliser la fonction summary pour decrire les donnees numeriques.
La distribution des échantillons est conséquente avec ce a quoi on s’attenderais, alors que pour chaque groupe la mediane et le range sont semblables. Pour tous les patients des valeurs d’environ: Min : 7.1 1st Qu.:7.35 Median : 7.5 Mean : 7.8 3rd Qu.:7.7 Max :15.0
Ces valeurs suggèrent aussi qu’il y a beaucoup de gènes qui s’expriment avec des lectures autour de 7.8 et quelques veleurs sont significativement plus élevées, à des valeurs allant jusqu’à 15.
Les rownames sont ils uniques?
Les noms de colonnes sont uniques et correspondent à 44 patients répartis en 3 groupes. Les noms de colonnes donnent aussi des informations sur l’état de santé des patients, par rapport au DG/nDG et au type de traitement (insuline, diète, NA):
DG : diabete gestationnel
nDG: non diabete gestationnel
DIETE/INSULIN: traitement applique chez les patients “diabete gestationnel”
NA : No Applicable (pas de traitement pour le diabete gestationnel)
Le nom des lignes est en deux parties, un entrezID numérique (numérique, jusqu’à 6 chiffres) et un Illumin_ID (alphanumérique du type: 01LVKJbtyL14UlF95Q). Le entrezID n’est pas unique, mais le Illumin_ID
A quoi vous font penser les rownames?
Ce sont des transcrits de gènes entrezID associés à des sondes Illumin_ID
2. Sélectionner au hasard 10000 sondes (lignes) en utilisant la fonction sample.
Cette approche nous permet d’optimiser le temps de calcul et de travailler sur une selection non biaisee.
exprDataSample = exprData[sample(nrow(exprData), 10000), ]
class(exprDataSample)## [1] "matrix"dim(exprDataSample)## [1] 10000 44Ce qui est attendu: une matrice de 10000 lignes appelee exprDataSample
3. former 3 groupes de donnees a partir de label (nDG, DG_DIETE, DG_INSULIN) appelee “mat_nDG”, “mat_DG_DIETE” et “mat_INSULIN” en utilisant la fonction grep
Ce qui est attendu: 3 matrices correspondant aux labels d’interet (nDG, DG_DIETE, DG_INSULIN)
4. boxplot
Presenter les donnees dans UN boxplot representant les groupes par 3 couleurs differentes. De plus, les labels dans les axes des X doivent etre representer verticalement (quel option pouvait vous utiliser) et comment faire pour que les labels soient visibles dans leur totalite
Ce qui est attendu : un boxplot avec les criteres demandes.
Ces boxplots démontrent un bon QC, alors que les données semblent constantes entre les patients. Cela confirme aussi ce qui a été souligné plus tôt, soit que la majorité des valeurs se situent sous 8 et quelques gènes ont valeurs beaucoup plus élevées (les points au haut des boites à moustaches).
5. dendogramme
Representer les echantillons dans un dendrogramme apres un clustering (et non les genes)
Dendrogramme avec un groupe surligné
Conclusions sur ce dendogramme a propos de l’organisation (la structure) des donnees en fonction des labels:
Les groupes ne sont pas fortement classifiés, alors que les feuilles voisines proviennent surtout de groupes séparés. Il y a une petit portion où les traitements sont regroupés, annoté par le carré jaune sur la deuxième figure. Les résultats des patients DG_INSULIN_2, DG_INSULIN_6, DG_INSULIN_9, DG_INSULIN_25, DG_INSULIN_36 sont tous adjacents et pourraient être considérée comme classés correctement.
6. PCA plot
QUESTIONS RDA:
A la vue du graphique, que pensez-vous de l’ordination des echantillons.
Les échantillons semblent encore ici (comme avec le Dendrogramme) difficile à regrouper sur la base de ce graphique. Les données ne semblent pas se séparer latéralement en groupes, selon l’axe de la PC1. Elles ne semblent pas non plus se séparer horizontalement, selon l’axe de la PC2. Il serait possible de distinguer un groupe plus ou moins clair en haut à droite, mais avec des ééchantillons associés à différents labels regroupés, il n’est pas possible de faire de conclusions.
Quelle est la valeur de PC1? qu’est ce que cela signifie?
La valeur d’inertie associée à la PC1 est d’environ 13,3% (32.43/243.3 : inertie PC1 / inertie totale). Cela veut dire que la PC1 explique 13,3% de la variation des données.
Quelle est la valeur de PC2? Est ce que PC1 ou PC2 semble etre explique par les conditions experimentales?
La valeur d’inertie associée à la PC2 est d’environ 12,3% (29.96/243.3 : inertie PC2 / inertie totale). Cela veut dire que la PC2 explique 12,3% de la variation des données. En somme, même avec les deux premières composantes principales, on ne peut qu’expliquer 25,6% de la variance totale. Les données ne semblent pas suivre les axes tracées par les conditions expérimentales sur le graphique de RDA présenté ci-dessous.
RDA: une extension de la PCA
Quelques explications sur RDA:
Ici, on structure les inputs pour l’analyse RDA.Pour cette analyse, on a besoin de la matrice de donnees (exprDataSample) et les vecteurs decrivant les niveaux des facteurs. C’est ce que fait ce code.
L’analyse proprement dite
Quel est la valeur de la variance totale?
243.326 (Total_variance)
Quelle est la valeur de la variance expliquee par les conditions experimentales.
Il y a environ 6% de la variance qui est expliquée par le modèle (Constrained_Proportion: 0.061) et donc 94% de non expliqué.
Le biplot de RDA:
La fonction goodness de rda:
Analyse differentielle des donnees d’expression: etude de l’effet de l’insuline sur le diabete gestationnel.
A partir de l’object exprData,
1. creez 2 matrices de donnees:
-une contenant les donnees diabete gestionnel-diete (label = “DG_DIETE”). Cette matrice s’appelera exprData_DGD.
-une contenant les donnees diabete gestionnel-insuline (label = “DG_INSULIN”). Cette matrice s’appelera exprData_DGI.
Au final: associez les 2 matrices en une matrice. Cette matrice s’appellera exprData_DG.
- creer la matrice de design: a partir des information precedentes.
Ce qui est attendu: object designLIMMA
- appliquer le modele lineaire de LIMMA
fit <- lmFit(exprData_DG,designLIMMA)
fit2= eBayes(fit)- extraire le top table et selectionnez des genes avec une Pvalue inferieur a 0.01.
res_TopTable = topTable(fit2 , coef=1 , adjust="fdr", p.value=0.01, number = 46889)
Il y a 46 889 gènes avec une variance significativement différente entre les groupes DG_INSULIN et DG_DIETE. Ils sont contenus dans l’objet “res_TopTable”
- A partir de l’input , recuperer le geneSymbol et le ensembl transcript ID. Utiliser BioMart le output s’appellera geneID.annotated
L’objet “geneID.annotated” est une liste des gènes sélectionnés par la TopTable de l’analyse différentielle, accompagnés de leur annotations respectives (“hgnc_symbol”,“entrezgene”,“ensembl_transcript_id”, “description”,“start_position”,“end_position”). C’est une liste de 88 774 gènes.
- A partir de l’objectgeneID.annotated, il faut identifier les CpGs qui sont associe aux transcripts selectionnees.
Decrivez l’object ensGene.
REPONSE: GRanges object with 63 280 ranges and 5 metadata columns, C’est une liste de tous les positions des gènes (transcrits et leurs variants) du génome humain Hg18. C’est une base de donnée de bioconductor qui utilise les informations sur les transcrits de Ensembl.
Decrivez ce que represente l’object mask.
REPONSE: Ce mask est l’intercept en commun de EnsGene avec les GeneID sélectionnés pus tôt. Cela permet de créer l’objet “ensGene_select”, qui contient les gènes sélectionnée avec toutes les informations associées au transcrit sélectionné. length(ensGene_select) = 18 457 indique le nombre de gènes sélectionnés et bien annotés.
acceder aux annotations.
Quels types d’information sont fournies par la librarie “FDb.InfiniumMethylation.hg18”
REPONSE (en quelques mots):
C’est une liste de 487 117 “GRanges object”" (avec 14 colonnes de metadata) qui correspondent à autant de sites de méthylation connus dans le génome humain.
Compare les 2 objects avec la fonction findoverlap
Comment interpretation le output de la fonction findoverlap?
REPONSE: La liste overlap_info contient 278 574 entrées, qui correspondent aux CpG répertoriés qui ont été retrouvés dans l’environnement immédiat (+/- 2000pb) des gènes sélectionnés plus tôt.
A ce stade nous avons selectionne des CpGs d’interet en accord avec les donnees transcriptomiques. Ces CpGs pourraient etre interessants pour continuer l’analyse.
FIN DE L’EVALUATION PRATIQUE 2.
Comment interpreter ces valeurs numeriques?
Le résultat de “goodness(rda.vegan)” est une liste des gènes les mieux expliqués par les composantes principales de la RDA. Ainsi, entrezID = “650058” est expliqué à 2.25% par la composante DRA1 et à 8.0% par la RDA2. En fait, dans cette liste de head, c’est la RDA2 qui semble le plus expliquer la variance observée pour ces gènes.