EBAIIn2_RNASeq_apprenants.Rmd

---
title: "EBAII niveau 2 - Atelier RNA-Seq"
author: "Hugo Varet$^{1,2}$ - hugo.varet@pasteur.fr"
date: "Mai 2021"
output:
    html_document:
      toc: true
      toc_float: true
      number_sections: true
---

$^1$ Hub de Bioinformatique et Biostatistique, Département de Biologie Computationnelle, Institut Pasteur

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

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.width=6, fig.height=5, cache=TRUE,
                      message=FALSE, out.width="60%", fig.align="center")
options(width=120)
```

```{r logo, echo=FALSE, fig.cap="", fig.align="center", out.width = '100%',include=TRUE}
knitr::include_graphics("logo.png")
```

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 :

```{r, warning=FALSE, message=FALSE}
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(EnrichmentBrowser)  # 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(ggVennDiagram)      # 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`.

# 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 :

```{r load_data_1factor}

```

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

## 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 :

 (1) la normalisation,
 (2) l'estimation de la dispersion,
 (3) l'ajustement du modèle linéaire généralisé.

```{r deseq2_1factor}

```

**Remarque** : nous utilisons ici DESeq2 avec le paramétrage par défaut, mais il existe de nombreuses subtilités mathématiques/statistiques. La [vignette](https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html) du package est extrêmement bien faite et explique les différentes options possibles.

## 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.

```{r homoscedasticity_1factor}

```

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

```{r pca_1factor}

```

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

```{r pca_1factor_factominer}

```

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

## 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 :

```{r results_1factor}

```

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é :

```{r maplot_interactif_1factor, warning=FALSE}

```

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

```{r maplot_1factor, warning=FALSE}

```

## 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` : 

```{r check_ids_1factor}

```

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 :

```{r compare_results_venn_1factor}

```

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 :

```{r compare_results_upset_1factor}

```

**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 :

```{r inconvenient_upset_1factor}

```

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,27$ et $-0,29$), 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 :

```{r compare_results_pairwise_1factor}

```

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 :

```{r filters_1factor}

```

## 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.

```{r heatmap_expression_1factor}

```

**Question** : quel message renvoie cette heatmap ? Pourquoi ?

```{r heatmap_expression_ok_1factor}

```

**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 :

```{r heatmap_log2FC_1factor}

```

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 :

```{r heatmap_log2FC_1factor_forced}

```

# 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.

```{r load_data_batch}

```

## Analyse en composantes principales

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

```{r deseq2_batch}

```

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

```{r pca_batch}

```

## 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 :

```{r results_batch}

```

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) :

```{r illustration_batch}

```

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

```{r results_ok_batch}

```

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

## 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% :

```{r heatmap_batch}

```

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()` :

```{r heatmap_ok_batch}

```

**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](https://bioinfo-fr.net/creer-des-heatmaps-a-partir-de-grosses-matrices-en-r) et [là](https://jokergoo.github.io/2020/06/30/rasterization-in-complexheatmap/).
 - Les heatmaps sont de formidables outils de visualisation et on peut imaginer une personnalisation à l'infini, notamment avec le package [ComplexHeatmap](https://jokergoo.github.io/ComplexHeatmap-reference/book/).

## 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 :

```{r ensembl_bm}

```

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* :

```{r ensembl_bm_dup_names}

```

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

```{r ensembl_bm_dup_ens_ids}

```

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

```{r heatmap_names}

```

**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](https://www.bioconductor.org/help/course-materials/2019/CSAMA/materials/lectures/lecture-08a-annotation.html#1) recense différents outils disponibles et leurs fonctionnements.

# 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 EnrichmentBrowser :

```{r import_KEGG}

```

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 :

```{r check_KEGG}

```

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 :

```{r ensembl2entrez}

```

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

```{r ensembl2entrez_dupentrez}

```

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 :

```{r merge_res_entrez}

```

## *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 :

```{r enrichment}

```

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/81}{1152/11488} = 3.07 $$
est significativement supérieur à 1 (P-valeur = $8.86 \times 10^{-6}$). 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+81)}{1152/(1152+11488)} = \frac{23,6}{9,1} = 2,59$$
On peut maintenant répéter l'opération grâce à une boucle sur tous les *pathways* KEGG :

```{r enrichment_loop}

```

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.

```{r enrichment_fdr}

```

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

```{r plot_enrichment, fig.width=8, fig.height=6}

```

**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. 

## 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}} $$

```{r camera}

```

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 :

```{r camera_illustration}

```

## 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 :

```{r fgsea}

```

On peut à nouveau illustrer l'enrichissement du *pathway* `hsa04060_Cytokine-cytokine_receptor_interaction` :

```{r fgsea_illustration}

```

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

```{r fgsea_illustration2}

```

**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](https://academic.oup.com/nar/article/40/17/e133/2411151)]
 - 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](https://www.pnas.org/content/102/43/15545)]
 - Nguyen, TM., Shafi, A., Nguyen, T. et al. *Identifying significantly impacted pathways: a comprehensive review and assessment*. Genome Biol 20, 203 (2019) [[Lien](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1790-4)]
 - 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](https://academic.oup.com/bib/article/22/1/545/5722384)]

## 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 EnrichmentBrowser. Il est également possible de les importer en utilisant le package msigdbr, lui même basé sur ces [données](https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) :

```{r kegg_msigdr}

```

On constate que 186 *pathways* KEGG sont répertoriés dans le tableau alors qu'EnrichmentBrower a permis d'en importer 343.

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

# 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 ?

## Ecriture de l'interaction dans le modèle

```{r load_data_2factors}

```

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é) :

```{r analysis_2factors}

```

## Visualisation de l'interaction

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

```{r acp_2factors}

```

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) :

```{r interaction_1gene}

```

## 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 :

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

```{r results_2factors}

```

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 :

```{r results_2factors2}

```

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.

## 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) :

```{r heatmap_2factors}

```

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) :

```{r heatmap_log2FC_2factors}

```

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

# Liens et références intéressants{.unlisted .unnumbered}

### RNA-Seq{.unlisted .unnumbered}

 - Chaine YouTube [StatQuest](https://www.youtube.com/c/joshstarmer/videos)
 - Site [DoItYourself Transcriptomics](https://diytranscriptomics.com/)
 - Site [RNA-Seqlopedia](https://rnaseq.uoregon.edu/)
 - 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](https://academic.oup.com/bib/article/19/5/776/3056951)]

### R{.unlisted .unnumbered}

 - Livre [*R for Data Science*](https://r4ds.had.co.nz/)
 - Livre [*R Cookbook*](https://rc2e.com/)
 - Livre [*Modern Data Science with R*](https://mdsr-book.github.io/mdsr2e/)
 - Livre [*Computational Genomics with R*](https://compgenomr.github.io/book/)

 - Livre [*ggplot2: elegant graphics for data analysis*](https://ggplot2-book.org/index.html)
 - Site de [référence ggplot2](https://ggplot2.tidyverse.org/reference/)
 - Site [extensions ggplot2](https://exts.ggplot2.tidyverse.org/gallery/)
 - Site [*A ggplot2 Tutorial for Beautiful Plotting in R*](https://www.cedricscherer.com/2019/08/05/a-ggplot2-tutorial-for-beautiful-plotting-in-r/)
 - Package R [patchwork](https://github.com/thomasp85/patchwork)
 - Livre [*Circular Visualization in R*](https://jokergoo.github.io/circlize_book/book/)
 - Livre [ComplexHeatmap](https://jokergoo.github.io/ComplexHeatmap-reference/book/)
 - Palettes de [couleurs](https://github.com/EmilHvitfeldt/r-color-palettes)

# Informations sur la session R{.unlisted .unnumbered}

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

```{r}
devtools::session_info()
```

**Remarque** : la sortie est plus lisible qu'avec le traditionnel `sessionInfo()`.