```{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"
url.main <- "http://jvanheld.github.io/stat1"
dir.current <- file.path(dir.main, "practicals", "02_yeast_annotations")
setwd(dir.current)

```


## But de ce TP

Durant ce TP, vous serez amenés à effectuer les tâches suivantes: 

1. Manipuler une table de données génomique (les annotations du génome de la levure).
2. Sélectionner un sous-ensemble des données en filtrant les lignes sur base d'un critère déterminé (type d'annotation, chromosome).
3. Générer des graphiques pour représenter différents aspects liés à ces données.
4. Calculer les estimateurs de tendance centrale et dispersion.
5. Calculer un intervalle de confiance autour de la moyenne.


## Rendu

A la fin du TP, vous déposerez deux fichiers sur Ametice. 

1. Votre **code R**.
2. Un **rapport synthétique** qui inclura une présentation des principaux résultats (figures, statistiques descriptives) et votre interprétation. 

### Attendus pour le code

1. Le code doit être **lisible et compréhensible**: donnez à vos variables des noms indiquant explicitement ce qu'elles contiennent.

2. Le code devra être **correctement documenté** (le symbole `#` en début ou en milieu de ligne indique que le reste de cette ligne est un commentaire).

- avant chaque bloc de code, expliquer ce que vous comptez faire, à quoi sert ce bloc de code; 

- si c'est utile, ajoutez quelques mots de commentaires pour justifier l'approche choisi;

- chaque fois que vous déinifissez une variable, ajoutez sur la même ligne un commentaire indiquant ce que cette variable représente. 

3. Le code doit être **tansportable**: après l'avoir téléchargé, on doit pouvoir l'exécuter sur une autre machine. Je testerai systématiquement si les fichiers de code peuvent être exécutés sur ma machine. Evitez donc tout recours à des chemins absolus (nous indiquons ci-dessous comment définir des chemins relatifs par rapport à la racine de votre compte). 

### Attendus pour le rapport d'interprétation

Le rapport doit être synthétique (1 page de texte maximum + autant de figures et tables que vous le désirez). 

Chaque question doit être exprimée explicitement avant de présenter les résultats qui y répondent et de fournir l'interprétation de ces résultats. 

Chaque figure ou table doit être documentée par une légende permettant à un lecteur naïf de comprendre ce qu'elle représente. L'interprétation des résultats affichés sur une figure ou table se trouvera dans le texte principal (avec une référence au numéro de figure ou table). 

## Exemple historique: génome de la levure

```{r load_yeast_cds_lengths, echo=FALSE}
## Read a GTF file
## Format specification: http://www.ensembl.org/info/website/upload/gff.html

## Download tab-delimited file with chromosome sizes (unless already there)
gtf.url <- "http://jvanheld.github.io/stat1/data/Saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.37.gtf.gz"
gtf.file <- file.path(dir.main, "data", "Saccharomyces_cerevisiae", "Saccharomyces_cerevisiae.R64-1-1.37.gtf.gz")

if (file.exists(gtf.file)) {
  message("GTF annotation file: ", gtf.file)
} else {
  download.file(gtf.URL, destfile = gtf.file)
}

feature.table <- read.delim(gtf.file, comment.char = "#", sep="\t", header=FALSE, row.names=NULL)
names(feature.table) <- c("seqname", "source", "feature", "start", "end", "score", "strand", "frame", "attribute")
# nrow(feature.table) ## Count feature number

## Select subset of features having "CDS" as "feature" attribute
cds <- subset(feature.table, feature=="CDS")
## nrow(cds) ## Count the number of "coding sequences", i.e. open reading frames annotated as cds

```

- 1992: publication du premier chromosome eucaryote complet, le 3ème chromosome de la levure. 
- 1996: publication du génome complet.

Sur base des gènes dU 3ème chromosome (échantillon) on peut estimer la taille moyenne d'un gène de levure. 

**Questions: ** 

(a) La moyenne d'échantillon (chromosome III) permettait-elle de prédire la moyenne de la population (génome complet) ? 

    Pour répondre à cette quesiton, nous imaginerons que nous sommes revenus en 1992, et utiliserons l'ensemble des gènes du chromosome III (considérés ici comme un échantillon du génome) pour estimer la taille moyenne des gènes pour l'ensemble du génome (la "population" de gènes").

(b) Cet échantillon peut-il être qualifié de "simple et indépendant" ?

## Analyse de la longueur des gènes de la levure du boulanger

```{r cds_length_histo, out.width="90%", fig.width=10, fig.height=4, fig.cap="Distribution of cds lengths for Saccharomyces cerevisiae. ", echo=FALSE}
par.ori <- par(no.readonly = TRUE)
par(mar = c(4.1,4.1,2.1,1.1))

## Add a column to the table with cds lengths
cds$length <- cds$end - cds$start +1
max.len <- max(cds$length)

## Select cds on the third chromosome
cds.III <- subset(cds, seqname == "III")
# View(cds.III)
par(mfrow=c(1, 2))

## Plot an histogram with cds lengths
hist(cds.III$length, 
     breaks=seq(from=0, to=max.len+100, by=100), 
     main="Chromosome III",
     xlab=NA, ylab="Number of cds", 
     col="#BBDDFF")

## Plot an histogram with cds lengths
hist(cds$length, 
     breaks=seq(from=0, to=max.len+100, by=100), 
     main="Full genome",
     xlab="cds length (bp)", ylab="Number of cds", 
     col="#BBFFDD")
par(mfrow=c(1,1))

par <- par.ori

```



## Tutoriel

Avant de passer aux exercices, nous vous montrons ici quelques éléments de base concernant la lecture, la manipulation et l'écriture des tableaux de données avec R.


### Le chemin de la maison (manuel)

Nous allons créer un dossier pour ce tuto, en partant de la racine de nore compte. 

Première possibilité (**rapide mais peu élégant**):  entrer (manuellement) le chemin de la racine de votre compte dans une variable

`dir.home <- /le/chemin/de/ma/maison`

- Avantage: rapide et pratique
- Désavantage: non portable, ne fonctionnera que sur votre ordinateur

### Le chemin de la maison (automatique)

Une solution plus générale: utiliser la commande **R** `Sys.getenv()`.

- Invoquée sans paramètre, cette commande liste toutes les variables d'environnement (votre configuration système). 
- On peut restreindre l'output à une variable d'environnement donnée, par exemple `Sys.getenv("HOME")` retourne le chemin de la racine de votre compte.

**Note:** écriture équivalente sous Linux: le symbole tilde `~` indique également le chemin de la rachine de votre compte. 

```{r get_home}
## Identify the home directory 
## by getting the environment variable HOME 
dir.home <- Sys.getenv("HOME")
print(dir.home)
```


### Création d'un dossier pour le TP

```{r create_practical_dir}
## Define a variable containing the path of the results for this tutorial
dir.tuto <- file.path(dir.home, "stat1", "TP2")

print(dir.tuto)

## Create the directory for this tutorial
dir.create(path = dir.tuto, showWarnings = FALSE, recursive = TRUE)

## Go to the tutorial directory
setwd(dir.tuto)

## List the files already present in the folder (if any)
list.files()

```

### Téléchargement du fichier GTF à partir d'EnsemblGenomes

**Astuce: ** avant de télécharger le fichier d'annotations (GTF) depuis EnsemblGenomes vers notre ordinateur, nous allons vérifier s'il est déjà présent (et dans ce cas on ne le re-télécharge pas). 

```{r download_GTF}
## Define the URL of the annotation file (GTF-formatted)
gtf.URL <- "ftp://ftp.ensemblgenomes.org/pub/release-37/fungi/gtf/saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.37.gtf.gz"

## Define the path where the local copy will be stored
local.GTF <- file.path(dir.tuto, "Saccharomyces_cerevisiae.R64-1-1.37.gtf.gz")

## If the local file file laready exists, skip the download
if (file.exists(local.GTF)) {
  message("GTF file already exists in the tutorial folder: ", local.GTF)
} else {
  ## Download annotation table in GTF format
  download.file(url = gtf.URL, destfile = local.GTF)
}
```

### Chargement d'un tableau de données

R comporte plusieurs types de structures tabulaires (matrix, data.frame, table). 

La structure la plus couramment utilisée est le `data.frame`, qui consiste en un tableau de valeurs (numériques ou chaînes de caractères) dont les lignes et les colonnes sont associées à des noms. 

La fonction `read.table()` permet de lire un fichier texte contenant un tableau de données, et de stocker le contenu dans une variable.

Plusieurs fonctions dérivées de `read.table()` facilitent la lecture de différents types de formats: 

- `read.delim()` pour les fichiers dont les colonnes sont délimitées par un caractère  particulier (généralement la tabulation, représentée par "\t").
- `read.csv()` pour les fichiers "comma-searated values".


1. Téléchargez le fichier suivant sur votre ordinateur: 

- [Saccharomyces_cerevisiae.R64-1-1.37.gtf](../../data/Saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.37.gtf)

2. Chargez-le au moyen de la fonction read.table (pour cela vous devez remplacer le chemin ci-dessous par celui de votre ordinateur). 

```{r read.table, echo=TRUE, result=TRUE}
## Read a GTF file with yeast genome annotations

## Load the feature table
feature.table <- read.table(
  local.GTF, 
  comment.char = "#", 
  sep="\t", 
  header=FALSE, 
  row.names=NULL)

## The bed format does not contain any column header, 
## so we set it manually based on the description of the format, 
## found here: 
##     http://www.ensembl.org/info/website/upload/gff.html
names(feature.table) <- c("seqname", "source", "feature", "start", "end", "score", "strand", "frame", "attribute")


```

### Exploration du contenu d'un tableau de données

La première chose à faire après avoir chargé un tableau de données est de vérifier ses dimensions

```{r echo=TRUE, result=TRUE, warning=FALSE}
dim(feature.table) ## Dimensions of the tbale
nrow(feature.table) ## Number of rows
ncol(feature.table) ## Number of columns
```

L'affichage du tableau d'annotations complet ne serait pas très lisible, puisqu'il comporte des dizaines de milliers de lignes. 

Nous pouvons afficher les premières lignes avec la fonction `head()`. 

**Note: ** la dernière colonne est particulièrement lourde (elle contient un tas d'information). Nous verrons plus loin comment sélectionner un sous-ensemble des colonnes pour alléger l'affichage.  

```{r echo=TRUE, eval=TRUE, result=TRUE}
## Display the 5 first rows of the feature table
head(feature.table, n = 5) 
```

La fonction `tail()` affiche les dernières lignes:

```{r echo=TRUE, eval=TRUE, result=TRUE}
## Display the 5 last rows of the feature table
tail(feature.table, n = 5) 
```

If you are using the **RStudio** environment, you can display the table in a dynamic viewer pane with the function  `View()`.

```{r echo=TRUE, eval=FALSE}
## In RStudio, display the table in a separate tab
View(feature.table) 
```


### Sélection de sous-ensembles d'un tableau

Sélection d'une ligne spécifiée par son indice.

```{r echo=TRUE, eval=TRUE, result=TRUE}
feature.table[12,]
```

Sélection d'une colonne spécifiée par son indice (affichage des premières valeurs seulement.

```{r echo=TRUE, eval=TRUE, result=TRUE}
head(feature.table[,3])
```

Selection d'une cellule en combinant indices de ligne et colonne.

```{r echo=TRUE, eval=TRUE, result=TRUE}
feature.table[12, 3]
```

Sélection d'un bloc de colonnes et/ou de lignes.

```{r echo=TRUE, eval=TRUE, result=TRUE}
feature.table[100:105, 1:6]
```
Sélection de colonnes "à la carte" (ici, les coordonnées génomiques de chaque "feature"):  chromosome, début, fin, brin. 

```{r echo=TRUE, eval=TRUE, result=TRUE}
feature.table[100:105, c(1,4,5,7)]
```


Sélectionner une colonne sur base de son nom.

```{r echo=TRUE, eval=TRUE, result=TRUE}
## Select the "start" column and print the 100 first results
head(feature.table$start, n=100)

## Print the 20 first values of the "feature" field, which indicates the feature type
head(feature.table$feature, n=20)
```

Sélection de plusieurs colonnes sur base de leurs noms. 


```{r echo=TRUE, eval=TRUE, result=TRUE}
## Select the "start" column and print the 100 first results
feature.table[100:106, c("seqname", "start", "end", "strand")]
```

**Note**: il est également possible de nommer les lignes d'un data.frame mais le tableau GTF ne se prête pas à cela. Nous verrons d'autres exemples ultérieurement. 



### Sélection d'un sous-ensemble de lignes sur base du contenu d'une colonne

La fonction `subset()` permet de sélectionner un sous-ensemble des lignes d'un data.frame sur base d'une condition appliquée à une ou plusieurs colonnes. 

Nous pouvons l'appliquer pour sélectionner le sous-ensemble des lignes du tableau d'annotations correspondant à des séquences codantes (CDS).

```{r echo=TRUE, eval=TRUE, result=TRUE}
## Select subset of features having "cds" as "feature" attribute
cds <- subset(feature.table, feature=="cds")

nrow(feature.table) ## Count the number of features
nrow(cds) ## Count the number of cds
```

### Décompte par valeur

La fonction `table()` permet de compter le nombre d'occurrences de chaque valeur dans un vecteur ou un tableau. Quelques exemples d'utilisation ci-dessous. 

```{r echo=TRUE, eval=TRUE, results=TRUE}
## Count the number of featues per chromosome
table(feature.table$seqname)

## Count the number of features per type
table(feature.table$feature)
```

On peut calculer des tables de contingence en comptant le nombre de combinaisons entre 2 vecteurs (ou 2 colonnes d'un tableau). 


```{r echo=TRUE, eval=TRUE, results=TRUE}
##  Table with two vectors
table(feature.table$feature, feature.table$seqname)

## Same result with a 2-column data frame
table(feature.table[, c("feature", "seqname")])
```

## Exercices

### 1. Spécifications du format GTF

Lisez les spécifications du format GTF.

- Ensembl (<http://www.ensembl.org/info/website/upload/gff.html>)
- UCSC (<https://genome.ucsc.edu/FAQ/FAQformat.html#format4>)

### 2. Création d'un dossier local pour le TP

Créez un dossier local (par exemple: `stat1/TP_levure` à partir de la racine de votre compte). Nous vous suggérons d'utiliser les fonctions suivantes: 

- `Sys.getenv("HOME")` (Linnux et Mac OS X), pour obtenir la racine de votre compte utilisateur;
- `file.path()` pour construire un chemin;
- `dir.create()` pour créer le dossier de ce TP. Lisez attendivement les options de cette fonction avec `help(dir.create)` 

### 3. Localisation du fichier d'annotations

Localisez le fichier d'annotations du génome de la levure en format GTF dans ce dossier local. 

- Site Ensembl Fungi: <http://fungi.ensembl.org/>
- Cliquez "Downloads" pour accéder au [site ftp](http://fungi.ensembl.org/info/website/ftp/)
- Dans la bo^îte de recherche, tapez *"saccharomyces cerevisiae"* et suivez le lien "[GTF](ftp://ftp.ensemblgenomes.org/pub/release-37/fungi/gtf/saccharomyces_cerevisiae)"
- COpiez l'adresse (URL) du ichier [Saccharomyces_cerevisiae.R64-1-1.37.gtf.gz](ftp://ftp.ensemblgenomes.org/pub/release-37/fungi/gtf/saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.37.gtf.gz)


### 4. Téléchargement d'un fichier à partir d'un site ftp

Fonctions suggérées: 

- `download.file()` (lisez l'aide pour conna^tre les arguments)

### 5. Chargement d'une table de données en R

Ecrivez un script qui charge la table de données dans une variable nommée `feature.table`, en utilisant la fonction R `read.delim()`. 

Veillez à ignorer les lignes de commentaires (qui commencent par un caractère `#`). 


### 6. Calcul de la longueur des gènes codants


- Ajoutez à la table d'annotations (`feature.table`) une colonne intitulée "length" qui indique la longueur de chaque élément génomique annoté. 

```{r}
## Add a colmn with feature lengths
feature.table[, "length"] <- feature.table[, "end"] - feature.table[, "start"] + 1

## Add a colmn with feature lengths: equivalent result with simpler notation
feature.table$length <- feature.table$end - feature.table$start + 1
```


- Comptez le nombre de lignes de la table correspondant à chaque type d'annotation (3ème colonne du GTF, "feature"). 

    - fonction `table()`

```{r}
~table(feature.table$feature)
```


- Sélectionnez les lignes correspondant à des régions codantes ("CDS") 

    - fonction `subset()`

```{r}
cds <- subset(feature.table, feature=="CDS")
```


- Comptez le nombre de CDS par chromosome.

    - fonction `table()`


```{r}
table(cds$seqname)
```

- Chargez la table des tailles de chromosomes [chrom_sizes.tsv](../../data/Saccharomyces_cerevisiae/chrom_sizes.tsv), et calculez la densité de gènes pour chaque chromosome (nombre de gènes par Mb). 

```{r}
## Download tab-delimited file with chromosome sizes (unless already there)
annot.url <- "http://jvanheld.github.io/stat1/data/Saccharomyces_cerevisiae/chrom_sizes.tsv"
chrom.size.file <- file.path(dir.tuto, "chrom_sizes.tsv")

if (!file.exists(chrom.size.file)) {
    download.file(annot.url, destfile = chrom.size.file)
}

## Read chromosome sizes
chrom.size <- read.delim(
  file = chrom.size.file,
  header = FALSE, row.names = 1)

## Assign a name to the columns
names(chrom.size) <- c("chromID", "size")
# View(chrom.size)

## print the size of hte third chromosome
chrom.size["III", "size"]

```


### 6. Histogramme de la longueur des gènes

Au moyen de la fonction `hist()`, dessinez un histogramme représentant la distribution de longueur des CDS. 

```{r}
hist(cds$length)
```

Choisissez les intervalles de classe de façon à ce que l'histogramme soit informaatif (ni trop ni trop peu de classes).

```{r}
## Take more or less 100 bins
h <- hist(cds$length, breaks=100)


```

Récupérez le résultat de `hist()` dans une variable nommée `cds.length.hist`.

```{r}
## Define breaks exactly in the way you wish
cds.length.hist <- hist(cds$length, breaks=seq(from=0, to=max(cds$length)+100, by=100))

```


Imprimez le résultat à l'écran (`print()`) et analysez la structure de la variable `cds.length.hist` (il s'agit d'une variable de type liste). 

Fonctions utiles: 
```{r}
## Display the values used to draw the histogram
print(cds.length.hist )

```

- `class(cds.length.hist)`
- `attributes(cds.length.hist)`


D'autres types de graphiques permettent d'explorer la distribution d'un ensemble des données. En particulier, les boîtes à moustaches (box plots) affichent, pour une série de données, la médiane, l'écart interquartile, un intervalle de confiance et les valeurs aberrantes. 

```{r fig.width=7, fig.height=5, fig.cap="Boîte à moustache indiquant la distribution de longueur des gènes par chromosome. "}
boxplot(length ~ seqname, data = cds, col="palegreen", horizontal=TRUE, las=1, xlab="Gene length", ylab="Chromosome")

```

### 7. Paramètres descriptifs

Calculez les paramètres de tendance centrale (moyenne, médiane, mode) et de dispersion (variance, écart-type, écart inter-quartile) 

- pour les gènes du chromosome III;
- pour l'ensemble des gènes de la levure.


```{r}
x <- subset(cds, seqname == "III", select = "length")
dim(x) 
class(x) ## Ah ah, this is not a vector but a data.frame

## Convert the data frame into a vector
x <- unlist(x)
class(x)
head(x)

## Compute the mean, either manually or with the ad hoc R function
n <- length(x) 
print(paste("Chromosome III contains", n, "CDS"))
message("Chromosome III contains ", n, " CDS")

m <- mean(x)
print(m)
message("mean(x) = ", round(m, digits = 1))

## Compute the mean manually to compare the result
m.recalc <- sum(x)/n
message("Manually computed sample mean: ", round(digits=1, m.recalc))

## Compute manually standard dev of the sample
sample.var <- sum((x - m)^2)/ n
sample.sd <- sqrt(sample.var)
message("Sample standard dev =", round(digits=1, sample.sd))

## Compute an estimate of the population standard deviation
pop.sd.est <- sqrt(sum((x - m)^2) / (n-1))
message("Sample-based estimate of population standard dev =", round(digits=1, pop.sd.est))

## Compute the standard deviation with R function sd()
R.sd <- sd(x)
message("Result of R sd() function =", round(digits=1, R.sd))


```

**Ah ah! (skeptical tone)** The R function `sd()` does **not** compute the standard deviation of the input numbers ($s$), but the **estimate of the standard deviaiton of the population** ($\hat{\sigma}$)

Affichez ces paramètres sur l'histogramme de la longueur des gènes, en utilisant la fonction `arrows()`

### 8. Intervalle de confiance

A partir des gènes du chromosome III (considérés comme l'échantillon disponible en 1992), calculez un intervalle de confiance autour de la moyenne, et formulez l'interprétation de cet intervalle de confiance. 
Evaluez ensuite si cet intervalle de confiance recouvrait ou non la moyenne de la population (tous les gènes du génome de la levure, qui devint disponible 4 ans après le chromosome III). 

$$ \bar{x} \pm \frac{\hat{\sigma}}{\sqrt(n)} \cdot t_{1-\alpha/2}^{n-1}$$

```{r}
## Define alpha, the risk
alpha <- 0.05

## Let us get the critical value for the t distribution
help("TDist")

## Which value corresponds to alpha/2 

## Beware ! by default the qt() function return the lower tail
qt(p = alpha/2, df =  n - 1)

## For confidence intervals we need a positive t value, we thus take the upper tail
t.value <- qt(p = alpha/2, df =  n - 1, lower.tail = FALSE)

IC.min <- m - pop.sd.est * t.value /sqrt(n)
IC.max <- m + pop.sd.est * t.value /sqrt(n)

message("Confidence interval: [", 
        round(digits=1, IC.min), 
        ", ", 
        round(digits=1, IC.max), "]")

```

Dessinez un polygnone des fréquences indiquant le nombre de gènes par classe (milieux de classe). 

```{r}

## Draw a frequency polygon
plot(cds.length.hist $mids, cds.length.hist $counts, 
     main="Frequency polygon", 
     xlab="Gene length", ylab="Number of genes", 
     type="l", col="darkgreen", lwd=2, xlim=c(0, 5000))
grid()

arrows(x0 = IC.min, y0 = 100, x1 = IC.max, y1=100, length = 0, angle = 00, col="violet", lwd=3)
abline(v = m, col="violet")

pop.mean <- mean(cds$length)
abline(v = pop.mean, col="blue")

legend("topright", legend = c("genome", "chr3"), col = c("blue", "violet"), lwd=1)
```


### 9. Distribution de la longueur des gènes

- A partir du résultt de `hist()`, récupérez un tableau (dans une variable de type `data.frame`) indiquant les fréquences absolues (`count`) en fonction de la taille médiane des classes (`mids`),

- Ajoutez à ce tableau une colonne indiquant la fréquence relative de chaque classe de longueurs de gènes.

- Ajoutez à ce tableau des colonnes indiquant la **fonction de répartition empirique** des longueurs de gènes (nombre de gènes d'une taille inférieure ou égale à chaque valeur $x$ observée, et fréquence relative de ce nombre). 

    - fonction de base: `cumsum()`
    - fonction avancée;`ecdf()`
    
- Au moyen des fonctions `plot()` et `lines()`, dessinez un graphe représentant la fréquence absolue par classe (médianes de classes en $X$, comptages en $Y$), et la fonction de répartition empirique. 

    - suggestion: superposez les  utilisez le type de lignes "h" pour les fréquences de classe, et "l" ou "s" pour la fonction de répartition. 


### 10. Distribution attendue au hasard pour la longueur des gènes

Sur base de la taille du génome (12.156.679 bp) et des fréquences génomiques de codons définies ci-dessous, calculez la distribution de longueurs de gènes attendue au hasard, et ajoutez-là au graphique. 

Vous pouvez télécharger les fréquences génomiques de tous les  trinucléotides ici:  [3nt_genomic_Saccharomyces_cerevisiae-ovlp-1str.tab](http://jvanheld.github.io/stat1/data/Saccharomyces_cerevisiae/oligo_freq/3nt_genomic_Saccharomyces_cerevisiae-ovlp-1str.tab)

Alternative: créez une variable `freq.3nt` et assignez-y manuellement les valeurs pour les 4 nucléotides nécessaires, à partir de la table ci-dessous.


```{r echo=FALSE}
## Download tab-delimited file with chromosome sizes (unless already there)
freq.3nt.url <- "http://jvanheld.github.io/stat1/data/Saccharomyces_cerevisiae/oligo_freq/3nt_genomic_Saccharomyces_cerevisiae-ovlp-1str.tab"
freq.3nt.file <- file.path(dir.tuto, "3nt_genomic_Saccharomyces_cerevisiae-ovlp-1str.tab")

if (file.exists(freq.3nt.file)) {
  message("Trinucleotide frequency file: ", freq.3nt.file)
} else {
    download.file(freq.3nt.url, destfile = freq.3nt.file)
}


## Download oligonucleotide frequencies
oligo.freq <- read.table(freq.3nt.file, sep="\t", row.names = NULL)
colnames(oligo.freq) <- c("sequence", "frequency", "occurrences")
oligo.freq$frequency <- signif(digits=3, oligo.freq$frequency)

row.names(oligo.freq) <- oligo.freq$sequence

#oligo.freq["...",] <- c("...", "...", "...")

selected.oligos <- c("AAA", "ATG", "TAA", "TAG", "TGA")

kable(data.frame(oligo.freq[selected.oligos, ]), row.names = FALSE)
#print(data.frame(oligo.freq[selected.oligos, ]))
```


### 11. 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()
```

