Commit f3331084 authored by Magali Monnoye's avatar Magali Monnoye
Browse files

fichier rmd

parent 3c919345
Pipeline #44522 passed with stage
in 3 minutes and 37 seconds
---
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
Du microbiote d'enfants atteints de troubles du spectre autistique a été donné à des souris axéniques, afin de voir l'impact de ce transfert sur leur comportement mais aussi sur divers paramètres biochimiques
Chez les donneurs, il y avait 4 groupes composés chacun de 4 enfants
Un pool des 4 échantillons d'un même groupe a été fait pour le donner aux souris
# 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>
<div class="alert alert-success" role="alert">En résumé, à part pour l'échantillon 2.1, les séquences sont de très bonnes qualités</div>
```{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"
Job 4046552
```
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"
Job 4046554
```
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
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-10-19.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"
Job 3983975
```
Lien vers le [rapport affiliations_stats_cleaned](FROGS/affiliations_stats_cleaned.html)
<div class="alert comment">Suite aux différents outils FROGS utilisés pour le nettoyage des séquences, l'échantillon 113pq2 n'a plus que 1250 reads alors que tous les autres échantillons ont un nombre de reads > 15000, cet échantillon devra être supprimé lors de la rarefaction pour éviter de rarefier sur trop peu de séquences</div>
### 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
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-10-19.biom --out-tree FROGS/tree.nhx --html FROGS/tree.html --log-file FROGS/tree.log && conda deactivate"
Job 3983983
```
Lien vers l'arbre [Tree](FROGS/tree.html)
# Metagenomic phyloseq analysis
## 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{
frogsBiom <- "cleaned_biom-2021-10-19.biom"
frogs.data.biom <- import_frogs(frogsBiom, taxMethod = "blast")
metadata<-read.table("metadatasBALBc.txt",row.names=1, header=T)
sample_data(frogs.data.biom)<-metadata
sample_names(frogs.data.biom) <- get_variable(frogs.data.biom, "Name")
frogs.data<-frogs.data.biom
treefile<- read.tree("tree.nhx")
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)
```
Composition des groupes:
* AG : enfants atteints de TSA et ayant aussi des troubles gastro-intestinaux
* S-AG: leurs frères et sœurs comme contrôles
* A: Enfants atteints de TSA sans troubles gastro-intestinaux
* S-A: leurs frères et sœurs comme contrôles
Pour vérifier l'implantation du microbiote et voir les éventuelles différences entre les groupes, les fèces ont été prélévées à 3 et 6 semaines après le transfert et une partie du contenu cæcal a été gardée pour faire l'analyse 16S.
* Nombre d'échantillons par groupe:
```{r variable-1, eval = TRUE, echo=TRUE}
table(get_variable(frogs.data, "Group"))
```
* Nombre d'échantillons par statut:
```{r variable-2, eval = TRUE, echo=TRUE}
table(get_variable(frogs.data, "Statut"))
```
* Nombre d'échantillons par cohorte:
```{r variable-3, eval = TRUE, echo=TRUE}
table(get_variable(frogs.data, "Cohorte"))
```
## Sélection des données
```{r sample selection 1,fig.width=10,fig.height=4, eval=TRUE, echo = TRUE}
frogs.data <- subset_samples(frogs.data, Cohorte == "Cohorte_1")
frogs.data
```
## 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
```
#### Tableaux relatives abondances Phylum
* Relatives Abondances par échantillons:
```{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',
pageLength = 20,
buttons = list('copy', 'csv', 'excel')))
```
#### Plot composition Phylum
```{r 05-composition-phylum, eval = TRUE, fig.height=15, fig.width=15}
correct.order <- c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec" )
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")
p<- p + facet_wrap(~Group, scales = "free_x", nrow = 3) +
theme(strip.text.x = element_text(size = 14, color = "black")) +
scale_y_continuous(label = scales::percent)
plot(p)
```
```{r 06-composition-phylum-merged, eval = TRUE, fig.height=10, fig.width=10}
frogs.data.merged<-merge_samples(frogs.data,group="Group")
sample_data(frogs.data.merged)$Group <- sample_names(frogs.data.merged)
correct.order <- c("S-A_F3w","A_F3w","S-AG_F3w","AG_F3w","S-A_F6w","A_F6w","S-AG_F6w","AG_F6w", "S-A_Caec", "A_Caec", "S-AG_Caec", "AG_Caec")
sample_data(frogs.data.merged)$Group <- factor(sample_data(frogs.data.merged)$Group,