\(^1\) Hub de Bioinformatique et Biostatistique, Centre de Ressources et Recherche en Informatique, Institut Pasteur

\(^2\) Plate-forme Métabolomique, Centre de Ressources et Recherches Technologiques, Institut Pasteur

Ce fichier Rmarkdown contient tout le matériel de l’atelier RNA-Seq de l’école de Bioinformatique Aviesan/IFB/Inserm niveau 2. Diverses notions seront abordées à travers plusieurs types de designs expérimentaux :

design simple avec 3 conditions,
heatmaps,
prise en compte d’un effet “batch”,
conversion d’identifiants de gènes,
analyse de pathways,
design complexe avec 2 facteurs en interaction,

chacune d’entre-elles pouvant nécessiter plusieurs heures d’approfondissement d’un point de vue statistique/méthodologique. Nous les aborderons donc dans les grandes lignes tout en insistant sur certains concepts et détails clés.

Le code ci-dessous permet de charger les packages R dont nous allons avoir besoin :

library(DESeq2)             # analyse differentielle
library(limma)              # removeBatchEffect et CAMERA
library(fgsea)              # gene-set enrichment analysis
library(biomaRt)            # faire la correspondance ensembl id - gene name
library(gage)               # recuperer gene-sets KEGG
library(msigdbr)            # recuperer gene-sets KEGG
library(ggplot2)            # plots
library(plotly)             # plots interactifs
library(ggrepel)            # ajouter des labels sur les plots
library(ggbeeswarm)         # plots avancés
library(scales)             # transformation log2 dans ggplot
library(dplyr)              # manipulation de donnees
library(ggvenn)             # diagrammes de Venn
library(UpSetR)             # UpSet plot
library(FactoMineR)         # ACP
library(factoextra)         # plot de l'ACP
library(pheatmap)           # pretty heatmap
library(RColorBrewer)       # palette pour les heatmaps
library(ExploreModelMatrix) # application shiny interpretation design
library(devtools)           # session info

Les données utilisées proviennent de diverses expériences RNA-Seq et sont disponibles dans des fichiers au format .RData que nous chargerons avec la fonction load(). Chaque fichier .RData contient deux objets :

counts : matrice de comptages avec les gènes en lignes et les échantillons en colonnes (après nettoyage, mapping et comptage des reads sur les gènes),
target : data.frame décrivant le design expérimental.

On suppose que les données sont de bonne qualité et qu’il n’y a pas eu de problème durant l’extraction des ARN, la préparation des librairies et le séquençage. Aussi, il n’y a pas d’outlier ni d’inversion d’échantillon. En résumé, tous les contrôles ont été faits grâce à des statistiques/figures descriptives simples. Enfin, les échantillons sont nommés et ordonnés de la même manière dans le tableau de design et dans les comptages.

Remarque : certaines données ont été anonymisées en modifiant le contexte du design expérimental et en utilisant des identifiants de gènes du type gene1, gene2, …, geneN.

1 Analyse différentielle simple : étude d’un facteur à 3 niveaux

Le jeu de données utilisé ici provient d’une expérience sur un champignon dont l’objectif était de suivre l’évolution transcriptomique au cours du temps. Pour cela, nous avons séquencé l’ARN de 3 échantillons à T=0h, T=4h et à T=8h. On commence par afficher le design de l’expérience ainsi que la matrice de comptages obtenue :

load("data_3_levels.RData")
ls()

## [1] "counts" "target"

print(target)

##      label time replicate
## T0-1  T0-1   T0        r1
## T0-2  T0-2   T0        r2
## T0-3  T0-3   T0        r3
## T4-1  T4-1   T4        r1
## T4-2  T4-2   T4        r2
## T4-3  T4-3   T4        r3
## T8-1  T8-1   T8        r1
## T8-2  T8-2   T8        r2
## T8-3  T8-3   T8        r3

print(head(counts))

##       T0-1 T0-2 T0-3 T4-1 T4-2 T4-3 T8-1 T8-2 T8-3
## gene1 2331 3079 3513 2250 3268 3141 1785 3488 2107
## gene2  772 1677  270 1023 1973  479 1165 1318  700
## gene3 2467 3655 3204 2442 3998 2595 1786 3032 1640
## gene4 4177 4632 5667 3798 5428 5543 3159 4880 3907
## gene5 4844 5501 6222 3731 4857 4522 3356 4906 2578
## gene6 3522 5206 3438 2643 5209 3239 1962 3163 1834

Question : que peut-on faire à partir de ces données et comment pouvons-nous les analyser ?

1.1 Utilisation de DESeq2

Les commandes suivantes permettent de créer un objet dds (qui est central et spécifique à l’analyse différentielle avec DESeq2) puis de réaliser les étapes suivantes :

la normalisation,
l’estimation de la dispersion,
l’ajustement du modèle linéaire généralisé.

dds <- DESeqDataSetFromMatrix(countData=counts,
                              colData=target,
                              design=~time)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in design formula are characters, converting
## to factors

dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)

Remarque : nous utilisons ici DESeq2 avec le paramétrage par défaut, mais il existe de nombreuses subtilités mathématiques/statistiques. La vignette du package est extrêmement bien faite et explique les différentes options possibles.

1.2 Analyse en composantes principales et homoscédasticité

L’analyse en composantes principales (ACP) est une méthode de réduction de la dimension qui permet d’explorer la structure du jeu de données et de le visualiser facilement sur un sous-espace bien choisi. Les figures générées permettent de détecter des échantillons outliers, des inversions, des effets “batch”, etc…

Les deux figures ci-dessous permettent de visualiser le caractère hétéroscédastique des données de comptages, et l’homoscédasticité après transformation rlog. Ce sont ces données transformées qui sont utilisées par les méthodes exploratoires telles que l’ACP car le calcul des distances entre les échantillons ne doit pas être drivé par seulement quelques gènes fortement exprimés.

data.frame(m=apply(counts, 1, mean) + 1, v=apply(counts, 1, var) + 1) %>%
  ggplot(aes(m, v)) +
  scale_x_log10() +
  scale_y_log10() +
  geom_point(alpha=0.2) +
  geom_abline(slope=1, intercept=0, linetype="dashed", col="red") +
  xlab("Mean") +
  ylab("Variance") +
  ggtitle("Comptages bruts : un point = un gène") +
  theme_light()

rld <- rlog(dds)
data.frame(m=apply(assay(rld), 1, mean), v=apply(assay(rld), 1, var)) %>%
  ggplot(aes(m, v)) +
  geom_point(alpha=0.2) +
  xlab("Mean") +
  ylab("Variance") +
  ggtitle("Comptages transformés : un point = un gène") +
  theme_light()

On utilise donc les données rendues homoscédastiques grâce à la méthode rlog() puis la fonction plotPCA() disponible dans DESeq2 :

plotPCA(rld, intgroup="time")

Une alternative est d’utiliser les packages R FactoMineR et factoextra afin d’avoir plus de flexibilité :

pca <- PCA(t(assay(rld)), scale.unit=FALSE, ncp=3, graph=FALSE)
fviz_pca_ind(pca, axes=c(1, 2), repel=TRUE, habillage=factor(target$time), mean.point=FALSE)

fviz_pca_ind(pca, axes=c(1, 3), repel=TRUE, habillage=factor(target$time), mean.point=FALSE)

Question : pourquoi y a-t-il une différence entre les deux approches (plotPCA() de DESeq2 et PCA() de FactoMineR) ?

1.3 Tests des 3 comparaisons

Le code ci-dessous permet d’extraire toutes les comparaisons 2 à 2 :

4h vs 0h,
8h vs 0h,
8h vs 4h,

et d’afficher un MA-plot :

resultsNames(dds)

## [1] "Intercept"     "time_T4_vs_T0" "time_T8_vs_T0"

# interpretation des coefficients
# ExploreModelMatrix() pour lancer l'application Shiny
res_4vs0 <- lfcShrink(dds, contrast=c("time", "T4", "T0"), type="ashr", quiet=TRUE)
res_8vs0 <- lfcShrink(dds, contrast=c("time", "T8", "T0"), type="ashr", quiet=TRUE)
res_8vs4 <- lfcShrink(dds, contrast=c("time", "T8", "T4"), type="ashr", quiet=TRUE)
summary(res_4vs0, alpha=0.05)

## 
## out of 2169 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 348, 16%
## LFC < 0 (down)     : 186, 8.6%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

summary(res_8vs0, alpha=0.05)

## 
## out of 2169 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 378, 17%
## LFC < 0 (down)     : 319, 15%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

summary(res_8vs4, alpha=0.05)

## 
## out of 2169 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 36, 1.7%
## LFC < 0 (down)     : 57, 2.6%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

# MA-plot de la comparaison 4h vs 0h
DESeq2::plotMA(res_4vs0, alpha=0.05)

Ce MA-plot est produit grâce à une fonction du package DESeq2. Il est possible de le reproduire manuellement et en ajoutant une couche d’interactivité :

p <- res_4vs0 %>%
  as.data.frame() %>%
  mutate(DE = ifelse(!is.na(padj) & padj <= 0.05, "DE", "Not DE")) %>%
  mutate(id=rownames(res_4vs0)) %>%
  ggplot(aes(x=baseMean, y=log2FoldChange, color=DE, label=id)) +
  geom_point(alpha=0.4) +
  labs(color="4h vs 0h") +
  geom_hline(yintercept = 0, linetype="dashed", color="grey") +
  scale_color_manual(values=c("DE"="red", "Not DE"="black")) +
  scale_x_log10() +
  theme_light()
toWebGL(ggplotly(p, tooltip="label"))

On peut aussi ajouter les noms de certains gènes bien choisis :

res_4vs0 %>%
  as.data.frame() %>%
  mutate(DE = ifelse(!is.na(padj) & padj <= 0.05, "DE", "Not DE")) %>%
  mutate(id=ifelse(log2FoldChange > 3, rownames(res_4vs0), "")) %>%
  ggplot(aes(x=baseMean, y=log2FoldChange, color=DE, label=id)) +
  geom_point(alpha=0.4) +
  labs(color="4h vs 0h") +
  geom_text_repel(show.legend=FALSE, size=3, point.padding=0.2, max.overlaps=20) +
  geom_hline(yintercept = 0, linetype="dashed", color="grey") +
  scale_color_manual(values=c("DE"="red", "Not DE"="black")) +
  scale_x_log10() +
  theme_light()

1.4 Croisement des résultats

Etant donné le design expérimental, nous avons obtenus les résultats d’une analyse différentielle pour 3 comparaisons. L’objectif est ici de croiser ces résultats sous la meilleure forme possible, mais nous devons pour cela nous assurer que les gènes sont ordonnés de la même manière dans les 3 objets res_4vs0, res_8vs0 et res_8vs4 :

all(rownames(res_4vs0) == rownames(res_8vs0)) & all(rownames(res_4vs0) == rownames(res_8vs4))

## [1] TRUE

Une manière assez naturelle est de croiser les listes de gènes différentiellement exprimés, en fixant par exemple un seuil de 5% sur la p-valeur ajustée. Ainsi, nous pouvons détecter des gènes dérégulés spécifiquement à 8h mais qui ne l’étaient pas à 4h. Les diagrammes de Venn sont souvent utilisés pour faire cela :

d <- data.frame(gene = rownames(res_4vs0),
                log2FC_4vs0 = res_4vs0$log2FoldChange, padj_4vs0 = res_4vs0$padj,
                log2FC_8vs0 = res_8vs0$log2FoldChange, padj_8vs0 = res_8vs0$padj,
                log2FC_8vs4 = res_8vs4$log2FoldChange, padj_8vs4 = res_8vs4$padj)
print(head(d), digits=3)

##    gene log2FC_4vs0 padj_4vs0 log2FC_8vs0 padj_8vs0 log2FC_8vs4 padj_8vs4
## 1 gene1     0.01005    0.9649      0.0438    0.8398     0.01630     0.952
## 2 gene2     0.07783    0.7625      0.1605    0.5879     0.01461     0.951
## 3 gene3    -0.00112    0.9950     -0.1656    0.3724    -0.09092     0.621
## 4 gene4     0.05119    0.8084      0.0681    0.7501     0.00289     0.998
## 5 gene5    -0.22933    0.0608     -0.2644    0.0366    -0.01051     0.972
## 6 gene6    -0.06300    0.7622     -0.3609    0.0805    -0.14105     0.411

x <- list("4h vs 0h"=d %>% filter(padj_4vs0 <= 0.05) %>% pull(gene),
          "8h vs 0h"=d %>% filter(padj_8vs0 <= 0.05) %>% pull(gene),
          "8h vs 4h"=d %>% filter(padj_8vs4 <= 0.05) %>% pull(gene))
ggvenn(x)

Le diagramme généré est plutôt lisible car nous avons croisé seulement 3 listes. Pour croiser davantage de données, les UpSet plots ont récemment vu le jour et sont souvent plus adaptés :

upset(fromList(x), order.by="freq")

Question : comment interpréter les chiffres issus du diagramme de Venn et de l’UpSet plot ?

Les UpSet plots sont souvent plus lisibles que les diagrammes de Venn mais ont toujours l’inconvénient d’avoir des “effets de seuil”. L’exemple ci-dessous permet de comprendre ce qu’il peut se passer :

d %>% 
  filter(abs(log2FC_8vs0 - log2FC_4vs0) < 0.1) %>%
  filter(padj_4vs0 > 0.05 & padj_8vs0 < 0.05) %>%
  dplyr::select(contains(c("gene", "log2FC", "padj"))) %>%
  head()

##      gene log2FC_4vs0 log2FC_8vs0 log2FC_8vs4  padj_4vs0  padj_8vs0 padj_8vs4
## 1   gene5  -0.2293309  -0.2644350 -0.01050506 0.06083904 0.03657322 0.9722079
## 2 gene374  -0.2786602  -0.3524606 -0.01992609 0.06849556 0.03167935 0.9484213
## 3 gene467  -0.2538975  -0.3431442 -0.03324774 0.07968930 0.02544719 0.8991029
## 4 gene473  -0.2855411  -0.3302746 -0.00588014 0.05550870 0.03794243 0.9912358
## 5 gene523  -0.2069653  -0.3004383 -0.04890325 0.10319406 0.02013790 0.8363462
## 6 gene716  -0.2702907  -0.3253369 -0.01213579 0.06908139 0.04138005 0.9705380

Par exemple, au seuil de 5%, le gène 5 est différentiellement exprimé à 8h mais pas à 4h. Pourtant, les log2(Fold-Change) à 4h et à 8h sont très proches (\(-0,23\) et \(-0,26\)), tout comme les P-valeurs ajustées (\(0,06\) et \(0,04\)). On constate d’ailleurs que la différence entre 4h et 8h n’est pas significative (P-valeur ajustée = \(0,97\)). Nous allons donc construire une figure qui permet de visualiser au mieux tous les résultats simultanément :

d <- d %>%
  mutate(DE_8vs4=ifelse(!is.na(padj_8vs4) & padj_8vs4 <= 0.05, "DE", "Not DE"),
         DE_both_times=ifelse(!is.na(padj_4vs0) & padj_4vs0 <= 0.05 & !is.na(padj_8vs0) & padj_8vs0 <= 0.05, "Yes", "No"))
ggplot(d, aes(x=log2FC_4vs0, y=log2FC_8vs0, color=DE_8vs4, size=DE_both_times)) +
  geom_point(alpha=0.5) +
  scale_color_manual(values=c("DE"="red", "Not DE"="black")) +
  scale_size_manual(values=c("Yes"=1, "No"=0.4)) +
  geom_hline(yintercept=0, lty=2, color="grey") +
  geom_vline(xintercept=0, lty=2, color="grey") +
  geom_abline(slope=1, intercept=0, lty=2, color="grey") +
  labs(color="8h vs 4h", size="DE at 4h and 8h") +
  xlab(bquote(log[2](FC)~"4h vs 0h")) +
  ylab(bquote(log[2](FC)~"8h vs 0h")) +
  ggtitle("Comparison of the 3 conditions") +
  theme_light()

Cette figure permet de repérer certains gènes pour lesquels l’effet à 4h est différent de l’effet à 8h, chacun de ces deux effets pouvant être significatifs ou non. On peut alors construire des filtres en fonction de la question biologique, par exemple :

d %>%
  filter(padj_8vs4 <= 0.05 & abs(log2FC_8vs4) > 2) %>%
  dplyr::select(contains(c("gene", "log2FC", "padj")))

##      gene log2FC_4vs0 log2FC_8vs0 log2FC_8vs4    padj_4vs0    padj_8vs0    padj_8vs4
## 1  gene27   3.9889546    7.994236    3.674804 9.962520e-21 3.678433e-76 3.074497e-16
## 2 gene322   0.4810958  -11.417758  -12.140662 1.529902e-02 1.448189e-61 5.244541e-68
## 3 gene898   1.6960953    5.228053    3.407727 1.314652e-09 2.950421e-83 1.085497e-35
## 4 gene970   4.2622347    8.106160    3.574516 2.723993e-25 6.628782e-86 9.362211e-17

1.5 Heatmaps

Les heatmaps sont également souvent utilisées pour visualiser les résultats. Pour autant, celles-ci sont plus complexes qu’il n’y parait et il existe de nombreux paramètres à définir :

quels gènes inclure ?
afficher leurs identifiants ou noms ?
clusteriser les lignes ? colonnes ? selon quelle méthode ?
centrer/réduire les lignes ? colonnes ?
quelle échelle de couleur choisir ?
prendre en compte un éventuel effet réplicat/batch ?

Le code R ci-dessous permet de sélectionner les gènes différentiellement exprimés dans les 3 comparaisons (avec un seuil de 0.001 sur la P-valeur ajustée) et d’afficher une heatmap (i) sans paramétrage spécifique et (ii) en spécifiant certains paramètres.

DE_genes <- d %>%
  filter(padj_4vs0 <= 0.001 | padj_8vs0 <= 0.001 | padj_8vs4 <= 0.001) %>%
  pull(gene)
mat_hm <- assay(rlog(dds))
pheatmap(mat_hm[DE_genes, ])

Question : quel message renvoie cette heatmap ? Pourquoi ?

pheatmap(mat_hm[DE_genes, ],
         scale="row",
         show_rownames=FALSE,
         clustering_distance_rows="euclidean",
         clustering_distance_cols="euclidean",
         clustering_method="ward.D",
         color=colorRampPalette(rev(brewer.pal(9, "RdBu")))(200))

Questions :

pourquoi cette nouvelle heatmap est-elle si différente de la première ? En quoi est-elle plus adaptée ?
le clustering est-il réalisé à partir des données brutes ou centrées-réduites ?

On peut également visualiser les log2(Fold-Change) plutôt que les données d’expression. Le code R ci-dessous permet de faire cela pour un certain nombre de gènes pré-sélectionnés :

rownames(d) <- d$gene
d <- d %>%
  filter(padj_8vs4 <= 0.05 & (padj_4vs0 <= 0.05 | padj_8vs0 <= 0.05)) %>% 
  filter(abs(log2FC_8vs0 - log2FC_4vs0) > 1) %>%
  dplyr::select(log2FC_4vs0, log2FC_8vs0) %>%
  as.matrix()
print(head(d))

##         log2FC_4vs0  log2FC_8vs0
## gene27    3.9889546  7.994236337
## gene135   3.9349746  6.211718231
## gene213   1.4518787  0.215970665
## gene220   1.1212949 -0.009835288
## gene237   1.7449548  0.733371069
## gene287   0.1331511 -1.122335229

pheatmap(d, 
         scale="none", 
         cluster_cols=FALSE,
         clustering_distance_rows="euclidean",
         angle_col=45,
         fontsize_row=5,
         color=colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         breaks=seq(-max(abs(d)), max(abs(d)), length=101))

La heatmap produite est influencée par quelques gènes ayant des log2(Fold-Change) très élevés ou très faibles. Dans ce cas, on peut borner l’échelle de couleur de manière à faire ressortir les autres gènes plus facilement. Toutefois, il faut faire attention à ne pas déformer la figure :

pheatmap(d, 
         scale="none", 
         cluster_cols=FALSE,
         clustering_distance_rows="euclidean",
         angle_col=45,
         fontsize_row=5,
         color=colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         breaks=seq(-3, 3, length=101))

2 Prise en compte d’un effet batch

Nous utilisons ici un jeu de données RNA-Seq humain : un traitement expérimental T et un traitement contrôle C ont été appliqués sur des cultures cellulaires provenant de 3 donneurs indépendants, l’objectif étant de détecter les gènes dérégulés par l’effet du traitement.

rm(list=ls())
load("data_batch_effect.RData")
ls()

## [1] "counts" "target"

print(target)

##       label condition donor
## C-d03 C-d03         C   d03
## C-d08 C-d08         C   d08
## C-d12 C-d12         C   d12
## T-d03 T-d03         T   d03
## T-d08 T-d08         T   d08
## T-d12 T-d12         T   d12

print(head(counts))

##                 C-d03 C-d08 C-d12 T-d03 T-d08 T-d12
## ENSG00000000003    64    55    37    62    45    50
## ENSG00000000005     0     0     0     0     0     0
## ENSG00000000419  8370  5420  6154  7823  5283  5849
## ENSG00000000457   970   811   567   950   669   545
## ENSG00000000460  1689  1113   981  2065  1082  1264
## ENSG00000000938    44    25    10    24     5     1

2.1 Analyse en composantes principales

Commençons par créer l’objet dds à partir des comptages et du tableau de design :

dds <- DESeqDataSetFromMatrix(countData=counts,
                              colData=target,
                              design=~condition)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in design formula are characters, converting
## to factors

Question : que révèle l’analyse en composante principale ?

rld <- rlog(dds)
pca <- PCA(t(assay(rld)), scale.unit=FALSE, ncp=3, graph=FALSE)
fviz_pca_ind(pca, axes=c(1, 2), repel=TRUE, habillage=factor(target$condition), mean.point=FALSE)

fviz_pca_ind(pca, axes=c(1, 3), repel=TRUE, habillage=factor(target$condition), mean.point=FALSE)

2.2 Importance du design dans DESeq2

Le code ci-dessous permet de réaliser l’analyse différentielle et visualiser la distribution des p-valeurs brutes :

dds <- DESeq(dds)
resultsNames(dds)

## [1] "Intercept"        "condition_T_vs_C"

res <- lfcShrink(dds, coef="condition_T_vs_C", type="ashr", quiet=TRUE)
summary(res, alpha=0.05)

## 
## out of 32453 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 17, 0.052%
## LFC < 0 (down)     : 69, 0.21%
## outliers [1]       : 106, 0.33%
## low counts [2]     : 10723, 33%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

ggplot(as.data.frame(res), aes(x=pvalue)) + 
  geom_histogram(color = "white", breaks = seq(0, 1, by = 0.025)) +
  xlab("Raw p-value") +
  ylab("Frequency") + 
  theme_light()

## Warning: Removed 32870 rows containing non-finite values (`stat_bin()`).

La distribution n’est pas celle attendue car la variabilité (information) générée par l’effet “donneur” n’est pas prise en compte dans le modèle. Nous pouvons visualiser ci-dessous cette variabilité pour un seul gène (bien choisi) :

ncounts <- counts(dds, norm=TRUE)
d <- data.frame(target, ncounts=ncounts["ENSG00000066279",])
print(d)

##       label condition donor   ncounts
## C-d03 C-d03         C   d03  575.2294
## C-d08 C-d08         C   d08 3654.0259
## C-d12 C-d12         C   d12  378.6506
## T-d03 T-d03         T   d03 1095.8348
## T-d08 T-d08         T   d08 4831.8431
## T-d12 T-d12         T   d12  598.5342

ggplot(d, aes(x=condition, y=ncounts, color=donor, shape=donor)) +
  geom_point(size=4) +
  scale_y_continuous(trans = log2_trans(),
    breaks = trans_breaks("log2", function(x) 2^x),
    labels = trans_format("log2", math_format(2^.x))) +
  xlab("Condition") +
  ylab("Normalized counts (log2 scale)") +
  labs(shape="Donor", color="Donor") +
  ggtitle("Illustration of the donor effect for gene ENSG00000066279") +
  theme_light()

Il est nécessaire de spécifier dans le design du modèle cet effet “donneur” :

dds <- DESeqDataSetFromMatrix(countData=counts,
                              colData=target,
                              design=~donor+condition)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in design formula are characters, converting
## to factors

dds <- DESeq(dds)
resultsNames(dds)

## [1] "Intercept"        "donor_d08_vs_d03" "donor_d12_vs_d03" "condition_T_vs_C"

res <- lfcShrink(dds, coef="condition_T_vs_C", type="ashr", quiet=TRUE)
summary(res, alpha=0.05)

## 
## out of 32453 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 577, 1.8%
## LFC < 0 (down)     : 729, 2.2%
## outliers [1]       : 0, 0%
## low counts [2]     : 13702, 42%
## (mean count < 10)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

ggplot(as.data.frame(res), aes(x=pvalue)) + 
  geom_histogram(color = "white", breaks = seq(0, 1, by = 0.025)) +
  xlab("Raw p-value") +
  ylab("Frequency") + 
  theme_light()

## Warning: Removed 32764 rows containing non-finite values (`stat_bin()`).

Question : quelles sont les améliorations apportées par la prise en compte de cet effet “donneur” ?

2.3 Elimination de l’effet batch pour la visualisation

Maintenant que l’analyse différentielle est correctement réalisée, nous pouvons visualiser les résultats sous la forme d’une heatmap, en sélectionnant les gènes différentiellement exprimés au seuil de 5% :

DE_genes <- res %>%
  as.data.frame() %>%
  filter(padj <= 0.05) %>%
  rownames()
mat_hm <- assay(rlog(dds))
pheatmap(mat_hm[DE_genes,],
         scale="row",
         show_rownames=FALSE,
         clustering_distance_rows="euclidean",
         clustering_distance_cols="euclidean",
         clustering_method="ward.D",
         color=colorRampPalette(rev(brewer.pal(9, "RdBu")))(200))

Dans cette heatmap, les échantillons clusterisent par donneur et pas par traitement alors que nous avons pourtant sélectionné les gènes selon ce critère. Une solution est d’ajuster la matrice avant de produire la figure grâce à la fonction removeBatchEffect() :

mat_hm_adj <- removeBatchEffect(mat_hm, batch=target$donor)
pheatmap(mat_hm_adj[DE_genes,], 
         scale="row",
         show_rownames=FALSE,
         clustering_distance_rows="euclidean",
         clustering_distance_cols="euclidean",
         clustering_method="ward.D",
         color=colorRampPalette(rev(brewer.pal(9, "RdBu")))(200),
         main="Heatmap based on the *adjusted* variance-stabilized counts")

Remarques :

Il ne faut pas confondre l’ajustement dans le design du modèle linéaire généralisé de DESeq2 (les données restent inchangées) et l’ajustement directement sur la matrice d’expression (les valeurs sont alors modifiées) dans une optique de visualisation.
Si la heatmap contient trop de cellules (en lignes et/ou colonnes) alors il se peut que l’affichage et le rendu soient contraints par le nombre de pixels de l’écran. Ces deux articles décrivent ce “problème” : ici et là.
Les heatmaps sont de formidables outils de visualisation et on peut imaginer une personnalisation à l’infini, notamment avec le package ComplexHeatmap.

2.4 Identifiants des gènes

Les identifiants utilisés pour l’analyse différentielle sont issus d’ensembl (et du fichier d’annotation au format GFF ou GTF) et ne sont pas toujours pratiques pour l’interprétation. Le code ci-dessous est une manière parmi tant d’autres de récupérer les noms des gènes correspondants :

ensembl <- useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl", version=92)
# attr <- listAttributes(ensembl)$name
# attr[grepl("ensembl_gene_id", attr)]
# attr[grepl("entrezgene", attr)]
# attr[grepl("external_gene_name", attr)]
bm <- getBM(attributes=c("ensembl_gene_id", "external_gene_name"), mart=ensembl)
head(bm)

##   ensembl_gene_id external_gene_name
## 1 ENSG00000276626            RF00100
## 2 ENSG00000201317           RNU4-59P
## 3 ENSG00000200823         SNORD114-2
## 4 ENSG00000221598            MIR1249
## 5 ENSG00000199595            RF00019
## 6 ENSG00000200839           RNA5SP48

Le tableau bm permet de faire le lien entre les identifiants ensembl sous la forme ENSG000... et les noms humainement interprétables. Ici, il est nécessaire d’explorer ce tableau afin de vérifier la présence de doublons dans chacune des colonnes. Commençons par les gene names :

bm %>%
  filter(!is.na(external_gene_name)) %>%
  pull(external_gene_name) %>%
  duplicated() %>%
  any()

## [1] TRUE

dup_gene_names <- bm %>%
  filter(!is.na(external_gene_name)) %>%
  filter(duplicated(external_gene_name)) %>%
  pull(external_gene_name)
bm %>%
  filter(external_gene_name %in% dup_gene_names) %>%
  arrange(external_gene_name) %>%
  head(n=10)

##    ensembl_gene_id external_gene_name
## 1  ENSG00000197953            AADACL2
## 2  ENSG00000261846            AADACL2
## 3  ENSG00000242908        AADACL2-AS1
## 4  ENSG00000262466        AADACL2-AS1
## 5  ENSG00000276072               AATF
## 6  ENSG00000275700               AATF
## 7  ENSG00000274099           ABCB10P1
## 8  ENSG00000280461           ABCB10P1
## 9  ENSG00000281173           ABCB10P1
## 10 ENSG00000261524           ABCB10P3

Plusieurs identifiants ensembl correspondent à un même gene name ! Regardons ensuite les identifiants ensembl :

bm %>%
  filter(!is.na(ensembl_gene_id)) %>%
  pull(ensembl_gene_id) %>%
  duplicated() %>%
  any()

## [1] FALSE

Les identifiants ensembl sont uniques et nous avons un seul gene name par identifiant ensembl. Nous allons donc pouvoir les remplacer facilement :

DE_genes <- res %>%
  as.data.frame() %>%
  filter(padj <= 0.0001 & (log2FoldChange > 2.5 | log2FoldChange < -3.5)) %>%
  rownames() %>%
  sort()
DE_gene_names <- bm %>%
  arrange(ensembl_gene_id) %>%
  filter(ensembl_gene_id %in% DE_genes) %>%
  pull(external_gene_name)
print(cbind(DE_genes, DE_gene_names))

##      DE_genes          DE_gene_names
## [1,] "ENSG00000111537" "IFNG"       
## [2,] "ENSG00000112116" "IL17F"      
## [3,] "ENSG00000115008" "IL1A"       
## [4,] "ENSG00000115523" "GNLY"       
## [5,] "ENSG00000124191" "TOX2"       
## [6,] "ENSG00000124766" "SOX4"       
## [7,] "ENSG00000136634" "IL10"       
## [8,] "ENSG00000169860" "P2RY1"

pheatmap(mat_hm_adj[DE_genes,], 
         scale="row",
         labels_row=DE_gene_names,
         clustering_distance_rows="euclidean",
         clustering_distance_cols="euclidean",
         clustering_method="ward.D",
         color=colorRampPalette(rev(brewer.pal(9, "RdBu")))(200),
         main="Heatmap based on the *adjusted* variance-stabilized counts")

Remarque : nous avons utilisé le package R biomaRt pour associer les identifiants ensembl aux gene names, mais les packages AnnotationHub et AnnotationDbi (entre autres) peuvent être intéressants. Ce document recense différents outils disponibles et leurs fonctionnements.

3 Du gène au gene-set

L’analyse conduite ci-dessus nous a permis de détecter un certain nombre de gènes dérégulés par l’effet du traitement. Nous pouvons alors nous demander si ces gènes appartiennent préférentiellement à des pathways/voies/catégories (i.e. ensembles de gènes cohérents biologiquement) connus dans la littérature. Autrement dit, nous allons rechercher des gene-sets enrichis en gènes différentiellement exprimés, i.e. des gene-sets perturbés par l’effet du traitement.

De nombreuses méthodes et plusieurs packages R et applications web existent pour réaliser cette analyse fonctionnelle de manière automatisée. Ici, nous allons la construire pas-à-pas en utilisant trois méthodes différentes et en illustrant les points importants.

La première étape est d’importer dans R les gene-sets que l’on souhaite tester : par exemple les pathways KEGG, les GO terms de la catégorie Molecular Function, etc… Ici, nous importons uniquement les pathways KEGG en utilisant le package R gage :

kegg <- kegg.gsets(species="human", id.type="kegg")$kg.sets
print(head(kegg, 3))

## $`hsa00010 Glycolysis / Gluconeogenesis`
##  [1] "10327"  "124"    "125"    "126"    "127"    "128"    "130"    "130589" "131"    "160287" "1737"   "1738"  
## [13] "2023"   "2026"   "2027"   "217"    "218"    "219"    "2203"   "221"    "222"    "223"    "224"    "226"   
## [25] "229"    "230"    "2538"   "2597"   "26330"  "2645"   "2821"   "3098"   "3099"   "3101"   "387712" "3939"  
## [37] "3945"   "3948"   "441531" "501"    "5105"   "5106"   "5160"   "5161"   "5162"   "5211"   "5213"   "5214"  
## [49] "5223"   "5224"   "5230"   "5232"   "5236"   "5313"   "5315"   "55276"  "55902"  "57818"  "669"    "7167"  
## [61] "80201"  "83440"  "84532"  "8789"   "92483"  "92579"  "9562"  
## 
## $`hsa00020 Citrate cycle (TCA cycle)`
##  [1] "1431"  "1737"  "1738"  "1743"  "2271"  "3417"  "3418"  "3419"  "3420"  "3421"  "4190"  "4191"  "47"    "48"   
## [15] "4967"  "50"    "5091"  "5105"  "5106"  "5160"  "5161"  "5162"  "55753" "6389"  "6390"  "6391"  "6392"  "8801" 
## [29] "8802"  "8803" 
## 
## $`hsa00030 Pentose phosphate pathway`
##  [1] "132158" "2203"   "221823" "226"    "229"    "22934"  "230"    "2539"   "25796"  "2821"   "414328" "51071" 
## [13] "5211"   "5213"   "5214"   "5226"   "5236"   "55276"  "5631"   "5634"   "6120"   "64080"  "6888"   "7086"  
## [25] "729020" "8277"   "84076"  "8789"   "9104"   "9563"

length(kegg)

## [1] 347

Les pathways KEGG sont importés sous la forme d’une liste : chaque élément est un vecteur (ensemble) de gènes représentant un pathway. Aussi, les gènes sont caractérisés par des numéros qui correspondent à leur identifiant Entrez (NCBI). On peut vérifier que chaque gene-set contient des identifiants uniques (i.e. pas de doublons) et qu’aucun gene-set n’est présent deux fois :

any(sapply(kegg, function(x) any(duplicated(x))))

## [1] FALSE

any(duplicated(lapply(kegg, function(x) sort(x))))

## [1] FALSE

L’étape suivante consiste à transformer les identifiants ensembl utilisés pour l’analyse différentielle en identifiants de type Entrez. Pour cela, nous pouvons utiliser la même méthode que précédemment et vérifier la présence de doublons :

bm <- getBM(attributes=c("ensembl_gene_id", "entrezgene"), mart=ensembl)
head(bm)

##   ensembl_gene_id entrezgene
## 1 ENSG00000210049         NA
## 2 ENSG00000211459       4549
## 3 ENSG00000210077         NA
## 4 ENSG00000210082       4550
## 5 ENSG00000209082         NA
## 6 ENSG00000198888       4535

bm %>%
  filter(!is.na(ensembl_gene_id)) %>%
  pull(ensembl_gene_id) %>%
  duplicated() %>%
  any()

## [1] TRUE

dup_ensembl <- bm %>%
  filter(!is.na(ensembl_gene_id)) %>%
  filter(duplicated(ensembl_gene_id)) %>%
  pull(ensembl_gene_id)
bm %>%
  filter(ensembl_gene_id %in% dup_ensembl) %>%
  arrange(ensembl_gene_id) %>%
  head(n=10)

##    ensembl_gene_id entrezgene
## 1  ENSG00000011454       2844
## 2  ENSG00000011454      23637
## 3  ENSG00000023171      57476
## 4  ENSG00000023171  100128242
## 5  ENSG00000058600      55718
## 6  ENSG00000058600  101060521
## 7  ENSG00000072195      10290
## 8  ENSG00000072195  100996693
## 9  ENSG00000079974      11158
## 10 ENSG00000079974      11159

On remarque qu’il est possible d’avoir différents identifiants Entrez qui pointent vers un même identifiant ensembl…

bm %>%
  filter(!is.na(entrezgene)) %>%
  pull(entrezgene) %>%
  duplicated() %>%
  any()

## [1] TRUE

dup_entrez <- bm %>%
  filter(!is.na(entrezgene)) %>%
  filter(duplicated(entrezgene)) %>%
  pull(entrezgene)
bm %>%
  filter(entrezgene %in% dup_entrez) %>%
  arrange(entrezgene) %>%
  head(n=10)

##    ensembl_gene_id entrezgene
## 1  ENSG00000236149         23
## 2  ENSG00000204574         23
## 3  ENSG00000232169         23
## 4  ENSG00000231129         23
## 5  ENSG00000236342         23
## 6  ENSG00000206490         23
## 7  ENSG00000225989         23
## 8  ENSG00000175164         28
## 9  ENSG00000281879         28
## 10 ENSG00000159842         29

On remarque que plusieurs identifiants ensembl (et donc plusieurs log2(Fold-Change), P-valeurs…) correspondent à un même identifiant Entrez. Les pathways étant définis par des identifiants Entrez, nous devons choisir quel identifiant ensembl associer à chaque Entrez. Ici, nous choisons l’identifiant ensembl ayant l’expression la plus forte, mais nous pourrions choisir tout autre critère :

res <- res %>%
  as.data.frame() %>%
  mutate(ensembl_gene_id=rownames(res)) %>%
  inner_join(bm, by="ensembl_gene_id", multiple="all") %>%
  arrange(desc(baseMean)) %>%
  filter(!duplicated(entrezgene))
print(head(res))

##   baseMean log2FoldChange      lfcSE      pvalue       padj ensembl_gene_id entrezgene
## 1 615460.0     0.16162487 0.09395108 0.009413525 0.10519249 ENSG00000075624         60
## 2 598624.2     0.05468944 0.07862636 0.319484933 0.75184011 ENSG00000100453       3002
## 3 428040.3    -0.09076801 0.08600251 0.113641579 0.48043456 ENSG00000166710        567
## 4 409505.7     0.07194170 0.07626481 0.189826115 0.60658308 ENSG00000074800       2023
## 5 407566.1     0.18814987 0.09673261 0.003371159 0.04877515 ENSG00000184009         71
## 6 400827.1     0.03958300 0.07516372 0.463933244 0.83995735 ENSG00000111640       2597

3.1 ORA : Over-Representation Analysis

Nous pouvons maintenant tester si le gene-set hsa05418 Fluid shear stress and atherosclerosis (par exemple) est enrichi en gènes différentiellement exprimés. Nous utilisons la méthode ORA (Over-Representation Analysis) avec un test de Fisher :

res <- res %>%
  mutate(geneset = ifelse(entrezgene %in% kegg$`hsa05418 Fluid shear stress and atherosclerosis`, "in", "out")) %>%
  mutate(DE = ifelse(padj <= 0.05 & !is.na(padj), "DE", "notDE"))
print(head(res), digits=3)

##   baseMean log2FoldChange  lfcSE  pvalue   padj ensembl_gene_id entrezgene geneset    DE
## 1   615460         0.1616 0.0940 0.00941 0.1052 ENSG00000075624         60      in notDE
## 2   598624         0.0547 0.0786 0.31948 0.7518 ENSG00000100453       3002     out notDE
## 3   428040        -0.0908 0.0860 0.11364 0.4804 ENSG00000166710        567     out notDE
## 4   409506         0.0719 0.0763 0.18983 0.6066 ENSG00000074800       2023     out notDE
## 5   407566         0.1881 0.0967 0.00337 0.0488 ENSG00000184009         71      in    DE
## 6   400827         0.0396 0.0752 0.46393 0.8400 ENSG00000111640       2597     out notDE

tab <- table(res$geneset, res$DE)
print(tab)

##      
##          DE notDE
##   in     25   113
##   out  1152 18047

print(fisher.test(tab), alternative="greater")

## 
##  Fisher's Exact Test for Count Data
## 
## data:  tab
## p-value = 8.253e-07
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  2.143326 5.408058
## sample estimates:
## odds ratio 
##   3.465486

Le gene-set hsa05418 Fluid shear stress and atherosclerosis est enrichi en gènes dérégulés car son Odds-Ratio défini comme :

\[ \text{OR} = \frac{\text{Odd in gene-set}}{\text{Odd out gene-set}} = \frac{25/113}{1152/18047} = 3.46 \] est significativement supérieur à 1 (P-valeur = \(8.25 \times 10^{-7}\)). On peut également calculer le Fold-Enrichment :

\[ \text{FE} = \frac{\text{Pct gene DE in gene-set}}{\text{Pct gene DE out gene-set }} = \frac{25/(25+113)}{1152/(1152+18047)} = \frac{0,181}{0,6} = 3,02\] On peut maintenant répéter l’opération grâce à une boucle sur tous les pathways KEGG :

res_kegg <- data.frame(kegg=names(kegg), ngenes=sapply(kegg, length), OR=NA, pvalue=NA)
for (name in names(kegg)){
  res <- res %>%
    mutate(geneset = ifelse(entrezgene %in% kegg[[name]], "in", "out")) %>%
    mutate(DE = ifelse(padj <= 0.05 & !is.na(padj), "DE", "notDE"))
  tab <- table(res$geneset, res$DE)
  if (nrow(tab) < 2 | ncol(tab) < 2) stop("tab is not a 2x2 table, can't perform the statistical test")
  f <- fisher.test(table(res$geneset, res$DE), alternative="greater")
  res_kegg[which(res_kegg$kegg == name), c("OR", "pvalue")] <- c(f$estimate, f$p.value)
}
print(head(res_kegg))

##                                                                                                kegg ngenes        OR
## hsa00010 Glycolysis / Gluconeogenesis                         hsa00010 Glycolysis / Gluconeogenesis     67 0.9795541
## hsa00020 Citrate cycle (TCA cycle)                               hsa00020 Citrate cycle (TCA cycle)     30 0.0000000
## hsa00030 Pentose phosphate pathway                               hsa00030 Pentose phosphate pathway     30 0.0000000
## hsa00040 Pentose and glucuronate interconversions hsa00040 Pentose and glucuronate interconversions     35 0.0000000
## hsa00051 Fructose and mannose metabolism                   hsa00051 Fructose and mannose metabolism     33 0.0000000
## hsa00052 Galactose metabolism                                         hsa00052 Galactose metabolism     31 0.5138860
##                                                      pvalue
## hsa00010 Glycolysis / Gluconeogenesis             0.5884936
## hsa00020 Citrate cycle (TCA cycle)                1.0000000
## hsa00030 Pentose phosphate pathway                1.0000000
## hsa00040 Pentose and glucuronate interconversions 1.0000000
## hsa00051 Fructose and mannose metabolism          1.0000000
## hsa00052 Galactose metabolism                     0.8574880

Etant donné le nombre important de pathways testés, il est probable que certaines p-valeurs soient significatives seulement par chance : on doit les ajuster de manière à prendre en compte la multiplicité des tests et ainsi contrôler le taux de faux positifs.

res_kegg <- res_kegg %>%
  mutate(FDR = p.adjust(pvalue)) %>%
  arrange(FDR)
print(head(res_kegg))

##                                                                                            kegg ngenes       OR
## hsa04060 Cytokine-cytokine receptor interaction hsa04060 Cytokine-cytokine receptor interaction    295 4.864888
## hsa05321 Inflammatory bowel disease                         hsa05321 Inflammatory bowel disease     65 9.196942
## hsa05323 Rheumatoid arthritis                                     hsa05323 Rheumatoid arthritis     93 6.782570
## hsa04659 Th17 cell differentiation                           hsa04659 Th17 cell differentiation    108 5.780550
## hsa04630 JAK-STAT signaling pathway                         hsa04630 JAK-STAT signaling pathway    166 4.217496
## hsa04062 Chemokine signaling pathway                       hsa04062 Chemokine signaling pathway    192 3.900735
##                                                       pvalue          FDR
## hsa04060 Cytokine-cytokine receptor interaction 3.742439e-22 1.298626e-19
## hsa05321 Inflammatory bowel disease             1.898498e-13 6.568804e-11
## hsa05323 Rheumatoid arthritis                   6.680686e-13 2.304837e-10
## hsa04659 Th17 cell differentiation              6.316017e-12 2.172710e-09
## hsa04630 JAK-STAT signaling pathway             7.727864e-11 2.650657e-08
## hsa04062 Chemokine signaling pathway            8.713272e-11 2.979939e-08

sum(res_kegg$FDR <= 0.05)

## [1] 42

Au seuil de 5%, on détecte 42 pathways KEGG enrichis en gènes dérégulés. On peut les représenter de cette manière :

res_kegg %>%
  filter(FDR <= 0.05) %>%
  mutate(kegg = gsub("_", " ", kegg)) %>%
  ggplot(aes(x=OR, y=reorder(kegg, OR), color=FDR, size=ngenes)) +
  geom_point() +
  labs(size="# genes", color="FDR") +
  xlab("Odds-ratio") +
  ylab("Gene-set") +
  ggtitle("Pathways KEGG enrichis") +
  theme_light()

Remarque : on aurait pu visualiser les Fold-Enrichment à la place des Odds-Ratios.

Cette analyse d’enrichissement utilise la méthode ORA (Over-Representation Analysis) mais il existe beaucoup d’autres méthodes dont les objectifs et hypothèses sont légèrement différents.

3.2 CAMERA : Competitive Gene Set Test Accounting for Inter-gene Correlation

La méthode CAMERA du package R limma utilise l’ensemble des résultats plutôt que de filtrer en utilisant un seuil sur la P-valeur ajustée (par exemple). Le code R ci-dessous permet de réaliser une telle analyse, en utilisant en entrée l’ensemble des statistiques de tests.

Remarque : la statistique de test d’un gène est défini comme :

\[ \text{stat} = \frac{\log_2(\text{Fold-Change})}{\text{standard-error}} \]

res <- res %>% mutate(stat=log2FoldChange/lfcSE)
statistic <- res %>% pull(stat)
names(statistic) <- res %>% pull(entrezgene)
statistic <- statistic[!is.na(statistic)]
print(head(statistic))

##         60       3002        567       2023         71       2597 
##  1.7203089  0.6955611 -1.0554111  0.9433145  1.9450511  0.5266237

res.cam <- cameraPR(statistic=statistic, index=kegg)
print(head(res.cam))

##                                                                        NGenes Direction       PValue          FDR
## hsa05323 Rheumatoid arthritis                                              93      Down 2.662690e-27 9.239535e-25
## hsa04060 Cytokine-cytokine receptor interaction                           294      Down 1.401363e-23 2.431364e-21
## hsa04061 Viral protein interaction with cytokine and cytokine receptor     99      Down 3.669337e-23 4.244199e-21
## hsa05332 Graft-versus-host disease                                         42      Down 1.859811e-22 1.613386e-20
## hsa04940 Type I diabetes mellitus                                          43      Down 2.179444e-17 1.512534e-15
## hsa05321 Inflammatory bowel disease                                        65      Down 4.155388e-16 2.403200e-14

Le gene-set hsa04060 Cytokine-cytokine receptor interaction est détecté comme étant plutôt down-régulé, et on peut visualiser les statistiques de test des gènes qui le composent :

res %>%
  filter(entrezgene %in% kegg[["hsa04060 Cytokine-cytokine receptor interaction"]]) %>%
  ggplot(aes(x=stat)) +
  geom_histogram(bins=50, color="white", fill="grey") +
  geom_vline(xintercept=0, linetype="dashed") +
  xlab("Statistique de test (par gène)") +
  ylab("Frequency") +
  ggtitle("Pathway hsa04060") +
  theme_light()

3.3 GSEA : Gene Set Enrichment Analysis

La méthode GSEA est également souvent utilisée et permet, comme CAMERA, de détecter si un gene-set est plutôt up-régulé ou down-régulé. Nous utilisons ici le package fgsea pour l’illustrer, avec en entrée l’ensemble des statistiques de tests :

res.gsea <- fgseaMultilevel(pathways=kegg, stats=statistic, eps=0) %>%
  arrange(padj)

## Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (7.91% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

print(head(res.gsea))

##                                                                   pathway         pval         padj   log2err
## 1:                        hsa04060 Cytokine-cytokine receptor interaction 1.139811e-20 3.955145e-18 1.1778933
## 2:                                          hsa05323 Rheumatoid arthritis 8.464193e-09 1.468537e-06 0.7477397
## 3:                                     hsa05332 Graft-versus-host disease 1.943499e-07 2.247981e-05 0.6901325
## 4:                                   hsa03013 Nucleocytoplasmic transport 1.065354e-06 7.493048e-05 0.6435518
## 5:                                      hsa04940 Type I diabetes mellitus 1.079690e-06 7.493048e-05 0.6435518
## 6: hsa04061 Viral protein interaction with cytokine and cytokine receptor 1.878431e-06 9.311652e-05 0.6435518
##            ES       NES size                             leadingEdge
## 1: -0.7604852 -2.343556  294 112744,3562,414062,6349,1437,388372,...
## 2: -0.7954374 -2.148027   93     414062,6349,1437,6348,6364,3552,...
## 3: -0.8788751 -2.083855   42        3558,3552,3458,941,3824,3122,...
## 4:  0.7058413  1.978587  108    9883,55706,3843,23279,3839,30000,...
## 5: -0.8599076 -2.045707   43        3558,3552,3458,941,3122,3123,...
## 6: -0.7447774 -2.048142   99   414062,6349,388372,9560,6348,6364,...

On peut à nouveau illustrer l’enrichissement du pathway hsa04060 Cytokine-cytokine receptor interaction :

plotEnrichment(kegg[["hsa04060 Cytokine-cytokine receptor interaction"]], statistic) + 
  labs(title="hsa04060 Cytokine-cytokine receptor interaction")

Enfin, la figure ci-dessous permet de visualiser les résultats pour les 10 pathways les plus up- et down-régulés :

topPathwaysUp <- res.gsea[ES > 0][head(order(pval), n=10), pathway]
topPathwaysDown <- res.gsea[ES < 0][head(order(pval), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))
plotGseaTable(kegg[topPathways], statistic, res.gsea, gseaParam=0.5)

Questions/remarques :

Quel background utiliser : tous les gènes ? uniquement ceux présents dans au moins un gene-set ?
Quels gene-sets tester (KEGG, GO, Reactome, WikiPathways…) ? Impact sur l’ajustement des p-valeurs ?
Quelle base de données choisir pour importer les gene-sets : packages R (EnrichmentBrowser, msigdbr, gage, …) ? sites web ? fichiers GMT ?
Quel est l’impact de la version de l’annotation ?
Quelle méthode choisir : ORA, CAMERA, GSEA ?
Avec la méthode ORA : tester l’enrichissement en gènes dérégulés ? up-régulés ? down-régulés ?
Comment faire lorsqu’on doit tester l’enrichissement pour plusieurs comparaisons (e.g. 0h - 4h - 8h) ?
Que faire avec les organismes plus ou moins bien annotés ?

Pour aller plus loin :

Di Wu, Gordon K. Smyth, Camera: a competitive gene set test accounting for inter-gene correlation, Nucleic Acids Research, Volume 40, Issue 17, 1 September 2012, Page e133 [Lien]
Aravind Subramanian, Pablo Tamayo, Vamsi K. Mootha et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences Oct 2005, 102 (43) 15545-15550 [Lien]
Nguyen, TM., Shafi, A., Nguyen, T. et al. Identifying significantly impacted pathways: a comprehensive review and assessment. Genome Biol 20, 203 (2019) [Lien]
Ludwig Geistlinger, Gergely Csaba, Mara Santarelli et al, Toward a gold standard for benchmarking gene set enrichment analysis, Briefings in Bioinformatics, Volume 22, Issue 1, January 2021, Pages 545–556 [Lien]

3.4 Bonus : import des pathways KEGG avec msigdbr

Dans les sections précédentes, nous avons utilisé des pathways KEGG importés à partir du package R gage. Il est également possible de les importer en utilisant le package msigdbr, lui même basé sur ces données :

kegg.msigdb <- msigdbr(species="Homo sapiens",
                       category="C2",
                       subcategory="KEGG")
print(head(as.data.frame(kegg.msigdb)))

##   gs_cat gs_subcat               gs_name gene_symbol entrez_gene    ensembl_gene human_gene_symbol human_entrez_gene
## 1     C2   CP:KEGG KEGG_ABC_TRANSPORTERS       ABCA1          19 ENSG00000165029             ABCA1                19
## 2     C2   CP:KEGG KEGG_ABC_TRANSPORTERS      ABCA10       10349 ENSG00000154263            ABCA10             10349
## 3     C2   CP:KEGG KEGG_ABC_TRANSPORTERS      ABCA12       26154 ENSG00000144452            ABCA12             26154
## 4     C2   CP:KEGG KEGG_ABC_TRANSPORTERS      ABCA13      154664 ENSG00000179869            ABCA13            154664
## 5     C2   CP:KEGG KEGG_ABC_TRANSPORTERS       ABCA2          20 ENSG00000107331             ABCA2                20
## 6     C2   CP:KEGG KEGG_ABC_TRANSPORTERS       ABCA3          21 ENSG00000167972             ABCA3                21
##   human_ensembl_gene  gs_id gs_pmid gs_geoid gs_exact_source                                              gs_url
## 1    ENSG00000165029 M11911                         hsa02010 http://www.genome.jp/kegg/pathway/hsa/hsa02010.html
## 2    ENSG00000154263 M11911                         hsa02010 http://www.genome.jp/kegg/pathway/hsa/hsa02010.html
## 3    ENSG00000144452 M11911                         hsa02010 http://www.genome.jp/kegg/pathway/hsa/hsa02010.html
## 4    ENSG00000179869 M11911                         hsa02010 http://www.genome.jp/kegg/pathway/hsa/hsa02010.html
## 5    ENSG00000107331 M11911                         hsa02010 http://www.genome.jp/kegg/pathway/hsa/hsa02010.html
## 6    ENSG00000167972 M11911                         hsa02010 http://www.genome.jp/kegg/pathway/hsa/hsa02010.html
##     gs_description
## 1 ABC transporters
## 2 ABC transporters
## 3 ABC transporters
## 4 ABC transporters
## 5 ABC transporters
## 6 ABC transporters

length(unique(kegg.msigdb$gs_name))

## [1] 186

On constate que 186 pathways KEGG sont répertoriés dans le tableau alors que gage a permis d’en importer 347.

Question : comment faire pour comparer le contenu des pathways importés grâce à msigdbr et gage ?

4 Analyse statistique d’un design complexe

Dans les sections ci-dessus, nous avons analysé des données provenant de designs expérimentaux simples. Or, pour répondre à certaines questions biologiques, il est nécessaire de créer des designs dits complexes, i.e. incluants plusieurs facteurs.

Nous allons étudier ici des données RNA-Seq provenant d’un tel design : l’objectif est d’évaluer l’effet d’un traitement (T vs C) sur deux souches d’une bactérie (WT et KO). Les deux facteurs/variables étudié(e)s sont donc le traitement et la souche, chacun prenant deux niveaux/modalités.

Question : que peut-on faire à partir de ces données et comment pouvons-nous les analyser ?

4.1 Ecriture de l’interaction dans le modèle

rm(list=ls())
load("data_interaction.RData")
ls()

## [1] "counts" "target"

print(target)

##         label strain treatment replicate
## WT-C-1 WT-C-1     WT         C        r1
## WT-C-2 WT-C-2     WT         C        r2
## WT-C-3 WT-C-3     WT         C        r3
## WT-T-1 WT-T-1     WT         T        r1
## WT-T-2 WT-T-2     WT         T        r2
## WT-T-3 WT-T-3     WT         T        r3
## KO-C-1 KO-C-1     KO         C        r1
## KO-C-2 KO-C-2     KO         C        r2
## KO-C-3 KO-C-3     KO         C        r3
## KO-T-1 KO-T-1     KO         T        r1
## KO-T-2 KO-T-2     KO         T        r2
## KO-T-3 KO-T-3     KO         T        r3

print(head(counts))

##       WT-C-1 WT-C-2 WT-C-3 WT-T-1 WT-T-2 WT-T-3 KO-C-1 KO-C-2 KO-C-3 KO-T-1 KO-T-2 KO-T-3
## gene1    412    590    656    698    537    737     35     47     36      9     21     19
## gene2   7483  11099   9825  20678  13621  19940  13170  14153  11198   8475   9479   8464
## gene3   7005  11450  10400  11932   9691  14038   8533  10837   7824   6814   8318   6668
## gene4   6774  12938   9603  13754  11554  16687   4405   4446   3661   4045   6133   4075
## gene5   2706   4018   3604   4927   3839   5663   4526   5084   4326   4435   4985   4502
## gene6  25276  42273  38634  49580  40575  56964  36755  40690  30480  28131  32439  27851

# mettre WT en reference de strain et C en reference de treatment
target$strain <- factor(target$strain, levels=c("WT", "KO"))
target$treatment <- factor(target$treatment, levels=c("C", "T"))

Le code R ci-dessous permet de créer l’objet dds avec le design adéquat et de réaliser les 3 étapes de l’analyse différentielle (normalisation, estimation de la dispersion et modèle linéaire généralisé) :

dds <- DESeqDataSetFromMatrix(countData=counts, 
                              colData=target, 
                              design=~replicate + strain + treatment + strain:treatment)

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in design formula are characters, converting
## to factors

dds <- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)
resultsNames(dds)

## [1] "Intercept"           "replicate_r2_vs_r1"  "replicate_r3_vs_r1"  "strain_KO_vs_WT"     "treatment_T_vs_C"   
## [6] "strainKO.treatmentT"

4.2 Visualisation de l’interaction

L’analyse en composante principale permet de visualiser la structure des données :

rld <- rlog(dds)
pca <- PCA(t(assay(rld)), scale.unit=FALSE, ncp=3, graph=FALSE)
fviz_pca_ind(pca, axes=c(1, 2), repel=TRUE, mean.point=FALSE, 
             habillage=factor(paste(target$strain, target$treatment, sep="-")))

fviz_pca_ind(pca, axes=c(1, 3), repel=TRUE, mean.point=FALSE, 
             habillage=factor(paste(target$strain, target$treatment, sep="-")))

Sur le premier plan factoriel formé par les axes 1 et 2 nous visualisons l’effet du traitement (axe 1) et l’effet de la souche (axe 2). Ensuite, l’axe 3 permet de visualiser l’interaction souche-traitement : l’effet du traitement sur les WT (haut vers le bas) est opposé à l’effet du traitement sur KO (bas vers le haut).

On peut également visualiser l’interaction pour un gène spécifique (bien choisi) :

ncounts <- counts(dds, norm=TRUE)
d <- data.frame(target, ncounts=ncounts["gene914",])
print(d)

##         label strain treatment replicate   ncounts
## WT-C-1 WT-C-1     WT         C        r1  57.68749
## WT-C-2 WT-C-2     WT         C        r2  43.65947
## WT-C-3 WT-C-3     WT         C        r3  44.75770
## WT-T-1 WT-T-1     WT         T        r1 491.84167
## WT-T-2 WT-T-2     WT         T        r2 276.67419
## WT-T-3 WT-T-3     WT         T        r3 338.98342
## KO-C-1 KO-C-1     KO         C        r1 275.29097
## KO-C-2 KO-C-2     KO         C        r2 116.71415
## KO-C-3 KO-C-3     KO         C        r3 266.30771
## KO-T-1 KO-T-1     KO         T        r1  44.58641
## KO-T-2 KO-T-2     KO         T        r2  23.61180
## KO-T-3 KO-T-3     KO         T        r3  27.25240

ggplot(d, aes(x=strain, y=ncounts, color=treatment, shape=replicate)) +
  geom_beeswarm(dodge.width=0.3, size=4) +
  scale_y_continuous(trans = log2_trans(),
    breaks = trans_breaks("log2", function(x) 2^x),
    labels = trans_format("log2", math_format(2^.x))) +
  xlab("Condition") +
  ylab("Normalized counts (log2 scale)") +
  labs(shape="Replicate", color="Treatment") +
  ggtitle("Illustration of the donor effect for gene gene914") +
  theme_light()

4.3 Ecriture des contrastes

Nous pouvons ensuite utiliser DESeq2 et des contrastes adaptés pour evaluer l’effet du traitement sur chaque souche et ensuite comparer ces deux effets :

T vs C pour les WT
T vs C pour les KO
différence entre (2) et (1)

# ExploreModelMatrix() pour lancer l'application Shiny
res_TvsC.WT <- lfcShrink(dds, contrast=c(0, 0, 0, 0, 1, 0), type="ashr", quiet=TRUE)
res_TvsC.KO <- lfcShrink(dds, contrast=c(0, 0, 0, 0, 1, 1), type="ashr", quiet=TRUE)
res_interaction <- lfcShrink(dds, contrast=c(0, 0, 0, 0, 0, 1), type="ashr", quiet=TRUE)
summary(res_TvsC.WT, alpha=0.05)

## 
## out of 1716 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 487, 28%
## LFC < 0 (down)     : 451, 26%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

summary(res_TvsC.KO, alpha=0.05)

## 
## out of 1716 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 463, 27%
## LFC < 0 (down)     : 506, 29%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

summary(res_interaction, alpha=0.05)

## 
## out of 1716 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 0 (up)       : 398, 23%
## LFC < 0 (down)     : 430, 25%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Afin de comparer les effets du traitement chez les WT et chez les KO, nous assemblons les résultats au sein d’un même tableau et nous créons de nouvelles variables indiquant d’une part si un gène est en interaction (DE_int) et d’autre part si les deux effets sont significatifs (DE_both_strains). Ensuite nous pouvons tracer la figure :

d <- data.frame(gene=rownames(res_TvsC.WT),
                log2FC_TvsC.WT=res_TvsC.WT$log2FoldChange,
                padj_TvsC.WT=res_TvsC.WT$padj,
                log2FC_TvsC.KO=res_TvsC.KO$log2FoldChange,
                padj_TvsC.KO=res_TvsC.KO$padj,
                log2FC_int=res_interaction$log2FoldChange,
                padj_int=res_interaction$padj)
d <- d %>%
  mutate(DE_int=ifelse(!is.na(padj_int) & padj_int <= 0.05, "Yes", "No"),
         DE_both_strains=ifelse(!is.na(padj_TvsC.WT) & padj_TvsC.WT <= 0.05 & !is.na(padj_TvsC.KO) & padj_TvsC.KO <= 0.05, "Yes", "No"))
print(head(d))

##    gene log2FC_TvsC.WT padj_TvsC.WT log2FC_TvsC.KO padj_TvsC.KO  log2FC_int    padj_int DE_int DE_both_strains
## 1 gene1    -0.09689838 5.178970e-01     -0.6354706 6.810186e-03 -0.47905538 0.052770555     No              No
## 2 gene2     0.49374402 8.312332e-06     -0.1073933 3.601136e-01 -0.58749593 0.000169735    Yes              No
## 3 gene3    -0.05797363 6.122993e-01      0.0907722 4.114685e-01  0.14584864 0.357576667     No              No
## 4 gene4     0.13107436 3.794633e-01      0.5116283 6.742476e-04  0.31723593 0.112832737     No              No
## 5 gene5     0.08540833 4.220296e-01      0.3899746 7.871748e-05  0.28317668 0.044938613    Yes              No
## 6 gene6     0.09137723 3.921496e-01      0.1185208 2.488123e-01  0.02705004 0.882849983     No              No

ggplot(d, aes(x=log2FC_TvsC.WT, y=log2FC_TvsC.KO, color=DE_int, size=DE_both_strains)) +
  geom_point(alpha=0.5) +
  scale_color_manual(values=c("Yes"="red", "No"="black")) +
  scale_size_manual(values=c("Yes"=1, "No"=0.4)) +
  geom_hline(yintercept=0, lty=2, color="grey") +
  geom_vline(xintercept=0, lty=2, color="grey") +
  geom_abline(slope=1, intercept=0, lty=2, color="grey") +
  labs(color="Strain-condition\ninteraction", size="DE for both\nWT and KO") +
  xlab(bquote(log[2](FC)~"T vs C for WT")) +
  ylab(bquote(log[2](FC)~"T vs C for KO")) +
  ggtitle("Illustration of the interaction") +
  theme_light()

Sur cette figure, les points éloignés de la diagonale (droite \(Y = X\)) correspondent aux gènes en interaction, i.e. pour lesquels l’effet du traitement diffère selon la souche. Autrement dit, le log2(Fold-Change) chez les WT est très différent du log2(Fold-Change) chez les KO. Cette différence d’effets est significative lorsque le point est rouge (au seuil de 5%), et les effets individuels (chez les WT et KO) sont significatifs simultanément si le point est agrandi.

4.4 Heatmap

On peut sélectionner ces gènes et les visualiser sous la forme d’une heatmap (en ajustant sur l’effet du réplicat) :

int_genes <- d %>%
  filter(padj_int <= 0.05) %>%
  pull(gene)
mat_hm <- assay(rlog(dds))
mat_hm_adj <- removeBatchEffect(mat_hm, batch=target$replicate)
pheatmap(mat_hm_adj[int_genes,],
         scale="row",
         show_rownames=FALSE,
         clustering_distance_rows="euclidean",
         clustering_distance_cols="euclidean",
         clustering_method="ward.D",
         color=colorRampPalette(rev(brewer.pal(9, "RdBu")))(200))

On voit sur cette heatmap que pour chaque ligne l’effet du traitement est opposé entre les deux souches (e.g. rouge vs bleu pour les WT contre bleu vs rouge pour les KO).

Le code ci-dessous permet de visualiser les log2(Fold-Change) :

rownames(d) <- d$gene
d <- d %>%
  filter(padj_int < 0.005 & sign(log2FC_TvsC.WT) != sign(log2FC_TvsC.KO)) %>%
  dplyr::select(log2FC_TvsC.WT, log2FC_TvsC.KO) %>%
  as.matrix()
pheatmap(d, 
         scale="none", 
         show_rownames=FALSE,
         cluster_cols=FALSE,
         clustering_distance_rows="euclidean",
         angle_col=45,
         fontsize_row=5,
         color=colorRampPalette(rev(brewer.pal(9, "RdBu")))(100),
         breaks=seq(-4, 4, length=101))

Question : comment analyser un design expérimental comportant plus de 3 facteurs d’intérêt ?

Liens et références intéressants

RNA-Seq

Chaine YouTube StatQuest
Site DoItYourself Transcriptomics
Site RNA-Seqlopedia
Article sur la normalisation : Ciaran Evans, Johanna Hardin, Daniel M Stoebel, Selecting between-sample RNA-Seq normalization methods from the perspective of their assumptions, Briefings in Bioinformatics, Volume 19, Issue 5, September 2018, Pages 776–792 [Lien]

R

Livre R for Data Science
Livre R Cookbook
Livre Modern Data Science with R
Livre Computational Genomics with R
Livre ggplot2: elegant graphics for data analysis
Site de référence ggplot2
Site extensions ggplot2
Site A ggplot2 Tutorial for Beautiful Plotting in R
Package R patchwork
Livre Circular Visualization in R
Livre ComplexHeatmap
Palettes de couleurs

Informations sur la session R

Dans un souci de reproductibilité, la commande ci-dessous permet d’afficher l’ensemble des packages utilisés ainsi que leurs versions :

devtools::session_info()

## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
##  setting  value
##  version  R version 4.2.1 (2022-06-23)
##  os       macOS Big Sur ... 10.16
##  system   x86_64, darwin17.0
##  ui       X11
##  language (EN)
##  collate  en_US.UTF-8
##  ctype    en_US.UTF-8
##  tz       Europe/Paris
##  date     2023-05-24
##  pandoc   2.19.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
## 
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
##  package              * version   date (UTC) lib source
##  abind                  1.4-5     2016-07-21 [1] CRAN (R 4.2.0)
##  annotate               1.74.0    2022-04-26 [1] Bioconductor
##  AnnotationDbi          1.58.0    2022-04-26 [1] Bioconductor
##  ashr                   2.2-54    2022-02-22 [1] CRAN (R 4.2.0)
##  assertthat             0.2.1     2019-03-21 [1] CRAN (R 4.2.0)
##  babelgene              22.9      2022-09-29 [1] CRAN (R 4.2.0)
##  backports              1.4.1     2021-12-13 [1] CRAN (R 4.2.0)
##  beeswarm               0.4.0     2021-06-01 [1] CRAN (R 4.2.0)
##  Biobase              * 2.56.0    2022-04-26 [1] Bioconductor
##  BiocFileCache          2.4.0     2022-04-26 [1] Bioconductor
##  BiocGenerics         * 0.42.0    2022-04-26 [1] Bioconductor
##  BiocParallel           1.30.4    2022-10-13 [1] Bioconductor
##  biomaRt              * 2.52.0    2022-04-26 [1] Bioconductor
##  Biostrings             2.64.1    2022-08-18 [1] Bioconductor
##  bit                    4.0.5     2022-11-15 [1] CRAN (R 4.2.0)
##  bit64                  4.0.5     2020-08-30 [1] CRAN (R 4.2.0)
##  bitops                 1.0-7     2021-04-24 [1] CRAN (R 4.2.0)
##  blob                   1.2.3     2022-04-10 [1] CRAN (R 4.2.0)
##  broom                  1.0.3     2023-01-25 [1] CRAN (R 4.2.0)
##  bslib                  0.4.2     2022-12-16 [1] CRAN (R 4.2.0)
##  cachem                 1.0.6     2021-08-19 [1] CRAN (R 4.2.0)
##  callr                  3.7.3     2022-11-02 [1] CRAN (R 4.2.0)
##  car                    3.1-1     2022-10-19 [1] CRAN (R 4.2.1)
##  carData                3.0-5     2022-01-06 [1] CRAN (R 4.2.0)
##  cli                    3.6.0     2023-01-09 [1] CRAN (R 4.2.0)
##  cluster                2.1.4     2022-08-22 [1] CRAN (R 4.2.0)
##  coda                   0.19-4    2020-09-30 [1] CRAN (R 4.2.0)
##  codetools              0.2-19    2023-02-01 [1] CRAN (R 4.2.0)
##  colorspace             2.1-0     2023-01-23 [1] CRAN (R 4.2.0)
##  cowplot                1.1.1     2020-12-30 [1] CRAN (R 4.2.0)
##  crayon                 1.5.2     2022-09-29 [1] CRAN (R 4.2.0)
##  curl                   5.0.0     2023-01-12 [1] CRAN (R 4.2.0)
##  data.table             1.14.6    2022-11-16 [1] CRAN (R 4.2.0)
##  DBI                    1.1.3     2022-06-18 [1] CRAN (R 4.2.0)
##  dbplyr                 2.3.0     2023-01-16 [1] CRAN (R 4.2.0)
##  DelayedArray           0.22.0    2022-04-26 [1] Bioconductor
##  DESeq2               * 1.36.0    2022-04-26 [1] Bioconductor
##  devtools             * 2.4.5     2022-10-11 [1] CRAN (R 4.2.0)
##  digest                 0.6.31    2022-12-11 [1] CRAN (R 4.2.0)
##  dplyr                * 1.1.0     2023-01-29 [1] CRAN (R 4.2.0)
##  DT                     0.27      2023-01-17 [1] CRAN (R 4.2.0)
##  ellipsis               0.3.2     2021-04-29 [1] CRAN (R 4.2.0)
##  emmeans                1.8.4-1   2023-01-17 [1] CRAN (R 4.2.0)
##  estimability           1.4.1     2022-08-05 [1] CRAN (R 4.2.0)
##  evaluate               0.20      2023-01-17 [1] CRAN (R 4.2.0)
##  ExploreModelMatrix   * 1.8.0     2022-04-26 [1] Bioconductor
##  factoextra           * 1.0.7     2020-04-01 [1] CRAN (R 4.2.0)
##  FactoMineR           * 2.7       2022-12-14 [1] CRAN (R 4.2.0)
##  fansi                  1.0.4     2023-01-22 [1] CRAN (R 4.2.0)
##  farver                 2.1.1     2022-07-06 [1] CRAN (R 4.2.0)
##  fastmap                1.1.0     2021-01-25 [1] CRAN (R 4.2.0)
##  fastmatch              1.1-3     2021-07-23 [1] CRAN (R 4.2.0)
##  fgsea                * 1.22.0    2022-04-26 [1] Bioconductor
##  filelock               1.0.2     2018-10-05 [1] CRAN (R 4.2.0)
##  flashClust             1.01-2    2012-08-21 [1] CRAN (R 4.2.0)
##  fs                     1.6.1     2023-02-06 [1] CRAN (R 4.2.0)
##  gage                 * 2.46.1    2022-08-28 [1] Bioconductor
##  genefilter             1.78.0    2022-04-26 [1] Bioconductor
##  geneplotter            1.74.0    2022-04-26 [1] Bioconductor
##  generics               0.1.3     2022-07-05 [1] CRAN (R 4.2.0)
##  GenomeInfoDb         * 1.32.4    2022-09-06 [1] Bioconductor
##  GenomeInfoDbData       1.2.8     2022-07-22 [1] Bioconductor
##  GenomicRanges        * 1.48.0    2022-04-26 [1] Bioconductor
##  ggbeeswarm           * 0.7.1     2022-12-16 [1] CRAN (R 4.2.0)
##  ggplot2              * 3.4.1     2023-02-10 [1] CRAN (R 4.2.0)
##  ggpubr                 0.6.0     2023-02-10 [1] CRAN (R 4.2.0)
##  ggrepel              * 0.9.3     2023-02-03 [1] CRAN (R 4.2.0)
##  ggsignif               0.6.4     2022-10-13 [1] CRAN (R 4.2.0)
##  ggVennDiagram        * 1.2.2     2022-09-08 [1] CRAN (R 4.2.0)
##  glue                   1.6.2     2022-02-24 [1] CRAN (R 4.2.0)
##  GO.db                  3.15.0    2023-04-06 [1] Bioconductor
##  graph                  1.74.0    2022-04-26 [1] Bioconductor
##  gridExtra              2.3       2017-09-09 [1] CRAN (R 4.2.0)
##  gtable                 0.3.1     2022-09-01 [1] CRAN (R 4.2.0)
##  highr                  0.10      2022-12-22 [1] CRAN (R 4.2.0)
##  hms                    1.1.2     2022-08-19 [1] CRAN (R 4.2.0)
##  htmltools              0.5.4     2022-12-07 [1] CRAN (R 4.2.0)
##  htmlwidgets            1.6.1     2023-01-07 [1] CRAN (R 4.2.0)
##  httpuv                 1.6.8     2023-01-12 [1] CRAN (R 4.2.0)
##  httr                   1.4.4     2022-08-17 [1] CRAN (R 4.2.0)
##  invgamma               1.1       2017-05-07 [1] CRAN (R 4.2.0)
##  IRanges              * 2.30.1    2022-08-18 [1] Bioconductor
##  irlba                  2.3.5.1   2022-10-03 [1] CRAN (R 4.2.0)
##  jquerylib              0.1.4     2021-04-26 [1] CRAN (R 4.2.0)
##  jsonlite               1.8.4     2022-12-06 [1] CRAN (R 4.2.0)
##  KEGGREST               1.36.3    2022-07-14 [1] Bioconductor
##  knitr                  1.42      2023-01-25 [1] CRAN (R 4.2.0)
##  labeling               0.4.2     2020-10-20 [1] CRAN (R 4.2.0)
##  later                  1.3.0     2021-08-18 [1] CRAN (R 4.2.0)
##  lattice                0.20-45   2021-09-22 [1] CRAN (R 4.2.1)
##  lazyeval               0.2.2     2019-03-15 [1] CRAN (R 4.2.0)
##  leaps                  3.1       2020-01-16 [1] CRAN (R 4.2.0)
##  lifecycle              1.0.3     2022-10-07 [1] CRAN (R 4.2.0)
##  limma                * 3.52.4    2022-09-27 [1] Bioconductor
##  locfit                 1.5-9.7   2023-01-02 [1] CRAN (R 4.2.0)
##  magrittr               2.0.3     2022-03-30 [1] CRAN (R 4.2.0)
##  MASS                   7.3-58.2  2023-01-23 [1] CRAN (R 4.2.0)
##  Matrix                 1.5-3     2022-11-11 [1] CRAN (R 4.2.0)
##  MatrixGenerics       * 1.8.1     2022-06-30 [1] Bioconductor
##  matrixStats          * 0.63.0    2022-11-18 [1] CRAN (R 4.2.0)
##  memoise                2.0.1     2021-11-26 [1] CRAN (R 4.2.0)
##  mime                   0.12      2021-09-28 [1] CRAN (R 4.2.0)
##  miniUI                 0.1.1.1   2018-05-18 [1] CRAN (R 4.2.0)
##  mixsqp                 0.3-48    2022-11-16 [1] CRAN (R 4.2.0)
##  msigdbr              * 7.5.1     2022-03-30 [1] CRAN (R 4.2.0)
##  multcomp               1.4-22    2023-02-10 [1] CRAN (R 4.2.0)
##  multcompView           0.1-8     2019-12-19 [1] CRAN (R 4.2.0)
##  munsell                0.5.0     2018-06-12 [1] CRAN (R 4.2.0)
##  mvtnorm                1.1-3     2021-10-08 [1] CRAN (R 4.2.0)
##  pheatmap             * 1.0.12    2019-01-04 [1] CRAN (R 4.2.0)
##  pillar                 1.8.1     2022-08-19 [1] CRAN (R 4.2.0)
##  pkgbuild               1.4.0     2022-11-27 [1] CRAN (R 4.2.0)
##  pkgconfig              2.0.3     2019-09-22 [1] CRAN (R 4.2.0)
##  pkgload                1.3.2     2022-11-16 [1] CRAN (R 4.2.0)
##  plotly               * 4.10.1    2022-11-07 [1] CRAN (R 4.2.0)
##  plyr                   1.8.8     2022-11-11 [1] CRAN (R 4.2.0)
##  png                    0.1-8     2022-11-29 [1] CRAN (R 4.2.0)
##  prettyunits            1.1.1     2020-01-24 [1] CRAN (R 4.2.0)
##  processx               3.8.0     2022-10-26 [1] CRAN (R 4.2.0)
##  profvis                0.3.7     2020-11-02 [1] CRAN (R 4.2.0)
##  progress               1.2.2     2019-05-16 [1] CRAN (R 4.2.0)
##  promises               1.2.0.1   2021-02-11 [1] CRAN (R 4.2.0)
##  ps                     1.7.2     2022-10-26 [1] CRAN (R 4.2.0)
##  purrr                  1.0.1     2023-01-10 [1] CRAN (R 4.2.0)
##  R6                     2.5.1     2021-08-19 [1] CRAN (R 4.2.0)
##  rappdirs               0.3.3     2021-01-31 [1] CRAN (R 4.2.0)
##  RColorBrewer         * 1.1-3     2022-04-03 [1] CRAN (R 4.2.0)
##  Rcpp                   1.0.10    2023-01-22 [1] CRAN (R 4.2.0)
##  RCurl                  1.98-1.10 2023-01-27 [1] CRAN (R 4.2.0)
##  remotes                2.4.2     2021-11-30 [1] CRAN (R 4.2.0)
##  rintrojs               0.3.2     2022-08-09 [1] CRAN (R 4.2.0)
##  rlang                  1.0.6     2022-09-24 [1] CRAN (R 4.2.0)
##  rmarkdown              2.20      2023-01-19 [1] CRAN (R 4.2.0)
##  RSQLite                2.2.20    2022-12-22 [1] CRAN (R 4.2.0)
##  rstatix                0.7.2     2023-02-01 [1] CRAN (R 4.2.0)
##  rstudioapi             0.14      2022-08-22 [1] CRAN (R 4.2.0)
##  RVenn                  1.1.0     2019-07-18 [1] CRAN (R 4.2.0)
##  S4Vectors            * 0.34.0    2022-04-26 [1] Bioconductor
##  sandwich               3.0-2     2022-06-15 [1] CRAN (R 4.2.0)
##  sass                   0.4.5     2023-01-24 [1] CRAN (R 4.2.0)
##  scales               * 1.2.1     2022-08-20 [1] CRAN (R 4.2.0)
##  scatterplot3d          0.3-42    2022-09-08 [1] CRAN (R 4.2.0)
##  sessioninfo            1.2.2     2021-12-06 [1] CRAN (R 4.2.0)
##  shiny                  1.7.4     2022-12-15 [1] CRAN (R 4.2.0)
##  shinydashboard         0.7.2     2021-09-30 [1] CRAN (R 4.2.0)
##  shinyjs                2.1.0     2021-12-23 [1] CRAN (R 4.2.0)
##  SQUAREM                2021.1    2021-01-13 [1] CRAN (R 4.2.0)
##  stringi                1.7.12    2023-01-11 [1] CRAN (R 4.2.0)
##  stringr                1.5.0     2022-12-02 [1] CRAN (R 4.2.0)
##  SummarizedExperiment * 1.26.1    2022-05-01 [1] Bioconductor
##  survival               3.5-3     2023-02-12 [1] CRAN (R 4.2.0)
##  TH.data                1.1-1     2022-04-26 [1] CRAN (R 4.2.0)
##  tibble                 3.1.8     2022-07-22 [1] CRAN (R 4.2.1)
##  tidyr                  1.3.0     2023-01-24 [1] CRAN (R 4.2.0)
##  tidyselect             1.2.0     2022-10-10 [1] CRAN (R 4.2.0)
##  truncnorm              1.0-8     2018-02-27 [1] CRAN (R 4.2.0)
##  UpSetR               * 1.4.0     2019-05-22 [1] CRAN (R 4.2.0)
##  urlchecker             1.0.1     2021-11-30 [1] CRAN (R 4.2.0)
##  usethis              * 2.1.6     2022-05-25 [1] CRAN (R 4.2.0)
##  utf8                   1.2.3     2023-01-31 [1] CRAN (R 4.2.0)
##  vctrs                  0.5.2     2023-01-23 [1] CRAN (R 4.2.0)
##  vipor                  0.4.5     2017-03-22 [1] CRAN (R 4.2.0)
##  viridisLite            0.4.1     2022-08-22 [1] CRAN (R 4.2.0)
##  withr                  2.5.0     2022-03-03 [1] CRAN (R 4.2.0)
##  xfun                   0.37      2023-01-31 [1] CRAN (R 4.2.0)
##  XML                    3.99-0.13 2022-12-04 [1] CRAN (R 4.2.0)
##  xml2                   1.3.3     2021-11-30 [1] CRAN (R 4.2.0)
##  xtable                 1.8-4     2019-04-21 [1] CRAN (R 4.2.0)
##  XVector                0.36.0    2022-04-26 [1] Bioconductor
##  yaml                   2.3.7     2023-01-23 [1] CRAN (R 4.2.0)
##  zlibbioc               1.42.0    2022-04-26 [1] Bioconductor
##  zoo                    1.8-11    2022-09-17 [1] CRAN (R 4.2.0)
## 
##  [1] /Library/Frameworks/R.framework/Versions/4.2/Resources/library
## 
## ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

Remarque : la sortie est plus lisible qu’avec le traditionnel sessionInfo().

EBAII niveau 2 - Atelier RNA-Seq

Hugo Varet\(^{1,2}\) - hugo.varet@pasteur.fr

Juin 2023