---
title: "Evaluation of an adapted semi-automated DNA extraction for human salivary shotgun metagenomics"
author:
- Victoria Meslier*
- Elisa Menozzi*
- Aymeric David*
- Christian Morabito
- Sara Lucas Del Pozo
- Alexandre Famechon
- Janet North
- Benoit Quinquis
- Sofia Koletsi
- Jane MacNaughtan
- Roxana Mezabrovschi
- S. Dusko Ehrlich
- Anthony HV. Schapira$
- Mathieu Almeida$
output:
  html_document: 
    fig_caption: yes
    highlight: kate
    number_sections: yes
    theme: journal
    toc: yes
    toc_depth: 4
    toc_float: yes
    df_print: kable
subtitle: Plots for paper 2023
---

This document presents the codes and plots used for the analysis related to the manuscript **Evaluation of an adapted semi-automated DNA extraction for human salivary shotgun metagenomics**. In this project, we evaluated the impact of two DNA extraction procedures on the overall DNA quantity and quality, the sequencing performances and the overall structure and composition of the oral microbiome. Salivary samples were obtained from 14 patients recruited in the RAPSODI study and the PD FRONTLINE platform, using the OMNIgene•ORAL OM-501 sample collection kit (DNAGenotek), and processed using P1 and P2 protocols as described in **Figure 1** on the manuscript. The libraries construction and shotgun sequencing using Ion Torrent technology was fully described in **PMID: 36369227**, and the raw read processing in **PMID: 32075887**. 


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning=FALSE, message=FALSE)
```

# Load
- Load necessary R packages :

```{r, message=FALSE}
library(ggpubr)
suppressMessages(library(stringr))
suppressMessages(library(tidyverse))
library(reshape2)
require(ggplot2)
require(pheatmap)
library(data.table)
library(Hmisc)
library(kableExtra)
library(scales)
require(pheatmap)
library(data.table)
library(vegan)
library(FactoMineR)
library(factoextra)
library(ComplexHeatmap, lib.loc = "/usr/lib64/R/library")
library(circlize)
# Install the momr package available forgemia.inra.fr qand published in Le Chatelier E, Nielsen T, Qin J, Prifti E, Hildebrand F, Falony G, et al. Richness of human gut microbiome correlates with metabolic markers. Nature. 2013;500(7464):541–6
# remotes::install_gitlab("metagenopolis/momr", host = "https://forgemia.inra.fr/") 
library(momr)
```

- Load data and metadata : 

```{r}
load("./data/microbial_data_ASAL_paper_2023.RData") #RData containing microbial table abundances and taxonomy file
load(file ="./data/metadata_ASAL_2023.RData") # Associated ASAL metadata
```


# Impact on DNA quantity and DNA quality
## DNA quantity and quality 
- protocol:

```{r, fig.height=3, fig.width=8}
Qdna_plot <- ggplot(metadata.asal, aes(x=protocol, y = DNA_MTQ_ug, color = protocol))+
  geom_boxplot()+
  theme_classic() +
  stat_compare_means(method = "wilcox.test", label.x =0.8, label.y= 1.5) +
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  ylab("DNA quantity (µg)")+
  theme(legend.text = element_text(size=12), axis.text=element_text(size=10), axis.title.x = element_blank())
genomic_plot <- ggplot(metadata.asal, aes(x=protocol, y = genomic_peak_size_bp, color = protocol))+
  geom_boxplot()+
  theme_classic() +
  stat_compare_means(method = "wilcox.test", label.x =1, label.y = 30000) +
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  ylab("Genomic peak size(bp)")+
  scale_y_continuous(labels = function(x) format(x, scientific = TRUE))+
  theme(legend.text = element_text(size=12), axis.text=element_text(size=10), axis.title.x = element_blank())
# tiff("./figures/DNAquantity_DNA_quality_protocol.tiff", width = 8, height = 3, units = "in", res=400)
ggarrange(Qdna_plot, genomic_plot, ncol=2, nrow=1, legend = "none")
# dev.off()
```

> No significant difference for DNA quantity. The size of the genomic peak is significantly higher in the P1 protocol. 


- protocol and sample type :

```{r, fig.height=5, fig.width=8}
Qdna_plot_pt <- ggplot(metadata.asal, aes(x=protocol_type, y = DNA_MTQ_ug, color = protocol_type))+
  geom_boxplot()+
  theme_classic() +
  stat_compare_means(method = "wilcox.test", comparisons = list(c("P1_raw", "P1_pooled"), c("P2_raw", "P2_pooled"), c("P1_raw", "P2_raw"), c("P1_pooled", "P2_pooled")),exact = FALSE)+
  scale_color_manual(values = c("P1_raw" = "cyan4", "P1_pooled" = "cyan3", "P2_raw" = "hotpink4", "P2_pooled" = "hotpink3"))+
  ylab("DNA quantity (µg)")+
  rotate_x_text(angle=45)+
  theme(legend.text = element_text(size=12), axis.text=element_text(size=10), axis.title.x = element_blank())
# tiff("./figures/DNAquantity_protocol_type.tiff", width = 8, height = 5, units = "in", res=400)
Qdna_plot_pt
# dev.off()
```

> There is no impact on the DNA quantity 


```{r, fig.height=5, fig.width=8}
genomic_plot_pt <- ggplot(metadata.asal, aes(x=protocol_type, y = genomic_peak_size_bp, color = protocol_type))+
  geom_boxplot()+
  theme_classic() +
  stat_compare_means(method = "wilcox.test", comparisons = list(c("P1_raw", "P1_pooled"), c("P2_raw", "P2_pooled"), c("P1_raw", "P2_raw"), c("P1_pooled", "P2_pooled")), exact = FALSE)+
  scale_color_manual(values = c("P1_raw" = "cyan4", "P1_pooled" = "cyan3", "P2_raw" = "hotpink4", "P2_pooled" = "hotpink3"))+
  ylab("Mean size genomic peak (bp)")+
  rotate_x_text(angle=45)+
  theme(legend.text = element_text(size=12), axis.text=element_text(size=10), axis.title.x = element_blank())+
  scale_y_continuous(labels = function(x) format(x, scientific = TRUE))
# tiff("./figures/DNAgenomic_quality_protocol_type.tiff", width = 8, height = 5, units = "in", res=400)
genomic_plot_pt
# dev.off()
```

> P1_raw and P1_pooled are equivalent. The DNA quality of P2_raw is the lowest. 


## DNA degradation profiles (protocol)

```{r,fig.height=3, fig.width=12}
f1_plot <- ggplot(metadata.asal, aes(x=protocol, y = as.numeric(fraction1_perc), color = protocol))+
  geom_boxplot()+
  theme_classic() +
  stat_compare_means(label.x =1, label.y = 55) +
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  ylab("fraction1_perc")+
  rotate_x_text(angle=45)+
  theme(legend.text = element_text(size=12), axis.text=element_text(size=10), axis.title.x = element_blank())+
  theme(legend.position= "none")
f2_plot <- ggplot(metadata.asal, aes(x=protocol, y = as.numeric(fraction2_perc), color = protocol))+
  geom_boxplot()+
  theme_classic() +
  stat_compare_means(label.x =1, label.y = 28) +
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  ylab("fraction2_perc")+
  rotate_x_text(angle=45)+
  theme(legend.text = element_text(size=12), axis.text=element_text(size=10), axis.title.x = element_blank())+
  theme(legend.position= "none")
f3_plot <- ggplot(metadata.asal, aes(x=protocol, y = as.numeric(fraction3_perc), color = protocol))+
  geom_boxplot()+
  theme_classic() +
  stat_compare_means(label.x =1, label.y = 100) +
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  ylab("fraction3_perc")+
  rotate_x_text(angle=45)+
  theme(legend.text = element_text(size=12), axis.text=element_text(size=10), axis.title.x = element_blank())
# tiff("./figures/DNA_degradation_profiles_per_protocol.tiff", width = 12, height = 3, units = "in", res=400)
ggarrange(f1_plot, f2_plot, f3_plot, ncol=3, nrow=1, widths = c(1,1,1.4))
# dev.off()
```


> No significance for fraction1 (the smaller DNA fraction size) and fraction3 (containing the genomic peak) but a significantly higher proportion of the DNA fraction 2 (degraded fraction) in the P2 protocol.  


# Shotgun Metagenomic Sequencing Validation 
## Total read count and removed host read count (protocol)
```{r, fig.height=3.5, fig.width=7}
data <- reshape2::melt(metadata.asal[,c("protocol", "total_read_count", "removed_host_reads_count")], id.vars = c("protocol"))
# tiff("./figures/total_read_host_read_protocol.tiff", width = 7, height = 3.5, units = "in", res=400)
ggboxplot(data, x="protocol", y="value", color="protocol")+
  stat_compare_means()+
  theme(axis.text=element_text(size=10))+
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  ylab("read values")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme(legend.text = element_text(size=7), legend.position = "none")+
  facet_wrap(~variable, scales = "free_y")
# dev.off()
```

> No significant differences. 


## High Quality read count

```{r, fig.height=3, fig.width=4}
# tiff("./figures/High_quality_read_protocol.tiff", width = 4, height = 3, units = "in", res=400)
HQ_read_plot <- ggboxplot(metadata.asal, x="protocol", y="HQ_clean_read_count", color="protocol")+
  theme_classic()+
  theme(axis.text=element_text(size=10))+
  ylab("High Quality read count")+
  stat_compare_means(label.x = 1, label.y = 1.5e7)+
  theme(legend.text = element_text(size=7))+
  theme(legend.position = "none")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))
HQ_read_plot
# dev.off()
```

> No significant difference. 


## Mapping rates onto the 8.4M human oral microbial genes catalog

```{r, fig.height=3, fig.width=4}
# tiff("./figures/perc_mapped_reads_procotol.tiff", width = 4, height = 3, units = "in", res=400)
perc_mapped <- ggboxplot(metadata.asal, x="protocol", y="perc_mapped_read_8.4oral", color="protocol")+
  theme_classic()+
  theme(axis.text=element_text(size=10))+
  ylab("% Mapped Reads")+
  stat_compare_means(label.x = 1, label.y = 65)+
  theme(legend.text = element_text(size=7))+
  theme(legend.position = "none")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))
perc_mapped
# dev.off()
```

> No significant difference. 

# Evaluation of the impact on alpha diversity metrics (gene and MSP species richness)
## Correlation of gene and MSP species richness with the DNA quality 
```{r fig.height=3, fig.width=8, message=FALSE, warning=FALSE}
correl_gene_plot <- ggscatter(metadata.asal,
          x="genomic_peak_size_bp",
          y="gene.richness",
          color = "protocol",
          add="reg.line",
          conf.int=TRUE,
          shape = "type",
          add.params=list(color="black", fill="lightgray")) +
  geom_smooth(formula = 'y ~ x', method = "lm")+
  stat_cor(method="spearman")+
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  theme(legend.position = "right")
correl_msp_plot <- ggscatter(metadata.asal,
          x="genomic_peak_size_bp",
          y="msp.richness",
          color = "protocol",
          add="reg.line",
          conf.int=TRUE,
          shape = "type",
          add.params=list(color="black", fill="lightgray")) +
  geom_smooth(formula = 'y ~ x', method = "lm")+
  stat_cor(method="spearman")+
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  theme(legend.position = "right")

# tiff("./figures/correlation_gene_species_richness_genomic_peak_protocol.tiff", width = 8, height = 3, units = "in", res=400)
ggarrange(correl_gene_plot,correl_msp_plot, ncol=2,nrow=1, common.legend = TRUE, legend = "right")
# dev.off()
```


## Impact of the protocol on richness data of the pooled
-all samples (raw and pooled):

```{r, fig.height=3, fig.width=6}
gene_plot <- ggboxplot(metadata.asal, x="protocol", y="gene.richness", color="protocol")+
  stat_compare_means()+
  theme(axis.text=element_text(size=10))+
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  ylab("gene richness")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme(legend.text = element_text(size=7), legend.position = "none")+
  ggtitle("all samples")

msp_plot <- ggboxplot(metadata.asal, x="protocol", y="msp.richness", color="protocol")+
  stat_compare_means()+
  theme(axis.text=element_text(size=10))+
  scale_color_manual(values = c("P1" = "cyan4", "P2" = "hotpink4"))+
  ylab("MSP richness")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme(legend.text = element_text(size=7), legend.position = "none")+
  ggtitle("all samples")

# tiff("./figures/gene_msp_richness_protocol.tiff", width = 6, height = 3, units = "in", res=400)
ggarrange(gene_plot,msp_plot, ncol=2, nrow=1,  legend = "none")
# dev.off()
```


> No significant difference. 


-pooled samples: 

```{r, fig.height=3, fig.width=6}
pooled_gene_plot <- ggboxplot(metadata.asal[which(metadata.asal$type == "pooled"),], x="protocol_type", y="gene.richness", color="protocol_type")+
  stat_compare_means()+
  theme(axis.text=element_text(size=10))+
  scale_color_manual(values = c("P1_pooled" = "cyan3", "P2_pooled" = "hotpink3"))+
  ylab("gene richness")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme(legend.text = element_text(size=7), legend.position = "none")
pooled_msp_plot <- ggboxplot(metadata.asal[which(metadata.asal$type == "pooled"),], x="protocol_type", y="msp.richness", color="protocol_type")+
  stat_compare_means()+
  theme(axis.text=element_text(size=10))+
  scale_color_manual(values = c("P1_pooled" = "cyan3", "P2_pooled" = "hotpink3"))+
  ylab("species richness")+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))+
  theme(legend.text = element_text(size=7), legend.position = "none")

# tiff("./figures/gene_species_richness_protocol_pooled_samples.tiff", width = 6, height = 3, units = "in", res=400)
ggarrange(pooled_gene_plot,pooled_msp_plot, ncol=2,nrow=1, legend = "none")
# dev.off()
```


> No significant difference. 


# Impact on the oral microbiome structure
## Spearman correlation based on gene composition

```{r, fig.height=5, fig.width=8}
ann <- data.frame(metadata.asal$protocol, metadata.asal$Pool_sample)
colnames(ann) <- c('protocol', 'Pool_sample')

colours <- list('protocol' = c("P1" = "cyan4", "P2" = "hotpink4"),'Pool_sample' = c("ind1"="azure2", "ind2"="firebrick1", "ind3"="grey", "ind4"="orange", "ind5"="darkorchid1", "ind6"="blue", "ind7"="darkgreen", "ind8"="dodgerblue2", "ind9/ind12"="gold1", "ind10/ind13"="chartreuse3", "ind11/ind14"="hotpink"))

colAnn <- HeatmapAnnotation(df = ann, which = 'col', col = colours, annotation_width = unit(c(1, 4), 'cm'), gap = unit(1, 'mm'))

# tiff(file = "./figures/Heatmap_gene_correlation.tiff", height = 5, width = 8, units = "in", res = 300)
ht <- Heatmap(as.matrix(gene.correlation.mat),
        clustering_method_rows = "ward.D",clustering_method_columns = "ward.D",
        clustering_distance_rows = "spearman", clustering_distance_columns = "spearman",
        name = "Sp correlation",
        col = colorRamp2(seq(min(gene.correlation.mat), max(gene.correlation.mat), length = 3), c("red", "yellow", "white")),
        column_names_gp = gpar(fontsize = 7),
        row_names_gp = gpar(fontsize = 7),
        top_annotation=colAnn)
ht
# dev.off()
```

> Pooled samples of 2 individuals extracted with P1 and P2 protocols clustered together.

## PCA based on MSP species composition
```{r}
#Ordering msp.mat as metadata
msp.mat<-msp.oral.asal[,match(rownames(metadata.asal), colnames(msp.oral.asal))]
#setting rownames with taxonomy
rownames(msp.mat)<-paste(rownames(msp.mat), taxo[match(rownames(msp.mat), rownames(taxo)), "species"], sep="_")
# PCA analysis
data_for_pca <- t(log10(msp.mat + 1e-10) + 10)
data<- cbind(data_for_pca, metadata.asal[match(rownames(data_for_pca), rownames(metadata.asal)),c(3:5)])
res<-PCA(data, quali.sup = c(854:856), graph=FALSE)
```

```{r, fig.height=4.5, fig.width=5.5}
# tiff(file = "./figures/PCA_MSP_species_protocol.tiff", height = 4.5, width = 5.5, units = "in", res = 300)
pca_plot <- fviz_pca_ind(res,	axes= c(1,2),
             # label = "all",
             label = "none",
             col.ind = data$protocol,
             geom.ind = "point", pointshape = 18, pointsize = 2,
             palette = c("P1" = "cyan4", "P2" = "hotpink4"),
             addEllipses = TRUE, ellipse.level=0.95,
             invisible="quali",
             # title="PCA on MSP species composition",
             repel = TRUE)
pca_plot
# dev.off()
```

> Pooled samples extracted from P1 and P2 protocols clustered together. No specific clustering by P1 and P2 protocols. 


## Speaman correlations based on MSP composition
```{r}
msp.mat = msp.oral.asal[,metadata.asal$sampleID]
data = Hmisc::rcorr(as.matrix(msp.mat), type = "spearman")
data_to_plot = data$r
# tiff(file = "./figures/Heatmap_MSP_correlation.tiff",width = 8, height = 6,  units = "in", res = 300)
Heatmap(as.matrix(data_to_plot), 
        clustering_method_rows = "ward.D",clustering_method_columns = "ward.D",
        clustering_distance_rows = "spearman", clustering_distance_columns = "spearman",
        name = "Sp correlation",
        col = colorRamp2(seq(min(data_to_plot), max(data_to_plot), length = 3), c("red", "yellow", "white")),
        column_names_gp = gpar(fontsize = 7), 
        row_names_gp = gpar(fontsize = 6),
        top_annotation=colAnn)
# dev.off()
```

## Bray curtis dissimilarity matrix using MSP species composition
```{r, fig.height=6, fig.width=8}
data = t(log10(msp.oral.asal + 1e-10) + 10)
BC = as.matrix(vegdist(data, method = "bray", diag = TRUE, upper= FALSE, na.rm = FALSE))

ann_colors = list(
  protocol = c(P1="cyan4", P2="hotpink4"),
  'Pool_sample' = c("ind1"="azure2", "ind2"="firebrick1", "ind3"="grey", "ind4"="orange", "ind5"="darkorchid1", "ind6"="blue", "ind7"="darkgreen", "ind8"="dodgerblue2", "ind9/ind12"="gold1", "ind10/ind13"="chartreuse3", "ind11/ind14"="hotpink")
  )

# tiff(file = "./figures/Bray-curtis_MSP_composition.tiff",width = 8, height = 6,  units = "in", res = 300)
pheatmap(as.matrix(BC), 
         clustering_distance_rows = "euclidean",
         clustering_distance_cols = "euclidean",
         # annotation_col=metadata.asal[-which(metadata.asal$fraction == "degraded"),c(3:4,9)],
         annotation_col=metadata.asal[,c("protocol", "Pool_sample")],
         annotation_colors=ann_colors,
         clustering_method = "ward.D2",         
         fontsize_row=5, fontsize_col=5,
         show_colnames = F, show_rownames = F,
         dendrogram="both", key=FALSE, treeheight_row=40,
         treeheight_col=40
         ) 
# dev.off()
```



## Analysis at different taxonomical ranks on the pooled samples
### phylum on pooled samples
```{r, fig.height=10, fig.width = 5.5}
phylum.t <- as.data.frame(t(phylum)) #15 phyla
phylum.t <- phylum.t[,colSums(phylum.t)>0] #12 phyla ok
phylum.t$protocol_type <- metadata.asal[match(rownames(phylum.t), rownames(metadata.asal)), "protocol_type"]
data <- melt(phylum.t, id.vars = c("protocol_type"))
data$protocol_type <- factor(data$protocol_type, levels = c("P2_raw", "P2_pooled","P1_raw", "P1_pooled"))

# tiff(file = "./figures/Phylum_composition_pooled_samples_protocol.tiff", height = 10, width = 5.5, units = "in", res = 300)
phylum_plot <- ggplot(data[which(data$protocol_type == "P1_pooled" | data$protocol_type == "P2_pooled"),], aes(x=reorder(variable, value, na.rm =TRUE), y=value)) +
  geom_boxplot(aes(fill=protocol_type), outlier.size = 0.1)+
  labs(x = "Phylum rank")+
  labs(y="Relative abundance")+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x),labels = trans_format("log10", math_format(10^.x))) +
  scale_fill_manual(values = c("P1_pooled" = "cyan3", "P2_pooled" = "hotpink3"))+
  coord_flip()+
  theme_classic()
phylum_plot
# dev.off()
```

```{r}
# impact of the protocol on the phylum relative abundance of the pooled samples
res.phylum <- testRelations2(x = as.matrix(phylum[,which(metadata.asal$type == "pooled")]), y = as.factor(metadata.asal[which(metadata.asal$type == "pooled"), "protocol"]), type = "nonparametric", multiple.adjust = "BH", paired = F)
res.phylum[,c(10:13)]%>%
  kbl(caption = "pvalues Wilcoxon tests") %>%
  kable_classic(full_width = F, html_font = "Cambria")
```

> No significant differences between pooled samples extracted with P1 or P2 protocols (p> 0.05).  


### family on pooled samples
```{r, fig.height=10, fig.width = 5.5}
fam.t <- as.data.frame(t(family)) #80
fam.t <- fam.t[rownames(metadata.asal),colSums(fam.t)>0] #55 fam remains
fam.t$protocol_type <- metadata.asal[match(rownames(fam.t), rownames(metadata.asal)), "protocol_type"]
data <- melt(fam.t, id.vars = c("protocol_type"))
data$protocol_type <- factor(data$protocol_type, levels = c("P2_raw", "P2_pooled","P1_raw", "P1_pooled"))
data2 <- data[which(data$protocol_type == "P1_pooled" | data$protocol_type == "P2_pooled"),]
# tiff(file = "./figures/Family_composition_pooled_samples_protocol.tiff", height = 10, width = 5.5, units = "in", res = 300)
family_plot <- ggplot(data2[data2$value>0,], aes(x=reorder(variable, value, na.rm =TRUE), y=value)) +
  geom_boxplot(aes(fill=protocol_type), outlier.size = 0.1)+
  labs(x = "Family rank")+
  labs(y="Relative abundance")+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), labels = trans_format("log10", math_format(10^.x)))+
  scale_fill_manual(values = c("P1_pooled" = "cyan3", "P2_pooled" = "hotpink3"))+
  coord_flip()+
  theme_classic()
family_plot
# dev.off()
```

```{r}
# impact of the protocol on the family relative abundance of the pooled samples
res.fam <- testRelations2(x = as.matrix(family[,which(metadata.asal$type == "pooled")]), y = as.factor(metadata.asal[which(metadata.asal$type == "pooled"), "protocol"]), type = "nonparametric", multiple.adjust = "BH", paired = F)
res.fam[,c(10:13)]%>%
  kbl(caption = "pvalues Wilcoxon tests") %>%
  kable_classic(full_width = F, html_font = "Cambria")
```

> No significant differences between pooled samples extracted with P1 or P2 protocols (p> 0.05). 


### genus on pooled sample
```{r, fig.height=10, fig.width = 5.5}
gen.t <- as.data.frame(t(genus)) #163
gen.t <- gen.t[rownames(metadata.asal),colSums(gen.t)>0] #86 gen remains
gen.t$protocol_type <- metadata.asal[match(rownames(gen.t), rownames(metadata.asal)), "protocol_type"]
data <- melt(gen.t, id.vars = c("protocol_type"))
data$protocol_type <- factor(data$protocol_type, levels = c("P2_raw", "P2_pooled","P1_raw", "P1_pooled"))
data3 <- data[which(data$protocol_type == "P1_pooled" | data$protocol_type == "P2_pooled"),]
# tiff(file = "./figures/Genus_composition_pooled_samples_protocol.tiff", height = 10, width = 5.5, units = "in", res = 300)
genus_plot <- ggplot(data3[data3$value>0,], aes(x=reorder(variable, value, na.rm =TRUE), y=value)) +
  geom_boxplot(aes(fill=protocol_type), outlier.size = 0.1)+
  labs(x = "Genus rank")+
  labs(y="Relative abundance")+
  scale_y_log10(breaks = trans_breaks("log10", function(x) 10^x), labels = trans_format("log10", math_format(10^.x)))+
  scale_fill_manual(values = c("P1_pooled" = "cyan3", "P2_pooled" = "hotpink3"))+
  coord_flip()+
  theme_classic()
genus_plot
# dev.off()
```

```{r}
# impact of the protocol on the genus relative abundance of the pooled samples
res.genus <- testRelations2(x = as.matrix(genus[,which(metadata.asal$type == "pooled")]), y = as.factor(metadata.asal[which(metadata.asal$type == "pooled"), "protocol"]), type = "nonparametric", multiple.adjust = "BH", paired = F)
res.genus[,c(10:13)]%>%
  kbl(caption = "pvalues Wilcoxon tests") %>%
  kable_classic(full_width = F, html_font = "Cambria")
```


> No significant differences between pooled samples extracted with P1 or P2 protocols (p> 0.05). 


### MSP species on pooled samples
```{r}
# impact of the protocol on the MSP species relative abundance of the pooled samples
res.msp <- testRelations2(x = as.matrix(msp.oral.asal[,which(metadata.asal$type == "pooled")]), y = as.factor(metadata.asal[which(metadata.asal$type == "pooled"), "protocol"]), type = "nonparametric", multiple.adjust = "BH", paired = F)
res.msp[,c(10:13)]%>%
  kbl(caption = "pvalues Wilcoxon tests") %>%
  kable_classic(full_width = F, html_font = "Cambria")
```

> No significant differences. Trends observed for two MSP species (p<0.1).

```{r, echo = FALSE, include=FALSE}
sessionInfo()
```


