report_Rotenone.Rmd 69 KB
Newer Older
Magali Monnoye's avatar
Magali Monnoye committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
---
title: Analyses 16s projet Rotenone Paula Perez Pardo
description: |
  Analyses 16S 
author:
  - name: Magali Monnoye
    affiliation: INRAE - MICALIS
date: "2021-11-16 (Last updated: `r Sys.Date()`)"
bibliography: resources/biblio.bib # don't change
csl: resources/biomed-central.csl # don't change
output:
  html_document:
    lib_dir: site_libs
    self_contained: false
    number_sections: FALSE
    code_folding: "hide"
    toc: true
    toc_depth: 5
    toc_float: true
    df_print: paged
    css: ['css/styles.css', 'https://use.fontawesome.com/releases/v5.0.9/css/all.css']
  pdf_document:
    toc: yes
  word_document:
    toc: yes  
always_allow_html: yes
    
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(eval=FALSE, echo =FALSE,  #echo=TRUE pour afficher le code dans rapport
                      cache = FALSE, message = FALSE, 
                      warning = FALSE, cache.lazy = FALSE,
                      fig.height = 5, fig.width = 10.5)
```


```{r datatable global options, echo=FALSE, eval=TRUE}
options(DT.options = list(pageLength = 10,
                          scrollX = TRUE,
                          language = list(search = 'Filter:'),
                          dom = 'lftipB',
                          buttons = c('copy', 'csv', 'excel')))
```


- Espace de travail collaboratif [GitLab](https://forgemia.inra.fr/magali.monnoye/analyses_16s/202111_Rotenone)  

- Adresse du dépot: git@forgemia.inra.fr:magali.monnoye/analyses_16s.git
  
- Repository sur la plateforme migale: <strong class="tool">ssh -X mmonnoye@migale.jouy.inrae.fr</strong> dans save/analyses_16S/202111_Rotenone/


```{bash}
#Pour faire des git clone / git add.....:
#Via le terminal, se positionner dans le repertoire ou se trouve le depot donc ici:
#~/work/Analyses_16S/202109_Utrech/
#et faire les git souhaités...
#Tout ce qu'on veut pusher
git add FROGS*.tsv
git add FROGS*.biom
git add *.html
git add *affiliation.tsv

#Enusite faire les commits
git commit -m "fichiers tsv"

#Puis faire les push
git push

#Si gitlab fatal: pack exceeds maximum allowed size
#Pour supprimer dernier commit:
git reset HEAD~

#Si gros conflits entre les git/Si impossible de faire de git pull/push (quelques pertes d'infos !!!!)
git fetch --all
git reset --hard origin/master

#Pour supprimer les fichiers venndiagram
rm -f VennDiagram*.log
#Pour supprimer de force un dossier et tout ce qu'il contient
rm -rf ... 

# Si j'ai ce message d'erreur: "The following untracked working tree files would be overwritten by merge"
# Faire ceci dan le terminal du rstudio:
git add * 
git stash
git pull

```

# Rappel du contexte de l'étude

magali's avatar
magali committed
94
Deux groupes de 10 souris ont été traitées ou non avec de la rotenone (modèle de Parkinson)
Magali Monnoye's avatar
Magali Monnoye committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205


# Metabarcoding analysis

## Import des séquences

En premier, créer un nouveau dossier dans GitLab / Aller dans Rstudio / Faire un pull / Ouvrir Project

Sur le terminal, je copie les dossiers nécéssaires pour l'export html

```{bash}

#sur le terminal, je me mets dans le repertoire cd 202110_BALBc/
#puis je copie les fichiers ci dessous dans le nouveau dossier Rotenone:

cp report_BALBc.Rmd ~/save/analyses_16s/202111_Rotenone/
cp gitlab-ci.yml ~/save/analyses_16s/202111_Rotenone/

cp -r resources/ ~/save/analyses_16s/202111_Rotenone/
cp -r site_libs/ ~/save/analyses_16s/202111_Rotenone/
cp -r css/ ~/save/analyses_16s/202111_Rotenone/

```

Puis j'importe les séquences

```{bash}
Import des séquences (fichiers FASTQ) sur migale : 
Je crée un zip de mon fichier séquences (Selection des séquences / clic droit / 7-zip / Ajouter à l'archive / Format de l'archive .zip) et j'importe avec Upload

Ou si trop de séquences, copier le dossier de séquences de l'ordi sur le serveur comme ceci: 
Ouvrir le terminal et après C:\Users\magali> copier cette ligne:
scp -r .\Documents\1-MANIPS\2021\2021-09_Lea\Sequences_BALBc\ mmonnoye@migale.jouy.inrae.fr:~/save/analyses_16s/202110_BALBc/

```

## Contrôle qualité

Je me place dans le repertoire 202111_Rotenone/

Je lance l'outil <strong class="tool">fastqc</strong> dédié à l'analyse de la qualité des basemodules issues d'un séquençage haut-débit

J'utilise aussi <strong class="tool">multiqc</strong>


```{bash}
#Creation dossier FASTQC
mkdir FASTQC
mkdir LOGS

#Je lance fastq
for i in *.fastq.gz ; do echo "conda activate fastqc-0.11.8 && fastqc $i -o FASTQC && conda deactivate" >> fastqc.sh ; done
qarray -cwd -V -N fastqc -o LOGS -e LOGS fastqc.sh

#Je lance multiqc (-o pour specifier le dossier de destination)
qsub -cwd -V -N multiqc -o LOGS -e LOGS -b y "conda activate multiqc-1.8 && multiqc FASTQC -o MULTIQC && conda deactivate"

```

Numéro du job **fastqc** sur le cluster de migale: 4045739.1-40:1

Numéro du job **multiqc** sur le cluster de migale: 4045740

Lien vers le [rapport multiqc](MULTIQC/multiqc_report.html)

<div class="alert alert-success" role="alert">**Sequence Counts**

L'échantillon 2.1 a un nombre de reads très faible.

Cet échantillon devra être supprimé lors de l'analyse biostatistique.

**Sequence Quality Histograms**

Les scores de qualité moyens sont plutôt bons. 

La courbe diminue un peu plus pour certains (zone orange) mais le chevauchement de R1 et R2 peut surmonter cela.

Touts les reads ont une longueur de 250 bases. Il indique qu'aucun rognage n'a été effectué précédemment.
     
**Per Sequence Quality Scores**

Une large majorité des reads à une qualité moyenne > 30 (99.9 % of confidence)

**Per Base Sequence Content**

L'outils est noté "failed" car la proportion des bases est biaisée en début de lecture ce qui est normal pour ce type de librairie.

Cela ne nuiera pas aux analyses en aval.

**Per Sequence GC Content**

Not informative for amplicon data
    
**Per Base N Content**

Il y a très peu de N bases présentes. 

**Sequence Duplication Levels**

Not informative for amplicon data

**Overrepresented sequences**

Not informative for amplicon data

**Adapter Content** 

Les adaptateurs Illumina sont présents à différents niveaux pour certains échantillons. C'est représentatif de petits fragments qui ont été séquencés. 

Ces séquences seront supprimées plus tard avec FROGS.</div> 

magali's avatar
magali committed
206
207

<div class="alert comment">En résumé, à part pour l'échantillon 2.1, les séquences sont de très bonnes qualités</div>
Magali Monnoye's avatar
Magali Monnoye committed
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379

```{bash}
#Creer archive 
tar zcvf data.tar.gz *.fastq.gz
#suppression des fichiers fastq car ils sont dans l'archive
rm -f *.fastq.gz
#suppression des fichiers zip dans FASTQC
rm -f FASTQC/*.zip

#Je copie l'archive data.tar.gz dans 202109_Utrech
cp data.tar.gz ~/save/analyses_16s/202110_BALBc/
cp -r MULTIQC ~/save/analyses_16s/202110_BALBc/

#Suppression dossier Sequences_BALBc
rm -rf Sequences_BALBc

```

## FROGS 

### Nettoyage des reads

Je vais d'abord utiliser l'outil <strong class="tool">preprocess</strong> 

```{bash}
#Je me déplace dans 202110_BALBc et je cree fichier FROGS
mkdir FROGS
mkdir LOGS

qsub -cwd -V -N preprocess -o LOGS -e LOGS -pe thread 8 -R y -b y "conda activate frogs-3.2.2 && preprocess.py illumina --input-archive data.tar.gz --min-amplicon-size 200 --max-amplicon-size 490 --merge-software pear --five-prim-primer ACGGGAGGCAGCAG --three-prim-primer GGATTAGATACCCTGGTA --R1-size 250 --R2-size 250 --nb-cpus 8 --output-dereplicated FROGS/preprocess.fasta --output-count FROGS/preprocess.tsv --summary FROGS/preprocess.html --log-file FROGS/preprocess.log && conda deactivate"

#Pour voir avancé d'un job
qstat
#Pour voir infos quand job terminé
qacct -j 3982109

```

Numéro du job **preprocess** sur le cluster de migale: 4046130

Les paramètres suivants ont été choisis : 

|Parametre | Valeur | 
|----------|--------|
|--min-amplicon-size | 200 | 
|--max-amplicon-size | 490 | 
|--already-contiged |Séquences non contiguées |

Lien vers le [rapport preprocess](FROGS/preprocess.html)

- Le nombre de séquences par échantillon

|Preprocess |Min |Max |
|---------|--------|--------|
|Avant|27598|49717|
|Après|26228|47574|
 
<div class="alert comment">Hormis pour l'échantillon 2.1, le pourcentrage de séquences gardées est de 92 à 96 %, ce qui est très bien et normal à ce stade</div>

### Clustering

Clustering à l'aide de <strong class="tool">swarm</strong> et d=1 

```{bash}
qsub -cwd -V -N clustering -o LOGS -e LOGS -pe thread 8 -R y -b y "conda activate frogs-3.2.2 && clustering.py --input-fasta FROGS/preprocess.fasta --input-count FROGS/preprocess.tsv --distance 1 --fastidious --nb-cpus 8 --log-file FROGS/clustering.log --output-biom FROGS/clustering.biom --output-fasta FROGS/clustering.fasta --output-compo FROGS/clustering_otu_compositions.tsv && conda deactivate"

#Pour voir infos quand job en cours:
qstat
#Pour avoir infos quand job fini:
qacct -j 4046530

```

Numéro du job **clustering** sur le cluster de migale: 4046530

### Remove Chimera

On détecte les chimeres avec <strong class="tool">vsearch</strong> 

```{bash}
qsub -cwd -V -N chimera -o LOGS -e LOGS -pe thread 8 -R y -b y "conda activate frogs-3.2.2 && remove_chimera.py --input-fasta FROGS/clustering.fasta --input-biom FROGS/clustering.biom --non-chimera FROGS/remove_chimera.fasta --nb-cpus 8 --log-file FROGS/remove_chimera.log --out-abundance FROGS/remove_chimera.biom --summary FROGS/remove_chimera.html && conda deactivate"

qacct -j 4046537

```

Numéro du job **chimera** sur le cluster de migale: 4046537

Lien vers le [rapport remove_chimera](FROGS/remove_chimera.html)

- Voici le nombre de séquences gardées et retirées  

|Remove Chimera|Kept|Removed
|---------|--------|--------|
|Clusters|48711|55712|
|Abundance|536743|158730|

<div class="alert comment">Le nombre de clusters gardés par échantillon varie lui de 2325 (3.6) à 3973(4.3), ceci est normal à ce stade car on retire les séquences chimériques</div>

### Filtres sur l'abondance et la prevalence

On retire les OTUs les moins abondants, limite 0,005%

On retire egalement les séquences phiX

J'utlise l'outil <strong class="tool">otu_filters</strong>

```{bash}
qsub -cwd -V -N filters -o LOGS -e LOGS -pe thread 8 -R y -b y "conda activate frogs-3.2.2 && otu_filters.py --input-fasta FROGS/remove_chimera.fasta --input-biom FROGS/remove_chimera.biom --output-fasta FROGS/filters.fasta --nb-cpus 8 --log-file FROGS/filters.log --output-biom FROGS/filters.biom --summary FROGS/filters.html --excluded FROGS/filters_excluded.tsv --contaminant /db/frogs_databanks/contaminants/phi.fa --min-sample-presence 1 --min-abundance 0.00005 && conda deactivate"

qacct -j 4046541
```

Numéro du job **filters** sur le cluster de migale: 4046541

Lien vers le [rapport filters](FROGS/filters.html)

- Voici le nombre d'OTUs et de séquences gardés et retirés  

|Filters|Kept|Removed |
|---------|--------|--------|
|OTUs|620|48091|
|Abundance|473260|63483|

### Affiliation 

Maintenant, je vais affillier les séquences sur la database <strong class="tool">silva</strong> v138

Création fichier **affiliation.biom** pour import **affiliationExplorer**

```{bash}
qsub -cwd -V -N affiliation -o LOGS -e LOGS -pe thread 8 -R y -b y "conda activate frogs-3.2.2 && affiliation_OTU.py --input-fasta FROGS/filters.fasta --input-biom FROGS/filters.biom --nb-cpus 8 --log-file FROGS/affiliation.log --output-biom FROGS/affiliation.biom --summary FROGS/affiliation.html --reference /db/frogs_databanks/assignation/silva_138_16S/silva_138_16S.fasta && conda deactivate"

Job 4046545
qstat

#Rapport de l'affiliation des OTUs
qsub -cwd -V -N affiliations_stats -o LOGS -e LOGS -b y "conda activate frogs-3.2.2 && affiliations_stat.py --input-biom FROGS/affiliation.biom --output-file FROGS/affiliations_stats.html --log-file FROGS/affiliations_stats.log --multiple-tag blast_affiliations --tax-consensus-tag blast_taxonomy --identity-tag perc_identity --coverage-tag perc_query_coverage  && conda deactivate"

```

Numéro du job **affiliation** sur le cluster de migale: 4046545

Lien vers le [rapport affiliation_OTU](FROGS/affiliation.html)

Numéro du job **affiliation_stats** sur le cluster de migale: 4046552

Lien vers le [rapport affiliations_stats](FROGS/affiliations_stats.html)

<div class="alert comment">La quasi totalité des échantillons ont un pourcentage d'identité >95% et un pourcentage de couverture de >99%</div>

### Multi affiliation

Création fichiers **multi_aff.tsv** pour import **affiliationExplorer** et **affiliation.tsv** qui sera la table d'abondance 

```{bash}
qsub -cwd -V -N biom_to_tsv -o LOGS -e LOGS -b y "conda activate frogs-3.2.2 && biom_to_tsv.py --input-biom FROGS/affiliation.biom --input-fasta FROGS/filters.fasta --output-tsv FROGS/affiliation.tsv --output-multi-affi FROGS/multi_aff.tsv --log-file FROGS/biom_to_tsv.log  && conda deactivate"

```

Numéro du job **biom_to_tsv** sur le cluster de migale: 4046554

Je vais "rectifier" toutes les multis affiliations sur **affiliationExplorer** en important **multi_aff.tsv** et **affiliation.biom**

* Lien pour accéder au module: [affiliationExplorer](https://shiny.migale.inrae.fr/app/affiliationexplorer)

Je refais un  <strong class="tool">affiliations_stats</strong> sur le biom nettoyé des multis affiliations

```{bash}
#Il faut importer le fichiers "cleaned_biom-2021-06-17.biom" généré par affiliationExplorer

#Affiliation stat pour comparer avant et après correction avec affiliationExplorer
Magali Monnoye's avatar
Magali Monnoye committed
380
qsub -cwd -V -N affiliations_stats_cleaned -o LOGS -e LOGS -b y "conda activate frogs-3.2.2 && affiliations_stat.py --input-biom cleaned_biom-2021-11-17.biom --output-file FROGS/affiliations_stats_cleaned.html --log-file FROGS/affiliations_stats_cleaned.log --multiple-tag blast_affiliations --tax-consensus-tag blast_taxonomy --identity-tag perc_identity --coverage-tag perc_query_coverage  && conda deactivate"
Magali Monnoye's avatar
Magali Monnoye committed
381
382
383

```

Magali Monnoye's avatar
Magali Monnoye committed
384
385
Numéro du job **affiliations_stats_cleaned** sur le cluster de migale: 4048514

Magali Monnoye's avatar
Magali Monnoye committed
386
387
Lien vers le [rapport affiliations_stats_cleaned](FROGS/affiliations_stats_cleaned.html)

Magali Monnoye's avatar
Magali Monnoye committed
388
<div class="alert comment">Suite aux différents outils FROGS utilisés pour le nettoyage des séquences, l'échantillon 2.1 n'a plus que 42 séquences alors que tous les autres échantillons ont un nombre de séquences > 19000, cet échantillon devra être supprimé pour les analyses statistiques.</div>  
Magali Monnoye's avatar
Magali Monnoye committed
389
390
391
392
393
394
395
396
397
398
399

### Tree

Je crée l'arbre phylogénique avec <strong class="tool">tree</strong>

```{bash}

#Pour avoir aide sur outils tree.py
conda activate frogs-3.2.2 && tree.py --help

#Creation arbre phylogenique
Magali Monnoye's avatar
Magali Monnoye committed
400
qsub -cwd -V -N tree -o LOGS -e LOGS -b y "conda activate frogs-3.2.2 && tree.py --input-sequences  FROGS/filters.fasta --biom-file cleaned_biom-2021-11-17.biom --out-tree FROGS/tree.nhx --html FROGS/tree.html --log-file FROGS/tree.log  && conda deactivate"
Magali Monnoye's avatar
Magali Monnoye committed
401
402
403

```

Magali Monnoye's avatar
Magali Monnoye committed
404
405
Numéro du job **tree** sur le cluster de migale: 4048515

Magali Monnoye's avatar
Magali Monnoye committed
406
407
Lien vers l'arbre [Tree](FROGS/tree.html)

Magali Monnoye's avatar
Magali Monnoye committed
408
409
410
411
412
413
414
415
416
417
```{bash}
Copier fichier du serveur sur mon ordi: 
1. Sur le terminal, je me place dans le dossier de destination des fichiers
PS C:\Users\magali> cd .\Documents\1-MANIPS\2021\2021-11_Rotenone\

2. Puis je copie cette ligne
scp mmonnoye@migale.jouy.inrae.fr:~/save/analyses_16s/202111_Rotenone/data.tar.gz .

```

Magali Monnoye's avatar
Magali Monnoye committed
418
419
# Metagenomic phyloseq analysis 

magali's avatar
magali committed
420

Magali Monnoye's avatar
Magali Monnoye committed
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
## Import packages et données

```{r load-packages, eval=TRUE, cache = FALSE}

library(kableExtra)
library(DT)
library(tidyverse)
library(phyloseq)
library(phyloseq.extended)
library(data.table)
library(plotly)
library(cowplot)
library(vegan)
```

J'importe les données 

```{r import-data, eval=TRUE}
#si impossible de trouver la fonction "read.tree"
library(ape) 

 if(file.exists("frogs.data.rds")){
   frogs.data <- readRDS("frogs.data.rds")
 }else{
magali's avatar
magali committed
445
  frogsBiom <- "cleaned_biom-2021-11-17.biom"
Magali Monnoye's avatar
Magali Monnoye committed
446
447
  frogs.data.biom <- import_frogs(frogsBiom, taxMethod = "blast")
  
magali's avatar
magali committed
448
  metadata<-read.table("metadatasRotenone.txt",row.names=1, header=T)
Magali Monnoye's avatar
Magali Monnoye committed
449
450
451
452
453
  
  sample_data(frogs.data.biom)<-metadata
  sample_names(frogs.data.biom) <- get_variable(frogs.data.biom, "Name")
  frogs.data<-frogs.data.biom
  
magali's avatar
magali committed
454
  treefile<- read.tree("FROGS/tree.nhx") 
Magali Monnoye's avatar
Magali Monnoye committed
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
  phy_tree(frogs.data) <- treefile
  
  tax_table(frogs.data) <- gsub("\"", "", tax_table(frogs.data))
  tax_table(frogs.data) <- gsub("unknown.*", "Unknown", tax_table(frogs.data))
  
  #sauvegarde de l'objet forgs.data et .rds
  saveRDS(frogs.data, file = "frogs.data.rds")
  
}

```

Je dispose de l'object phyloseq suivant:

```{r data-preview, eval=TRUE, echo=TRUE}
frogs.data

```

qui contient les métadonnées suivantes:

```{r metadata, eval=TRUE, echo=TRUE}
sample_variables(frogs.data)
```

magali's avatar
magali committed
480
## Sélection des données
Magali Monnoye's avatar
Magali Monnoye committed
481

magali's avatar
magali committed
482
<div class="alert alert-danger" role="alert">Je supprime l'échantillon 2.1 car il possède trop peu de reads</div>
Magali Monnoye's avatar
Magali Monnoye committed
483
484


magali's avatar
magali committed
485
486
487
```{r sample selection,fig.width=10,fig.height=4, eval=TRUE}
#Je supprime les échantillons TF51_2, T026, TF50_2 et T023_2 qui ont un nombre de séquences < 5000
frogs.data <- subset_samples(frogs.data, !(sample_names(frogs.data) %in% c("2.1")) )
Magali Monnoye's avatar
Magali Monnoye committed
488

magali's avatar
magali committed
489
490
frogs.data 
```
Magali Monnoye's avatar
Magali Monnoye committed
491
492
493
494
495
496

* Nombre d'échantillons par groupe:

```{r variable-1, eval = TRUE, echo=TRUE}
table(get_variable(frogs.data, "Group"))
```
magali's avatar
magali committed
497
Le groupe **Vehicule** sans rotenone
Magali Monnoye's avatar
Magali Monnoye committed
498

magali's avatar
magali committed
499
Le groupe **Rotenone** avec rotenone
Magali Monnoye's avatar
Magali Monnoye committed
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528


## Visualisation de l'abondance par échantillon {.tabset}

### Abondances brutes

```{r 01-plot abondances, eval = TRUE}
p<-plot_bar(frogs.data)
plot(p)
```

### Abondances niveau Phylum 

```{r 02-plot abondances phylum, eval = TRUE}
p <- plot_bar(frogs.data, fill = "Phylum")
plot(p)
```


## Composition des communautés 

### Représentation des communautés au niveau Phylum {.tabset}

* Nombre d'OTUs par Phylum:

```{r count-phylum, eval = TRUE}
#Get count of phyla
table<-table(phyloseq::tax_table(frogs.data)[, "Phylum"])
table
magali's avatar
magali committed
529

Magali Monnoye's avatar
Magali Monnoye committed
530
531
532
533
534
535
```

#### Tableaux relatives abondances Phylum 

* Relatives Abondances par échantillons:

magali's avatar
magali committed
536
537
<div class="alert alert-info" role="alert">Abandance Relative = nb total des individus d'une éspèce par rapport au nb total des individus de toutes les espèces présentes</div>

Magali Monnoye's avatar
Magali Monnoye committed
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
```{r relative abundance-phylum, eval = TRUE}

#Convert to relative abundance
data_rel_abund = phyloseq::transform_sample_counts(frogs.data, function(x) 100 * x/sum(x)) 

#Agglomerate to phylum-level and rename
data_phylum <- phyloseq::tax_glom(data_rel_abund, "Phylum")

phyloseq::taxa_names(data_phylum) <- phyloseq::tax_table(data_phylum)[, "Phylum"]
rel_abun<-phyloseq::otu_table(data_phylum)

rel_abun %>%
 DT::datatable(extensions = 'Buttons', 
                                        options = list(dom = 'Bfrtip', 
                                                       pageLength = 20,
                                                       buttons = list('copy', 'csv', 'excel'))) 

```

* Relatives abondances par groupes: 

```{r 9-relative abundance-phylum, eval = TRUE}
#Convert to relative abundance
data_rel_abund = phyloseq::transform_sample_counts(frogs.data, function(x) 100 * x/sum(x))

#Agglomerate to phylum-level and rename
data_phylum <- phyloseq::tax_glom(data_rel_abund, "Phylum")

ps1 <- merge_samples(data_phylum, "Group")
data_phylum_2 <- transform_sample_counts(ps1, function(x) 100 * x/sum(x))

phyloseq::taxa_names(data_phylum_2) <- phyloseq::tax_table(data_phylum_2)[, "Phylum"]
rel_abun2<-phyloseq::otu_table(data_phylum_2)
rel_abun2 %>%
 DT::datatable(extensions = 'Buttons', 
                                        options = list(dom = 'Bfrtip', 
magali's avatar
magali committed
574
                                                       pageLength = 10,
Magali Monnoye's avatar
Magali Monnoye committed
575
576
577
578
579
580
                                                       buttons = list('copy', 'csv', 'excel')))
```


#### Plot composition Phylum

magali's avatar
magali committed
581
582
```{r 03-composition-phylum, eval = TRUE, fig.height=5, fig.width=10}
correct.order <- c("Vehicule","Rotenone")
Magali Monnoye's avatar
Magali Monnoye committed
583
584
585
586
sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group,
levels = correct.order)

p <- plot_composition(frogs.data, "Kingdom", "Bacteria", "Phylum",numberOfTaxa = 6, fill = "Phylum")
magali's avatar
magali committed
587
p<- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +
Magali Monnoye's avatar
Magali Monnoye committed
588
589
590
591
592
593
594
    theme(strip.text.x = element_text(size = 14, color = "black"))  +
      scale_y_continuous(label = scales::percent)

plot(p)
```


magali's avatar
magali committed
595
```{r 04-composition-phylum-merged, eval = TRUE, fig.height=5, fig.width=6}
Magali Monnoye's avatar
Magali Monnoye committed
596
597
598
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)

magali's avatar
magali committed
599
correct.order <- c("Vehicule","Rotenone")
Magali Monnoye's avatar
Magali Monnoye committed
600
601
602
603
sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group,
levels = correct.order)

p <- plot_composition(frogs.data.merged, "Kingdom", "Bacteria", "Phylum",numberOfTaxa = 6, fill = "Phylum")
magali's avatar
magali committed
604
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
Magali Monnoye's avatar
Magali Monnoye committed
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
p <- p + ggtitle("Phylum Composition (12 top Phylum)")+
      scale_y_continuous(label = scales::percent) +  theme(axis.text.x = element_blank())

plot(p)

```

#### Stats Phylum

##### Boxplot Phylum avec stats

On va se focaliser sur les 6 Phylums et représenter leurs abondances dans les différents groupe sous forme de boxplot. 

Le nombre d'étoiles codent le niveau de significativité, les comparaisons non-significatives ne sont pas indiqués

magali's avatar
magali committed
620
<div class="alert alert-info" role="alert">Le test wilcox compare Vehicule à Rotenone</div>
Magali Monnoye's avatar
Magali Monnoye committed
621

magali's avatar
magali committed
622
```{r 05-Phylum stats,fig.width=6,fig.height=8, eval = TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
623
library(ggpubr)
magali's avatar
magali committed
624
625
626
627
628
library(cowplot)

#Select groups
frogs.data_bis <- subset_samples(frogs.data, Group %in% c("Vehicule","Rotenone"))

Magali Monnoye's avatar
Magali Monnoye committed
629
## Select only some families  
magali's avatar
magali committed
630
631
632
families <- c("Actinobacteriota","Bacteroidota","Deferribacterota", "Desulfobacterota" ,"Firmicutes", "Proteobacteria")
phy <- frogs.data_bis %>% subset_taxa(Phylum %in% families) %>% tax_glom(taxrank="Phylum") #agglomerate at family level

Magali Monnoye's avatar
Magali Monnoye committed
633
634
635
636
 ## Transform count to relative abundances 
depth <- sample_sums(frogs.data)[1] 
plotdata<-psmelt(phy) %>% 
  mutate(Abundance  = Abundance / depth, 
magali's avatar
magali committed
637
         Group = factor(Group, labels = c("Vehicule","Rotenone" )))  
Magali Monnoye's avatar
Magali Monnoye committed
638
639
640
641
642
  
p <- ggplot(plotdata,aes(x = Group,y=Abundance, color = Group, Group = Group)) + 
  stat_boxplot(geom = "errorbar", width = 0.5) + 
  geom_boxplot(outlier.alpha = 1, 
               outlier.size = 0.8) + 
magali's avatar
magali committed
643
644
  facet_wrap(~Phylum, scales = "free_y", ncol = 3) + 
  scale_color_manual(values = c("Vehicule" = "#227432", "Rotenone" = "#8C67A9"), guide = "none") + 
Magali Monnoye's avatar
Magali Monnoye committed
645
646
647
648
649
650
  
  # theme_classic() + ## fond blanc
  labs(x = NULL) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10))

## Compare reference level T0 to all other levels using a wilcoxon test and adjust p-values using the holm correction. 
magali's avatar
magali committed
651
p <- p + stat_compare_means(aes(label = ..p.signif..), method = "wilcox.test", p.adjust.method = "holm", ref.group = "Vehicule", hide.ns = T, 
Magali Monnoye's avatar
Magali Monnoye committed
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
                       label.y.npc = c(0.90),
                       size = 7, 
                       fontface = "bold")
plot(p)

```

##### Kruskal test

Test de Kruskal (Anova non paramétrique) pour chaque Phylum pour tester si les abondances sont similaires (pour ce Phylum) entre les différents groupes. 

```{r data.-phylum, eval=TRUE}
# depth <- sample_sums(frogs.data)[1] 
data.phylum <- frogs.data %>% 
  transform_sample_counts(function(x) { x / sum(x)}) %>% ## transform counts to proportions
  fast_tax_glom(taxrank = "Phylum") %>% 
  psmelt() %>% 
magali's avatar
magali committed
669
  mutate(Group = factor(Group, labels = c("Vehicule","Rotenone")))
Magali Monnoye's avatar
Magali Monnoye committed
670
671
672
673
674
675
676
677
678
679
680
681
```

```{r test kruskal-phylum, eval=TRUE}
data.test <- compare_means(Abundance ~ Group, data = data.phylum, method = "kruskal", group.by = "Phylum")
data.test
```


### Représentation des communautés au niveau Family {.tabset} 

* Nombre d'OTUs par Family:

magali's avatar
magali committed
682
```{r count-family, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
683
#Get count of phyla
magali's avatar
magali committed
684
685
 table(phyloseq::tax_table(frogs.data)[, "Family"])

Magali Monnoye's avatar
Magali Monnoye committed
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
```

#### Tableaux relatives abondances Family 

* Relative Abondance par échantillon:

```{r relative abundance-phylum-family-1, eval = TRUE}
library(openxlsx)

#Convert to relative abundance
data_rel_abund = phyloseq::transform_sample_counts(frogs.data, function(x) 100 * x/sum(x)) 

#Agglomerate to phylum-level and rename
data_phylum <- phyloseq::tax_glom(data_rel_abund, "Family")

phyloseq::taxa_names(data_phylum) <- phyloseq::tax_table(data_phylum)[, "Family"]
rel_abun_3<-phyloseq::otu_table(data_phylum)

rel_abun_3 %>%
 DT::datatable(extensions = 'Buttons', 
                                        options = list(dom = 'Bfrtip', 
                                                       pageLength = 10,
                                                       buttons = list('copy', 'csv', 'excel')))
```

* Relative Abondance par groupe:

magali's avatar
magali committed
713
```{r 2.1-relative abundance-family-2, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
#Convert to relative abundance
data_rel_abund = phyloseq::transform_sample_counts(frogs.data, function(x) 100 * x/sum(x))

#Agglomerate to phylum-level and rename
data_phylum <- phyloseq::tax_glom(data_rel_abund, "Family")

ps1 <- merge_samples(data_phylum, "Group")
data_phylum_2 <- transform_sample_counts(ps1, function(x) 100 * x/sum(x))


phyloseq::taxa_names(data_phylum_2) <- phyloseq::tax_table(data_phylum_2)[, "Family"]
rel_abun_4<-phyloseq::otu_table(data_phylum_2)

rel_abun_4 %>%
 DT::datatable(extensions = 'Buttons', 
                                        options = list(dom = 'Bfrtip', 
                                                       pageLength = 10,
                                                       buttons = list('copy', 'csv', 'excel')))
```


#### Plot composition Family 

magali's avatar
magali committed
737
738
```{r 06-composition-phylum,fig.width=10,fig.height=5, eval=TRUE}
correct.order <- c("Vehicule","Rotenone")
Magali Monnoye's avatar
Magali Monnoye committed
739
740
741
742
743
sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group,
levels = correct.order)


p <- plot_composition(frogs.data, "Kingdom", "Bacteria", "Family",numberOfTaxa = 10, fill = "Family")
magali's avatar
magali committed
744
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
Magali Monnoye's avatar
Magali Monnoye committed
745
746
747
748
749

plot(p)

```

magali's avatar
magali committed
750
```{r 07-composition-family merged,fig.width=6,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
751
752
753
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)

magali's avatar
magali committed
754
correct.order <- c("Vehicule","Rotenone")
Magali Monnoye's avatar
Magali Monnoye committed
755
756
757
758
sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group,
levels = correct.order)

p <- plot_composition(frogs.data.merged, "Kingdom", "Bacteria", "Family",numberOfTaxa = 10, fill = "Family")
magali's avatar
magali committed
759
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +  theme(axis.text.x = element_blank())   
Magali Monnoye's avatar
Magali Monnoye committed
760
761
762
763
764
765
766
767
768
769
770
771
plot(p)

```

#### Stats Family

##### Boxplot Family avec stats

On va se focaliser sur les 10 Family et représenter leurs abondances dans les différents groupe sous forme de boxplot. 

Le nombre d'étoiles codent le niveau de significativité, les comparaisons non-significatives ne sont pas indiqués

magali's avatar
magali committed
772
773
774
775
776
777
778
<div class="alert alert-info" role="alert">Le test wilcox compare Vehicule à Rotenone </div>

```{r 08-Family stats,fig.width=10,fig.height=8, eval = TRUE}

#Select groups
frogs.data_quat <- subset_samples(frogs.data, Group %in% c("Vehicule","Rotenone"))

Magali Monnoye's avatar
Magali Monnoye committed
779
## Select only some families  
magali's avatar
magali committed
780
781
782
families <- c("Muribaculaceae","Marinifilaceae", "Rikenellaceae", "Tannerellaceae","Prevotellaceae" , "Bacteroidaceae", "Lachnospiraceae", "Ruminococcaceae","Anaerovoracaceae", "Peptostreptococcaceae")
phy <- frogs.data_quat %>% subset_taxa(Family %in% families) %>% tax_glom(taxrank="Family") #agglomerate at family level

Magali Monnoye's avatar
Magali Monnoye committed
783
 ## Transform count to relative abundances 
magali's avatar
magali committed
784
depth <- sample_sums(frogs.data_quat)[1] 
Magali Monnoye's avatar
Magali Monnoye committed
785
786
plotdata<-psmelt(phy) %>% 
  mutate(Abundance  = Abundance / depth, 
magali's avatar
magali committed
787
         Group = factor(Group, labels = c("Vehicule","Rotenone")))  
Magali Monnoye's avatar
Magali Monnoye committed
788
789
790
791
792
  
p <- ggplot(plotdata,aes(x = Group,y=Abundance, color = Group, Group = Group)) + 
  stat_boxplot(geom = "errorbar", width = 0.5) + 
  geom_boxplot(outlier.alpha = 1, 
               outlier.size = 0.8) + 
magali's avatar
magali committed
793
794
  facet_wrap(~Family, scales = "free_y", ncol = 5) + 
  scale_color_manual(values = c("Vehicule" = "#227432", "Rotenone" = "#8C67A9"), guide = "none") + 
Magali Monnoye's avatar
Magali Monnoye committed
795
796
797
798
799
800
  
  # theme_classic() + ## fond blanc
  labs(x = NULL) + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 10))

## Compare reference level T0 to all other levels using a wilcoxon test and adjust p-values using the holm correction. 
magali's avatar
magali committed
801
p <- p + stat_compare_means(aes(label = ..p.signif..), method = "wilcox.test", p.adjust.method = "holm", ref.group = "Vehicule", hide.ns = T, 
Magali Monnoye's avatar
Magali Monnoye committed
802
803
804
805
806
807
808
809
                       label.y.npc = c(0.90),
                       size = 7, 
                       fontface = "bold")
plot(p)
```

##### Kruskal test

magali's avatar
magali committed
810
Test de Kruskal (Anova non paramétrique) pour chaque Family pour tester si les abondances sont similaires (pour la Family) entre les différents groupes. 
Magali Monnoye's avatar
Magali Monnoye committed
811
812
813
814
815
816
817

```{r data.-family, eval=TRUE}
# depth <- sample_sums(frogs.data)[1] 
data.phylum <- frogs.data %>% 
  transform_sample_counts(function(x) { x / sum(x)}) %>% ## transform counts to proportions
  fast_tax_glom(taxrank = "Family") %>% 
  psmelt() %>% 
magali's avatar
magali committed
818
  mutate(Group = factor(Group, labels = c("Vehicule","Rotenone")))
Magali Monnoye's avatar
Magali Monnoye committed
819
820
821
822
823
824
825
826
827
828
829
830
831
```

```{r test kruskal-family, eval=TRUE}
data.test <- compare_means(Abundance ~ Group, data = data.phylum, method = "kruskal", group.by = "Family")
data.test
```

### Zoom au niveau Family dans chacun des Phylums {.tabset}  

Les 6 phylums d'intérêts sont les Actinobacteriota, les Bacteroidota, les Desulfobacterota, les Fibrobacterota, les Firmicutes et les Proteobacteria.

#### Actinobacteriota

magali's avatar
magali committed
832
833
```{r 09-composition-Actinobacteriota,fig.width=8,fig.height=5, eval=TRUE}
correct.order <- c("Vehicule","Rotenone")
Magali Monnoye's avatar
Magali Monnoye committed
834
835
836
837
sample_data(frogs.data)$Group <- factor(sample_data(frogs.data)$Group,
levels = correct.order)

p <- plot_composition(frogs.data, "Phylum", "Actinobacteriota", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
838
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
Magali Monnoye's avatar
Magali Monnoye committed
839
840
841
842

plot(p)
```

magali's avatar
magali committed
843
```{r 10-composition-Actinobacteriota merged,fig.width=6,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
844
845
846
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)

magali's avatar
magali committed
847
correct.order <- c("Vehicule","Rotenone")
Magali Monnoye's avatar
Magali Monnoye committed
848
849
850
851
sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group,
levels = correct.order)

p <- plot_composition(frogs.data.merged, "Phylum", "Actinobacteriota", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
852
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +  theme(axis.text.x = element_blank())   
Magali Monnoye's avatar
Magali Monnoye committed
853
854
855
856
857
858
plot(p)

```

#### Bacteroidota

magali's avatar
magali committed
859
```{r 11-composition-bacteroidetes,fig.width=8,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
860
861

p <- plot_composition(frogs.data, "Phylum", "Bacteroidota", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
862
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
Magali Monnoye's avatar
Magali Monnoye committed
863
864
865
866

plot(p)
```

magali's avatar
magali committed
867
```{r 12-composition-Bacteroidota merged,fig.width=6,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
868
869
870
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)

magali's avatar
magali committed
871
correct.order <- c("Vehicule","Rotenone")
Magali Monnoye's avatar
Magali Monnoye committed
872
873
874
875
sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group,
levels = correct.order)

p <- plot_composition(frogs.data.merged, "Phylum", "Bacteroidota", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
876
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +  theme(axis.text.x = element_blank())   
Magali Monnoye's avatar
Magali Monnoye committed
877
878
879
880
881
882
plot(p)

```

#### Desulfobacterota

magali's avatar
magali committed
883
```{r 13-composition-tenericutes,fig.width=8,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
884
885

p <- plot_composition(frogs.data, "Phylum", "Desulfobacterota", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
886
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
Magali Monnoye's avatar
Magali Monnoye committed
887
888
889
890

plot(p)
```

magali's avatar
magali committed
891
```{r 14-composition-Desulfobacterota merged,fig.width=6,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
892
893
894
895
896
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)


p <- plot_composition(frogs.data.merged, "Phylum", "Desulfobacterota", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
897
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +  theme(axis.text.x = element_blank())   
Magali Monnoye's avatar
Magali Monnoye committed
898
899
900
901
plot(p)

```

magali's avatar
magali committed
902
#### Deferribacterota
Magali Monnoye's avatar
Magali Monnoye committed
903

magali's avatar
magali committed
904
```{r 15-composition-Deferribacterota,fig.width=8,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
905

magali's avatar
magali committed
906
907
p <- plot_composition(frogs.data, "Phylum", "Deferribacterota", "Family",numberOfTaxa = 5, fill = "Family")
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
Magali Monnoye's avatar
Magali Monnoye committed
908
909
910
911

plot(p)
```

magali's avatar
magali committed
912
```{r 16-composition-Deferribacterota merged,fig.width=6,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
913
914
915
916
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)


magali's avatar
magali committed
917
918
p <- plot_composition(frogs.data.merged, "Phylum", "Deferribacterota", "Family",numberOfTaxa = 5, fill = "Family")
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +  theme(axis.text.x = element_blank())   
Magali Monnoye's avatar
Magali Monnoye committed
919
920
921
922
923
924
plot(p)

```

#### Firmicutes

magali's avatar
magali committed
925
```{r 17-composition-firmicutes,fig.width=8,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
926
927

p <- plot_composition(frogs.data, "Phylum", "Firmicutes", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
928
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
Magali Monnoye's avatar
Magali Monnoye committed
929
930
931
932

plot(p)
```

magali's avatar
magali committed
933
```{r 18-composition-Firmicutes merged,fig.width=6,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
934
935
936
937
938
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)


p <- plot_composition(frogs.data.merged, "Phylum", "Firmicutes", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
939
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +  theme(axis.text.x = element_blank())   
Magali Monnoye's avatar
Magali Monnoye committed
940
941
942
943
944
945
plot(p)

```

#### Proteobacteria

magali's avatar
magali committed
946
```{r 19-composition-proteobacteria,fig.width=8,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
947
948

p <- plot_composition(frogs.data, "Phylum", "Proteobacteria", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
949
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1)
Magali Monnoye's avatar
Magali Monnoye committed
950
951
952
953

plot(p)
```

magali's avatar
magali committed
954
```{r 20-composition-Proteobacteria merged,fig.width=6,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
955
956
957
958
959
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)


p <- plot_composition(frogs.data.merged, "Phylum", "Proteobacteria", "Family",numberOfTaxa = 5, fill = "Family")
magali's avatar
magali committed
960
p <- p + facet_wrap(~Group, scales = "free_x", nrow = 1) +  theme(axis.text.x = element_blank())   
Magali Monnoye's avatar
Magali Monnoye committed
961
962
963
964
965
966
967
968
969
970
971
972
973
974
plot(p)

```

## Diversité $\alpha$

### Courbes de raréfaction

Comme ces distances sont sensibles à la profondeur d'échantillonnage, on va raréfier les échantillons avant de calculer les distances. 

Normalisation par rarefaction: même nombre de séquences pour chaque échantillon.



magali's avatar
magali committed
975
976
977
978
```{r rarefaction,echo=TRUE, eval = TRUE}
#On peut rajouter le filtre sur l'abondance >= 10000
#abundant.samples <- (sample_sums(frogs.data) >= 10000) 
#frogs.data.rare <- rarefy_even_depth(prune_samples(abundant.samples, frogs.data), rngseed = 20170317)
Magali Monnoye's avatar
Magali Monnoye committed
979
980
981
982
983
984
985
986

frogs.data.rare<-rarefy_even_depth(frogs.data,rngseed=20170329)
sample_sums(frogs.data.rare)[1:5]
```

Avant de calculer les diversités $\alpha$, on va faire des courbes de raréfaction pour vérifier si on a saturé la richesse sous-dominante (i.e. celle qui passe les filtres d'abondances). 


magali's avatar
magali committed
987
```{r 21-ggrare, message=FALSE,fig.width=10,fig.height=6, eval = TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
988
989


magali's avatar
magali committed
990
991
if(file.exists("frogs.data.courbes.rds")){
   p <- readRDS("frogs.data.courbes.rds")
Magali Monnoye's avatar
Magali Monnoye committed
992
993
994
 }else{
  p <- ggrare(frogs.data, step = 100, color = "Group", plot = FALSE)

magali's avatar
magali committed
995
    saveRDS(p, file = "frogs.data.courbes.rds")
Magali Monnoye's avatar
Magali Monnoye committed
996
997
998
999
  }

rare.level <- min(sample_sums(frogs.data))

magali's avatar
magali committed
1000
1001
family.palette<- c("Vehicule" = "#227432",
                   "Rotenone" = "#8C67A9")  
Magali Monnoye's avatar
Magali Monnoye committed
1002
1003

p <- p + facet_wrap(~Group) + geom_vline(xintercept = rare.level, color = "Gray60")+
magali's avatar
magali committed
1004
1005
  xlab("Sample Size (nb de séquences)") + ylab("Species Richness (nb d'OTU)")+ scale_color_manual(values = family.palette)+ scale_fill_manual(breaks = c("Vehicule","Rotenone"), 
                       values=c("#227432","#8C67A9"))
Magali Monnoye's avatar
Magali Monnoye committed
1006
1007
1008
1009
1010
1011

plot(p)
```

### Effet du régime sur la diversité

magali's avatar
magali committed
1012
alpha diversité: Diversité au sein d'un échantillon / Une valeur par échantillon
Magali Monnoye's avatar
Magali Monnoye committed
1013

magali's avatar
magali committed
1014
1015
<div class="alert alert-info" role="alert">
  * **OBSERVED**: compte le nombre d'OTU observés dans chaque echantillon.
Magali Monnoye's avatar
Magali Monnoye committed
1016

magali's avatar
magali committed
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
  * **CHAO1**: Nb d'OTU + Nb d'OTU que l'on a pas vu, applique des corrections pour prendre en compte la diversité.
  
  * **SHANNON**: sensible aux variations d'importance des espèces les plus RARES / la valeur de la fonction de Shannon augmente avec la diversité
	
      + Shannon minimal -> Diversité la plus faible, tous les individus de la population sont de la même espèce (ou une seule espèce est sur-représentée alors que toutes les autres ne contiennent qu'un seul individu).
    
      + Shannon maximal -> Tous les individus sont équitablement répartis entre les espèces. 
	
  * **InvSIMPSON** sensible aux variations d'importance des espèces les plus abondantes (Indice de dominance). 
  
    L'inverse de Simpson permet de faire varier l'indice dans le même sens que la diversité :plus la diversité spécifique est élevée plus l'indice est fort.
</div>    
  
```{r 22-plot-richness,fig.width=10,fig.height=7, eval = TRUE}
#Ne pas utiliser frogs.data.rare car cela altère le resultats sur la diversité
p <- plot_richness(frogs.data, color = "Group",measures = c("Observed", "Chao1", "Shannon", "InvSimpson"),x = "Group") +theme_bw() + geom_boxplot(aes(fill = Group)) +  geom_point() + theme(axis.text.x = element_blank())+ scale_color_manual(values = family.palette)+ scale_fill_manual(breaks = c("Vehicule","Rotenone"), 
  values=c("#227432","#8C67A9"))
Magali Monnoye's avatar
Magali Monnoye committed
1034
1035
1036
1037
1038
1039
1040

plot(p)

```

### Richness Table

magali's avatar
magali committed
1041
1042
1043
1044
```{r create-richness-table, message = FALSE, eval=TRUE}
#On va tester si les différences de richess observées sont significatives en prenant soin de corriger pour les profondeurs de séquençage #inégales (on pourrait aussi raréfier et ne pas corriger)

richness.table <- estimate_richness(frogs.data, measures = c("Observed", "Chao1", "Shannon", "InvSimpson"))
Magali Monnoye's avatar
Magali Monnoye committed
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
richness.table <- cbind(richness.table, sample_data(frogs.data))
richness.table$Depth <- sample_sums(frogs.data)


richness.table %>%
 DT::datatable(extensions = 'Buttons', 
                                        options = list(dom = 'Bfrtip', 
                                                       pageLength = 10,
                                                       buttons = list('copy', 'csv', 'excel')))
```


## Statistiques sur la diversité {.tabset}

### Statistiques sur Richness Observed

#### Anova sur la Richness Observed 

Analyse de la variance, compare les moyennes d'échantillons.

magali's avatar
magali committed
1065
1066
1067
```{r 33-anova, message = FALSE, eval=TRUE}
div_data <- cbind(estimate_richness(frogs.data, measures = "Observed"),  
                  sample_data(frogs.data))
Magali Monnoye's avatar
Magali Monnoye committed
1068
1069
1070
1071
1072
1073
1074
1075
1076
model <- aov(Observed ~ 0 + Group, data = div_data)
anova(model)

## impact of Groupe and Time_point on richness after correction for depth / seulement possible si T0 et Tf
#observed.aov <- aov(Observed ~ Depth + Groupe + Time_point, data = richness.table)
#anova(observed.aov)                  
```

Df: nombre de degrés de liberté
magali's avatar
magali committed
1077

Magali Monnoye's avatar
Magali Monnoye committed
1078
Sum Sq: inertie inter et intra respectivement
magali's avatar
magali committed
1079

Magali Monnoye's avatar
Magali Monnoye committed
1080
Mean Sq: c'est la variance, obtenue en faisant le quotient de la 2eme colonne par la 1ere (moyenne des inerties)
magali's avatar
magali committed
1081

Magali Monnoye's avatar
Magali Monnoye committed
1082
F value: obtenue en faisant quotient des 2 valeurs trouvées dans la 3eme colonne
magali's avatar
magali committed
1083

Magali Monnoye's avatar
Magali Monnoye committed
1084
1085
1086
1087
1088
Pr(>F): pvalue

#### Coefficient Richness Observed

coef function: permet d'extraire que les coefficients estimés, mesure l'écart des groupes à la moyenne.
magali's avatar
magali committed
1089

Magali Monnoye's avatar
Magali Monnoye committed
1090
1091
Groupes équilibrés si les valeurs sont du même ordre de grandeur.

magali's avatar
magali committed
1092
```{r 32-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1093
1094
1095
1096
1097
1098
coef(model)
```

#### Test de comparaisons multiples sur Richness Observed 

Utilisé pour déterminer les différences significatives entre les moyennes de groupes dans une analyse de variance.
magali's avatar
magali committed
1099

Magali Monnoye's avatar
Magali Monnoye committed
1100
1101
Pour étudier les différences inter-groupes, permet de distinguer parmi les échantillons s’il y en a qui diffèrent significativement des autres.

magali's avatar
magali committed
1102
```{r 34-tukey, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1103
1104
1105
TukeyHSD(model)     
```
diff: différences entre les moyennes obséervées
magali's avatar
magali committed
1106

Magali Monnoye's avatar
Magali Monnoye committed
1107
lwr et upr: donnent les bornes inférieure et supérieur de l'intervalle
magali's avatar
magali committed
1108

Magali Monnoye's avatar
Magali Monnoye committed
1109
1110
1111
p adj: pvalue après ajustement pour les comparaisons multiples

```{r wilcox test, eval=TRUE}
magali's avatar
magali committed
1112
#   comparison_data <- compare_means(
Magali Monnoye's avatar
Magali Monnoye committed
1113
#   Observed ~ Group,
magali's avatar
magali committed
1114
#   data = div_data,
Magali Monnoye's avatar
Magali Monnoye committed
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
#   method = "wilcox.test"
#   ) %>% filter(p.adj <= 0.05)
# comparison_data
```


### Statistiques sur Richness Chao1

#### Anova sur la Richness Chao1

magali's avatar
magali committed
1125
```{r 35-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1126

magali's avatar
magali committed
1127
1128
div_data <- cbind(estimate_richness(frogs.data, measures = "Chao1"),  
                  sample_data(frogs.data))
Magali Monnoye's avatar
Magali Monnoye committed
1129
1130
1131
1132
1133
1134
1135
model <- aov(Chao1 ~ 0 + Group, data = div_data)
anova(model)
                  
```

#### Coefficient Chao1

magali's avatar
magali committed
1136
```{r 37-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1137
1138
1139
1140
1141
coef(model)
```

#### Test de comparaisons multiples sur Richness Chao1

magali's avatar
magali committed
1142
```{r 38-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
TukeyHSD(model)     
```

```{r wilcox test chao, eval=TRUE}
# comparison_data <- compare_means(
#   Chao1  ~ Group,
#   data = div_data, 
#   method = "wilcox.test"
#   ) %>% filter(p.adj <= 0.05)
# comparison_data
```


### Statistiques sur Richness shannon

#### Anova sur la Richness Shannon

magali's avatar
magali committed
1160
```{r 43-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1161

magali's avatar
magali committed
1162
1163
div_data <- cbind(estimate_richness(frogs.data, measures = "Shannon"),  
                  sample_data(frogs.data))
Magali Monnoye's avatar
Magali Monnoye committed
1164
1165
1166
1167
1168
1169
1170
model <- aov(Shannon ~ 0 + Group, data = div_data)
anova(model)
                  
```

#### Coefficient Shannon

magali's avatar
magali committed
1171
```{r 44-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1172
1173
1174
1175
1176
coef(model)
```

#### Test de comparaisons multiples sur Richness Shannon

magali's avatar
magali committed
1177
```{r 45-create-richness-table,, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
TukeyHSD(model)     
```


```{r wilcox test Shannon, eval=TRUE }
# comparison_data <- compare_means(
#   Shannon   ~ Group,
#   data = div_data, 
#   method = "wilcox.test"
#   ) %>% filter(p.adj <= 0.05)
# comparison_data
```

### Statistiques sur Richness InvSimpson

#### Anova sur la Richness InvSimpson

magali's avatar
magali committed
1195
```{r 49-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1196

magali's avatar
magali committed
1197
1198
div_data <- cbind(estimate_richness(frogs.data, measures = "InvSimpson"),  
                  sample_data(frogs.data))
Magali Monnoye's avatar
Magali Monnoye committed
1199
1200
1201
1202
1203
1204
1205
model <- aov(InvSimpson ~ 0 + Group, data = div_data)
anova(model)
                  
```

#### Coefficient InvSimpson

magali's avatar
magali committed
1206
```{r 50-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1207
1208
1209
1210
1211
coef(model)
```

#### Test de comparaisons multiples sur Richness InvSimpson

magali's avatar
magali committed
1212
```{r 51-create-richness-table, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
TukeyHSD(model)     
```

```{r wilcox test InvSimpson, eval=TRUE }
# comparison_data <- compare_means(
#   InvSimpson   ~ Group,
#   data = div_data, 
#   method = "wilcox.test"
#   ) %>% filter(p.adj <= 0.05)
# comparison_data
```

magali's avatar
magali committed
1225
## Diversité $\beta$
Magali Monnoye's avatar
Magali Monnoye committed
1226

magali's avatar
magali committed
1227
Diversité entre les communautés, comparaison des échantillons entre eux / les communautés entre elles.
Magali Monnoye's avatar
Magali Monnoye committed
1228

magali's avatar
magali committed
1229
1230
1231
1232
1233
<div class="alert alert-info" role="alert">
  * **Bray-Curtis** (dissimilarité) / Info abondance
	
Bray curtis nous permet de visualiser la composition, l'abondance, la dominance, etc... de communautés dans l'écosystème et peut être utilisé pour obtenir un résultat du degré de dissemblance entre les espèces.
L'indice de Bray-Curtis donne le même poids aux différences d'abondances observées pour les espèces rares que pour les espèces importantes.
Magali Monnoye's avatar
Magali Monnoye committed
1234

magali's avatar
magali committed
1235
1236
1237
Si =1 cela signifie que les communautés sont totalement différentes.
	
  * **Jaccard** (similarité) / Info présence - absence
Magali Monnoye's avatar
Magali Monnoye committed
1238

magali's avatar
magali committed
1239
Cet indice va de 0  les deux échantillons partagent les mêmes OTUs à 1  les deux échantillons n'ont aucunOTU en commun
Magali Monnoye's avatar
Magali Monnoye committed
1240

magali's avatar
magali committed
1241
  * **Unifrac** (phylogénie) prend en compte arbre phylogénique et l'abondance : la distance UniFrac entre deux échantillons se base sur les longueurs des branches de l'arbre partagées entre ces deux échantillons.
Magali Monnoye's avatar
Magali Monnoye committed
1242

magali's avatar
magali committed
1243
1244
1245
      + unifrac = 0 pas de branche partagée.
      
      + unifrac = 1 pas de branche partagée.
Magali Monnoye's avatar
Magali Monnoye committed
1246

magali's avatar
magali committed
1247
  * **Weighted UniFrac** (pondéré) intègre les abondances lors du calcul des longueurs de branches partagées / non partagées pour calculer la distance.
Magali Monnoye's avatar
Magali Monnoye committed
1248

magali's avatar
magali committed
1249
1250
UniFrac non pondéré est plus sensible aux différences de caractéristiques de faible abondance. 
</div>
Magali Monnoye's avatar
Magali Monnoye committed
1251
1252
1253
1254

### Calcul des distances 


magali's avatar
magali committed
1255
```{r 23-dist-as-heatmap,fig.width=12,fig.height=11, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278

SampleOrder <- levels(reorder(sample_names(frogs.data.rare), as.numeric(get_variable(frogs.data.rare, "Group"))))
dist.bc <- distance(frogs.data.rare, "bray")
p.bc<-plot_dist_as_heatmap(dist.bc, order = SampleOrder,title="Bray Distances", show.names = T)

SampleOrder <- levels(reorder(sample_names(frogs.data.rare), as.numeric(get_variable(frogs.data.rare, "Group"))))
dist.jac <- distance(frogs.data.rare, "cc")
p.jac<-plot_dist_as_heatmap(dist.jac, order = SampleOrder,title="Jaccard Distances", show.names = T)

SampleOrder <- levels(reorder(sample_names(frogs.data.rare), as.numeric(get_variable(frogs.data.rare, "Group"))))
dist.uf <- distance(frogs.data.rare, "unifrac")
p.uf<-plot_dist_as_heatmap(dist.uf, order = SampleOrder,title="Unifrac Distances", show.names = T)

SampleOrder <- levels(reorder(sample_names(frogs.data.rare), as.numeric(get_variable(frogs.data.rare, "Group"))))
dist.wuf <- distance(frogs.data.rare, "wunifrac")
p.wuf<- plot_dist_as_heatmap(dist.wuf, order = SampleOrder,title="Weighted Unifrac Distances", show.names = T)

plot_grid(
  p.bc, p.jac, p.uf, p.wuf, ncol = 2)
```

### Ordination des échantillons {.tabset} 

magali's avatar
magali committed
1279
<div class="alert alert-info" role="alert">On va faire une MDS (non-parametric multi-dimensional scaling) aussi appelé PCoA pour avoir une représentation en 2D de la distance entre tous les échantillons. 
Magali Monnoye's avatar
Magali Monnoye committed
1280
1281
Permet de visusaliser les similitudes entre les groupes.

magali's avatar
magali committed
1282
La mise à l'échelle multidimensionnelle (MDS) est une approche populaire pour représenter graphiquement les relations entre des objets (par exemple, des tracés ou des échantillons) dans un espace multidimensionnel.</div> 
Magali Monnoye's avatar
Magali Monnoye committed
1283
1284
1285

#### Plot ordination Bray Curtis

magali's avatar
magali committed
1286
```{r 24-ord-plot-bray,fig.width=10,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1287
1288
ord <- ordinate(frogs.data.rare, method = "MDS", distance = "bray")

magali's avatar
magali committed
1289
1290
1291
1292
1293
family.palette<- c("Vehicule" = "#227432",
                   "Rotenone" = "#8C67A9")


p <- plot_ordination(frogs.data.rare, ord, color = "Group",shape="Group")
Magali Monnoye's avatar
Magali Monnoye committed
1294
p <- p + theme_bw() + ggtitle("MDS + Bray Curtis")+ scale_color_manual(values = family.palette)
magali's avatar
magali committed
1295
p <- p +geom_point(size=3)+ stat_ellipse(aes(group = Group), size = 1)
Magali Monnoye's avatar
Magali Monnoye committed
1296
1297

plot(p)
magali's avatar
magali committed
1298

Magali Monnoye's avatar
Magali Monnoye committed
1299
1300
```

magali's avatar
magali committed
1301

Magali Monnoye's avatar
Magali Monnoye committed
1302
<div class="alert alert-info" role="alert">Pour confirmer l'effet du régime sur les distances, on va faire une permanova. 
magali's avatar
magali committed
1303

Magali Monnoye's avatar
Magali Monnoye committed
1304
1305
Analyse multivariée de la variance par permutations basée sur les matrices de distances</div> 

magali's avatar
magali committed
1306
* Adonis sur la matrice Bray Curtis 
Magali Monnoye's avatar
Magali Monnoye committed
1307
1308


magali's avatar
magali committed
1309
1310
```{r 35-adonis-bc-1, eval=TRUE}
adonis(dist.bc ~ Group ,data = as(sample_data(frogs.data.rare), 'data.frame'))
Magali Monnoye's avatar
Magali Monnoye committed
1311
1312
```
Df: degree of freedom. 
magali's avatar
magali committed
1313
1314
1315

Sums Of Sqs: sum of squares.

Magali Monnoye's avatar
Magali Monnoye committed
1316
MeanSqs: sum of squares by degree of freedom. 
magali's avatar
magali committed
1317

Magali Monnoye's avatar
Magali Monnoye committed
1318
F: F statistics. 
magali's avatar
magali committed
1319

Magali Monnoye's avatar
Magali Monnoye committed
1320
R2: partial R-squared. 
magali's avatar
magali committed
1321

Magali Monnoye's avatar
Magali Monnoye committed
1322
1323
1324
1325
1326
Pr(>F) p-value based


#### Plot ordination Jaccard

magali's avatar
magali committed
1327
1328

```{r 28-ord-plot-jac,fig.width=10,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1329
1330
1331
ord <- ordinate(frogs.data.rare, method = "MDS", distance = "cc")


magali's avatar
magali committed
1332
1333
1334
p <- plot_ordination(frogs.data.rare, ord, color = "Group",shape="Group")
p <- p + theme_bw() + ggtitle("MDS + Jaccard")+ scale_color_manual(values = family.palette)
p <- p +geom_point(size=3)+ stat_ellipse(aes(group = Group), size = 1)
Magali Monnoye's avatar
Magali Monnoye committed
1335
1336
1337
1338
plot(p)

```

magali's avatar
magali committed
1339
* Adonis sur la matrice Jaccard 
Magali Monnoye's avatar
Magali Monnoye committed
1340

magali's avatar
magali committed
1341
```{r 35-adonis-jac-1, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1342
1343
1344
1345
1346
adonis(dist.jac ~ Group,data = as(sample_data(frogs.data.rare), 'data.frame'))
```

#### Plot ordination Unifrac

magali's avatar
magali committed
1347
1348

```{r 30-ord-plot-unifrac,fig.width=10,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1349
1350
ord <- ordinate(frogs.data.rare, method = "MDS", distance = "unifrac")

magali's avatar
magali committed
1351
1352
1353
p <- plot_ordination(frogs.data.rare, ord, color = "Group", shape="Group")
p <- p + theme_bw() + ggtitle("MDS + Unifrac")+ scale_color_manual(values = family.palette)
p <- p +geom_point(size=3)+ stat_ellipse(aes(group = Group), size = 1)
Magali Monnoye's avatar
Magali Monnoye committed
1354
1355
1356
1357
1358

plot(p)

```

magali's avatar
magali committed
1359
* Adonis sur la matrice Unifrac
Magali Monnoye's avatar
Magali Monnoye committed
1360

magali's avatar
magali committed
1361
```{r 36-adonis-uni-1, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1362
1363
1364
1365
1366
adonis(dist.uf ~ Group,data = as(sample_data(frogs.data.rare), 'data.frame'))
```

#### Plot ordination Weighted Unifrac

magali's avatar
magali committed
1367
```{r 31-ord-plot-Wunifrac,fig.width=10,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1368
1369
ord <- ordinate(frogs.data.rare, method = "MDS", distance = "wUnifrac")

magali's avatar
magali committed
1370
1371
1372
p <- plot_ordination(frogs.data.rare, ord, color = "Group", shape="Group")
p <- p + theme_bw() + ggtitle("MDS + Weighted Unifrac")+ scale_color_manual(values = family.palette)
p <- p +geom_point(size=3)+ stat_ellipse(aes(group = Group), size = 1)
Magali Monnoye's avatar
Magali Monnoye committed
1373
1374
1375
1376
1377

plot(p)

```

magali's avatar
magali committed
1378
* Adonis sur la matrice wUnifrac
Magali Monnoye's avatar
Magali Monnoye committed
1379

magali's avatar
magali committed
1380
```{r 37-adonis-wuni-1, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1381
1382
1383
adonis(dist.wuf ~ Group,data = as(sample_data(frogs.data.rare), 'data.frame'))
```

magali's avatar
magali committed
1384
### Clustering des échantillons {.tabset} 
Magali Monnoye's avatar
Magali Monnoye committed
1385
1386
1387

Clustering des échantillons pour vérifier s'ils se regroupent par groupe (ou autre). 

magali's avatar
magali committed
1388
1389
<div class="alert alert-info" role="alert">
* **Ward** : cest la plus courante. Elle consiste à réunir les deux clusters dont le regroupement fera le moins baisser linertie interclasse. Cest la distance de Ward qui est utilisée : la distance entre deux classes est celle de leurs barycentres au carré, pondérée par les effectifs des deux clusters. 
Magali Monnoye's avatar
Magali Monnoye committed
1390
1391
Cette technique tend à regrouper les petites classes entre elles.

magali's avatar
magali committed
1392
* **Median** : moyenne de tous les liens, quils soient entre observations de deux clusters différents ou intraclasses. Cette méthode est la seule qui sattache directement au cluster obtenu et non aux caractéristiques des clusters candidats.
Magali Monnoye's avatar
Magali Monnoye committed
1393

magali's avatar
magali committed
1394
1395
* **Average** : le logiciel mesure tous les liens entre chaque observation du cluster A et chaque observation du cluster B et en fait une moyenne. Cest une des méthodes les plus efficaces. Elle tend à réunir des clusters aux inerties faibles.
</div>
Magali Monnoye's avatar
Magali Monnoye committed
1396

magali's avatar
magali committed
1397
#### Clustering sur la matrice de Bray Curtis
Magali Monnoye's avatar
Magali Monnoye committed
1398

magali's avatar
magali committed
1399
```{r 33-plot-clust-bray curtis,fig.width=10,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1400

magali's avatar
magali committed
1401
plot_clust(frogs.data.rare, dist = dist.bc, method = "ward.D2", color ="Group", palette = family.palette)
Magali Monnoye's avatar
Magali Monnoye committed
1402
1403
1404
1405
1406
plot_clust(frogs.data.rare, dist = dist.bc, method = "ward.D", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.bc, method = "median", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.bc, method = "average", color = "Group", palette = family.palette)
```

magali's avatar
magali committed
1407
#### Clustering sur la matrice de Jaccard
Magali Monnoye's avatar
Magali Monnoye committed
1408

magali's avatar
magali committed
1409
```{r 34-plot-clust-jaccard,fig.width=10,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1410
1411
1412
1413
1414
1415
plot_clust(frogs.data.rare, dist = dist.jac, method = "ward.D2", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.jac, method = "ward.D", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.jac, method = "median", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.jac, method = "average", color = "Group", palette = family.palette)
```

magali's avatar
magali committed
1416
#### Clustering sur la matrice Unifrac 
Magali Monnoye's avatar
Magali Monnoye committed
1417

magali's avatar
magali committed
1418
```{r 35-plot-clust-unifrac,fig.width=10,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1419
1420
1421
1422
1423
1424
plot_clust(frogs.data.rare, dist = dist.uf, method = "ward.D2", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.uf, method = "ward.D", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.uf, method = "median", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.uf, method = "average", color = "Group", palette = family.palette)
```

magali's avatar
magali committed
1425
#### Clustering sur la matrice Weighted Unifrac 
Magali Monnoye's avatar
Magali Monnoye committed
1426

magali's avatar
magali committed
1427
```{r 36-plot-clust-wunifrac,fig.width=10,fig.height=5, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1428
1429
1430
1431
1432
1433
plot_clust(frogs.data.rare, dist = dist.wuf, method = "ward.D2", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.wuf, method = "ward.D", color="Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.wuf, method = "median", color = "Group", palette = family.palette)
plot_clust(frogs.data.rare, dist = dist.wuf, method = "average", color = "Group", palette = family.palette)
```

magali's avatar
magali committed
1434
### Network {.tabset}
Magali Monnoye's avatar
Magali Monnoye committed
1435

magali's avatar
magali committed
1436
#### Network Bray curtis
Magali Monnoye's avatar
Magali Monnoye committed
1437

magali's avatar
magali committed
1438
<div class="alert alert-success" role="alert">Network representation of microbiomes, sample-wise or taxa-wise, based on a user-defined ecological distance and (potentially arbitrary) threshold.
Magali Monnoye's avatar
Magali Monnoye committed
1439

magali's avatar
magali committed
1440
Les réseaux sont représentés avec une **distance max = 0.4** 
Magali Monnoye's avatar
Magali Monnoye committed
1441

magali's avatar
magali committed
1442
Default 0.4. The maximum ecological distance (as defined by distance) allowed between two samples to still consider them connected by an edge in the graphical model</div>
Magali Monnoye's avatar
Magali Monnoye committed
1443

magali's avatar
magali committed
1444
```{r 37-network bray,fig.width=8,fig.height=8, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1445

magali's avatar
magali committed
1446
1447
1448
1449
g <- make_network(frogs.data.rare, distance = dist.bc, max.dist = 0.4, keep.isolates = TRUE)
p <- plot_network(g, physeq = frogs.data.rare, color = "Group", label = "value", hjust = 2, title = "Beta diversity network")
p <- p+ scale_color_manual(values = family.palette)
p
Magali Monnoye's avatar
Magali Monnoye committed
1450
1451
```

magali's avatar
magali committed
1452
#### Network Jaccard
Magali Monnoye's avatar
Magali Monnoye committed
1453

magali's avatar
magali committed
1454
```{r 38-network jaccard,fig.width=8,fig.height=8, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1455

magali's avatar
magali committed
1456
1457
1458
1459
g <- make_network(frogs.data.rare, distance = dist.jac, max.dist = 0.4, keep.isolates = TRUE)
p <- plot_network(g, physeq = frogs.data.rare, color = "Group", label = "value", hjust = 2, title = "Jaccard diversity network")
p <- p+ scale_color_manual(values = family.palette)
p
Magali Monnoye's avatar
Magali Monnoye committed
1460
1461
```

magali's avatar
magali committed
1462
#### Network Unifrac
Magali Monnoye's avatar
Magali Monnoye committed
1463

magali's avatar
magali committed
1464
```{r 39-network unifrac,fig.width=8,fig.height=8, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1465

magali's avatar
magali committed
1466
1467
1468
1469
g <- make_network(frogs.data.rare, distance = dist.uf, max.dist = 0.4, keep.isolates = TRUE)
p <- plot_network(g, physeq = frogs.data.rare, color = "Group", label = "value", hjust = 2, title = "Unifrac diversity network")
p <- p+ scale_color_manual(values = family.palette)
p
Magali Monnoye's avatar
Magali Monnoye committed
1470
1471
```

magali's avatar
magali committed
1472
#### Network Weighted Unifrac
Magali Monnoye's avatar
Magali Monnoye committed
1473

magali's avatar
magali committed
1474
```{r 40-network Wunifrac,fig.width=8,fig.height=8, eval=TRUE}
Magali Monnoye's avatar
Magali Monnoye committed
1475