```{r setup, include=FALSE, size="huge"}
library(knitr)
## Default parameters for displaying the slides
knitr::opts_chunk$set(
  echo = TRUE, 
  eval = TRUE, 
  fig.width = 7, 
  fig.height = 5, 
  fig.align = "center", 
  fig.path = "figures/",
  size = "tiny", 
  warning = FALSE, 
  results = TRUE, 
  message = FALSE, 
  comment = "")

dir.main <- "~/stat1"
dir.current <- file.path(dir.main, "practicals", "tests_multiples_kmers")
setwd(dir.current)

```


## But de ce tutoriel

- Apprendre à manipuler des tableaux de données (objet R data.frame)
- Travailler sur les colonnes de tableaux de données
- Se familiariser par la pratique avec les différentes corrections de tests multples. 
- Représentations graphiques associées aux tests multiples (histogramme des P-valeurs, MA plot, Volcano plot)

## Jeu de données : contexte

Nous avons utilisé le logiciel *oligo-analysis* de la suite Regulatory Sequence Analysis Tools (RSAT; <http://rsat.eu/>) pour compter les occurrences de chaque octanucléotide dans les pics de liaison du facteur transcriptionnel Sox2 détectés par une expérience de ChIP-seq.

### Conditions d'analyse

- Nombre de séquences: $S = 8014$
- Taille totale des séquences: $L = 2.776.509$

- On dénombre uniquement les occurrences renouvelantes (non-chevauchantes). 

- Les octanucléotides ont été regroupés par paires réverse complémentaires (analyse insensible au brin). 

- Le programme a également calculé la probabilité a priori (colonne `prior` du tableau) de chaque octanucléotide sur base d'un modèle de Markov d'ordre 5. Autrement dit, la probabilité a priori d'un octanucléotide ($p$) est estimée sur bases des probabilités d'hexanucléotides observées dans l'ensemble des séquences. 


## Questions à traiter

Durant ce tutoriel, nous allons explorer ce jeu de données, et calculer une série de statistiques pour détecter les octanucléotides (abbréviation: **8nt**) significativement sur-représentés.


## Téléchargement du tableau de données


1. A la racine de votre compte, créez un dossier de travail nommé `TP_stat1_multiple_tests_kmers`.

```{r workdir}

## Create work directory
work.dir <- '~/TP_stat1_multiple_tests_kmers'
dir.create(work.dir, recursive = TRUE, showWarnings = FALSE)
setwd(work.dir)

```


2. Téléchargez dans ce dossier le fichier de comptage des k-mères, qui se trouve sur le site Wbe du cours. 

<https://jvanheld.github.io/stat1/data/Sox2_oligos/Sox2-peaks_8nt_mkv5.tsv.gz>


```{r download_data}
kmer.file <- 'Sox2-peaks_8nt_mkv5.tsv.gz'
kmer.url <- file.path("https://jvanheld.github.io/stat1/data/Sox2_oligos/", kmer.file)


## Download the kmer occurrence file in the work directory
# download.file(url = kmer.url, destfile = kmer.file)
```

## Chargement d'un fichier tabulaire de données

Nous allons charger le tableau de comptage des k-mères dans une variable nomée `kmer.table`, de type `data.frame`.

Les data frames sont des tableaux à deux dimensions (lignes/colonnes) qui permettent d'associer 
- un nom à chaque ligne, 
- un nom à chaque colonne
- un type de données (nombres, chaînes de caractères, ...) est défini automatiquement pour chaque colonne sur base de son contenu

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

## Check the files
list.files()

## Load the file
kmer.table <- read.delim(
  sep = "\t",
  file = kmer.file, 
  header = 1, 
  row.names = 1)

## Check the content of the original data table
names(kmer.table) ## Column content
head(kmer.table)  ## First rows
tail(kmer.table) ## Last rows
dim(kmer.table)   ## Dimensions of the data frame


```

## Sélection de lignes et de colonnes

Il existe plusieurs façons de sélectionner des lignes et colonnes d'un tableau. 


```{r show_columns, echo=TRUE, eval=TRUE}
# Afficher le contenu de la cellule en 10ème ligne, 2ème colonne
kmer.table[10,2] 

## Afficher la 10ème ligne
kmer.table[10,] 

## Afficher la 2ème colonne
## Pour éviter d'afficher des dizaines de milliers de valeurs, nous affichonss eulement les 200 premières valeurs en englobant le résultat dans la fonction head()
head(kmer.table[,2] , n = 200)

## Afficher les lignes 10 à 15
kmer.table[10:15,] 

## Afficher les colonnes 1 et 2 pour les lignes 10 à 20
kmer.table[10:15, 1:2] 

## Sélection d'une ligne par son nom
kmer.table["gattacca|tggtaatc",]
kmer.table[c("gattacca|tggtaatc", "aaaaaaaa|tttttttt"),]

## Sélection d'une colonne par son nom
## (ici, 50 premières lignes seulement )
head(n = 50, kmer.table[, "occurrences"])

## Notation rapide
head(n = 50, kmer.table$occurrences)

## Notation équivalente
kmer.table$occurrences[1:50]

```



### Exercice

1. Utilisez la fonction `plot()` pour dessiner un nuage de points indiquant le nombre d'occurrences ($Y$) en fonction de la probabilité a priori de chaque 8nt ($X$). Soignez le graphique comme au TP précédent. 

```{r occ_prior_scatter, eval=FALSE, echo=FALSE}
plot(kmer.table$prior, kmer.table$occurrences,
     xlab = "Prior k-mer probability",
     ylab = "Observed occurrences",
     col = "darkgrey", 
     las = 1)
grid(lty = "solid")
```

2. Dessinez un histogramme des occurrences au moyen de la fonction `hist()`.

## Calcul de statistiques sules colonnes

Nous allons calculer le nombre attendu d'occurrences pour chaque oligonucléotide, en tenant compte de sa probabilité a priori. 

```{r exp_occ}
kmer.sum <- sum(kmer.table$occurrences)

## Check the sum of prior probabilities 
## (~1 but there are some rounding errors)
sum(kmer.table$prior)


kmer.table$exp.occ <- kmer.table$prior * kmer.sum

plot(kmer.table$exp.occ, kmer.table$occurrences, las = 1,
     col = densCols(kmer.table$exp.occ, kmer.table$occurrences),
     xlab = "Expected occurrences", 
     ylab = "Observed occurrences")
grid(lty = "solid", col = "gray")
abline(a = 0, b = 1)
```


## Exercice

Calculez les statistiques suivantes pour l'ensemble des oligos. 

A chaque étape, ajoutez à la table `kmer.table` une colonne avec la nouvelle statistique calculée. Ceci permettra d'exporter ensuite l'ensemble des résultats.  

Dessinez pour chaque statistique calculée un graphique qui représente le mieux possible les résultats, en faisant ressortir les points pertinents. 

1. Rapport entre occurrences observées et attendues (obs/exp)
2. Le logarithme en base 2 de ce rapport (log-ratio)
3. $P(X= x)$: probabilité binomiale d'observer exactement $X$ occurrences 
4. P-valeur de sur-représentation: $P(X \ge x)$ au moyen de la binomiale.
5. P-valeur de sur-représentation: $P(X \ge x)$ au moyen de la Poisson. 

6. Dessinez un nuage de points comparant les P-valeurs binomiale et Poisson. Utilisez des axes logarithmiques pour mettre en évidence les faibles p-valeurs. 

7. Calculez la E-valeur binomiale.

8. Calculez la significativité de cette E-valeur, définie ici comme $sig_{E} = -log_{10}(E)$. 

9. Dessinez un graphique avec la significativité ($Y$) en fonction du log-ratio ($X$). 

10. Dessinez un histogramme des p-valeurs avec des intervalles de 5%. 

11. Générez une colonne de nombre d'occurrences aléatoires tirés selon une loi de Poisson dont l'espérance est prise dans la colonne `prior`. 

12. Refaites les mêmes calculs et dessins (log-ratio, significativité, volcano plot, histogramme des p-valeurs) avec ces occurrences aléatoires. 

13. Sauvegardez le tableau de résultat dans un fichier nommé `Sox2-peaks_8nt_mkv5_overstats.tsv`.

## Solutions


## 1. Rapport entre occurrences observées et attendues (obs/exp)

```{r exp-occ_ratio, echo=TRUE}
## Observed / expected ratio
kmer.table$ratio <- kmer.table$occurrences / kmer.table$exp.occ

## Plot ratio as a function of the expected occurrences
plot(kmer.table$exp.occ, kmer.table$ratio, 
     col = densCols(kmer.table$exp.occ, kmer.table$ratio),
     xlab = "expected occurrences",
     ylab = "obs/exp ratio"
     )
grid(lty = "solid")
abline(h = 1)

```

Méfiez-vous du rapport obs/exp. 

## 2. Le logarithme en base 2 de ce rapport (log-ratio)

```{r log-ratio, echo=TRUE}
## log-ratio
kmer.table$log.ratio <- log(kmer.table$ratio, base = 2)

## Plot log-ratio as a function of the expected occurrences
plot(kmer.table$exp.occ, kmer.table$log.ratio, 
     col = densCols(kmer.table$exp.occ, kmer.table$log.ratio),
     xlab = "expected occurrences",
     ylab = "log2(obs/exp ratio)"
     )
grid(lty = "solid")
abline(h = 0)

```


## 3. $P(X= x)$: probabilité binomiale d'observer exactement $X$ occurrences 

```{r dbinom, fig.width=10, fig.height=4, out.width = "90%"}
## Number of trials is the sum of all k-mer occurrences
N <- sum(kmer.table$occurrences)

## Binomial PMF
kmer.table$dbinom <- dbinom(x = kmer.table$occurrences, size = N, prob = kmer.table$prior)

par(mfrow = c(1,2))
plot(kmer.table$exp.occ, kmer.table$dbinom, 
     col = densCols(kmer.table$exp.occ, kmer.table$dbinom),
     xlab = "Expected occurrences",
     ylab = "Binomial density"
     )
grid(lty = "solid")
abline(h = c(0,1))

## log Y axis
plot(kmer.table$exp.occ, kmer.table$dbinom, 
     col = densCols(kmer.table$exp.occ, kmer.table$dbinom),
     log = "y",
     xlab = "Expected occurrences",
     ylab = "Binomial density"
     )
grid(lty = "solid")
abline(h = 1)


```


## 4. P-valeur de sur-représentation: $P(X \ge x)$ au moyen de la binomiale.

```{r Pbinom, fig.width=10, fig.height=4, out.width = "90%"}
## Binomial P-value
kmer.table$Pbinom <- pbinom(
  q = kmer.table$occurrences - 1, 
  size = N, 
  prob = kmer.table$prior, 
  lower.tail = FALSE)

par(mfrow = c(1,2))
plot(kmer.table$exp.occ, kmer.table$Pbinom, 
     col = densCols(kmer.table$exp.occ, kmer.table$Pbinom),
     xlab = "Expected occurrences",
     ylab = "Binomial P-value"
     )
grid(lty = "solid")
abline(h = c(0,1))

## log Y axis
plot(kmer.table$exp.occ, kmer.table$Pbinom, 
     col = densCols(kmer.table$exp.occ, kmer.table$Pbinom),
     log = "y",
     xlab = "Expected occurrences",
     ylab = "Binomial P-value"
     )
grid(lty = "solid")
abline(h = 1)

par(mfrow = c(1,1))

```

## 5. P-valeur de sur-représentation: $P(X \ge x)$ au moyen de la Poisson. 

```{r Ppois, fig.width=10, fig.height=4, out.width = "90%"}
## Binomial P-value
kmer.table$Ppois <- ppois(
  q = kmer.table$occurrences - 1, 
  lambda = kmer.table$exp.occ, 
  lower.tail = FALSE)


par(mfrow = c(1,2))
plot(kmer.table$exp.occ, kmer.table$Ppois, 
     col = densCols(kmer.table$exp.occ, kmer.table$Pois),
     xlab = "Expected occurrences",
     ylab = "Poisson P-value"
     )
grid(lty = "solid")
abline(h = c(0,1))

## log Y axis
plot(kmer.table$exp.occ, kmer.table$Ppois, 
     col = densCols(kmer.table$exp.occ, kmer.table$Pois),
     log = "y",
     xlab = "Expected occurrences",
     ylab = "Poisson P-value"
     )
grid(lty = "solid")
abline(h = 1)
par(mfrow = c(1,1))

```

## 6. Dessinez un nuage de points comparant les P-valeurs binomiale et Poisson. Utilisez des axes logarithmiques pour mettre en évidence les faibles p-valeurs. 

```{r Pbinom_vs_Ppois}
plot(kmer.table$Ppois, kmer.table$Pbinom,
     col = densCols(kmer.table$Ppois, kmer.table$Pbinom),
     xlab = "Poisson P-value", ylab = "Binomial P-value",
     log = "xy",
     las = 1)
grid(lty = "solid")
abline(a = 0, b = 1)
abline(h = 1)
abline(v = 1)

```


## 7. Calculez la E-valeur binomiale.

```{r Ebinom, fig.width=10, fig.height=4, out.width = "90%"}
nb.tests <- nrow(kmer.table) # Number of tests
kmer.table$Ebinom <- kmer.table$Pbinom * nb.tests


par(mfrow = c(1,2))
plot(kmer.table$exp.occ, kmer.table$Ebinom, 
     col = densCols(kmer.table$exp.occ, kmer.table$Ebinom),
     las = 1,
     xlab = "Expected occurrences",
     ylab = "Binomial E-value"
     )
grid(lty = "solid")
abline(h = c(0, N))
abline(h = 0, col = "brown")

## log Y axis
plot(kmer.table$exp.occ, kmer.table$Ebinom, 
     col = densCols(kmer.table$exp.occ, kmer.table$Ebinom),
     las = 1,
     log = "y",
     xlab = "Expected occurrences",
     ylab = "Binomial E-value"
     )
grid(lty = "solid")
abline(h = c(0, N))
abline(h = 0, col = "brown")
par(mfrow = c(1,1))

```


## 8. Calculez la significativité de cette E-valeur, définie ici comme $sig_{E} = -log_{10}(E)$. 

```{r binom_sig}
kmer.table$sig <- -log(kmer.table$Ebinom, base = 10)

plot(kmer.table$exp.occ, kmer.table$sig, 
     col = densCols(kmer.table$exp.occ, kmer.table$sig),
     las = 1,
     xlab = "Expected occurrences",
     ylab = "Binomial significance"
     )
grid(lty = "solid")
abline(h = 0, col = "brown")

```


## 9. Dessinez un graphique avec la significativité ($Y$) en fonction du log-ratio ($X$). 

```{r volcano_plot}
plot(kmer.table$log.ratio, kmer.table$sig,
     col = densCols(kmer.table$log.ratio, kmer.table$sig),
     las = 1,
     panel.first = grid(lty="solid"))
abline(v=0)
abline(h=0, col="brown")
abline(h=2, col="red")

```


## 10. Dessinez un histogramme des p-valeurs avec des intervalles de 5%. 

```{r pval_hist_kmers}
hist(kmer.table$Pbinom, 
     breaks=seq(from = 0, to = 1, by = 0.05), las = 1,
     xlab = "Binomial P-value", ylab = "Number of k-mers",
     col = "lightgray",
     main = "P-value histogram")
```


## 11. Générez une colonne de nombre d'occurrences aléatoires tirés selon une loi de Poisson dont l'espérance est prise dans la colonne `prior`. 

```{r}
kmer.table$rand.occ <- rpois(n = nrow(kmer.table), lambda = kmer.table$exp.occ)
```


##  12. Refaites les calculs et dessins ci-dessus avec ces occurrences aléatoires. 


```{r pval_hist_rand-occ}
kmer.table$rand.occ <- rpois(n = nrow(kmer.table), lambda = kmer.table$exp.occ)

kmer.table$rand.Pbinom <- pbinom(q = kmer.table$rand.occ -1, size = N, prob = kmer.table$prior, lower.tail = FALSE)

kmer.table$rand.Ppois <- ppois(q = kmer.table$rand.occ -1, lambda = kmer.table$exp.occ, lower.tail = FALSE)


hist(kmer.table$rand.Ppois, 
     breaks = seq(from = 0, to = 1, by = 0.05), las = 1,
     xlab = "Binomial P-value", ylab = "Number of k-mers",
     col = "lightgray",
     main = "P-value histogram - Random Poisson")

```

## 13. Sauvegardez le tableau de résultat dans un fichier nommé `Sox2-peaks_8nt_mkv5_overstats.tsv`.

```{r write_table}
write.table(x = kmer.table, 
            file = "Sox2-peaks_8nt_mkv5_overstats.tsv", 
            sep = "\t", 
            row.names = TRUE, 
            col.names = NA, quote = FALSE
            )

```



## Avant de terminer : conservez la trace de votre session

La traçabilité constitue un enjeu essentiel en sciences. 
La fonction ***R***  `sessionInfo()` fournit un résumé des conditions d'une session de travail: version de R, système opérateur, bibliothèques de fonctions utilisées. 


```{r}
sessionInfo()
```

