MAD Plate AR

library(DspikeIn)

## Warning: replacing previous import 'ggplot2::alpha' by 'microbiome::alpha' when
## loading 'microbiomeutilities'

## 
## Thank you for using **DspikeIn**! 🎉
## For support,📧 contact: Mitra Ghotbi (mitra.ghotbi@gmail.com)
##  GitHub Repository: https://github.com/mghotbi/DspikeIn 
## 🐛 Found a bug or have a suggestion? Open an issue on GitHub: 
##    https://github.com/mghotbi/DspikeIn/issues

data("metadata_full", package="DspikeIn")

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.3

library(ggrepel)

## Warning: package 'ggrepel' was built under R version 4.3.3

library(patchwork)

## Warning: package 'patchwork' was built under R version 4.3.3

library(DspikeIn)

color_palette <- DspikeIn::color_palette$MG

filtered_data <- metadata_full %>%
  filter(plate.ID == "AlexRurik_MAD_plate1",
         Animal.type %in% c("Turtle", "Snake")) %>%
  group_by(Animal.type) %>%
  slice_head(n = 1)  

Lets say 10000 reads can be min and 2.5 good shanon index

acceptable_reads <- 10000  
acceptable_shannon <- 2.5   

Lets plot it

p1 <- ggplot(filtered_data, aes(x = Animal.type, y = Total_Reads_total, fill = Animal.type)) +
  geom_boxplot(alpha = 0.6, outlier.shape = NA) +
  geom_jitter(aes(size = Total_Reads_spiked), width = 0.2, alpha = 0.7, color = "black") +
  scale_y_log10() + 
  scale_fill_manual(values = color_palette) +
  geom_hline(yintercept = acceptable_reads, linetype = "dashed", color = "red2", linewidth = 1) +
  annotate("text", x = 1.5, y = acceptable_reads * 1.2, 
           label = "Minimum Acceptable Reads", color = "red2", fontface = "italic") +
  labs(title = "Sample Quality: Total Reads",
       y = "Total Reads (log scale)",
       x = "Animal Type",
       fill = "Animal Type") +
  theme_minimal() +
  theme(legend.position = "none")

#  Richness vs. Diversity
p2 <- ggplot(filtered_data, aes(x = Observed, y = Shannon, color = Animal.type)) +
  geom_point(aes(size = Total_Reads_total), alpha = 0.8) +
  geom_text_repel(aes(label = sample.id), size = 4, box.padding = 0.5) +
  scale_color_manual(values = color_palette) +
  scale_size_continuous(range = c(3, 8)) +
  geom_hline(yintercept = acceptable_shannon, linetype = "dashed", color = "blue2", linewidth = 1) +
  annotate("text", x = max(filtered_data$Observed) * 0.8, y = acceptable_shannon + 0.2, 
           label = "Minimum Diversity Threshold", color = "blue2", fontface = "italic") +
  labs(title = "Richness vs. Diversity",
       x = "Observed Richness",
       y = "Shannon Diversity Index",
       color = "Animal Type",
       size = "Total Reads") +
  theme_minimal() +
  theme(legend.position = "bottom")

Patch it

final_plot <- p1 + p2 + plot_annotation(title = "Assessing Sample Quality for Re-Sequencing")

print(final_plot)

Or you can plot it this way

library(circlize)

## Warning: package 'circlize' was built under R version 4.3.3

## ========================================
## circlize version 0.4.16
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
## 
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
##   in R. Bioinformatics 2014.
## 
## This message can be suppressed by:
##   suppressPackageStartupMessages(library(circlize))
## ========================================

sector_colors <- setNames(DspikeIn::color_palette$MG, 
                          unique(c(filtered_data$Animal.type, 
                                   "Total Reads", "Spiked Reads")))

# chord diagram ( Total Reads & Spiked Reads)
chord_data <- data.frame(
  from = rep(filtered_data$Animal.type, 2),
  to = c(rep("Total Reads", nrow(filtered_data)), 
         rep("Spiked Reads", nrow(filtered_data))),
  value = c(filtered_data$Total_Reads_total, 
            filtered_data$Total_Reads_spiked)
)
circos.clear()
circos.par(start.degree = 90, gap.degree = 5, track.margin = c(0.05, 0.05))

# chord diagram
chordDiagram(chord_data, 
             grid.col = sector_colors, 
             transparency = 0.3, 
             annotationTrack = "grid", 
             preAllocateTracks = list(track.height = 0.1))

circos.track(track.index = 1, panel.fun = function(x, y) {
  sector.name <- get.cell.meta.data("sector.index")
  xlim <- get.cell.meta.data("xlim")
  ylim <- get.cell.meta.data("ylim")

  circos.text(x = mean(xlim), 
              y = ylim[1] + 0.1,  
              labels = sector.name, 
              facing = "inside", 
              cex = 1.4,  
              niceFacing = TRUE, 
              font = 2)  
}, bg.border = NA)

circos.clear()