Comparing bioinformatics pipelines

16S short reads

1 Import dataset

asv_table <- read.csv("../data/asv_table.csv",sep=";")
asv_table[is.na(asv_table)] <- 0
asv_table[,-1] <- apply(asv_table[,-1],2,function(x) x/sum(x))

taxa_table <- read.csv("../data/taxa_table.csv",sep=";")
colors <- taxa_table$Color
colors[colors==""] <- NA
asv_table$Color <- colors
taxa_table <- taxa_table[,-9]

metadata <- read.csv("../data/metadata.csv",sep=";")
metadata$SettingID <- factor(metadata$SettingID,
                             levels=(c("Mock_reference",
                                      paste0("DA",1:48),
                                      paste0("VU",1:32),
                                      paste0("DE",1:32))))

source("custom_functions.R")

2 Stacked-barlot

df_p <- melt(asv_table) %>% 
  `colnames<-`(c("SeqID","color", "sample","value")) %>%
  dplyr::mutate(SettingID=sapply(strsplit(as.character(sample),"_"),function (x) x[1]))

Using SeqID, Color as id variables

df_p$SettingID <- gsub("^Mock$","Reference",df_p$SettingID)
df_p$SettingID <-  factor(df_p$SettingID,
                             levels= (c("Reference",
                                      paste0("DA",1:48),
                                      paste0("VU",1:32),
                                      paste0("DE",1:32))))


df_p %<>% dplyr::group_by(SettingID,SeqID) %>%
  summarize(
    value=mean(value),
    color=unique(color)
  )

`summarise()` has grouped output by 'SettingID'. You can override using the
`.groups` argument.

colors <- df_p$color
names(colors) <- df_p$SeqID

p <- ggplot(df_p, aes(fill=SeqID, y=value, x=SettingID)) + 
    geom_bar(position="stack", stat="identity",width=0.8) +
  theme_minimal() + 
  theme(panel.border = element_rect(color = "black", fill = NA, size = 0),
        axis.text.x = element_text(size = 8,color="black"),
        axis.text.y = element_text(size = 10,color="black"),
        axis.ticks.x = element_line(size=0.3,color="black"),
        axis.ticks.y = element_line(size=0.3,color="black"),
        axis.ticks.length = unit(4,"pt"),
        panel.grid = element_blank()) + 
  scale_fill_manual(values=colors) + 
  ylab("Relative abundance") + 
  xlab("Pipeline") + 
  theme(legend.position = "none") +
  scale_x_discrete(guide = guide_axis(angle = 90)) 

Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.

Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.

p

3 PCoA

# Taxonomic classifier
p_pcoa_classifier <- pca_plot_custom(asv_table,taxa_table,metadata)


p_pcoa_replicate <- pca_plot_custom(asv_table,taxa_table,metadata,variable="Replicate")


p_pcoa_prep <- pca_plot_custom(asv_table,taxa_table,metadata,variable="Preprocessing")


p_pcoa_tool <- pca_plot_custom(asv_table,taxa_table,metadata,variable="Tool")


p_pcoa <- ggarrange(p_pcoa_classifier,p_pcoa_replicate,p_pcoa_tool,p_pcoa_prep, ncol=2,nrow=2)

p_pcoa

4 Bray-Curtis

p_bray <- bray_custom(asv_table,metadata) 

`summarise()` has grouped output by 'prep_and_tool'. You can override using the
`.groups` argument.

p_bray 

5 Performance metrics

metadata <- precision_recall(asv_table,metadata)

p_performance <- performance_plot(metadata)

Using SettingID as id variables

p_performance