Préparation des données pour pouvoir faire l’analyse d’expression différentielle

DESeq à besoin d’un object dans un format particulier pour fonctionner. Celui-ci contient les données d’expression ainsi que les informations sur les échantillons nécessaires à la comparaison.
Afin de créer cet objet il faut les informations suivantes:
- Le dossier dans lequel il y a les données d’expression (avec un fichier d’expression par échantillon)
- Le fichier contenant les informations des échantillons (génotype, condition environnementale)

• Nous créons cet objet avec le code suivant:

# Définition du dossier contenant les données d'expression pour tous 
# les échantillons de la comparaison analysée 

directory_0.5mM <- "../../../data/0.5mM_WT_vs_arp6/"
sampleFiles_0.5mM <- list.files(directory_0.5mM)

# Chargement du fichier contenant les information sur les échantillons
sampleCondition_0.5mM <- read_tsv("../../../data/Data_info_WT_arp6_0.5mM.txt")

# Création d'un objet combinant le nom des fichiers des échantillons et 
# les informations sur ces échantillons
sampleTable_0.5mM <- data.frame(sampleName = sampleFiles_0.5mM,
                          fileName = sampleFiles_0.5mM,
                          condition = sampleCondition_0.5mM)
sampleTable_0.5mM$condition <- factor(sampleTable_0.5mM$condition.genotype)

# Vérifiez que tout est bon à cette étape avec la commande:
View(sampleTable_0.5mM)


# Création de l'objet spécifique à DESeq:
ddsHTSeq_0.5mM <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable_0.5mM,
                                       directory = directory_0.5mM,
                                       design= ~ condition)
ddsHTSeq_0.5mM

• Nous devons ensuite enlever les gènes non exprimés ou très faiblement exprimés:

keep <- rowSums(counts(ddsHTSeq_0.5mM)) >= 6
ddsHTSeq_0.5mM <- ddsHTSeq_0.5mM[keep,]

• DESeq doit savoir quelle est la condition de référence, et quelle est la condition que nous comparons à cette référence. Ceci se fait avec la fonction suivante:

ddsHTSeq_0.5mM$condition <- relevel(ddsHTSeq_0.5mM$condition, ref = "WT")

Analyse d’expression différentielle

Maintenant que tout est préparé, nous pouvons enfin effectuer l’analyse d’expression différentielle. Cette analyse intègre une normalisation des données pour corriger pour la taille de la librarie (raison pour laquelle nous n’avons pas besoin de le faire comme à la séance précédente).
Ensuite DESeq va effectuer un test statistique afin de définir si chaque gène est différentiellement exprimé. Ce test est fait en partant du principe que la majorité des gènes ne changent pas d’expression entre les deux conditions. Il va attribuer une p-value ajustée (en utilisant la méthode de benjamini-hochberg) pour tenir compte du fait qu’un nombre important de tests statistiques sont effectués (un test par gène sachant qu’il y a environ 20000 gènes).

• L’analyse d’expression différentielle est effectuée avec la fonction suivante:

dds_0.5mM <- DESeq(ddsHTSeq_0.5mM)

• Maintenant nous pouvons ajouter un seuil pour la p-value ajustée à 5%. Les gènes avec une p-value ajustée inférieure à cette valeur sont différentiellement exprimés.

res_0.5mM <- results(dds_0.5mM, alpha=0.05)
res_0.5mM

## log2 fold change (MLE): condition arp6 vs WT 
## Wald test p-value: condition arp6 vs WT 
## DataFrame with 24942 rows and 6 columns
##             baseMean log2FoldChange     lfcSE         stat    pvalue      padj
##            <numeric>      <numeric> <numeric>    <numeric> <numeric> <numeric>
## AT1G01010   137.1266     -0.6392509  0.267640    -2.388472 0.0169186  0.121908
## AT1G01020   190.2127     -0.1222662  0.226461    -0.539901 0.5892656  0.845239
## AT1G01030    42.0664      0.2243782  0.436388     0.514171 0.6071322  0.853653
## AT1G01040   630.0558     -0.0141657  0.133220    -0.106334 0.9153177  0.978841
## AT1G01050   801.5885     -0.0430895  0.113406    -0.379957 0.7039774  0.900702
## ...              ...            ...       ...          ...       ...       ...
## ATMG00650    2.83986    0.183766174  1.845123  0.099595633  0.920665  0.979897
## ATMG00700    1.07625    1.407237803  3.276264  0.429525108  0.667541        NA
## ATMG01275    8.32380    0.000787451  1.153750  0.000682515  0.999455  0.999789
## ATMG01380    8.99937   -0.044668361  0.819902 -0.054480143  0.956553  0.989149
## ATMG01390 4380.02306    0.329337764  0.335331  0.982128562  0.326037  0.669278

• Afin de connaitre le nombre de gènes différentiellement exprimés pour cette comparaison nous faisons:

sum(res_0.5mM$padj < 0.05, na.rm=TRUE)

## [1] 2347

• Pour visualiser l’analyse d’expression différentielle nous utilisons un graphique MA:

plotMA(res_0.5mM, ylim=c(-2,2))

Ce graphique représente la différence d’expression entre les deux conditions en fonction de l’expression moyenne pour chaque gène. Les gènes différentiellement exprimés sont en bleu.

Que pensez-vous de ce graphique?
Est-ce que la détection de gènes différentiellement exprimés est la même peut importe le niveau d’expression?

• Nous allons maintenant récupérer le résultat de DESeq dans un objet de format data.frame, qui nous permettra de faire des analyses sur ces gènes:

resOrdered_0.5mM <- res_0.5mM[order(res_0.5mM$padj),]
resOrdered_0.5mM

## log2 fold change (MLE): condition arp6 vs WT 
## Wald test p-value: condition arp6 vs WT 
## DataFrame with 24942 rows and 6 columns
##            baseMean log2FoldChange     lfcSE      stat       pvalue
##           <numeric>      <numeric> <numeric> <numeric>    <numeric>
## AT2G25450  3002.159       -3.53804 0.1464012  -24.1667 4.98057e-129
## AT2G04700   895.935       -2.18926 0.1128378  -19.4019  7.44167e-84
## AT2G33220   621.495       -2.91022 0.1515776  -19.1995  3.73624e-82
## AT1G02500  3495.873       -1.59592 0.0842435  -18.9441  4.93711e-80
## AT2G17420  1592.868       -2.21170 0.1168775  -18.9232  7.34475e-80
## ...             ...            ...       ...       ...          ...
## ATCG01200   1.42129        1.74277   2.13309  0.817015     0.413920
## ATCG01240   1.17072        2.59761   2.46154  1.055280     0.291297
## ATCG01310   1.06519       -2.29352   2.52609 -0.907933     0.363914
## ATMG00180   1.71439       -1.79036   2.44552 -0.732100     0.464107
## ATMG00700   1.07625        1.40724   3.27626  0.429525     0.667541
##                   padj
##              <numeric>
## AT2G25450 1.19399e-124
## AT2G04700  8.91996e-80
## AT2G33220  2.98563e-78
## AT1G02500  2.95894e-76
## AT2G17420  3.52152e-76
## ...                ...
## ATCG01200           NA
## ATCG01240           NA
## ATCG01310           NA
## ATMG00180           NA
## ATMG00700           NA

DESeq_result_0.5mM <- as.data.frame(resOrdered_0.5mM) %>% 
  rownames_to_column(var="gene")

• Regardez le tableau avec les données d’analyse d’expression différentielle avec la fonction suivante.

View(DESeq_result_0.5mM)

Est-ce que le gène ARP6 est bien différentiellement exprimé?

• Finalement, exportez le résultat afin de le conserver sans avoir à refaire toute l’analyse:

write_tsv(DESeq_result_0.5mM, file="../DEseq_0.5mM_WT_vs_arp6.txt")