# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .r
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: R
#     language: R
#     name: ir
# ---

# # M1 MEG - UE1 
# # TP Transcriptomique 1 - TP5 - 2021-2022

#
# ## Analyse de données de puces (MicroArrays)
#
# ### Adaptation du TP MeV Testis Vs Ovary, sous forme de Jupyter Notebook en R
#
# Selon les explications du site https://wiki.bits.vib.be/index.php/Analyze_your_own_microarray_data_in_R/Bioconductor

# ---
# **Avant toute chose**
#
# <mark>Avant de démarrer cette analyse, sauvegardez une copie de ce notebook pour backup. <mark>
#
# Vous pouvez aussi faire des copies de ce notebook au fur et à mesure de la progression
#
# ---
#
# __*Rappels sur l'utilisation des notebooks*__
#
# - Pour ajouter une nouvelle cellule, cliquez sur l'icône "+" dans la barre des menus 
# - Vous pouvez "cliquer-glisser" pour bouger une cellule 
# - Vous pouvez choisir le type de cellule dans le petit menu déroulant dans la barre des menus : 
#     - 'Code' pour entrer des lignes de commandes à executer 
#     - 'Markdown cells' pour ajouter simplement du texte, formatable avec quelques signes 
# - Pour executer une cellule 'Code', pressez SHIFT+ENTER ou cliquez sur l'icône "play" 
# - Pour afficher une cellule 'Markdown', pressez SHIFT+ENTER ou cliquez sur l'icône "play" 
# - Pour modifier une cellule 'Markdown', double-cliquez dessus
#
# ---

# ## I - Préparation de l'espace de travail
#
# Ceci est un jupyter notebook en R, ce qui signifie que les commandes seront des commandes R directement interprétables par Plasma. 
# Comme pour tout projet en R, nous devons d'abord choisir le répertoire de travail. Par défaut, ce sera celui de l'environnement déployé lors de l'ouverture du TP.
# Puis nous devons indiquer à R les paquets à charger qui seront nécessaires pour l'analyse. 
#

# +
# Le répertoire de travail est d'office attribué au répertoire de l'environnement par Jupyter Hub
# La commande suivante n'a donc pas besoin d'être exécutée.
# setwd("~/m1meg-ue1-tp5-r/")

# +
## Cellule de Code n°1 ##

# Chargement des packages dans l'environnement présent.
# Vous allez voir un grand nombre de lignes en rose s'afficher, c'est normal !

library(Matrix)
library(lattice)
library(fdrtool)
library(rpart)
library(Biobase)
library(Biostrings)
library(mouse4302cdf)
library(ggplot2)
library(limma)
library(affy)
library(affyPLM)
library(simpleaffy)
library(genefilter)
# -

# Nous pouvons vérifier ces éléments avec les commandes suivantes : 

# +
## Cellule de Code n°2 ##

getwd()
sessionInfo()
# -

# ---
# ## II - Analyse des puces 
#
# Nous allons analyser un ensemble de données d’expression, obtenues sur puces d’expression de type Affymetrix. Il s’agit d’une portion de la série **GSE12989**, accessible par une des deux bases de données d’expression génique (“repository”) : Gene Expression Omnibus (GEO) http://www.ncbi.nlm.nih.gov/geo/  
# Sur GEO, un ensemble d'échantillons composent une **Série**, avec un identifiant commençant par GSE.
#
# Ouvrez le lien https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12989
#
# La description de la série indique :
#
# Status 			Public on Jun 24, 2009
# Title 			Foxl2 functions throughout mouse ovary development
# Organism 		Mus musculus
# Experiment type 	Expression profiling by array
# Summary 		This SuperSeries is composed of the following subset Series:
# - GSE12905: Foxl2 functions in sex determination and histogenesis throughout mouse ovary development, analyzed by Affymetrix arrays
# - GSE12942: Foxl2 functions in sex determination and histogenesis throughout mouse ovary development, analyzed by Agilent arrays
# Citation(s) 
# •	<NAME>, <NAME>, <NAME>, <NAME> et al. Foxl2 functions in sex determination and histogenesis throughout mouse ovary development. BMC Dev Biol 2009 Jun 18;9:36. PMID: 19538736
#
# Les puces Affymetrix utilisées dans le subset GSE12905 suivent la plateforme (design de la puce) GPL1261, c’est à dire le Affymetrix Mouse Genome 430 2.0 Array (nom commercial de la puce).
#

# > <mark>**Question 1 :**<mark>  
# > Combien d’échantillons analysés avec cette puce sont présents dans GEO ? 

# *Double-cliquez ici pour entrer votre réponse ou des notes*  

# ### 1 - Récupération des données
#
# Les fichiers CEL qui nous intéressent ont déjà été télécharchés sur Plasmabio. Il suffit donc d'indiquer à R où aller les chercher, et stocker cette information dans la variable *celpath*
#

# +
## Cellule de Code n°3 ##

# Chemin d'accès aux données sur plasmabio :
celpath = "/srv/data/meg-m1-ue1/DataTP5/"
# -

# Les puces commerciales Affymetrix étant très répandues, des "packages" R dédiés ont été créés pour analyser les données issues de ces puces. C'est le cas des librairies affy, affyPLM et simpleaffy que nous avons chargé plus haut.    
# Dans ces packages, il existe donc des fonctions spécialement créées pour prendre en charge ces données. Nous allons en utiliser quelques unes. 
# Nous allons maintenant charger les données de ces fichiers dans un objet R approprié : AffyBatch

# +
## Cellule de Code n°4 ##

# Import des données d'intensité de fluorescence brutes des sondes depuis les fichiers CEL, 
# et ajout des données dans un objet R AffyBatch
data = ReadAffy(celfile.path=celpath)
# -

# Cet objet AffyBatch est accessible via des méthodes contenues dans le package affy.

# +
## Cellule de Code n°5 ##

# Récupération de l'annotation des échantillons, qu'on met dans un objet appelé ph (pour phénotype)
ph = <EMAIL>@<EMAIL>
ph
ph@data
# -

# > <mark>**Question 2 :**<mark>  
# > Combien d’échantillons de testicules et d’ovaires sont présents ? 
# Combien y a t-il de réplicats par stade de développement ? 

# Remplir votre réponse

# *Double-cliquez ici pour entrer votre réponse ou des notes. Par la suite, il vous suffira de sélectionner la cellule de la question, puis d'ajouter une cellule markdown avec le "+" pour entrez vos réponses*   
#

# +
## Cellule de Code n°6 ##

#Récupération de l'annotation des sondes
feat = data@featureData


# +
## Cellule de Code n°7 ##

# Récupération du nom du fichier CDF associé aux puces
cdfName(data)

# +
## Cellule de Code n°8 ##

# Récupération du nombre de probe sets représentés sur les puces
# (La fonction length() compte le nombre d'items dans un vecteur)
length(featureNames(data))
# -

# > <mark>**Question 3 :**<mark>  
# > Combien de probes sets sont analysés grâce à cette puce ? 

# ####   
#
# ### 2 - Contrôle qualité des données de puces
#
# Nous commençons par attribuer des noms un peu plus gérables aux échantillons :
#

# +
## Cellule de Code n°9 ##

# Attribution de noms informatifs aux échantillons
ph@data[ ,1] = c("Testis_E13_1","Testis_E13_2","Testis_E13_3","Testis_P0_1","Testis_P0_2","Testis_P0_3","Ovary_E13_1","Ovary_E13_2","Ovary_E13_3", "Ovary_P0_1", "Ovary_P0_2", "Ovary_P0_3")
ph@data
# -

# Une première étape importante d'une analyse de données est de se rendre compte de l'aspect global des données, en termes de répartition et de distribution.   
# Pour celà, nous allons utiliser différents modes de représentation graphique de nos données.

# #### 2.1 - Histogrammes de distribution des données 
# Puis nous allons créer 2 histogrammes de la distribution des données (avec ggplot2).
# Le premier affichera la répartition des données brutes, le second celle des valeurs passées en log2.

# +
## Cellule de Code n°10 ##

# Récupérations des intensités PM (sondes Perfect Match dans les probe sets)
pmexp = pm(data)

# La méthode dim() fournit le nombre d'éléments dans la matrice, 
# c'est à dire le nomdre de sondes (lignes) et le nombre d'échantillons (colonnes)
dim(pmexp)
# -

# On crée des vecteurs qui vont nous servir comme colonnes de la dataframe pour les histogrammes
#   - un pour y mettre les noms des échantillons, appelé **sampleNames**
#   - un pour y stocker les intensités de fluorescences brutes, applelé **nologs**
#   - un pour y stocker les intensités passées en log2, appelé **logs**

# +
## Cellule de Code n°11 ##

# initialisation de 3 vecteurs
sampleNames = vector()
nologs = vector()
logs = vector()

# remplissage des vecteurs avec les données (la boucle for permet de parcourir les 12 colonnes)
for (i in 1:12)
{
  sampleNames = c(sampleNames,rep(ph@data[i,1],dim(pmexp)[1]))
  nologs = c(nologs, pmexp[,i])
  logs = c(logs,log2(pmexp[,i]))
}

# Vérification avec les premières valeurs
nologs[1:10]
logs[1:10]

# +
## Cellule de Code n°12 ##

# Après avoir rempli les vecteurs, nous combinons sampleNames et nologs ou logs dans 2 dataframes
nologData = data.frame(sampleName=sampleNames, nologInt=nologs)
logData = data.frame(sampleName=sampleNames, logInt=logs)

# Vérification avec les premières valeurs
nologData[1:5,]
logData[1:5,]

# On vérifie le nombre d'éléments dans nos 2 dataframes
dim(nologData)
dim(logData)

# -

# Maintenant nous pouvons créer les 2 histogrammes: 

# +
## Cellule de Code n°13 ##

# Nous décidons de la taille des graphs 
options(repr.plot.width=18, repr.plot.height= 7)

# Histogramme avec les valeurs brutes
ggplot(nologData, aes(nologInt, colour = sampleName)) +
 geom_density() +
 theme_bw(base_size = 20)


# +
## Cellule de Code n°14 ##

# Histogramme avec les valeurs en log2
ggplot(logData, aes(logInt, colour = sampleName)) +
 geom_density() +
 theme_bw(base_size = 20)
# -

# > <mark>**Question 4 :**<mark>  
# > Que pouvez-vous dire de la distribution des données ? Vous semble t-elle logique pour des transcriptomes ?

# #### 2.2 - Boxplots de distribution des données 
# Nous utilisons les mêmes dataframes pour créer 2 boxplots, 
# de façon à comparer la répartition des données entre échantillons
#

# +
## Cellule de Code n°15 ##

# On choisit une taille de graph moins importante (les valeurs choisies précédemment restent valables tant qu'elles ne sont pas modifiées)
options(repr.plot.width=12, repr.plot.height= 7)

# Création du boxplot
ggplot(nologData,aes(sampleName,nologInt)) + 
 geom_boxplot() +
 theme_grey(base_size = 18) +
 theme(axis.text.x = element_text(size = 10)) +
 ggtitle("Données noLog - Avant normalisation")

# +
## Cellule de Code n°15 ##

# Même grap pour les valeurs en log2
ggplot(logData,aes(sampleName,logInt)) + 
 geom_boxplot() +
 theme_grey(base_size = 18) +
 theme(axis.text.x = element_text(size = 10)) +
 ggtitle("Données Log - Avant normalisation")
# -

# > <mark>**Question 5 :**<mark>  
# > Pouvez-vous déterminer si les données ont été normalisées entre puces ou pas ?  
# > Quel est l'intérêt de la transformation log2 pour ces 2 types de représentation ?    

# #### 2.3 - MA plots
# Les plots MA permettent de visualiser la variablité de l'expression des gènes, et sa répartition en fonction du niveau d'expression.   
# Ici, chaque puce est compré à une *pseudo-puce moyenne*, pour laquelle chaque sonde est associée à la valeur médiane de son intensité dans tous les échantillons.
#
# Dans un plot MA, M est figuré en fonction de A :
#
#     - M est la différence d'intensité d'une sonde sur la puce testée et l'intensité de cette même sonde sur la pseudo-puce
#     M = logPMInt_array - logPMInt_pseudoarray
#     - A est la moyenne des intensités de cette sonde sur la puce testée et sur la pseudo-puce
#     A = (logPMInt_array + logPMInt_pseudoarray)/2   
#     
# (les plots se créent dans des fichier séparés, colonne de gauche)

# +
## Cellule de Code n°16 ##

# Création des MA plots pour chaque échantillon.
for (i in 1:12)
{
  name = paste("MAplot",i,".jpg",sep="")
  jpeg(name)
  MAplot(data,which=i)
  dev.off()
}

# -

# Ouvrez les MA plots et examinez-les.
#
# > <mark>**Question 6 :**<mark>  
# > Que pouvez-vous constater ? Pensez-vous qu'il faille normaliser les donnéees ?

# ####   
#
#
# ### 3 - Normalisation des données
# Comme indiqué en cours, la méthode RMA est la méthode standard pour analyser des données de microarray.   
#
# RMA est une des rares méthodes de normalisation qui utilise seulement les sondes PM.  
# - Correction du bruit de fond, pour corriger la variation spatiale au sein de chaque puce : une intensité corrigée est calculée pour chaque sonde PM de façon à ce que toutes les intensités corrigées soient positives.  
# - Transformation log pour améliorer la distribution des données : le logarithme en base 2 est calculé pour chaque sonde sur les intensités corrigées du bruit de fond.  
# La transformation en log2 va permettre que les données soient moins "tordues" et suivent une distribution plus normale, et va faire en sorte que les ratios d'up-régulation et de down-régulation soient répartis de façon égale.     
# - Normalisation par la méthode des quantiles pour corriger les variations entre les puces : égalise la distribution des données entre les échantillons pour les rendre complètement comparables.   
# - Normalisation des sondes pour corriger les variations au sein des probesets : égalise le comportement des sondes entre les puces, et combine les valeurs normalisées des sondes d'un même probeset en une seule valeur pour le probeset.  
#
# La méthode GCRMA est basée sur RMA, avec tous ses aspects positifs.    
# La différence porte sur la correction du bruit de fond, tous les autres étapes de l'analyse sont identiques.  
# Pour le bruit de fond, GCRMA corrige l'hybridation non-spécifique, contrairement à RMA qui ignore ce point.  
# GCRMA prend en compte l'information de séquence des sondes (d'où le "GC") pour estimer la susceptibilité de chaque sonde à l'hybridation non-spécifique.  

# *Normalement, nous utiliserions la méthode gcrma directement dans R, grâce à cette commande :   
# library("gcrma")  
# data.gcrma = gcrma(data)*
#
# *Pour des raisons techniques en cours de résolution, cela ne marche pas actuellement sur Plasma.*
#
# *Nous allons donc charger directement le résultat de la normalisation pour la suite de l'analyse*
#

# +
## Cellule de Code n°17 ##

# Chargement de la library, nécessaire pour la suite
# library("gcrma")

# Chargement de l'objet "data_gcrma.RData", pour le réintégrer dans la procédure d'analyse
load("data_gcrma.RData")
# -

# #### 
#
# #### 3.1 - Effet de la normalisation GCRMA : box-plots
#
# Nous allons maintenant vérifier l'effet de la normalisation, avec un boxplot des données normalisées.

# +
## Cellule de Code n°18 ##

# Box plots des valeurs normalisées
# création d'une data frame avec seulement les valeurs d'expression normalisées
normexpr.gcrma = data.frame(exprs(data.gcrma))
normexpr.gcrma[1:5,]


# initialisation des vecteurs
sampleNames = vector()
normlogs = vector()

# remplissage des vecteurs
for (i in 1:12)
{
    sampleNames = c(sampleNames,rep(ph@data[i,1],dim(data.gcrma)[1]))
    normlogs = c(normlogs,normexpr.gcrma[,i])
}

normlogs[1:10]

# création de la dataframe avec les 2 vecteurs
normData = data.frame(sampleName=sampleNames, norm_logInt=normlogs)
normData[1:5,]

# création des graphs
ggplot(normData,aes(sampleName,norm_logInt)) + 
 geom_boxplot() +
 theme_grey(base_size = 18) +
 theme(axis.text.x = element_text(size = 10)) +
 ggtitle("Données Log - Après normalisation")

# -

# > <mark>**Question 7 :**<mark>
# Que pouvez-vous constater ? Pensez-vous que la normalisation entre échantillons a fonctionné ?

# ####   
#
# #### 3.2 - Effet de la normalisation GCRMA : MA plots
#
# Pour voir ce que la normalisation a fait aux données de chaque échantillon, nous allons refaire les MA plots.

# +
## Cellule de Code n°19 ##
# MAplots des données normalisées 

for (i in 1:12)
{
  name = paste("MAplot-normData",i,".jpg",sep="")
  jpeg(name)
  MAplot(data.gcrma,which=i)
  dev.off()
}

# -

# > <mark>**Question 8 :**<mark>
# Qu'en pensez-vous ?

# ####   
#
#
# ### 4 - Groupement des échantillons : PCA
#
# Pour voir si les échantillons se regroupent de façon cohérente après normalisation, nous allons tenter une analyse PCA (Principal Componant Analysis).  
# Cette analyse vous sera expliquée plus en détail lors du TP PCA de l'UE5 (avec <NAME>).  
#

# +
## Cellule de Code n°20 ##

# Nous utilisons la méthode prcomp() pour lire la matrice de données issues de la normalisation (i.e normexpr.gcrma)
# (Plus d'explications là : https://www.youtube.com/watch?v=0Jp4gsfOLMs)

data.PCA = prcomp(t(normexpr.gcrma),scale.=FALSE)

# la part de la variation expliquée par chaque axe est le carré de sdev (standard deviation), un des éléments fournis par la fonction prcomp().
pca.var = data.PCA$sdev^2
# et on l'exprime en % du total :
pca.var.per = round(pca.var/sum(pca.var)*100, 1)


# On peut maintenant afficher la PCA avec ggplot2, qui permet de générer un plot avec les infos calculées :
# formatage des données comme ggplot2 les aime : une dataframe avec les coordonnées de chaque échantillon
pca.datagg = data.frame(Sample = ph@data[ ,1],
                        X = data.PCA$x[,1],
                        Y = data.PCA$x[,2])


# PCA plot selon les 2premiers axes avec ggplot2, les différents éléments indiqués permettent de jouer sur l'apparence du plot
ggplot(data = pca.datagg, aes (x=X, y=Y, label = Sample, colour = Sample)) +
  xlab(paste ("PC1 - ", pca.var.per[1], " %", sep = "")) +
  ylab(paste ("PC2 - ", pca.var.per[2], " %", sep = "")) +
  geom_text(size = 5, nudge_y = -8) +
  geom_point(size=4) +
  theme_classic() +
  ggtitle("Données normalisées avec GCRMA, plot PCA") +
  theme(text=element_text(size=20)) 

# -

# > <mark>**Question 9 :**<mark>  
# > Que pensez-vous de la répartition des réplicats ?   
# > D'après vous, à quel paramètre biologique est dûe la variation présentée sur le 1er axe ?   
# > Sur le 2ème ?

# ####   
#
# ### 5 - Identification des gènes DE (Différentiellement Exprimés)
#
# Nous allons maintenant identifier les gènes dont l'expression varie entre les conditions.   
# Nous allons d'abord effectuer une comparaison simple entre 2 groupes, ici le sexe gonadique.  
#
#

# #### 5.1 - Comparaison simple entre 2 groupes   
#
# La méthode utilisée ici est la méthode limma, qui permet d'identifier les gènes différentiellement exprimés au moyen de régression linéaire et de modèle Bayésien.  
# Nous attribuons à chaque échantillon sa valeur selon le critère choisi ("Testis" ou "Ovary"), puis nous indiquons à R que c'est cet élément qui va nous servir de facteur de comparaison.  
# Puis une matrice de contraste est créée pour comparer les 2 groupes.
#

# +
## Cellule de Code n°21 ##

# Ajout d'une colonne à la dataframe data pour caractériser les échantillons
ph@data[ ,2] = c("Testis","Testis","Testis","Testis","Testis","Testis", "Ovary", "Ovary", "Ovary", "Ovary", "Ovary", "Ovary")
colnames(ph@data)[2]="SexGonad"

# On vérifie que le facteur a été correctement attribué
ph@data

# +
## Cellule de Code n°22 ##

# On forme les 2 groupes
groups = ph@data$SexGonad

# On utilise ces groupes pour indiquer le facteur de classement, ainsi que les niveaux nécessaires (ici seulement 2 niveaux)
f = factor(groups,levels=c("Testis","Ovary"))
f

# +
## Cellule de Code n°23 ##

# Création de la matrice de contraste, qui permet de comparer les groupes selon ce facteur
design = model.matrix(~0 + f)
colnames(design) = c("Testis","Ovary")

data.fit = lmFit(data.gcrma,design)

# montre les 10 premières valeurs dans les 2 groupes
data.fit$coefficients[1:10,]

contrast.matrix = makeContrasts(Ovary-Testis,levels=design)
data.fit.con = contrasts.fit(data.fit,contrast.matrix)
data.fit.eb = eBayes(data.fit.con)
# -

# On affiche ensuite les résultats pour les 10 gènes les plus significativement différents :

# +
## Cellule de Code n°24 ##

# indique le noms des colonnes dans la matrice de résultats
names(data.fit.eb)

# montre les p-values associées au test pour les 10 premiers gènes
data.fit.eb$p.value[1:10,]
# -

# ### 5.2 - Volcano plot   
#
# On peut visualiser l'ensemble des p-values obtenues avec un **Volcano Plot**  
# (se crée aussi dans un fichier séparé, colonne de gauche)

# +
## Cellule de Code n°25 ##

# Volcano plot
name = "Volcano1.jpg"
jpeg(name)
volcanoplot(data.fit.eb,coef=1,highlight=50, style = "p-value", xlab = "log2 Fold Change", ylab = "-log10 p", names = NULL, hl.col = "blue")
dev.off()

# -

# ><mark>**Question 10 :**<mark>  
# >    Pensez-vous que les p-values soient interprétables directement ? Quelle étape importante doit d'abord être effectuée ?

# ####   
#
# ### 5.3 - Correction des p-values pour les tests multiples
#
# La fonction topTable() renvoie une table avec les gènes ordonnés selon le critère de p-value.  
# De plus, la fonction topTable() va ajuster la p-value obtenue, selon la methode choisie avec l'argument adjust.method. Ici, cet ajustement est effectué en utilisant la méthode de Benjamini and Hochberg, la plus utilisée pour analyser les transcriptomes sur puces. 
#

# +
## Cellule de Code n°26 ##

options(digits=2)
tab = topTable(data.fit.eb,coef=1,number=5000,adjust.method="BH")

# on affiche les 10 premier gènes, et 10 gènes pris plus loin dans le tableau
tab[1:10,]
tab[300:310,]

# la colonne B indique le log odds scores, 
# c'est à dire le logarithme de la vraisemblance de ll'hypothèse "expression différentielle" par rapport à l'hypothèse nulle
# (ratio des vraisemblances, transformée en log).
# Souvenez-vous des cours et TP sur la méthode des lod-scores en génétique humaine de L3 !)

# -

# On peut se servir de la même fonction topTable() pour sélectionner les gènes selon des seuils, par exemple de p-value ajustée.

# +
## Cellule de Code n°27 ##

# Sélection des gènes selon la adjusted p-value, sous un certain seuil

topgenesAdjP = tab[tab[, "adj.P.Val"] < 0.001, ]

# nombre de gènes correspondant au critère :
dim(topgenesAdjP)[1]

topgenesAdjP[1:10,]


# -

# ><mark>**Question 11 :**<mark>  
# >    Indiquez le nombre de gènes qui sont en dessous des seuils de p-value suivants :  
# >     adj p < 0.05  
# >     adj p < 0.01  
# >     adj p < 0.001  

# +
## Cellule de Code n°28 ##

# pour séparer les gènes up- and down-regulated, on inclut un seuil sur le log fold change (colonne logFC) :
  
topUpAdjP = topgenesAdjP[topgenesAdjP[, "logFC"] > 0, ]
dim(topUpAdjP)[1]
topDownAdjP = topgenesAdjP[topgenesAdjP[, "logFC"] < 0, ]
dim(topDownAdjP)[1]

# -

# ><mark>**Question 12 :**<mark>  
# >    Pour chacun des seuils de p-value précédents, indiquez le nombre de gènes up et down-régulés

# ><mark>**Question 13 :**<mark>  
# >    Est-ce que l'identification de ces gènes DE nous permet de capturer la totalité des variations des transcriptomes dans ces échantillons ?  
# > Quel autre facteur n'a pas été pris en compte ?  
# > Pensez-vous qu'une analyse entre 2 groupes en factorisant par ce 2ème facteur suffise ?

# ### 5.4 - Analyse ANOVA
# En fait nous avons réellement 4 groupes, et il vaut mieux utiliser un test statistique plus fin pour bien voir les différences entre ces 4 groupes.   
# Une ANOVA à 2 facteurs est le bon choix ici. Cette méthode permet de tester l'effet de chaque facteur ainsi que de leur éventuelle interaction.  
# Nous commençons par redéfinir des groupes, mais ce coup-ci, en établissant 2 facteurs.   
#  - fs : facteur Sexe Gonadique
#  - fd : facteur de stades de développement  
# Chacun des 2 facteurs possèdent 2 niveaux (on pourrait donc en spécifier plus si nécessaire, par exemple pour les stades de développement).  
# On effectue ensuite l'analyse en disant que l'on veut l'effet des 2 facteurs et leur interaction (signe * dans la fonction model.matrix().  
#
#

# +
## Cellule de Code n°29 ##

ph@data[ ,3] = c("Testis","Testis","Testis","Testis","Testis","Testis", "Ovary", "Ovary", "Ovary", "Ovary", "Ovary", "Ovary")
colnames(ph@data)[3]="GonadSex"
ph@data[ ,4] = c("E13","E13","E13","P0","P0","P0", "E13", "E13", "E13", "P0", "P0", "P0")
colnames(ph@data)[4]="DevStage"
ph@data

# +
## Cellule de Code n°30 ##

# On forme les groupes
groupsS = ph@data$GonadSex 
groupsD = ph@data$DevStage
fs = factor(groupsS,levels=c("Testis","Ovary"))
fd = factor(groupsD,levels=c("E13","P0"))

factored.design = model.matrix(~fs*fd)
data.fit4 = lmFit(data.gcrma,factored.design)

# Affichage des résultats :
# data.fit4[1:10,]
# lmfit() will fit a linear model to the data. It will create a data frame called data.fit4 containing 4 columns:
#   the first column contains the intercept of the linear model: it's the mean log expression in the testis samples
#   the second column compares Ovary samples with Testis
#   the third column compares P0 with E13
#   the fourth column contains the interaction showing you if the genes that respond to the gonadal sex are influenced by the developmental stage

# Performing the two-factor ANOVA is now done using eBayes(). 
data.fit.eb4 = eBayes(data.fit4)

# Affichage des p-value pour chacun des facteurs et leur interactions, pour les 10 premiers gènes
data.fit.eb4$p.value[1:10,1:4]
# -

# ### Correction des tests multiples avec une comparaison multiple (i.e une 2-way ANOVA)

# La correction des tests multiples s'effectue ici avec une autre fonction : decideTests()  
# On peut là aussi spécifier des seuils pour ne garder que les gènes à considérer comme DE.  
# La fonction retourne alors des valeurs 0, 1, -1 selon que le gène remplit les critères de sélection.

# +
## Cellule de Code n°31 ##

# The decideTests() method will perform multiple testing adjustment on these p-values. 
# Additionally, it will evaluate for each gene whether the results data.fit.eb fulfill the criteria for differential expression that you specify.
# For each gene and each comparison it will generate the following output:
#   -1: significantly downregulated
#    0: no significant evidence of differential expression
#    1: significantly upregulated
# The p.value argument specifies the FDR and the lfc argument specifies the minimal fold change that is required to be considered DE.

DEresults = decideTests(data.fit.eb4,method='global',adjust.method="BH",p.value=0.05,lfc=1)
DEresults[1:10,]

DEresultsStrict = decideTests(data.fit.eb4,method='global',adjust.method="BH",p.value=0.01,lfc=2)
DEresultsStrict[1:10,]
# -

# ### 5.5- Diagramme de Venn
#
# Pour terminer, on peut afficher des diagrammes de Venn qui montrent comment se répartissent les gènes DE selon les facteurs considérés.

# +
## Cellule de Code n°32 ##

# On choisit une taille de graph plus adaptée (format carré pour ne pas déformer les cercles)
options(repr.plot.width= 12, repr.plot.height= 12)

# Diagramme de Venn des gènes DE
# "both" for all differentially expressed genes, "up" for up-regulated genes only or "down" for down-regulated genes only.

vennDiagram(DEresults[,2:4], include="down", names= c("Gonad Sex", "Dev Stage", "Interaction"), counts.col=c("black"), circle.col = c("red", "blue", "green3"))
vennDiagram(DEresultsStrict[,2:4], include="up", names= c("Gonad Sex", "Dev Stage", "Interaction"), counts.col="red", circle.col = c("red4", "blue4", "green4"))


# -

# ><mark>**Question 14 :**<mark>   
# > <mark>   Reprenez et modifiez le code des cellules précedentes de façon à produire un diagramme de Venn :  <mark>  
# > <mark>  - sur les gènes avec une adj pvalue < 0.01  <mark>  
# > <mark>  - et avec un log fold change supérieur à 1.5      <mark>                                     
# > <mark>  - les gènes up et down-régulés   <mark>  
# > <mark>  - avec l'écriture des comptages en bleu  <mark> 
#   

# +
# rentrez ici le code modifié pour répondre à la question 13


# -

# ><mark>**Question 15 :**<mark>   
# > Combien de gènes sont DE à la fois selon le sexe et le stade de développement ?  
# > Combien de gènes ont leur expression modifiée selon le stade de développement, mais pas différemment selon le sexe ?  
# > Combien de gènes ont leur expression modifiée par l'interaction entre les 2 facteurs ?  

# ###  
# ---
#
# <div class="alert alert-block alert-success"> <b> Bravo ! </b> <br> 
# Pensez à sauver votre notebook, et à en garder une copie en format <b>html</b> <br>
# - Ouvrez "File" dans le Menu <br>
# - Selectionnez "Export Notebook As" <br>
# - Exportez le notebook en HTML <br>
# - Vous pouvez maintenant l'ouvrir dans votre browser  sans même être connecté à adénine ! </div>

StarcoderdataJupyterScriptsDedupFiltered