Simulations with spiked E. coli strain analysed using Sylph Profile

Load data and prepare for simulations

Fit ZIB models

profile 95

save_path = "output_rds/spike_ecoli_zib_p_95.rds"
if(file.exists(save_path)){
  fit_zib_95 <- readRDS(save_path)
} else {
  system.time({
    fit_zib_95 <- glmZiBFit(d_p95, design, nthreads = 10)
  })
  saveRDS(fit_zib_95, save_path)
}

plot_manhattan(fit_zib_95, taxonomy = tax_95, method = "HMP", 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.

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_volcano(fit_zib_95, label = T)

Found 487 tophits for spikedTRUE at alpha = 1 using holm 

plot_ani_dist(d_p95, phenotype = 'spiked', contigs = top_hits(fit_zib_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 310 rows containing non-finite outside the scale range
(`stat_ydensity()`).

profile 98

save_path = "output_rds/spike_ecoli_zib_p_98.rds"
if(file.exists(save_path)){
  fit_zib_98 <- readRDS(save_path)
} else {
  system.time({
    fit_zib_98 <- glmZiBFit(d_p98, design, nthreads = 10)
  })
  saveRDS(fit_zib_98, save_path)
}

plot_manhattan(fit_zib_98, taxonomy = tax_98, method = "holm", 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_volcano(fit_zib_98, label = T)

Found 662 tophits for spikedTRUE at alpha = 1 using holm 

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

Found 3 tophits for spikedTRUE at alpha = 0.05 using holm 

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

profile 99

save_path = "output_rds/spike_ecoli_zib_p_99.rds"
if(file.exists(save_path)){
  fit_zib_99 <- readRDS(save_path)
} else {
  system.time({
    fit_zib_99 <- glmZiBFit(d_p99, design, nthreads = 10)
  })
  saveRDS(fit_zib_99, save_path)
}

plot_manhattan(fit_zib_99, taxonomy = tax_99, method = "HMP", 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_volcano(fit_zib_99, label = T)

Found 403 tophits for spikedTRUE at alpha = 1 using holm 

plot_ani_dist(d_p99, phenotype = 'spiked', contigs = top_hits(fit_zib_99, alpha = 0.05)$Contig_name)

Found 4 tophits for spikedTRUE at alpha = 0.05 using holm 

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

Fit OB models

profile 95

save_path = "output_rds/spike_ecoli_ob_p_95.rds"
if(file.exists(save_path)){
  fit_ob_95 <- readRDS(save_path)
} else {
  system.time({
    fit_ob_95 <- glmObFit(d_p95, design, nthreads = 10)
  })
  saveRDS(fit_ob_95, save_path)
}

plot_manhattan(fit_ob_95, taxonomy = tax_95, method = "HMP", 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_volcano(fit_ob_95, label = T)

Found 487 tophits for spikedTRUE at alpha = 1 using holm 

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

Found 3 tophits for spikedTRUE at alpha = 0.5 using holm 

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

profile 98

save_path = "output_rds/spike_ecoli_ob_p_98.rds"
if(file.exists(save_path)){
  fit_ob_98 <- readRDS(save_path)
} else {
  system.time({
    fit_ob_98 <- glmObFit(d_p98, design, nthreads = 10)
  })
  saveRDS(fit_ob_98, save_path)
}

plot_manhattan(fit_ob_98, taxonomy = tax_98, method = "HMP", 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_volcano(fit_ob_98, label = T)

Found 662 tophits for spikedTRUE at alpha = 1 using holm 

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

Found 2 tophits for spikedTRUE at alpha = 0.05 using holm 

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

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

Found 2 tophits for spikedTRUE at alpha = 0.05 using holm 

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

Warning: Cannot compute density for groups with fewer than two datapoints.

profile 99

save_path = “output_rds/spike_ecoli_ob_p_99.rds” if(file.exists(save_path)){ fit_ob_99 <- readRDS(save_path) } else { system.time({ fit_ob_99 <- glmObFit(d_p99, design, nthreads = 10) }) saveRDS(fit_ob_99, save_path) }

plot_manhattan(fit_ob_99, taxonomy = tax_99, method = “HMP”, tax_levels = c(“Phylum”, “Order”, “Class”, “Genus”, “Species”)) plot_volcano(fit_ob_99, label = T) plot_ani_dist(d_p99, phenotype = ‘spiked’, contigs = top_hits(fit_ob_99, alpha = 0.05)\(Contig_name, show_points = T, plot_type = 'box') plot_ani_dist(d_p99, phenotype = 'spiked', contigs = top_hits(fit_ob_99, alpha = 0.05)\)Contig_name, show_points = T, plot_type = ‘violin’)