Simulations with spiked E. coli strain analysed using Sylph Query

Load data and prepare for simulations

Fit ZIB models

query 95

save_path = "output_rds/spike_ecoli_zib_q_95.rds"
if(file.exists(save_path)){
  fit_zib_95 <- readRDS(save_path)
} else {
  system.time({
    fit_zib_95 <- glmZiBFit(d_q95, design, nthreads = parallel::detectCores())
  })
  saveRDS(fit_zib_95, save_path)
}

plot_manhattan(fit_zib_95, taxonomy = tax_95, method = "BH")

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
ℹ The deprecated feature was likely used in the strainspy package.
  Please report the issue at <https://github.com/gtonkinhill/strainspy/issues>.

plot_manhattan(fit_zib_95, taxonomy = tax_95, aggregate_by_taxa = F)

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_volcano(fit_zib_95, label = T)

Found 1737 tophits for spikedTRUE at alpha = 1 using holm 

plot_ani_dist(d_q95, phenotype = 'spiked', contigs = top_hits(fit_zib_95)$Contig_name, show_points = T, plot_type = 'violin')

Found 2 tophits for spikedTRUE at alpha = 0.05 using holm 

Warning: Removed 218 rows containing non-finite outside the scale range
(`stat_ydensity()`).

query 98

save_path = "output_rds/spike_ecoli_zib_q_98.rds"
if(file.exists(save_path)){
  fit_zib_98 <- readRDS(save_path)
} else {
  system.time({
    fit_zib_98 <- glmZiBFit(d_q98, design, nthreads = parallel::detectCores())
  })
  saveRDS(fit_zib_98, save_path)
}


plot_manhattan(fit_zib_98, taxonomy = tax_98, method = "BH")

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_volcano(fit_zib_98, label = T)

Found 6124 tophits for spikedTRUE at alpha = 1 using holm 

plot_ani_dist(d_q98, phenotype = 'spiked', contigs = top_hits(fit_zib_98, alpha = 0.05)$Contig_name)

Found 4 tophits for spikedTRUE at alpha = 0.05 using holm 

Warning: Removed 282 rows containing non-finite outside the scale range
(`stat_boxplot()`).

query 99

save_path = "output_rds/spike_ecoli_zib_q_99.rds"
if(file.exists(save_path)){
  fit_zib_99 <- readRDS(save_path)
} else {
  system.time({
    fit_zib_99 <- glmZiBFit(d_q99, design, nthreads = parallel::detectCores())
  })
  saveRDS(fit_zib_99, save_path)
}

plot_manhattan(fit_zib_99, taxonomy = tax_99, method = "BH", tax_levels = c("Order", "Phylum", "Genus", "Species"))

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_manhattan(fit_zib_99, taxonomy = tax_99, aggregate_by_taxa = F)

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_volcano(fit_zib_99, label = T)

Found 9320 tophits for spikedTRUE at alpha = 1 using holm 

plot_ani_dist(d_q99, phenotype = 'spiked', contigs = top_hits(fit_zib_99)$Contig_name)

Found 18 tophits for spikedTRUE at alpha = 0.05 using holm 

Warning: Removed 1258 rows containing non-finite outside the scale range
(`stat_boxplot()`).

contigs_to_plot = top_hits(fit_zib_99)$Contig_name

Found 18 tophits for spikedTRUE at alpha = 0.05 using holm 

plot_ani_dist(d_q99, phenotype = 'spiked', contigs = top_hits(fit_zib_99)$Contig_name, plot_type = 'box', show_points = T, contig_names = sub("NZ_", "", unname(sapply(contigs_to_plot, function(x) unlist(strsplit(x," "))[1])) ))

Found 18 tophits for spikedTRUE at alpha = 0.05 using holm 

Warning: Removed 1258 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Fit OB models

query 95

save_path = "output_rds/spike_ecoli_ob_q_95.rds"
if(file.exists(save_path)){
  fit_ob_95 <- readRDS(save_path)
} else {
  system.time({
    fit_ob_95 <- glmObFit(d_q95, design, nthreads = parallel::detectCores())
  })
  saveRDS(fit_ob_95, save_path)
}

plot_manhattan(fit_ob_95, taxonomy = tax_95, method = "BH", tax_levels = c("Phylum", "Order", "Class", "Genus"))

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_manhattan(fit_ob_95, taxonomy = tax_95, aggregate_by_taxa = F)

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_volcano(fit_ob_95, label = T)

Found 1737 tophits for spikedTRUE at alpha = 1 using holm 

plot_ani_dist(d_q95, phenotype = 'spiked', contigs = top_hits(fit_ob_95, alpha = 0.05)$Contig_name, show_points = T, plot_type = 'violin')

Found 2 tophits for spikedTRUE at alpha = 0.05 using holm 

Warning: Removed 218 rows containing non-finite outside the scale range
(`stat_ydensity()`).

query 98

save_path = "output_rds/spike_ecoli_ob_q_98.rds"
if(file.exists(save_path)){
  fit_ob_98 <- readRDS(save_path)
} else {
  system.time({
    fit_ob_98 <- glmObFit(d_q98, design, nthreads = parallel::detectCores())
  })
  saveRDS(fit_ob_98, save_path)
}


plot_manhattan(fit_ob_98, taxonomy = tax_98, method = "BH", tax_levels = c("Phylum", "Order", "Class", "Genus", "Species"))

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_manhattan(fit_ob_98, taxonomy = tax_98, aggregate_by_taxa = F)

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_volcano(fit_ob_98, label = T)

Found 6124 tophits for spikedTRUE at alpha = 1 using holm 

plot_ani_dist(d_q98, phenotype = 'spiked', contigs = top_hits(fit_ob_98, alpha = 0.05)$Contig_name, show_points = T, plot_type = 'violin')

Found 4 tophits for spikedTRUE at alpha = 0.05 using holm 

Warning: Removed 282 rows containing non-finite outside the scale range
(`stat_ydensity()`).

query 99

save_path = "output_rds/spike_ecoli_ob_q_99.rds"
if(file.exists(save_path)){
  fit_ob_99 <- readRDS(save_path)
} else {
  system.time({
    fit_ob_99 <- glmObFit(d_q99, design, nthreads = parallel::detectCores())
  })
  saveRDS(fit_ob_99, save_path)
}

plot_manhattan(fit_ob_99, taxonomy = tax_99, tax_levels = c("Phylum", "Order", "Class", "Genus", "Species"))

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_manhattan(fit_ob_99, taxonomy = tax_99, aggregate_by_taxa = F)

! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.
! # Invaild edge matrix for <phylo>. A <tbl_df> is returned.

plot_volcano(fit_ob_99, label = T)

Found 9320 tophits for spikedTRUE at alpha = 1 using holm 

plot_ani_dist(d_q99, phenotype = 'spiked', contigs = top_hits(fit_ob_99, alpha = 0.05)$Contig_name, show_points = T, plot_type = 'violin')

Found 18 tophits for spikedTRUE at alpha = 0.05 using holm 

Warning: Removed 1258 rows containing non-finite outside the scale range
(`stat_ydensity()`).