```{r include=FALSE, echo=FALSE, eval=TRUE}
library(knitr)
library(kableExtra)
# library(formattable)

options(width = 300)
# options(encoding = 'UTF-8')
knitr::opts_chunk$set(
  fig.width = 7, fig.height = 5, 
  fig.path = 'figures/07_tests_multiples',
  fig.align = "center", 
  size = "tiny", 
  echo = TRUE, eval = TRUE, 
  warning = FALSE, message = FALSE, 
  results = FALSE, comment = "")

options(scipen = 12) ## Max number of digits for non-scientific notation
# knitr::asis_output("\\footnotesize")

```


## Se connecter au serveur RStudio de Roscoff

<http://r.sb-roscoff.fr/auth-sign-in>

## Aller dans son dossier de travail

Définir une variable qui indique le chemin du dossier de travail

```{r}
work.dir <- "~/intro_R" 
```

S’il n’existe pas encore, créer le dossier de travail.
(Commande Unix équivalente: "mkdir -p ~/intro_R")

```{r}
dir.create(work.dir, recursive = TRUE, showWarnings = FALSE)
```

Aller dans ce dossier de travail.
(Commande Unix équivalente: "cd ~/intro_R")

```{r}
setwd(work.dir)	
```

Où suis-je ? 
(Commande Unix équivalente: "pwd")

```{r}
getwd()
```

Qu'y a-t-il par ici ? 
(Commande Unix équivalente: "ls")

```{r}
list.files()
```

## R vu comme une calculatrice

```{r}
2 + 3
4 * 5
6 / 4
```

## Notion de variable/objet


```{r}
a <- 2      ## Assigner une valeur à une variable
print(a)    ## Afficher la valeur de la variable a
```

```{r}
b <- 3      ## Assigner une valeur à une seconde variable
c <- a + b  ## Effectuer un calcul avec 2 variables
print(c)    ## Afficher le contenu de la variable c
```


```{r}
a <- 7      ## Changer la valeur de a
print(c)    ## Note: le contenu de c n'est pas modifié
```


## Télécharger un fichier

La commande `download()` permet de télécharger un fichier à partir d'un serveur. 

```{r}
## download.file(url = "https://goo.gl/9QVAg6", destfile = "expression.txt")
```

```{r}
## download.file(url = "https://goo.gl/NQWnHg", destfile = "annotation.csv")
```

## Chargement des données

Charger le contenu du fichier "expression.txt" dans une variable nommée "exprs".

```{r}
exprs <- read.table(file = "expression.txt", header = TRUE, sep = "\t")
```

Accéder à l'aide d'une fonction

```{r}
help(read.table)
```

Notation alternative

```{r}
?read.table
```


## Affichage de l'objet "exprs"

Imprimer toutes les valeurs.

```{r results = TRUE}
print(exprs)
```

## Affichage des premières lignes de l'objet

```{r results = TRUE}
head(exprs)
```

## Un peu plus de lignes

```{r results = TRUE}
head(exprs, n = 15)
```

## Caractéristiques d'un tableau

Dimensions

```{r}
dim(exprs)    ## Dimensions
ncol(exprs)   ## Nombre de colonnes
nrow(exprs)   ## Nombre de lignes
```

Noms des lignes et colonnes

```{r}
colnames(exprs)
rownames(exprs)
```


## Résumé rapide des données par colonne

```{r results = TRUE}
  summary(exprs)
```

## Sélection de colonnes d'un tableau

Valeurs stockées dans la colonne nommée "WT1"

```{r}
exprs$WT1
```

Notation alternative

```{r}
exprs[, "WT1"]  ## Sélection de la colonne WT1
```

Sélection de plusieurs colonnes. 

```{r}
exprs[, c("WT1", "WT2")]
```

Sélection de colonnes par leur indice

```{r}
exprs[, 2]
exprs[, c(2, 3)]
```

## Histogramme des valeurs d'expression pour WT1

```{r fig.width=7, fig.height=5, out.width="60%"}
hist(exprs$WT1)
```

## Histogramme du logarithme de ces valeurs

```{r fig.width=7, fig.height=5, out.width="60%"}
hist(log(exprs$WT1))
```

## Nuages de points -- Expressions KO1 vs WT1

```{r fig.width=7, fig.height=5, out.width="60%"}
plot(x = log(exprs$WT1), y = log(exprs$KO1))
```

## Personnalisation des paramètres graphiques

```{r fig.width=7, fig.height=5, out.width="60%"}
plot(x = log(exprs$WT1),   ## données pour l’abscisse
     y = log(exprs$KO1),   ## données pour l’ordonnée
     main = "Expression KO1 vs WT1",  ## Titre principal
     xlab = "WT1",  ## légende de l’axe X
     ylab = "KO1",  ## légende de l’axe Y
     pch = 16,      ## caractère pour marquer les points
     las = 1,       ## écrire les échelles horizontalement
     col = "red")   ## couleur des points
grid()  ## Ajout d’une grille
abline(a = 0, b = 1)   ## Ajouter la droite X = Y (intercept a = 0, pente b = 1).
```

## Sélection de lignes d'un tableau

Sélection des lignes 4 et 11 du tableau des expressions

```{r}
exprs[c(4, 11), ]
```

Indices des lignes correspondant aux IDs ENSG00000253991 et ENSG00000099958

```{r}
which(exprs$id %in% c("ENSG00000253991", "ENSG00000099958"))
```

Afficher les lignes correspondantes

```{r}
exprs[which(exprs$id %in% c("ENSG00000253991", "ENSG00000099958")), ]
```

## Calculs sur des colonnes

Calcul de moyennes par ligne (`rowMeans`) pour un sous-ensemble donné des colonnes (WT1 et WT2).

```{r}
rowMeans(exprs[,c("WT1","WT2")])
```

Ajout de colonnes avec les expressions moyennes des WT et des KO.
 
```{r}
exprs$meanWT <- rowMeans(exprs[,c("WT1","WT2")])
exprs$meanKO <- rowMeans(exprs[,c("KO1","KO2")])

head(exprs) ## Check the result
```

Fold-change KO vs WT

```{r}
exprs$FC <- exprs$meanKO / exprs$meanWT
head(exprs) ## Check the result
```

Moyenne de tous les échantillons

```{r}
exprs$mean <- rowMeans(exprs[,c("WT1","WT2","KO1","KO2")])
```

## MA-plot: log2FC vs intensité

$M$ est le logarithme en base 2 du rapport d'expression.

$$M = log_{2} (\text{FC}) = log_{2} \left( \frac{\text{KO}}{\text{WT}} \right) = log_2 (\text{KO}) - log_2(\text{WT})$$

```{r}
exprs$M <- log2(exprs$FC)
```

$A$ (average intensity) est la moyenne des logarithmes des valeurs d'expression.

$$A = \frac{1}{2} log_2 (\text{KO} \cdot \text{WT}) = \frac{1}{2} \left( log_2 (\text{KO}) + log_2(\text{WT}) \right)$$
```{r}
exprs$A <- rowMeans(log2(exprs[,c("meanWT", "meanKO")]))
```



## MA-plot: log2FC vs intensité

```{r  fig.width=7, fig.height=5, out.width="60%"}
plot(x = exprs$A, y = exprs$M, main = "MA plot", las = 1,
     col = "blue", pch = 16, xlab = "A = intensity", ylab = "M = log2FC")
grid(lty = "solid", col = "lightgray")
abline(h = 0)
```

## Charger les annotations des gènes

```{r}
annot <- read.table(file = "annotation.csv", header = TRUE, sep = ";")
dim(annot)   ## Vérifier les dimensions
head(annot)  ## Afficher quelques lignes
```

Combien de gènes par chromosome ?

```{r}
table(annot$chr)
```

Question: combien de gènes sur le chromosome 8 ? Et sur le X ?


## Diagramme en bâtons -- gènes par chromosomes

```{r fig.width = 5, fig.height = 7, out.width = "40%"}
barplot(sort(table(annot$chr)), horiz = TRUE, las = 1,
        main = "Genes per chromosome", ylab = "Chromosome", 
        col = "lightblue", xlab = "Number of genes")

```


## Sélectionner les données du chromosome 8

1ere étape: fusionner les deux tableaux exprs et annot

```{r}
exprs.annot <- merge(exprs, annot, by = "id")
head(exprs.annot)
```

2eme étape: sous-ensemble des lignes pour lesquelles chr vaut 8

```{r}
exprs8 <- exprs.annot[which(exprs.annot$chr == 8),]
print(exprs8)

```

## Exporter exprs8 dans un fichier

```{r}
write.table(x = exprs8, file = "exprs8.txt", sep = "\t",
            row.names = TRUE, col.names = NA)

```

## A ajouter: ALLER PLUS LOIN (optionnel)

- charger un tableau complet de données (RNA-seq ou ChIP-seq et refaire quelques plots)
- hist(breaks = 100)
- apply
- sort, order
- lm() ?

