# =============================================================================
# Unified Sampling Across European Food Processing Environments Reveals
# Complex Genetic and Contextual Factors Underlying Listeria monocytogenes
# Persistence
# =============================================================================
# Author: Felix Spiegel
# Affiliation: Austrian Competence Centre for Feed and Food Quality, Safety, and Innovation - FFoQSI GmbH
# Date: March 2026
# Description: Data visualisation and Statistical analysis supporting manuscript findings
# =============================================================================
## -------------------------------------------------------------------------------------------------------------------------------------
library(ape)
library(Biostrings)
library(circlize)
library(ComplexHeatmap)
library(dplyr)
library(ggplot2)
library(ggrepel)
library(ggnewscale)
library(gridExtra)
library(grid)
library(gtable)
library(ggtree)
library(kableExtra)
library(knitr)
library(patchwork)
library(phangorn)
library(pheatmap)
library(readxl)
library(tidyverse)
library(treeio)
library(viridis)
library(writexl)
## -------------------------------------------------------------------------------------------------------------------------------------
## directory
metadata <- read_excel("~/Desktop/AllPartners_FoodSafeR_WP2.3_Facility_Sampling_Metadata_FS.xlsx")
# Display basic info about the dataset
head(metadata)
nrow(metadata)
ncol(metadata)
colnames(metadata)
##### Count unique ST and CT ########
# Count unique ST and CT
unique_ST <- length(unique(na.omit(metadata$ST)))
unique_CT <- length(unique(na.omit(metadata$CT)))
unique_ST
unique_CT
# now build table
final_table <- metadata %>%
filter(!is.na(Sequence_ID_Lmono), Sequence_ID_Lmono != "", Sequence_ID_Lmono != "NA") %>%
group_by(
Meta_Facility_Code,
FoodSafeR_Sample_ID,
Sequence_ID_Lmono,
ST,
CT,
Sample_Location,
Sample_Zone,
Connected_Samples,
Sampling_Timepoint,
pre_post_Disinfection,
Sampling_Date
) %>%
summarise(
Sequence_ID_Lmono = paste(Sequence_ID_Lmono, collapse = ", "),
.groups = "drop"
)
# save result
write_xlsx(final_table, "~/Desktop/Listeria_Isolates_Summary_01.xlsx")
## -------------------------------------------------------------------------------------------------------------------------------------
sample_count <- nrow(metadata)
count_mono <- sum(metadata$PCR_Lmono_Y_N == "Y", na.rm = TRUE)
count_FC <- sum(metadata$`Food_Contact/Non-Contact_Surface` == "FC", na.rm = TRUE)
count_NFC <- sum(metadata$`Food_Contact/Non-Contact_Surface` == "NFC", na.rm = TRUE)
count_pre <- sum(metadata$pre_post_Disinfection == "pre", na.rm = TRUE)
count_post <- sum(metadata$pre_post_Disinfection == "post", na.rm = TRUE)
count_na <- sum(metadata$pre_post_Disinfection == "NA", na.rm = TRUE)
count_meat <- sum(metadata$Facility_Type == "Meat", na.rm = TRUE)
count_fish <- sum(metadata$Facility_Type == "Fish", na.rm = TRUE)
count_dairy <- sum(metadata$Facility_Type == "Dairy", na.rm = TRUE)
count_austria <- sum(metadata$Facility_Country == "Austria", na.rm = TRUE)
count_greece <- sum(metadata$Facility_Country == "Greece", na.rm = TRUE)
count_spain <- sum(metadata$Facility_Country == "Spain", na.rm = TRUE)
# Build new dataframe
final_df <- data.frame(
Subcategory = c(
"", # Total Samples (no subcategory name)
"Austria", "Greece", "Spain", # Facility Country
"Meat", "Fish", "Dairy", # Facility Type
"FC", "NFC", # Food Contact Surface
"pre", "post", "NA", # Pre/Post Disinfection
"Positive", "Negative" # L. monocytogenes
),
Count = c(
sample_count,
count_austria, count_greece, count_spain,
count_meat, count_fish, count_dairy,
count_FC, count_NFC,
count_pre, count_post, count_na,
count_mono, sample_count - count_mono
)
)
# calculate percentages and ad them to the dataframe
final_df$Percentage <- round((final_df$Count / sample_count) * 100, 1)
final_df$Percentage[1] <- "" # Total samples — no percentage
# make the table
final_table <- kable(
final_df,
format = "html",
align = "c",
row.names = FALSE,
col.names = c("Subcategory", "Count", "Percentage"),
caption = "General Sampling Information
"
) %>%
kable_styling(full_width = FALSE, position = "center", font_size = 14) %>%
group_rows("Total Samples", 1, 1) %>%
group_rows("Facility Country", 2, 4) %>%
group_rows("Facility Type", 5, 7) %>%
group_rows("Food Contact Surface", 8, 9) %>%
group_rows("Pre/Post Disinfection", 10, 12) %>%
group_rows("L. monocytogenes", 13, 14)
# look at the table
final_table
# Save the table
save_kable(final_table, file = "~/Desktop/general_sampling_info.html")
# Count co-occurence of L. monocytogenes with L. spp
count_both_isolates <- sum(
metadata$Lineage != "NA" &
metadata$Listeria_spp_Isolate_Ids != "NA",
na.rm = TRUE
)
count_both_isolates
# Subset rows where both Lineage and Listeria_spp_Isolate_Ids are not "NA"
valid_rows <- metadata[
metadata$Lineage != "NA" &
metadata$Listeria_spp_Isolate_Ids != "NA",
]
# Count Persistent vs Non-persistent
persistent_counts <- table(valid_rows$Persistent)
persistent_counts
## -------------------------------------------------------------------------------------------------------------------------------------
# Define helper function for stacked positive/negative plots
make_stacked_plot <- function(data, group_var, title, colors, order_levels = NULL) {
df <- data %>%
group_by({{ group_var }}) %>%
summarise(
Positive = sum(PCR_Lmono_Y_N == "Y", na.rm = TRUE),
Negative = sum(PCR_Lmono_Y_N == "N", na.rm = TRUE),
.groups = "drop"
) %>%
pivot_longer(cols = c("Positive","Negative"),
names_to = "Result", values_to = "Count") %>%
group_by({{ group_var }}) %>%
mutate(Percentage = round(Count / sum(Count) * 100, 1))
if (!is.null(order_levels)) {
df[[deparse(substitute(group_var))]] <- factor(df[[deparse(substitute(group_var))]], levels = order_levels)
}
ggplot(df, aes(x = {{ group_var }}, y = Percentage, fill = Result)) +
geom_bar(stat = "identity", colour = "black", width = 0.8) + # wider bars
geom_text(aes(label = Percentage), position = position_stack(vjust = 0.5), size = 5) +
scale_fill_manual(values = colors) +
scale_y_continuous(limits = c(0, 100)) + # y-axis 0–100
theme_minimal() +
labs(title = title, x = NULL, y = "% positive Swab Samples") +
theme(
legend.position = "bottom",
legend.title = element_blank(),
legend.text = element_text(size = 14, face = "bold"), # bigger legend text
plot.title = element_text(hjust = 0.5, face = "bold", size = 20), # bigger, bold title
axis.title.x = element_text(size = 16, face = "bold"), # bigger x-axis label
axis.title.y = element_text(size = 16, face = "bold"), # bigger y-axis label
axis.text = element_text(size = 14, face = "bold") # bigger axis ticks
)
}
# Overall L. monocytogenes
count_mono <- sum(metadata$PCR_Lmono_Y_N == "Y", na.rm = TRUE)
sample_count <- nrow(metadata)
df_lmono <- data.frame(
Category = "L. monocytogenes",
Result = c("Positive", "Negative"),
Count = c(count_mono, sample_count - count_mono)
) %>%
mutate(Percentage = round(Count / sum(Count) * 100, 1))
p_lmono <- ggplot(df_lmono, aes(x = Category, y = Percentage, fill = Result)) +
geom_bar(stat = "identity", colour = "black", width = 0.8) +
geom_text(aes(label = Percentage), position = position_stack(vjust = 0.5), size = 5) +
scale_fill_manual(values = c("Positive" = "#238A8D", "Negative" = "#DCE319")) +
scale_y_continuous(limits = c(0, 100)) +
theme_minimal() +
labs(title = "Overall", x = NULL, y = "% positive Swab Samples") +
theme(
legend.position = "bottom",
legend.title = element_blank(),
legend.text = element_text(size = 14, face = "bold"), # bigger legend text
plot.title = element_text(hjust = 0.5, face = "bold", size = 20), # bigger, bold title
axis.title.x = element_text(size = 16, face = "bold"), # bigger x-axis label
axis.title.y = element_text(size = 16, face = "bold"), # bigger y-axis label
axis.text = element_text(size = 14, face = "bold") # bigger axis ticks
)
# Facility Type (Meat, Fish, Dairy)
p_facility <- make_stacked_plot(metadata, Facility_Type, "Facility Type",
c("Positive" = "#238A8D", "Negative" = "#DCE319"),
order_levels = c("Meat", "Fish", "Dairy"))
# Food Contact Surface (FC, NFC, NA)
p_surface <- make_stacked_plot(metadata, `Food_Contact/Non-Contact_Surface`,
"Surface Type",
c("Positive" = "#238A8D", "Negative" = "#DCE319"),
order_levels = c("FC", "NFC", "NA"))
# Disinfection State (pre, post, NA)
p_disinf <- make_stacked_plot(metadata, pre_post_Disinfection,
"Disinfection State",
c("Positive" = "#238A8D", "Negative" = "#DCE319"),
order_levels = c("pre", "post", "NA"))
# Helper to extract legend
get_legend <- function(myplot) {
tmp <- ggplotGrob(myplot)
leg <- gtable::gtable_filter(tmp, "guide-box")
return(leg)
}
# Remove legends from individual plots
p_lmono_nolegend <- p_lmono + theme(legend.position = "none")
p_facility_nolegend <- p_facility + theme(legend.position = "none")
p_surface_nolegend <- p_surface + theme(legend.position = "none")
p_disinf_nolegend <- p_disinf + theme(legend.position = "none")
# Extract shared legend from one plot
shared_legend <- get_legend(p_lmono)
# Arrange plots with shared legend and add a header
final_plot <- arrangeGrob(
arrangeGrob(p_lmono_nolegend, p_facility_nolegend, p_surface_nolegend, p_disinf_nolegend, ncol = 2),
shared_legend,
nrow = 2,
heights = c(10, 1),
top = textGrob("L. monocytogenes Positive Swab Samples Rates", gp = gpar(fontsize = 16, fontface = "bold"))
)
# Specify SVG file path
svg_file <- "~/Desktop/lmono_plot.svg"
# Open SVG device
svg(svg_file, width = 12, height = 8) # width/height in inches
# Draw the plot
grid.draw(final_plot)
# Close device
dev.off()
## -------------------------------------------------------------------------------------------------------------------------------------
# Helper function to compute counts, percentages, and add sample info
make_barplot_table <- function(data, group_var, group_name, reorder_levels = NULL) {
data <- data %>%
filter(!is.na({{ group_var }}))
if(!is.null(reorder_levels)) {
data <- data %>%
mutate({{ group_var }} := factor({{ group_var }}, levels = reorder_levels))
}
data %>%
group_by({{ group_var }}) %>%
summarise(
Positive = sum(PCR_Lmono_Y_N == "Y", na.rm = TRUE),
Negative = sum(PCR_Lmono_Y_N == "N", na.rm = TRUE),
.groups = "drop"
) %>%
mutate(
Total = Positive + Negative,
Perc_Positive = round(Positive / Total * 100, 1),
Perc_Negative = round(Negative / Total * 100, 1),
Sample_Percentage = round(Total / nrow(data) * 100, 1), # % of total samples
Category = group_name,
Group = as.character({{ group_var }})
) %>%
select(Category, Group, Total, Sample_Percentage, Positive, Negative, Perc_Positive, Perc_Negative)
}
# Overall (single line)
overall_table <- data.frame(
Category = "Overall",
Group = "All samples",
Positive = sum(metadata$PCR_Lmono_Y_N == "Y", na.rm = TRUE),
Negative = sum(metadata$PCR_Lmono_Y_N == "N", na.rm = TRUE)
) %>%
mutate(
Total = Positive + Negative,
Sample_Percentage = 100, # overall is 100%
Perc_Positive = round(Positive / Total * 100, 1),
Perc_Negative = round(Negative / Total * 100, 1)
) %>%
select(Category, Group, Total, Sample_Percentage, Positive, Negative, Perc_Positive, Perc_Negative)
# Facility Country
facility_country_table <- make_barplot_table(metadata, Facility_Country, "Facility Country")
# Facility Type
facility_table <- make_barplot_table(metadata, Facility_Type, "Facility Type")
# Surface Type
surface_table <- make_barplot_table(metadata, `Food_Contact/Non-Contact_Surface`, "Surface Type")
# Disinfection State (reorder: pre, post, NA)
disinf_table <- make_barplot_table(metadata, pre_post_Disinfection, "Disinfection State",
reorder_levels = c("pre", "post", NA))
# Combine all tables
barplot_summary_table <- bind_rows(
overall_table,
facility_country_table,
facility_table,
surface_table,
disinf_table
)
# Blank first column for all rows (only headers shown via group_rows)
barplot_summary_table$Category <- ""
# Determine row indices for grouping
rows_overall <- 1
rows_country <- (rows_overall + 1):(rows_overall + nrow(facility_country_table))
rows_facility <- (max(rows_country) + 1):(max(rows_country) + nrow(facility_table))
rows_surface <- (max(rows_facility) + 1):(max(rows_facility) + nrow(surface_table))
rows_disinf <- (max(rows_surface) + 1):(max(rows_surface) + nrow(disinf_table))
# Create HTML table
final_table <- kable(
barplot_summary_table,
format = "html",
align = "c",
row.names = FALSE,
col.names = c("", "Group", "Total Samples", "% of Total Samples", "Positive", "Negative", "% Positive", "% Negative"),
caption = "General Sampling Information
"
) %>%
kable_styling(full_width = FALSE, position = "center", font_size = 14) %>%
group_rows("Overall", rows_overall, rows_overall) %>%
group_rows("Facility Country", min(rows_country), max(rows_country)) %>%
group_rows("Facility Type", min(rows_facility), max(rows_facility)) %>%
group_rows("Surface Type", min(rows_surface), max(rows_surface)) %>%
group_rows("Disinfection State", min(rows_disinf), max(rows_disinf))
# View table
final_table
# Save as HTML
save_kable(final_table, file = "~/Desktop/general_sampling_info_combined.html")
## -------------------------------------------------------------------------------------------------------------------------------------
summary_table <- metadata %>%
filter(
!is.na(Sequence_ID_Lmono), # Exclude real NA
Sequence_ID_Lmono != "NA", # Exclude the literal string "NA"
!(
(is.na(Lineage) | Lineage == "NA") &
(is.na(ST) | ST == "NA") &
(is.na(CT) | CT == "NA")) # Exclude rows where all 3 are "NA"
) %>%
group_by(Facility_Type, Lineage, ST, CT) %>%
summarise(Count = n(), .groups = 'drop') %>%
arrange(Facility_Type, desc(Count))
# Show result
summary_table
# Total count in summary_table for sanity check
total_detected <- sum(summary_table$Count)
total_detected
# Format the summary_table nicely
formatted_table <- kable(
summary_table,
format = "html",
align = "c",
row.names = FALSE,
col.names = c("Facility Type", "Lineage", "ST", "CT", "Count"),
caption = "L. monocytogenes Typing Results
"
) %>%
kable_styling(full_width = FALSE, position = "center", font_size = 14)
formatted_table
# Save the table
save_kable(formatted_table, file = "~/Desktop/Lin_ST_MLST_info.html")
summary_CT_STonly <- summary_table %>%
group_by(ST) %>%
summarise(
CTs = paste(sort(unique(CT)), collapse = ", "),
Total_Count = sum(Count),
.groups = 'drop'
) %>%
arrange(as.numeric(ST)) # Sort STs numerically
summary_CT_STonly
## -------------------------------------------------------------------------------------------------------------------------------------
# summarize counts
summary_lineage <- summary_table %>%
group_by(Facility_Type, Lineage) %>%
summarise(Total_Count = sum(Count), .groups = 'drop')
# Set factor order for Facility_Type
summary_lineage$Facility_Type <- factor(
summary_lineage$Facility_Type,
levels = c("Meat", "Dairy", "Fish")
)
summary_lineage_plot <- ggplot(summary_lineage, aes(x = Facility_Type, y = Total_Count, fill =Lineage)) +
geom_bar(stat = "identity", color = "black", size = 0.3, width = 0.3)+ # super thin bars
geom_text(aes(label = Total_Count),
position = position_stack(vjust = 0.5),
color = "black",
size = 4) +
labs(
title = "Lineages Isolated in Facility Types",
x = "Facility Type",
y = "Isolate Count",
fill = "Lineage"
) +
theme_minimal(base_size = 10) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold", hjust = 0.5)
) +
scale_fill_brewer(palette = "Viridis")
summary_lineage_plot
ggsave(
filename = "~/Desktop/summary_lineage_plot.svg",
plot = summary_lineage_plot,
width = 8, height = 6, units = "in")
## -------------------------------------------------------------------------------------------------------------------------------------
# summarize counts
summary_ST <- summary_table %>%
group_by(Facility_Type, ST) %>%
summarise(Total_Count = sum(Count), .groups = 'drop')
# Set factor order for Facility_Type
summary_ST$Facility_Type <- factor(
summary_ST$Facility_Type,
levels = c("Meat", "Dairy", "Fish")
)
# Sort ST levels numerically
summary_ST$ST <- factor(summary_ST$ST, levels = sort(as.numeric(unique(summary_ST$ST))))
# Plot
summary_ST_plot <- ggplot(summary_ST, aes(x = Facility_Type, y = Total_Count, fill =ST)) +
geom_bar(stat = "identity", color = "black", size = 0.3, width = 0.3)+ # super thin bars
labs(
title = "Sequence Types Isolated in Facility Types",
x = "Facility Type",
y = "Isolate Count",
fill = "ST"
) +
theme_minimal(base_size = 10) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold", hjust = 0.5)
) +
scale_fill_viridis_d(option = "D")
summary_ST_plot
ggsave( filename = "~/Desktop/summary_ST_plot.svg",
plot = summary_ST_plot,
width = 8, height = 6, units = "in")
## -------------------------------------------------------------------------------------------------------------------------------------
summary_ST_plot2 <- ggplot(summary_ST, aes(x = ST, y = Total_Count, fill = Facility_Type)) +
geom_col(color = "black", width = 0.7, size = 0.3) +
geom_text(
aes(label = Total_Count),
position = position_stack(vjust = 0.5), # center inside each stacked segment
size = 5, # bigger text
fontface = "bold", # bold labels
color = "white" # better contrast on colored bars
) +
labs(
title = "Sequence Types by Facility Type",
x = "Sequence Type",
y = "Isolate Count",
fill = "Facility"
) +
theme_minimal(base_size = 12) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
axis.text.y = element_text(size = 12),
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
legend.position = "top",
legend.title = element_text(face = "bold", size = 12),
legend.text = element_text(size = 11)
) +
scale_fill_manual(
values = c("Meat" = "#404788", "Dairy" = "#55C667", "Fish" = "#FDE725")
)
summary_ST_plot2
ggsave(
filename = "~/Desktop/summary_ST_plot_by_facility.svg",
plot = summary_ST_plot2,
width = 10, height = 6, units = "in"
)
## add lineages
ST_lineage_lookup <- summary_table %>%
distinct(ST, Lineage)
summary_ST2 <- summary_ST %>%
left_join(ST_lineage_lookup, by = "ST") ## add lineage from summary_ST
colnames(summary_ST2)
# "Facility_Type" "ST" "Total_Count" "Lineage"
# order STs nummerically
summary_ST2 <- summary_ST2 %>%
mutate(ST_num = as.numeric(ST)) %>% # numeric version for sorting
arrange(Lineage, ST_num) %>% # sort by Lineage, then numeric ST
mutate(ST = factor(ST, levels = unique(ST))) # keep original labels
# plot
summary_ST_plot2 <- ggplot(summary_ST2,
aes(x = ST, y = Total_Count, fill = Facility_Type)) +
geom_col(color = "black", width = 0.7, size = 0.3) +
geom_text(
aes(label = Total_Count),
position = position_stack(vjust = 0.5),
size = 3, # smaller numbers
fontface = "bold",
color = "white"
) +
facet_wrap(~ Lineage, scales = "free_x") +
labs(
title = "Sequence Types by Facility Type and Lineage",
x = "Sequence Type (ST)",
y = "Number of Isolates",
fill = "Facility Type"
) +
theme_minimal(base_size = 12) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 11),
axis.text.y = element_text(size = 11),
axis.title.x = element_text(face = "bold", size = 14),
axis.title.y = element_text(face = "bold", size = 14),
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
legend.position = "top",
strip.text = element_text(face = "bold", size = 12)
) +
scale_fill_manual(
values = c("Meat" = "#404788FF",
"Dairy" = "#55C667FF",
"Fish" = "#FDE725FF")
)
summary_ST_plot2
ggsave(
filename = "~/Desktop/summary_ST_plot_by_facility_lineage.svg",
plot = summary_ST_plot2,
width = 10, height = 6, units = "in"
)
## -------------------------------------------------------------------------------------------------------------------------------------
# summarize counts
summary_CT <- summary_table %>%
group_by(Facility_Type, CT) %>%
summarise(Total_Count = sum(Count), .groups = 'drop')
# Set factor order for Facility_Type
summary_CT$Facility_Type <- factor(
summary_CT$Facility_Type,
levels = c("Meat", "Dairy", "Fish"))
# Sort ST levels numerically
summary_CT$CT <- factor(summary_CT$CT, levels = sort(as.numeric(unique(summary_CT$CT))))
summary_CT
# Plot
summary_CT_plot <- ggplot(summary_CT, aes(x = Facility_Type, y = Total_Count, fill =CT)) +
geom_bar(stat = "identity", color = "black", size = 0.3, width = 0.3)+ # super thin bars
labs(
title = "cgMLSTs Isolated in Facility Types",
x = "Facility Type",
y = "Isolate Count",
fill = "CT"
) +
theme_minimal(base_size = 10) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold", hjust = 0.5)
) +
scale_fill_viridis_d(option = "D")
summary_CT_plot
ggsave( filename = "~/Desktop/summary_CT_plot.svg",
plot = summary_CT_plot,
width = 8, height = 6, units = "in")
## -------------------------------------------------------------------------------------------------------------------------------------
# Correct the data types to character
metadata <- metadata %>%
mutate(ST = as.character(ST), CT = as.character(CT), Meta_Facility_Code = as.character(Meta_Facility_Code), Sampling_Date =as.Date(Sampling_Date) )
# Iterate through each unique company, ST and CT type
# Initialize an empty list to store results
results_list <- list()
#Make an empty list to store the errors
error_log <- list()
# Define parent directory for all companies
parent_dir <- "~/Desktop/SNP_persister_classification/"
parent_dir
# List all company directories
company_dirs <- list.dirs(parent_dir, recursive = FALSE)
company_dirs
# Loop over each company directory
for (company_dir in company_dirs) {
# Extract company name from directory path
company <- basename(company_dir)
# List all SNP matrix and tree files in the directory
snp_files <- list.files(company_dir, pattern = "snp_distance.*\\.tsv$", full.names = TRUE)
#tree_files <- list.files(company_dir, pattern = "snpma_preserved.*\\.treefile$", full.names = TRUE)
# Process each SNP file
for (snp_file in snp_files) {
tryCatch({
# Extract ST and CT identifiers from file name
st_id <- str_extract(snp_file, "(?<=ST)\\d+")
ct_id <- str_extract(snp_file, "(?<=CT)\\d+")
# Subset metadata for the current ST, CT, and company
subset_metadata <- metadata %>%
filter(ST == st_id, CT == ct_id, Meta_Facility_Code == company)
# Skip if no metadata for this combination
if (nrow(subset_metadata) == 0) {
warning(paste("No matching metadata found for ST:", st_id, "CT:", ct_id, "Company:", company))
next
}
# Load the SNP distance matrix
#snp_matrix <- read_tsv(snp_file, show_col_types = FALSE)
snp_matrix <- read.delim(snp_file, row.names = 1, check.names = FALSE)
# Subset SNP matrix to isolates in the metadata subset
isolate_ids <- subset_metadata$Sequence_ID_Lmono
filtered_snp_matrix <- snp_matrix[isolate_ids, isolate_ids, drop = FALSE]
# Skip if filtered SNP matrix is empty
if (nrow(filtered_snp_matrix) == 0 || ncol(filtered_snp_matrix) == 0) {
warning(paste("Filtered SNP matrix is empty for ST:", st_id, "CT:", ct_id, "Company:", company))
next
}
# Perform clustering based on SNP distances (≤ 10 SNPs)
snp_clusters <- hclust(as.dist(filtered_snp_matrix)) # Hierarchical clustering
cluster_groups <- cutree(snp_clusters, h = 10) # Define clusters with ≤ 10 SNPs
# Assign clusters back to the subset metadata
subset_metadata <- subset_metadata %>%
mutate(Cluster = cluster_groups[match(Sequence_ID_Lmono, names(cluster_groups))])
persistent_clusters <- subset_metadata %>%
group_by(Cluster) %>%
mutate(sampling_date = n_distinct(Sampling_Date)) %>%
filter(sampling_date >= 3) %>%
summarize(Isolates = list(Sequence_ID_Lmono), .groups = "drop")
# Annotate Persistence
subset_metadata <- subset_metadata %>%
mutate(Persistence = ifelse(Sequence_ID_Lmono %in% unlist(persistent_clusters$Isolates),
"Persistent", "Non-Persistent"))
# Save Results for the Current Subset
results_list[[paste0(company, "_ST", st_id, "_CT", ct_id)]] <- subset_metadata
}, error = function(e) {
# Log errors and print them
error_log[[paste0(company, "_ST", st_id, "_CT", ct_id)]] <- e$message
cat("Error in file processing for Company:", company, "ST:", st_id, "CT:", ct_id, "\n")
cat("Error message:", e$message, "\n")
})
}
}
# Check the error log after the loop
if (length(error_log) > 0) {
cat("The following datasets caused errors:\n")
print(names(error_log))
}
#Combine results for all companies, STs, and CTs into a single dataframe
final_metadata <- bind_rows(results_list)
# Save the combined results
write_delim(final_metadata, "AllPartners_FoodSafeR_WP2.3_Facility_Sampling_Metadata_FS_persiter_10SNP_3Y.csv", delim = ";")
#After export of this file, singletons and associated metadata were added manually!
## Get table that summarizes only the persistent clusters
# Ensure relevant columns are character if needed
final_metadata <- final_metadata %>%
mutate(across(c(Lineage, ST, Clonal_Complex, CT, Facility_Type, Meta_Facility_Code, Cluster), as.character))
# Filter for persistent isolates only and summarize clusters
persistent_cluster_summary <- final_metadata %>%
filter(Persistence == "Persistent", !is.na(Cluster)) %>%
group_by(Lineage, ST, Clonal_Complex, CT, Facility_Type, Meta_Facility_Code, Cluster) %>%
summarize(
Num_Isolates = n(),
.groups = "drop"
) %>%
arrange(Lineage, ST, CT, Cluster)
# View the summary table
print(persistent_cluster_summary)
# Optionally save to a CSV file
write_delim(persistent_cluster_summary, "Persistent_Cluster_Summary_Table_FS.csv", delim = ";")
# Reorder and sort the summary
persistent_cluster_summary <- persistent_cluster_summary %>%
select(Meta_Facility_Code, Facility_Type, Cluster, Lineage, ST, Clonal_Complex, CT, Num_Isolates) %>%
arrange(Meta_Facility_Code)
# Create styled table with centered caption
persister_cluster <- persistent_cluster_summary %>%
kable(format = "html",
caption = "Summary of Persistent Clusters
",
align = "c",
col.names = c("Facility Code", "Facility Type", "Cluster", "Lineage", "ST", "Clonal Complex", "CT", "Count"),
position = "center") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE,
font_size = 14)
persister_cluster
# Save the HTML table
save_kable(persister_cluster, file = "~/Desktop/persister_cluster.html")
## -------------------------------------------------------------------------------------------------------------------------------------
# Filter for persistent isolates
persistent_metadata <- final_metadata %>%
filter(Persistence == "Persistent")
# Summarize counts of persistent isolates by ST and Facility_Type
summary_persistent_ST <- persistent_metadata %>%
group_by(ST, Facility_Type) %>%
summarize(Total_Count = n(), .groups = "drop")
# Sort ST numerically for plotting
summary_persistent_ST <- summary_persistent_ST %>%
mutate(ST = factor(ST, levels = sort(as.numeric(unique(ST)))))
# Plot persistent isolates only
persistent_ST_plot <- ggplot(summary_persistent_ST, aes(x = ST, y = Total_Count, fill = Facility_Type)) +
geom_col(color = "black", width = 0.7, size = 0.2) +
geom_text(
aes(label = Total_Count),
position = position_stack(vjust = 0.5),
size = 5,
fontface = "bold",
color = "white"
) +
labs(
title = "Persistent Sequence Types by Facility Type",
x = "Sequence Type (ST)",
y = "Number of Persistent Isolates",
fill = "Facility Type"
) +
theme_minimal(base_size = 12) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
axis.text.y = element_text(size = 12),
axis.title.x = element_text(face = "bold", size = 14), # bigger & bold x-axis label
axis.title.y = element_text(face = "bold", size = 14), # bigger & bold y-axis label
plot.title = element_text(face = "bold", hjust = 0.5, size = 18), # bigger & bolder title
legend.position = "top",
legend.title = element_text(face = "bold", size = 12),
legend.text = element_text(size = 11)
) +
scale_fill_manual(
values = c("Meat" = "#404788FF", "Dairy" = "#55C667FF", "Fish" = "#FDE725FF")
)
# Show the plot
persistent_ST_plot
# Optionally save as SVG
ggsave(
filename = "~/Desktop/persistent_ST_plot_by_facility.svg",
plot = persistent_ST_plot,
width = 10, height = 6, units = "in"
)
##### Add Lineage info to persistent summary
ST_lineage_lookup <- summary_table %>%
distinct(ST, Lineage)
summary_persistent_ST2 <- summary_persistent_ST %>%
left_join(ST_lineage_lookup, by = "ST")
# Order STs numerically within each Lineage
summary_persistent_ST2 <- summary_persistent_ST2 %>%
mutate(ST_num = as.numeric(ST)) %>%
arrange(Lineage, ST_num) %>%
mutate(ST = factor(ST, levels = unique(ST))) # keep ST labels
# Plot persistent isolates faceted by Lineage
persistent_ST_plot2 <- ggplot(summary_persistent_ST2,
aes(x = ST, y = Total_Count, fill = Facility_Type)) +
geom_col(color = "black", width = 0.7, size = 0.3) +
geom_text(
aes(label = Total_Count),
position = position_stack(vjust = 0.5),
size = 3, # smaller text
fontface = "bold",
color = "white"
) +
facet_wrap(~ Lineage, scales = "free_x") +
labs(
title = "Sequence Types of Persisters by Facility Type and Lineage",
x = "Sequence Type (ST)",
y = "Number of Persistent Isolates",
fill = "Facility Type"
) +
theme_minimal(base_size = 12) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1, size = 11),
axis.text.y = element_text(size = 11),
axis.title.x = element_text(face = "bold", size = 14),
axis.title.y = element_text(face = "bold", size = 14),
plot.title = element_text(face = "bold", hjust = 0.5, size = 16),
legend.position = "top",
strip.text = element_text(face = "bold", size = 12)
) +
scale_fill_manual(
values = c("Meat" = "#404788FF",
"Dairy" = "#55C667FF",
"Fish" = "#FDE725FF")
)
# Show the plot
persistent_ST_plot2
# Save as SVG
ggsave(
filename = "~/Desktop/persistent_ST_plot_by_facility_lineage.svg",
plot = persistent_ST_plot2,
width = 10, height = 6, units = "in"
)
## -------------------------------------------------------------------------------------------------------------------------------------
# Filter for confirmed L. monocytogenes isolates
lmono_metadata_PCR <- metadata %>%
filter(PCR_Lmono_Y_N == "Y")
nrow(lmono_metadata_PCR)
colnames(metadata)
# filter the ones sequenced
lmono_metadata <- metadata %>%
filter(
Sequence_ID_Lmono != "NA",
Lineage != "NA")
nrow(lmono_metadata)
#
ct_counts_facility <- lmono_metadata %>%
group_by(Facility_Type, CT) %>%
summarize(CT_Count = n(), .groups = "drop")
singletons_per_facility <- ct_counts_facility %>%
filter(CT_Count <= 2)
lmono_singletons_facility <- lmono_metadata %>%
inner_join(singletons_per_facility, by = c("Facility_Type", "CT"))
nrow(lmono_singletons_facility)
# Ensure relevant columns are character
lmono_metadata <- lmono_metadata %>%
mutate(across(c(Meta_Facility_Code, Facility_Type, Lineage, ST, Clonal_Complex, CT), as.character))
# Summarize persistence by facility, lineage, ST, clonal complex, and CT
persistent_summary <- lmono_metadata %>%
group_by(Meta_Facility_Code, Facility_Type, Lineage, ST, Clonal_Complex, CT) %>%
summarize(
Total_Isolates = n(),
Persistent_Isolates = sum(Persistent == "Persistent", na.rm = TRUE),
Percent_Persistent = Persistent_Isolates / Total_Isolates * 100,
.groups = "drop"
) %>%
# Exclude rows with no persistent isolates
filter(Persistent_Isolates > 0) %>%
arrange(Meta_Facility_Code, Lineage, ST, CT) %>%
mutate(Percent_Persistent = round(Percent_Persistent, 1))
# Calculate overall persistence
overall_summary <- lmono_metadata %>%
summarize(
Meta_Facility_Code = "All",
Facility_Type = "-",
Lineage = "-",
ST = "-",
Clonal_Complex = "-",
CT = "-",
Total_Isolates = n(),
Persistent_Isolates = sum(Persistent == "Persistent", na.rm = TRUE),
Percent_Persistent = round(Persistent_Isolates / Total_Isolates * 100, 1)
)
# Combine overall summary with detailed persistent table
final_table <- bind_rows(overall_summary, persistent_summary)
## add historical isolates
# Create an empty historical row with the same columns and types
historical_row <- final_table[0, ] # 0 rows, keeps columns & types
# Fill in the historical values
historical_row[1, "Meta_Facility_Code"] <- "D (Historical)"
historical_row[1, "Facility_Type"] <- "Dairy"
historical_row[1, "Lineage"] <- "1"
historical_row[1, "ST"] <- "4"
historical_row[1, "Clonal_Complex"] <- "4"
historical_row[1, "CT"] <- "4605"
historical_row[1, "Total_Isolates"] <- 37
historical_row[1, "Persistent_Isolates"]= 37
historical_row[1, "Percent_Persistent"] <- 100
# Append historical row
final_table_df <- bind_rows(final_table, historical_row)
# Now run kable() on final_table_df
persistent_summary_table <- final_table_df %>%
kable(
format = "html",
col.names = c(
"Facility Code", "Facility Type", "Lineage", "ST",
"Clonal Complex", "cgMLST", "Total Isolates",
"Persistent Isolates", "Percent Persistent (%)"
),
align = "c",
escape = FALSE,
na = ""
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed", "responsive"),
full_width = FALSE,
font_size = 14
) %>%
row_spec(nrow(final_table_df), bold = TRUE, color = "grey", background = "white")
#Show and save
persistent_summary_table
save_kable(persistent_summary_table, "~/Desktop/Lmono_Persistent_Metadata_Table.html")
# Check: Categorize isolates by persistence status
# Purpose: verify that all sequenced isolates are accounted for and correctly
# classified as persistent, singleton, isolated twice, or multi-isolate non-persistent
# Summary counts — confirm numbers add up to total sequenced isolates
lmono_metadata <- lmono_metadata %>%
add_count(Meta_Facility_Code, CT, name = "CT_Count") %>%
mutate(
Persistent_Flag = Persistent == "Persistent",
Isolation_Category = case_when(
Persistent_Flag ~ "Persistent",
CT_Count == 1 ~ "Singleton",
CT_Count == 2 ~ "Isolated Twice",
CT_Count >= 3 ~ "Multi-isolate, Non-Persistent",
TRUE ~ NA_character_
)
)
# Summary counts
lmono_metadata %>%
count(Isolation_Category) %>%
mutate(Percent = round(n / sum(n) * 100, 1)) %>%
bind_rows(summarise(., Isolation_Category = "Total", n = sum(n), Percent = sum(Percent)))
# Trace the non-persistent multi-isolate group
# Purpose: understand why some isolates share a CT with persistent isolates
# but are not classified as persistent themselves
# Detailed breakdown of multi-isolate non-persistent isolates
multi_nonpersistent_summary <- lmono_metadata %>%
filter(Isolation_Category == "Multi-isolate, Non-Persistent") %>%
group_by(Meta_Facility_Code, CT, CT_Count) %>%
summarize(n_non_persistent = n(), .groups = "drop") %>%
mutate(n_persistent_in_group = CT_Count - n_non_persistent) %>%
left_join(
persistent_summary %>% select(Meta_Facility_Code, CT, Persistent_Isolates, Total_Isolates),
by = c("Meta_Facility_Code", "CT")
) %>%
mutate(Subtype = ifelse(is.na(Persistent_Isolates),
"No persistent isolates in CT group",
"CT group has persistent isolates"))
multi_nonpersistent_summary
table(multi_nonpersistent_summary$Subtype)
# Isolate-level breakdown
multi_nonpersistent_summary %>%
group_by(Subtype) %>%
summarize(
n_CT_groups = n(),
n_isolates = sum(n_non_persistent),
.groups = "drop"
)
# The 21 non-persistent isolates visible in your table
table_nonpersistent <- persistent_summary %>%
mutate(n_non_persistent = Total_Isolates - Persistent_Isolates) %>%
filter(n_non_persistent > 0) %>%
summarize(total = sum(n_non_persistent))
cat("Non-persistent isolates visible in persistence table:", table_nonpersistent$total, "\n")
# Cross-check: confirm total non-persistent isolates matches expectation
multi_nonpersistent_summary %>%
group_by(Subtype) %>%
summarize(n_isolates = sum(n_non_persistent)) %>%
ungroup() %>%
mutate(total_83 = sum(n_isolates))
## -------------------------------------------------------------------------------------------------------------------------------------
# Load gene presence/absence data
panaroo_data <- read_csv("~/Desktop/gene_presence_absence.csv")
colnames(panaroo_data)
ncol(panaroo_data)
# See first few column names
colnames(panaroo_data)[1:10]
## this counts the first three column names which are not my isolate IDs. Hence its 356 - 3 = 353. Everything there.
## genes analysed
n_genes <- nrow(panaroo_data)
n_genes
#### Make heatmap of gene presence and absence
# Extract binary matrix: rows = genes, columns = genomes
# Adjust the range to remove metadata columns
binary_matrix <- panaroo_data %>%
select(matches("^[0-9]+Lm$|^MRL-")) %>%
mutate_all(~ ifelse(is.na(.), 0, 1)) %>%
as.matrix()
binary_matrix
# Add rownames as gene IDs
rownames(binary_matrix) <- panaroo_data$Gene
# Open PNG device
png("~/Desktop/pangenome_heatmap.png", width = 1200, height = 1000, res=150)
# plot the heatmap
pangenome_heatmap <- pheatmap(
binary_matrix,
cluster_rows = TRUE,
cluster_cols = TRUE,
color = c( "#FDE725FF","#440154FF" ),
breaks = c(-0.1, 0.5, 1.1),
show_rownames = FALSE,
show_colnames = FALSE,
border_color = NA,
fontsize = 6
)
# Close the device
dev.off()
## -------------------------------------------------------------------------------------------------------------------------------------
# make sure the presence column exists
gene_presence_long <- panaroo_data %>%
select(Gene, matches("^[0-9]+Lm$")) %>%
pivot_longer(-Gene, names_to = "Genome", values_to = "Presence") %>%
mutate(Presence = ifelse(is.na(Presence), 0, 1))
# categorize genes by frequency
gene_freq <- gene_presence_long %>%
group_by(Gene) %>%
summarise(PresenceCount = sum(Presence),
GenomeCount = n_distinct(Genome)) %>%
mutate(Frequency = PresenceCount / GenomeCount,
Category = case_when(
Frequency == 1 ~ "Core",
Frequency >= 0.95 ~ "Soft Core",
Frequency > 0.15 ~ "Shell",
TRUE ~ "Cloud"
))
# Ensure category is a factor in desired order
gene_freq$Category <- factor(gene_freq$Category, levels = c("Core", "Soft Core", "Shell", "Cloud"))
# get the numbers
# Total genes analyzed
total_genes <- nrow(gene_freq) # counts all unique genes
print(total_genes)
# calculate
gene_summary <- gene_freq %>%
group_by(Category) %>%
summarise(
Count = n()
) %>%
mutate(
Percentage = round(Count / total_genes * 100, 1)
)
gene_summary
# Barplot with counts
barplot_gene_categories <- ggplot(gene_freq, aes(x = Category, fill = Category)) +
geom_bar() +
geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.5, size = 4) +
scale_fill_viridis_d(option = "D") +
scale_y_continuous(limits = c(0, 3000)) +
theme_minimal(base_size = 14) +
theme(plot.title = element_text(hjust = 0.5)) +
labs(
title = "Pangenome Gene Categories",
x = "Category",
y = "Number of Genes"
)
barplot_gene_categories
# make sure background is white
barplot_gene_categories <- barplot_gene_categories +
theme(
plot.background = element_rect(fill = "white", color = NA),
panel.background = element_rect(fill = "white", color = NA)
)
# save Barplot
ggsave(
filename = "~/Desktop/barplot_gene_categories.svg",
plot = barplot_gene_categories,
width = 8,
height = 6,
dpi = 300
)
## -------------------------------------------------------------------------------------------------------------------------------------
# Select gene presence/absence (genome columns only)
presence_absence <- panaroo_data %>%
select(matches("^[0-9]+Lm$")) %>%
mutate_all(~ ifelse(is.na(.), 0, 1))
simulate_pangenome <- function(matrix, n_iter = 100) {
n_genomes <- ncol(matrix)
core_accum <- matrix(0, nrow = n_iter, ncol = n_genomes)
pan_accum <- matrix(0, nrow = n_iter, ncol = n_genomes)
for (i in 1:n_iter) {
sampled_order <- sample(ncol(matrix))
for (j in 1:n_genomes) {
subset <- matrix[, sampled_order[1:j], drop = FALSE] # fix here
gene_sums <- rowSums(subset)
pan_accum[i, j] <- sum(gene_sums > 0)
core_accum[i, j] <- sum(gene_sums == j)
}
}
tibble(
Genomes = 1:n_genomes,
Core = colMeans(core_accum),
Pangenome = colMeans(pan_accum)
)
}
set.seed(123)
curve_data <- simulate_pangenome(as.matrix(presence_absence), n_iter = 100)
# Pivot for ggplot
curve_long <- curve_data %>%
pivot_longer(cols = c(Core, Pangenome), names_to = "Type", values_to = "Genes")
pangenome_curves <- ggplot(curve_long, aes(x = Genomes, y = Genes, color = Type)) +
geom_line(size = 1.2) +
scale_color_viridis_d(option = "D") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5) # Center the title
) +
labs(
title = "Pangenome and Core Genome Curves",
x = "Number of Genomes",
y = "Number of Genes",
color = "Genome Part"
)
pangenome_curves
# save pangenome_curves
ggsave(
filename = "~/Desktop/pangenome_curves.svg",
plot = pangenome_curves,
width = 8,
height = 6,
dpi = 300
)
## -------------------------------------------------------------------------------------------------------------------------------------
# simulate rarefraction
simulate_rarefaction <- function(matrix, n_iter = 100) {
n_genomes <- ncol(matrix)
gene_accum <- matrix(0, nrow = n_iter, ncol = n_genomes)
for (i in 1:n_iter) {
sampled_order <- sample(ncol(matrix))
for (j in 1:n_genomes) {
subset <- matrix[, sampled_order[1:j], drop = FALSE]
gene_accum[i, j] <- sum(rowSums(subset) > 0)
}
}
tibble(
Genomes = 1:n_genomes,
GeneFamilies = colMeans(gene_accum)
)
}
set.seed(42)
rarefaction_data <- simulate_rarefaction(binary_matrix, n_iter = 100)
# Plot
rarefaction_curve <- ggplot(rarefaction_data, aes(x = Genomes, y = GeneFamilies)) +
geom_line(color = "steelblue", size = 1.2) +
theme_minimal() +
labs(title = "Rarefaction Curve: Gene Family Diversity",
x = "Number of Genomes",
y = "Number of Gene Families") +
theme(plot.title = element_text(hjust = 0.5))
rarefaction_curve
## save Plot
ggsave(
filename = "~/Desktop/rarefaction_curves_genefamilies.svg",
plot = rarefaction_curve,
width = 8,
height = 6,
dpi = 300)
## -------------------------------------------------------------------------------------------------------------------------------------
colnames(metadata)
## read metadata for historical isolates (files was created manually)
metadata_historical_isolates <- read_delim(
"~/Desktop/metadata_historical_pca.csv",
delim = ";",col_types = cols(.default = "c")) # treat all columns as character)
## filter metadata file for needed information
metadata_filtered_pca <- metadata %>%
filter(Sequence_ID_Lmono != "NA") %>%
select(
Sequence_ID_Lmono,
Meta_ID ,
Lineage,
Clonal_Complex,
ST,
CT,
Facility_Country,
Facility_Type
)
## combine both
combined_metadata <- bind_rows(metadata_filtered_pca, metadata_historical_isolates)
# PCA on binary matrix (samples = columns, so transpose)
pca_res <- prcomp(t(binary_matrix), scale. = FALSE)
# Create PCA data frame for plotting
pca_df <- as.data.frame(pca_res$x[,1:2])
pca_df$Genome <- colnames(binary_matrix)
# Join PCA results with filtered metadata by genome/sample ID
pca_annotated <- pca_df %>%
left_join(combined_metadata, by = c("Genome" = "Sequence_ID_Lmono"))
# Plot PCA with larger and bolder labels
p <- ggplot(pca_annotated, aes(x = PC1, y = PC2, color = Facility_Type)) +
geom_point(size = 3, alpha = 0.8) +
stat_ellipse(aes(group = Facility_Type), linetype = "dashed", size = 1) +
theme_minimal() +
labs(
title = "PCA of Gene Content by Facility Type",
x = "PC1",
y = "PC2",
color = "Facility Type"
) +
theme(
plot.title = element_text(hjust = 0.5, size = 18, face = "bold"), # Title bigger and bold
axis.title.x = element_text(size = 14, face = "bold"), # X axis title
axis.title.y = element_text(size = 14, face = "bold"), # Y axis title
axis.text = element_text(size = 12), # Axis labels
legend.title = element_text(size = 14, face = "bold"), # Legend title
legend.text = element_text(size = 12) # Legend labels
) +
scale_color_manual(
values = c(
"Meat" = "#404788", # dark blue
"Dairy" = "#55C667", # green
"Fish" = "#FDE725" # yellow
)
)
# If you want to add circles around points of each Facility_Type (optional):
pca_fin <- p + stat_ellipse(aes(group = Facility_Type), linetype = "dashed", size = 1)
pca_fin
## save Plot
ggsave(
filename = "~/Desktop/PCA_GeneContent_Facility_Type.svg",
plot = pca_fin,
width = 8,
height = 6,
dpi = 300)
## -------------------------------------------------------------------------------------------------------------------------------------
# Ensure ST is numeric and ordered
pca_annotated$ST <- factor(as.numeric(as.character(pca_annotated$ST)),
levels = sort(unique(as.numeric(as.character(pca_annotated$ST)))))
# Calculate centroids for Facility_Type
ellipse_labels <- pca_annotated %>%
group_by(Facility_Type) %>%
summarize(
PC1 = mean(PC1, na.rm = TRUE),
PC2 = mean(PC2, na.rm = TRUE),
.groups = "drop"
)
# Final PCA plot
pca_cgMLST <- ggplot(pca_annotated, aes(x = PC1, y = PC2)) +
geom_point(aes(color = ST), size = 3, alpha = 0.8) +
stat_ellipse(aes(group = Facility_Type), linetype = "dashed", size = 0.5, color = "black") +
geom_text_repel(
data = ellipse_labels,
aes(x = PC1, y = PC2, label = Facility_Type),
fontface = "bold", color = "black", size = 3,
box.padding = 0.4, point.padding = 0.5,
segment.color = "grey50"
) +
scale_color_viridis_d(option = "D") +
theme_minimal() +
labs(
title = "PCA of Gene Content",
x = "PC1",
y = "PC2",
color = "Sequence Type (ST)"
) +
theme(plot.title = element_text(hjust = 0.5))
pca_cgMLST
ggsave(
filename = "~/Desktop/PCA_cgMLST_Facility_Type.svg",
plot = pca_cgMLST,
width = 8,
height = 6,
dpi = 300)
## -------------------------------------------------------------------------------------------------------------------------------------
# Load data
scoary_results <- read.csv("~/Desktop/scoary_eggnog_merged.csv")
nrow(scoary_results)
# Prepare data
scoary_results <- scoary_results %>%
mutate(
log2OR = log2(Odds_ratio),
negLog10BH = -log10(Benjamini_H_p),
Association = case_when(
log2OR > 0 & Benjamini_H_p < 0.05 ~ "Persistence",
log2OR < 0 & Benjamini_H_p < 0.05 ~ "Non-persistence",
TRUE ~ "Not significant"
))
nrow(scoary_results)
# Count the number of genes per association category
gene_counts <- scoary_results %>%
summarise(
total_genes = n(),
persistence = sum(log2OR > 0 & Benjamini_H_p < 0.05, na.rm = TRUE),
non_persistence = sum(log2OR < 0 & Benjamini_H_p < 0.05, na.rm = TRUE),
not_significant = sum(!(log2OR > 0 & Benjamini_H_p < 0.05) &
!(log2OR < 0 & Benjamini_H_p < 0.05), na.rm = TRUE)
)
print(gene_counts)
# Filter out points where y > 5
scoary_filtered <- scoary_results %>%
filter(negLog10BH <= 5)
nrow(scoary_filtered)
# Volcano plot without points at y = 5
volcano_plot <- ggplot(scoary_filtered, aes(x = log2OR, y = negLog10BH, color = Association)) +
geom_point(alpha = 0.7, size = 2) +
scale_color_manual(values = c(
"Persistence" = "#3B528BFF",
"Non-persistence" = "#55C667FF",
"Not significant" = "grey70"
)) +
geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "black") +
labs(
title = "Genes Associated with Persistence",
x = "log2(Odds Ratio)",
y = "-log10(BH-adjusted p-value)",
color = "Association"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
legend.position = "right"
) +
coord_cartesian(xlim = c(-5, 5)) +
scale_x_continuous(breaks = seq(-5, 5, by = 1))
volcano_plot
# Save to PDF
ggsave(
filename = "~/Desktop/volcano_persistence_genes_clipped.svg",
plot = volcano_plot,
width = 10,
height = 6,
dpi = 300
)
## -------------------------------------------------------------------------------------------------------------------------------------
# log2 OR and confidence intervals
# calculate confidence intervals for the values
scoary_results <- scoary_results %>%
mutate(
log_OR = log(Odds_ratio),
se_log_or = sqrt(1/Number_pos_present_in + 1/Number_neg_present_in + 1/Number_pos_not_present_in + 1/Number_neg_not_present_in),
lower_ci = exp(log_OR - 1.96 * se_log_or),
upper_ci = exp(log_OR + 1.96 * se_log_or),
log2_OR = log2(Odds_ratio),
log2_lower = log2(lower_ci),
log2_upper = log2(upper_ci)
)
nrow(scoary_results)
# kick our OR of infinite
scoary_results <- scoary_results %>%
filter(
is.finite(log_OR) &
is.finite(log2_OR) &
is.finite(log2_lower) &
is.finite(log2_upper)
)
# Filter significant and strong associations
filtered_genes <- scoary_results %>%
mutate(
odds_ratio = as.numeric(Odds_ratio),
# sensitivity = as.numeric(Sensitivity),
# specificity = as.numeric(Specificity),
adjusted_p = as.numeric(Benjamini_H_p)
) %>%
filter(
adjusted_p < 0.05, # Statistically significant
odds_ratio >= 2 | odds_ratio <= 0.5, # Strong association
# sensitivity >= 50, # Good detection of phenotype-positive
# specificity >= 50 # Good discrimination of phenotype-negative
)
nrow(filtered_genes)
colnames(filtered_genes)
# convert COG category as character and replace missing values with NA
filtered_genes$COG_category_y <- as.character(filtered_genes$COG_category_y)
filtered_genes$COG_category_y[is.na(filtered_genes$COG_category_y)] <- "NA"
# Define COG categories and descriptions
unique(filtered_genes$COG_category_y)
cog_list <- c("NA", "L", "V", "S", "M", "-", "IU", "C", "D", "K", "NU", "U", "MT","T", "F", "KL", "O", "E", "P", "G", "N", "KT", "J", "I", "Q", "GM", "EQ")
COG_descriptions <- c(
"NA" = "Not assigned",
"L" = "Replication, recombination and repair",
"V" = "Defense mechanisms",
"S" = "Function unknown",
"M" = "Cell wall/membrane/envelope biogenesis",
"-" = "Unknown",
"IU" = "Lipid transport and metabolism, Intracellular trafficking, secretion, and vesicular transport",
"C" = "Energy production and conversion",
"D" = "Cell cycle control, cell division, chromosome partitioning",
"K" = "Transcription",
"NU" = "Cell motility, Intracellular trafficking, secretion, and vesicular transport",
"U" = "Intracellular trafficking, secretion, and vesicular transport",
"MT" = "Cell wall/membrane/envelope biogenesis, Signal transduction mechanisms",
"T" = "Signal transduction mechanisms",
"F" = "Nucleotide transport and metabolism",
"KL" = "Transcription, Replication, recombination and repair ",
"O" = "Post-translational modification, protein turnover, chaperones",
"E" = "Amino acid transport and metabolism",
"P" = "Inorganic ion transport and metabolism",
"G" = "Carbohydrate transport and metabolism",
"N" = "Cell motility",
"KT" = "Transcription, Signal transduction mechanism",
"J" = "Translation, ribosomal structure and biogenesis",
"I" = "Lipid transport and metabolism",
"Q" = "Secondary metabolites biosynthesis, transport and catabolism",
"GM" = "Carbohydrate transport and metabolis, Cell wall/membrane/envelope biogenes",
"EQ" = "Amino acid transport and metabolism, Secondary metabolites biosynthesis, transport and catabolism"
)
# Set y-axis limits for plotting, so that CIs are nicely displayed)
min_log2 <- floor(min(filtered_genes$log2_lower, na.rm = TRUE))
max_log2 <- ceiling(max(filtered_genes$log2_upper, na.rm = TRUE))
## type the genes to enriched or depleted
filtered_genes <- filtered_genes %>%
mutate(
Type = case_when(
log2_OR > 0 ~ "Enriched",
log2_OR < 0 ~ "Depleted",
TRUE ~ "Neutral"
)
)
## Replace NA PFAMs with Unknown
filtered_genes$PFAMs[is.na(filtered_genes$PFAMs)] <- "Unknown"
pdf("~/Desktop/PFAM_COG_all_categories.pdf", width = 12, height = 8)
# Replace missing PFAMs
filtered_genes$PFAMs[is.na(filtered_genes$PFAMs)] <- "Unknown"
for(cog in cog_list) {
df <- filtered_genes %>% filter(COG_category_y == cog)
if(nrow(df) == 0) next
cog_desc <- COG_descriptions[cog]
# Order PFAMs by log2_OR descending
df <- df %>%
mutate(PFAMs = fct_reorder(PFAMs, log2_OR, .desc = TRUE))
# Plot
p <- ggplot(df, aes(x = PFAMs, y = log2_OR, color = Type)) +
geom_point(size = 3.5, position = position_dodge(width = 0.7)) +
geom_errorbar(aes(ymin = log2_lower, ymax = log2_upper),
width = 0.2, position = position_dodge(width = 0.7)) +
coord_flip() +
scale_color_manual(values = c(
"Enriched" = "#404788FF",
"Depleted" = "#55C667FF",
"Neutral" = "grey50"
)) +
scale_y_continuous(limits = c(-6, 6)) +
labs(
title = paste0("COG Category ", cog, " (", cog_desc, ")"),
x = "PFAM",
y = "log2(Odds Ratio)",
color = "Type"
) +
theme_minimal(base_size = 14) +
theme(
axis.text.y = element_text(size = 12, face = "bold", margin = margin(r = 14)),
axis.title.y = element_text(size = 14, face = "bold"),
axis.title.x = element_text(size = 14, face = "bold"),
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
legend.position = "bottom",
legend.text = element_text(size = 12, face = "bold")
)
print(p)
}
dev.off()
## -------------------------------------------------------------------------------------------------------------------------------------
# Summarize — recode missing/dash COG to "Unassigned"
cog_summary <- filtered_genes %>%
mutate(COG_category_y = case_when(
is.na(COG_category_y) | COG_category_y == "" | COG_category_y == "-" ~ "Unassigned",
TRUE ~ COG_category_y
)) %>%
group_by(COG_category_y, Type) %>%
summarise(
Mean_log2OR = mean(log2OR, na.rm = TRUE),
Gene_Count = n(),
PFAMs = paste(unique(na.omit(PFAMs)), collapse = "; "),
PFAM_Count = length(unique(na.omit(PFAMs))),
.groups = "drop"
)
# Order COGs by total gene count descending (same logic as first code block), Unassigned always last
cog_order <- cog_summary %>%
group_by(COG_category_y) %>%
summarise(Total_Genes = sum(Gene_Count)) %>%
arrange(desc(Total_Genes)) %>%
pull(COG_category_y)
# Move Unassigned to end
cog_order <- c("Unassigned", cog_order[cog_order != "Unassigned"])
# Apply factor levels — left to right = highest to lowest, Unassigned last
cog_summary <- cog_summary %>%
mutate(COG_category_y = factor(COG_category_y, levels = cog_order))
# Save
write_csv(cog_summary %>%
select(COG_category_y, Type, Mean_log2OR, Gene_Count, PFAM_Count, PFAMs) %>%
left_join(
cog_summary %>%
group_by(COG_category_y) %>%
summarise(Overall_Mean_log2OR = mean(Mean_log2OR, na.rm = TRUE)),
by = "COG_category_y"
),
"~/Desktop/cog_summary_filtered.csv")
# Plot
p <- ggplot(cog_summary, aes(x = COG_category_y, y = Mean_log2OR, size = Gene_Count, color = Type)) +
geom_point(alpha = 0.8) +
geom_hline(yintercept = 0, color = "black", linetype = "solid", linewidth = 1) +
scale_color_manual(values = c("Enriched" = "#404788FF", "Depleted" = "#55C667FF")) +
scale_size_continuous(
breaks = c(1, 10, 50, 100, 1500),
labels = c("1", "10", "50", "100", "1500"),
range = c(1, 10)
) +
scale_y_continuous(limits = c(
floor(min(cog_summary$Mean_log2OR, na.rm = TRUE)),
ceiling(max(cog_summary$Mean_log2OR, na.rm = TRUE))
)) +
labs(
title = "log2(OR) and Gene Counts per COG",
x = "COG Category",
y = "Mean log2(OR)",
size = "Gene Count",
color = "Type"
) +
theme_bw(base_size = 16) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 20),
axis.text.x = element_text(angle = 45, hjust = 1, face = "bold", size = 14),
axis.text.y = element_text(face = "bold", size = 18),
axis.title = element_text(face = "bold", size = 18),
legend.title = element_text(face = "bold", size = 18),
legend.text = element_text(size = 18)
)
p
ggsave("~/Desktop/COG_summary_plot.svg", plot = p, width = 12, height = 8, device = "svg")
## -------------------------------------------------------------------------------------------------------------------------------------
# Filter enriched PFAMs for selected COG categories
selected_cogs <- c("KL", "MT", "N", "NU", "O", "Q")
enriched_selected <- filtered_genes%>%
filter(COG_category_y %in% selected_cogs & is.finite(log2_OR))
# Replace NA in COG_category_y with "NA" string if present
enriched_selected$COG_category_y <- as.character(enriched_selected$COG_category_y)
enriched_selected$COG_category_y[is.na(enriched_selected$COG_category_y)] <- "NA"
# Assign descriptive labels
COG_descriptions <- c(
"KL" = "KL: Transcription & Replication/Repair",
"MT" = "MT: Cell wall & Transport",
"N" = "N: Cell motility",
"NU" = "NU: Cell motility & Unknown",
"O" = "O: Post-translational modification, protein turnover, chaperones",
"Q" = "Q: Secondary metabolites biosynthesis, transport and catabolism"
)
enriched_selected$COG_label <- COG_descriptions[enriched_selected$COG_category_y]
# Optional: order PFAMs within each COG by log2_OR
enriched_selected <- enriched_selected %>%
group_by(COG_label) %>%
arrange(log2_OR) %>%
mutate(PFAMs_ordered = factor(PFAMs, levels = PFAMs)) %>%
ungroup()
# Plot with stacked facets and borders
p_enriched_only <- ggplot(enriched_selected, aes(x = PFAMs_ordered, y = log2_OR)) +
geom_point(color = "#440154FF", size = 3) +
geom_errorbar(aes(ymin = log2_lower, ymax = log2_upper), width = 0.2, color = "#440154FF") +
coord_flip() +
facet_wrap(~ COG_label, ncol = 1, scales = "free_y") + # stack vertically
scale_y_continuous(limits = c(0, 6)) +
labs(
title = "PFAMs Enriched in Persistent Isolates",
x = "PFAM",
y = "log2(Odds Ratio)"
) +
theme_minimal(base_size = 14) +
theme(
axis.text.y = element_text(size = 16, margin = margin(r = 14)),
axis.title.y = element_text(face = "bold"),
axis.title.x = element_text(face = "bold"),
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
strip.text = element_text(face = "bold", size = 9),
strip.background = element_rect(color = "black", fill = "lightgray", linewidth = 0.8),
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
panel.spacing = unit(0.5, "lines")
)
p_enriched_only
# Save as SVG
ggsave(
filename = "~/Desktop/PFAM_enriched_selected_COGs_stacked_borders.svg",
plot = p_enriched_only,
width = 10,
height = 8,
dpi = 300
)
# Count genes per COG category
genes_per_cog <- enriched_selected %>%
group_by(COG_label) %>%
summarise(n_genes = n()) %>%
arrange(desc(n_genes))
print(genes_per_cog)
# If you want the total across all selected COGs
total_genes <- nrow(enriched_selected)
cat("Total genes in enriched COG categories:", total_genes, "\n")
# Check the column names to see what you have
colnames(enriched_selected)
# Look at the first few rows
head(enriched_selected, 10)
# Check if there's a gene/locus identifier column
# (common names: gene_id, locus_tag, protein_id, gene, ID, etc.)
str(enriched_selected)
# If you have a gene ID column, check for duplicates
# Replace 'gene_id' with whatever your actual column name is
enriched_selected %>%
group_by(Gene) %>% # use your actual column name
filter(n() > 1) %>%
arrange(Gene)
# Compare: number of rows vs unique PFAMs
nrow(enriched_selected) # total rows
n_distinct(enriched_selected$PFAMs) # unique PFAMs
## -------------------------------------------------------------------------------------------------------------------------------------
# Filter depleted PFAMs for selected COG categories
selected_cogs_depleted <- c("F", "GM","IU")
depleted_selected <- filtered_genes %>%
filter(COG_category_y %in% selected_cogs_depleted & is.finite(log2_OR))
# Replace NA in COG_category_y with "NA" string if present
depleted_selected$COG_category_y <- as.character(depleted_selected$COG_category_y)
depleted_selected$COG_category_y[is.na(depleted_selected$COG_category_y)] <- "NA"
# Assign descriptive labels
COG_descriptions_depleted <- c(
"F" = "F: Nucleotide transport and metabolism",
"GM" = "GM: Carbohydrate transport and metabolism & Cell wall/membrane/envelope biogenesis",
"IU" = "IU: Lipid transport and metabolism, Intracellular trafficking, secretion, and vesicular transport"
)
depleted_selected$COG_label <- COG_descriptions_depleted[depleted_selected$COG_category_y]
# Optional: order PFAMs within each COG by log2_OR
depleted_selected <- depleted_selected %>%
group_by(COG_label) %>%
arrange(log2_OR) %>%
mutate(PFAMs_ordered = factor(PFAMs, levels = PFAMs)) %>%
ungroup()
# Plot with stacked facets and borders
p_depleted <- ggplot(depleted_selected, aes(x = PFAMs_ordered, y = log2_OR)) +
geom_point(color = "#55C667FF", size = 3) + # different color for depleted
geom_errorbar(aes(ymin = log2_lower, ymax = log2_upper), width = 0.2, color = "#55C667FF") +
coord_flip() +
facet_wrap(~ COG_label, ncol = 1, scales = "free_y") + # stack vertically
scale_y_continuous(limits = c(-6, 0)) +
labs(
title = "PFAMs Depleted in Persistent Isolates",
x = "PFAM",
y = "log2(Odds Ratio)"
) +
theme_minimal(base_size = 14) +
theme(
axis.text.y = element_text(size = 16, margin = margin(r = 14)),
axis.title.y = element_text(face = "bold"),
axis.title.x = element_text(face = "bold"),
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
strip.text = element_text(face = "bold", size = 12),
strip.background = element_rect(color = "black", fill = "lightgray", size = 0.8),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
panel.spacing = unit(0.5, "lines")
)
p_depleted
# Save as SVG
ggsave(
filename = "~/Desktop/PFAM_depleted_selected_COGs_stacked_borders.svg",
plot = p_depleted,
width = 10,
height = 8,
dpi = 300
)
## combine it nicely
combined_vertical <- p_enriched_only / p_depleted +
plot_layout(heights = c(1, 1)) # adjust if needed
ggsave(
"~/Desktop/PFAM_enriched_depleted_combined.svg",
plot = combined_vertical,
width = 380, # mm (≈ single-column width * 2 if needed)
height = 480, # taller for stacking
units = "mm"
)
## -------------------------------------------------------------------------------------------------------------------------------------
# ENRICHED GENES
# Filter enriched PFAMs for selected COG categories (now including E and S)
selected_cogs_enriched <- c("KL", "MT", "N", "NU", "O", "Q", "E", "S")
enriched_selected <- filtered_genes %>%
filter(COG_category_y %in% selected_cogs_enriched & is.finite(log2_OR) & Type == "Enriched")
# Count genes per enriched COG category
genes_enriched_cogs <- enriched_selected %>%
group_by(COG_category_y) %>%
summarise(n_genes = n(), .groups = "drop") %>%
arrange(desc(n_genes))
print(genes_enriched_cogs)
# Count total PFAM domains
total_pfams_enriched <- enriched_selected %>%
mutate(n_pfams = str_count(PFAMs, ",") + 1) %>%
pull(n_pfams) %>%
sum()
# DEPLETED GENES
# Filter depleted PFAMs for selected COG categories (now including E and S)
selected_cogs_depleted <- c("F", "GM", "IU", "E", "S")
depleted_selected <- filtered_genes %>%
filter(COG_category_y %in% selected_cogs_depleted & is.finite(log2_OR) & Type == "Depleted")
# Count genes per depleted COG category
genes_depleted_cogs <- depleted_selected %>%
group_by(COG_category_y) %>%
summarise(n_genes = n(), .groups = "drop") %>%
arrange(desc(n_genes))
print(genes_depleted_cogs)
# Count total PFAM domains
total_pfams_depleted <- depleted_selected %>%
mutate(n_pfams = str_count(PFAMs, ",") + 1) %>%
pull(n_pfams) %>%
sum()
# COMBINED SUMMARY WITH GENE NAMES
summary_table_with_genes <- bind_rows(
enriched_selected %>%
select(Gene, COG_category_y, PFAMs) %>%
mutate(Type = "Enriched"),
depleted_selected %>%
select(Gene, COG_category_y, PFAMs) %>%
mutate(Type = "Depleted")
) %>%
arrange(Type, COG_category_y, Gene)
print(summary_table_with_genes)
# Save to CSV for easy viewing
write.csv(summary_table_with_genes,
"~/Desktop/enriched_depleted_genes_summary_fin.csv",
row.names = FALSE)
## -------------------------------------------------------------------------------------------------------------------------------------
# Select COGs S, K, L for enriched and depleted
selected_cogs <- c("S", "K", "L")
combined_selected <- filtered_genes %>%
filter(COG_category_y %in% selected_cogs & is.finite(log2_OR)) %>%
mutate(
COG_category_y = as.character(COG_category_y),
COG_category_y = ifelse(is.na(COG_category_y), "NA", COG_category_y)
)
# Assign descriptive COG labels
COG_descriptions <- c(
"S" = "S: Function unknown",
"K" = "K: Transcription",
"L" = "L: Replication, recombination and repair"
)
combined_selected <- combined_selected %>%
mutate(COG_label = COG_descriptions[COG_category_y])
# Assign colors
combined_selected <- combined_selected %>%
mutate(color_fill = ifelse(Type == "Enriched", "#440154FF", "#55C667FF")) %>%
group_by(COG_label) %>%
arrange(log2_OR) %>%
mutate(PFAMs_ordered = factor(make.unique(as.character(PFAMs)), levels = make.unique(as.character(PFAMs)))) %>%
ungroup()
# Create separate plot for each COG with custom dimensions
for(cog in unique(combined_selected$COG_label)) {
df <- combined_selected %>% filter(COG_label == cog)
p <- ggplot(df, aes(x = PFAMs_ordered, y = log2_OR, color = color_fill)) +
geom_point(size = 3) +
geom_errorbar(aes(ymin = log2_lower, ymax = log2_upper), width = 0.2) +
scale_color_identity(name = "Type", labels = c("Enriched", "Depleted")) +
coord_flip() +
scale_y_continuous(limits = c(-8, 8),breaks = seq(-8, 8, by = 2)) +
labs(
title = paste("", cog),
x = "PFAM",
y = "log2(Odds Ratio)"
) +
theme_minimal(base_size = 14) +
theme(
axis.text.y = element_text(size = 12, face = "bold", margin = margin(r = 14)),
axis.title.y = element_text(size = 14, face = "bold"),
axis.title.x = element_text(size = 14, face = "bold"),
plot.title = element_text(hjust = 0.5, face = "bold", size = 16),
panel.border = element_rect(color = "black", fill = NA, size = 0.8),
legend.position = "bottom",
legend.title = element_blank()
)
# Define custom dimensions
width_plot <- ifelse(cog == COG_descriptions["S"], 10, 12)
height_plot <- ifelse(cog == COG_descriptions["S"], 20, 12)
# Save each COG plot as SVG
filename <- paste0("~/Desktop/PFAM_enriched_depleted_", gsub(": ", "_", cog), ".svg")
ggsave(filename, plot = p, width = width_plot, height = height_plot, dpi = 300)
}
# combine the tree figures for publication
selected_cogs <- c("S", "K", "L")
COG_descriptions <- c(
"S" = "S: Function unknown",
"K" = "K: Transcription",
"L" = "L: Replication, recombination and repair"
)
combined_selected <- filtered_genes %>%
filter(COG_category_y %in% selected_cogs & is.finite(log2_OR)) %>%
mutate(
COG_category_y = as.character(COG_category_y),
COG_category_y = ifelse(is.na(COG_category_y), "NA", COG_category_y),
COG_label = COG_descriptions[COG_category_y],
color_fill = ifelse(Type == "Enriched", "#440154FF", "#55C667FF")
) %>%
group_by(COG_label) %>%
arrange(log2_OR) %>%
mutate(
PFAMs_ordered = factor(make.unique(as.character(PFAMs)),
levels = make.unique(as.character(PFAMs)))
) %>%
ungroup()
# Plot function
make_cog_plot <- function(cog_label) {
df <- combined_selected %>% filter(COG_label == cog_label)
y_margin <- ifelse(cog_label == "S: Function unknown", 40, 10)
line_height <- ifelse(cog_label == "S: Function unknown", 2.5, 1.8)
ggplot(df, aes(x = PFAMs_ordered, y = log2_OR, color = color_fill)) +
geom_point(size = 5) + # increased
geom_errorbar(aes(ymin = log2_lower, ymax = log2_upper),
width = 0.4, linewidth = 0.8) + # increased
scale_color_identity(name = "Type", labels = c("Enriched", "Depleted")) +
coord_flip() +
scale_y_continuous(limits = c(-8, 8), breaks = seq(-8, 8, by = 2)) +
labs(
title = cog_label,
x = "PFAM",
y = "log2(Odds Ratio)"
) +
theme_minimal(base_size = 18) + # increased
theme(
axis.text.y = element_text(size = 16, face = "bold", # increased
margin = margin(r = y_margin),
lineheight = line_height),
axis.text.x = element_text(size = 16), # increased
axis.title.y = element_text(size = 18, face = "bold",
margin = margin(r = 10)),
axis.title.x = element_text(size = 18, face = "bold",
margin = margin(t = 10)),
plot.title = element_text(hjust = 0.5, face = "bold", size = 20),
panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
plot.margin = margin(10, 15, 10, 15),
legend.position = "bottom",
legend.title = element_blank(),
legend.text = element_text(size = 16),
legend.key.size = unit(0.8, "cm")
)
}
# Rebuild plots
p_S <- make_cog_plot("S: Function unknown")
p_K <- make_cog_plot("K: Transcription")
p_L <- make_cog_plot("L: Replication, recombination and repair")
# Combine
combined_plot <- ((p_K / p_L) | p_S) +
plot_layout(widths = c(1, 1))
ggsave(
"~/Desktop/PFAM_combined_COG_S_K_L.svg",
plot = combined_plot,
width = 570,
height = 720,
units = "mm",
device = "svg"
)
## -------------------------------------------------------------------------------------------------------------------------------------
# Read data
virulence_summary <- read.csv(
"~/Desktop/virulence_summary_output.csv",
sep = ";", header = TRUE, check.names = FALSE
)
colnames (virulence_summary)
# Gene groups
# source: https://www.mgc.ac.cn/cgi-bin/VFs/genus.cgi?Genus=Listeria
# Define gene groups for Listeria virulence factors
# Existing gene groups
gene_groups <- list(
adherence = c("ami", "lap", "lapB", "fbpA", "lpeA", "lplA1"),
invasion = c("inlA", "inlB", "inlC", "inlF", "inlJ", "aut", "iap"),
motility = c("actA", "mpl"),
hemolysins = c("hly", "plcA", "plcB"),
toxins = c("llsA", "llsB", "llsD", "llsG", "llsH", "llsP", "llsX", "llsY"),
immune_modulation = c("vip", "inlK", "lntA"),
cell_envelope = c("oatA", "pdgA", "gtcA", "lspA"),
stress = c("bsh", "clpC", "clpE", "clpP", "prsA2"),
regulation = c("prfA"),
metabolic = c("hpt")
)
# Columns that exist in the dataset
existing_genes <- colnames(virulence_summary)
existing_genes <- existing_genes[!existing_genes %in% c("Isolate", "Phenotype", "ST")]
# Filter gene_groups to keep only genes present in the dataset
gene_groups <- lapply(gene_groups, function(genes) {
genes[genes %in% existing_genes]
})
# Check result
gene_groups
# Make sure Phenotype & ST are factors
virulence_summary$Phenotype <- factor(
virulence_summary$Phenotype,
levels = c("Non-Persistent", "Persistent")
)
virulence_summary$ST <- factor(virulence_summary$ST)
# built matrix
vir_mat <- virulence_summary %>%
select(-Isolate, -Phenotype, -ST) %>%
as.matrix()
rownames(vir_mat) <- virulence_summary$Isolate
# order isolates by ST and phenotype
virulence_summary <- virulence_summary %>%
arrange(ST, desc(Phenotype))
vir_mat <- vir_mat[virulence_summary$Isolate, ]
# map genes to groups
gene_group_mapping <- unlist(lapply(names(gene_groups), function(group) {
setNames(rep(group, length(gene_groups[[group]])), gene_groups[[group]])
}))
gene_group_mapping <- gene_group_mapping[colnames(vir_mat)]
genes_ordered <- names(sort(gene_group_mapping))
vir_mat <- vir_mat[, genes_ordered]
virulence_summary$Phenotype <- factor(
virulence_summary$Phenotype,
levels = c("Non-Persistent", "Persistent")
)
# Prepare matrix
vir_mat <- virulence_summary %>%
select(-Isolate, -Phenotype) %>%
as.matrix()
rownames(vir_mat) <- virulence_summary$Isolate
# Arrange isolates by phenotype
virulence_summary <- virulence_summary %>%
arrange(desc(Phenotype))
vir_mat <- vir_mat[virulence_summary$Isolate, ]
# Map genes to groups
gene_group_mapping <- unlist(lapply(names(gene_groups), function(group) {
setNames(rep(group, length(gene_groups[[group]])), gene_groups[[group]])
}))
gene_group_mapping <- gene_group_mapping[colnames(vir_mat)]
# Order genes by group
genes_ordered <- names(sort(gene_group_mapping))
vir_mat <- vir_mat[, genes_ordered]
# Row annotation for Phenotype
row_ha <- rowAnnotation(
Phenotype = virulence_summary$Phenotype,
col = list(Phenotype = c("Persistent"="#440154FF", "Non-Persistent"="#55C667FF")),
annotation_legend_param = list(
title_gp = gpar(fontsize = 12, fontface = "bold"), # legend title
labels_gp = gpar(fontsize = 12) # legend labels
)
)
# Heatmap colors
heat_col <- c("0" = "#D3D3D3", "1" = "#87CEFA")
# Column split factor
col_split <- factor(gene_group_mapping[genes_ordered])
# Replace rownames with ST
rownames(vir_mat) <- virulence_summary$ST
ht <- Heatmap(
vir_mat,
name = "Presence",
col = heat_col,
cluster_rows = FALSE,
cluster_columns = FALSE,
show_row_names = TRUE, # display ST per row
row_names_side = "left", # left side
row_names_gp = gpar(fontsize=10), # horizontal by default
show_column_names = TRUE,
column_names_gp = gpar(fontsize=12),
column_title_rot = 45,
column_split = col_split,
heatmap_legend_param = list(
title_gp = gpar(fontsize=12, fontface="bold"),
labels_gp = gpar(fontsize=12)
)
)
# Draw heatmap centered with padding
svg("~/Desktop/virulence_heatmap_final_centered.svg", width=14, height=8)
draw(row_ha +
ht,
heatmap_legend_side = "right",
padding = unit(c(0.5, 1, 0.5, 1), "mm") # top, right, bottom, left
)
dev.off()
## -------------------------------------------------------------------------------------------------------------------------------------
# Key genes to check for absence (VIP included)
genes_absent <- c("actA","ami","inlF","inlJ","aut","vip")
# Toxins
toxins <- c("llsA","llsB","llsD","llsG","llsH","llsP","llsX","llsY")
# Load virulence summary
virulence_summary <- read.csv(
"~/Desktop/virulence_summary_output.csv",
sep = ";", header = TRUE, check.names = FALSE
)
# Load metadata
metadata <- read_excel(
"~/Desktop/AllPartners_FoodSafeR_WP2.3_Facility_Sampling_Metadata_FS_Amplicon_controls_historical_isolates.xlsx"
)
metadata$Sequence_ID_Lmono <- trimws(metadata$Sequence_ID_Lmono)
# Function to map isolate IDs to metadata (keep ST and CT)
map_to_metadata <- function(isolate_ids, metadata) {
merge(
data.frame(Sequence_ID_Lmono = isolate_ids),
metadata[, c("Sequence_ID_Lmono", "ST", "CT")],
by = "Sequence_ID_Lmono",
all.x = TRUE
)
}
# Filter persistent isolates
persistent <- virulence_summary %>%
filter(Phenotype == "Persistent")
# Find isolates missing genes
missing_per_gene <- lapply(genes_absent, function(gene) {
isolates <- persistent$Isolate[persistent[[gene]] == 0]
isolates <- isolates[!startsWith(isolates, "MS")] # exclude MS
map_to_metadata(isolates, metadata)
})
names(missing_per_gene) <- genes_absent
# Find isolates harboring toxin
# At least one toxin present
any_toxin_isolates <- persistent$Isolate[rowSums(persistent[, toxins]) >= 1]
any_toxin_isolates <- any_toxin_isolates[!startsWith(any_toxin_isolates, "MS")]
any_toxin_meta <- map_to_metadata(any_toxin_isolates, metadata)
# Individual toxins
harbor_toxins_list <- lapply(toxins, function(toxin) {
isolates <- persistent$Isolate[persistent[[toxin]] == 1] # present
isolates <- isolates[!startsWith(isolates, "MS")]
map_to_metadata(isolates, metadata)
})
names(harbor_toxins_list) <- toxins
# Combine into one list
final_result_list <- list(
missing_per_gene = missing_per_gene, # includes VIP
harbor_any_toxin = harbor_toxins_list
)
# Convert to long table
combined_long <- do.call(rbind, lapply(names(final_result_list), function(category) {
sublist <- final_result_list[[category]]
do.call(rbind, lapply(names(sublist), function(feature) {
df <- sublist[[feature]]
if(nrow(df) == 0) return(NULL)
df$Feature <- feature
# Status labeling
if(category == "missing_per_gene") {
df$Status <- "Absent"
} else { # toxins
df$Status <- "Present"
}
return(df)
}))
}))
# Reorder columns
combined_long <- combined_long[, c("Sequence_ID_Lmono","ST","CT","Feature","Status")]
# Export CSV
write.csv(
combined_long,
"~/Desktop/interesting_features_all_isolates.csv",
row.names = FALSE
)
## -------------------------------------------------------------------------------------------------------------------------------------
## source: http://bacmet.biomedicine.gu.se/
# Define gene categories
# Downloaded database path
bacmet_fasta <- "~/Desktop/BacMet2_predicted_database.fasta"
# Read fasta file
bacmet <- readAAStringSet(bacmet_fasta)
# Inspect
bacmet
names(bacmet)[1:10]
# Extract functional information from headers
# gene: last [ ... ] in the header
gene <- str_extract(names(bacmet), '\\[([^\\]]+)\\]$') %>%
str_remove_all('\\[|\\]') %>%
toupper()
# organism: [ ... ] just before last one
organism <- str_extract(names(bacmet), '\\[([^\\]]+)\\](?= \\[[^\\]]+\\]$)') %>%
str_remove_all('\\[|\\]')
# function: text between last | and first [
function_desc <- str_extract(names(bacmet), '(?<=\\| ).*?(?= \\[)') %>%
str_trim()
# Make dataframe
bacmet_df <- tibble(
header = names(bacmet),
sequence = as.character(bacmet),
gene = gene,
organism = organism,
`function` = function_desc
)
# Rename dataframe for clarity
bacmet_database_info <- bacmet_df
# Subset info needed for plotting: only gene and function
bacmet_info_needed <- bacmet_database_info %>%
select("gene", "function")
# Inspect first 10 rows
head(bacmet_info_needed, 10) ## looks good
# Unique genes in BacMet info
unique_genes <- unique(bacmet_info_needed$gene)
## get gene names from gene presence absence matrix
gene_columns <- setdiff(colnames(bacmet_df), c("Isolate", "Phenotype"))
gene_columns
colnames(bacmet_df)
#make it upper letters
gene_columns_upper <- toupper(c(
"acn", "amvA", "arsA", "arsB", "arsC", "arsD", "arsR", "arsT",
"baeR", "bcrA", "bcrB", "bcrC", "bfrA", "blt", "cadC", "cdeA",
"chrR", "chtR", "chtS", "copA", "copB", "copC", "copR", "copS",
"copY/tcrY", "corC", "corT/coaT", "cpxR", "csoR", "cueR/ybbI", "czcR", "dpr/dps",
"emmdR", "emrBsm", "emrD-3", "evgA", "fabK", "fabL/ygaA", "farR", "fecD",
"fecE", "fetA/ybbL", "fetB/ybbM", "fieF/yiip", "gadA", "gadB", "gadC/xasA", "galE",
"glpF", "hmrR", "klaB/telA/kilB", "lde", "lmrS", "mdeA", "mdrL/yfmO", "mepA",
"merA", "mgtA", "mntH/yfeP", "modA", "modB", "nfsA", "nikB", "nikC",
"nikD", "nikE", "norB", "pdrM", "perR", "pstA", "pstB", "qacG",
"qacJ", "qacZ", "recG", "rpoS", "ruvB", "sh-fabI", "sitA", "sitC",
"smdA", "smdB", "smfY", "smrA", "sodA", "sugE", "troB", "tupC",
"vcrM", "wtpC", "ybtP", "ychH", "ycnJ", "yfeA", "yfeC", "yfeD",
"ziaA", "ziaR", "zitB/ybgR", "zntR/yhdM", "znuB/yebI"
))
gene_columns_upper
# Filter to keep only genes in gene_columns_upper
bacmet_filtered <- bacmet_info_needed %>%
filter(gene %in% gene_columns_upper) %>%
distinct(`function`, .keep_all = TRUE) # use backticks for the function column
# View the summarized result
bacmet_summary <- bacmet_filtered %>%
group_by(gene) %>%
summarise(
# collapse the 'function' column with commas
function_summary = paste(`function`, collapse = ", "),
# If you want, you can summarize other columns the same way:
# organism_summary = paste(unique(organism), collapse = ", "),
.groups = "drop"
)
# View the result
bacmet_summary
# Export to CSV
write_csv(bacmet_summary, "~/Desktop/bacmet_summary.csv")
## start visualizing
# Read bacmed data
bacmed_summary <- read.csv(
"~/Desktop/bacmet_matrix.csv",
sep = ";", header = TRUE, check.names = FALSE
)
## make groups
bacmet_gene_groups <- list(
arsenic_resistance = c("arsA","arsB","arsC","arsD","arsR","arsT"),
multidrug_efflux = c("acn","amvA","emmdR","emrBsm","emrD-3","evgA","farR","smdA","smdB","smfY","smrA"),
copper_resistance = c("copA","copB","copC","copR","copS","copY/tcrY"),
zinc_resistance = c("zntR/yhdM","znuB/yebI","ziaA","ziaR","zitB/ybgR"),
iron_homeostasis = c("fetA/ybbL","fetB/ybbM","fieF/yiip"),
cadmium_resistance = c("cadC","cdeA"),
stress_response = c("baeR","chrR","chtR","chtS","cpxR","csoR","cueR/ybbI"),
`antibiotic_target/fatty_acid_biosynthesis` = c("fabK","fabL/ygaA")
)
# Convert columns to factors
bacmed_summary$Phenotype <- factor(bacmed_summary$Phenotype,
levels = c("Persistent", "Non-Persistent"))
bacmed_summary$ST <- factor(bacmed_summary$ST)
# Build matrix for heatmap
bac_mat <- bacmed_summary %>%
select(-Isolate, -Phenotype, -ST) %>%
as.matrix()
rownames(bac_mat) <- bacmed_summary$Isolate
# Order isolates first by Phenotype, then by ST
bacmed_summary <- bacmed_summary %>%
arrange(Phenotype, ST)
bac_mat <- bac_mat[bacmed_summary$Isolate, ]
# Map genes to groups
gene_group_mapping <- unlist(lapply(names(bacmet_gene_groups), function(group) {
setNames(rep(group, length(bacmet_gene_groups[[group]])), bacmet_gene_groups[[group]])
}))
gene_group_mapping <- gene_group_mapping[colnames(bac_mat)]
# Order genes by group
genes_ordered <- names(sort(gene_group_mapping))
bac_mat <- bac_mat[, genes_ordered]
# Row annotation for Phenotype
row_ha <- rowAnnotation(
Phenotype = bacmed_summary$Phenotype,
col = list(Phenotype = c("Persistent"="#440154FF", "Non-Persistent"="#55C667FF")),
annotation_legend_param = list(
title_gp = gpar(fontsize = 12, fontface = "bold"),
labels_gp = gpar(fontsize = 12)
)
)
# Heatmap colors
heat_col <- c("0" = "#D3D3D3", "1" = "#87CEFA")
# Column split by gene group
col_split <- factor(gene_group_mapping[genes_ordered])
# Replace rownames with ST for display
rownames(bac_mat) <- bacmed_summary$ST
# Create heatmap
ht <- Heatmap(
bac_mat,
name = "Presence",
col = heat_col,
cluster_rows = FALSE,
cluster_columns = FALSE,
show_row_names = TRUE,
row_names_side = "left",
row_names_gp = gpar(fontsize=10),
show_column_names = TRUE,
column_names_gp = gpar(fontsize=12),
column_title_rot = 45,
column_split = col_split,
heatmap_legend_param = list(
title_gp = gpar(fontsize=12, fontface="bold"),
labels_gp = gpar(fontsize=12)
)
)
# Draw and save heatmap
svg("~/Desktop/bacmed_heatmap_final_centered.svg", width=14, height=8)
draw(row_ha + ht,
heatmap_legend_side = "right",
padding = unit(c(0.5, 1, 0.5, 1), "mm"))
dev.off()
## -------------------------------------------------------------------------------------------------------------------------------------
# Genes to analyze (with correct stress response gene name)
env_genes <- c(
"arsA", "arsC", "arsD", "arsR",
"cadC",
"copB", "copC", "copY/tcrY",
"zntR/yhdM", "ziaA", "ziaR",
"fetA/ybbL", "fetB/ybbM",
"amvA", "evgA", "farR", "smrA",
"cueR/ybbI"
)
# Define the gene groups
bacmet_gene_groups <- list(
arsenic_resistance = c("arsA","arsC","arsD","arsR"),
multidrug_efflux = c("amvA","evgA","farR","smrA"),
copper_resistance = c("copB","copC","copY/tcrY"),
zinc_resistance = c("zntR/yhdM","ziaA","ziaR"),
iron_homeostasis = c("fetA/ybbL","fetB/ybbM"),
cadmium_resistance = c("cadC"),
stress_response = c("cueR/ybbI")
)
# Load bacmet summary
bacmet_summary <- read.csv(
"~/Desktop/bacmet_matrix.csv",
sep = ";", header = TRUE, check.names = FALSE
)
# Load metadata
metadata <- read_excel(
"~/Desktop/AllPartners_FoodSafeR_WP2.3_Facility_Sampling_Metadata_FS_Amplicon_controls_historical_isolates.xlsx")
metadata$Sequence_ID_Lmono <- trimws(metadata$Sequence_ID_Lmono)
# Map isolate IDs to metadata
map_to_metadata <- function(isolate_ids, metadata) {
merge(
data.frame(Sequence_ID_Lmono = isolate_ids),
metadata[, c("Sequence_ID_Lmono", "ST", "CT")],
by = "Sequence_ID_Lmono",
all.x = TRUE
)
}
# Filter persistent isolates
persistent <- bacmet_summary %>%
filter(Phenotype == "Persistent")
# Function to find isolates harboring a gene
find_present <- function(gene) {
gene_cols <- grep(paste0("^", gene, "$"), colnames(persistent), value = TRUE)
if(length(gene_cols) == 0) return(data.frame())
isolates <- persistent$Isolate[rowSums(persistent[, gene_cols, drop = FALSE]) >= 1]
isolates <- isolates[!startsWith(isolates, "MS")]
map_to_metadata(isolates, metadata)
}
# Run for all genes
harbor_env_genes <- lapply(env_genes, find_present)
names(harbor_env_genes) <- env_genes
# Combine into one long table
combined_env_long <- do.call(rbind, lapply(names(harbor_env_genes), function(feature) {
df <- harbor_env_genes[[feature]]
if(nrow(df) == 0) return(NULL)
df$Feature <- feature
df$Status <- "Present"
# Assign gene group
df$Gene_Group <- sapply(feature, function(f) {
group <- names(bacmet_gene_groups)[sapply(bacmet_gene_groups, function(g) f %in% g)]
if(length(group) == 0) return(NA) else return(group)
})
return(df)
}))
# Reorder columns
combined_env_long <- combined_env_long[, c("Sequence_ID_Lmono","ST","CT","Feature","Gene_Group","Status")]
# Export CSV
write.csv(
combined_env_long,
"~/Desktop/persistent_env_genes_present.csv",
row.names = FALSE
)
## -------------------------------------------------------------------------------------------------------------------------------------
## code based on tutorial https://arftrhmn.net/creating-a-publication-quality-phylogeny-using-ggtree/
# Load tree
tree <- read.iqtree("~/Desktop/listeria_tree.treefile")
tree@phylo <- midpoint(tree@phylo)
# Load and clean metadata
meta <- read_excel("~/Desktop/AllPartners_FoodSafeR_WP2.3_Facility_Sampling_Metadata_FS_Amplicon_controls_historical_isolates.xlsx") %>%
filter(!is.na(Sequence_ID_Lmono)) %>%
filter(Sequence_ID_Lmono != "NA") %>%
arrange(Sequence_ID_Lmono, !is.na(ST)) %>%
distinct(Sequence_ID_Lmono, .keep_all = TRUE) %>%
as.data.frame()
# Build base tree - circular layout
t1 <- ggtree(tree, layout = "circular")
# Join metadata via left_join
t1$data <- t1$data %>%
left_join(meta, by = c("label" = "Sequence_ID_Lmono"))
# Add bootstrap support
t1$data$bootstrap <- '0'
t1$data[which(t1$data$SH_aLRT >= 70 & t1$data$UFboot >= 70),]$bootstrap <- '1'
# Add tip points and bootstrap coloring
t1 <- t1 +
geom_tippoint(size = 0.75 ) +
new_scale_color() +
geom_tree(aes(color = bootstrap == '1')) +
scale_color_manual(name = 'Bootstrap',
values = setNames(c('black', 'grey'), c(T, F)),
guide = "none")
t1
# Prepare heatmap dataframes
meta.persistent <- as.data.frame(meta[, 'Persistent'])
colnames(meta.persistent) <- 'Persistent'
rownames(meta.persistent) <- meta$Sequence_ID_Lmono
meta.st <- as.data.frame(meta[, 'ST'])
colnames(meta.st) <- 'ST'
rownames(meta.st) <- meta$Sequence_ID_Lmono
meta.facility <- as.data.frame(meta[, 'Facility_Type'])
colnames(meta.facility) <- 'Facility_Type'
rownames(meta.facility) <- meta$Sequence_ID_Lmono
# Sort ST levels numerically for legend
st_levels <- meta.st$ST[!is.na(meta.st$ST)] %>%
unique() %>%
as.numeric() %>%
sort() %>%
as.character()
# Ring 1 - Persistent
t2 <- gheatmap(t1, meta.persistent, width = 0.06, offset = 0.005,
colnames = FALSE) +
scale_fill_manual(name = "Phenotype",
breaks = c("Persistent", "Non-persistent"),
labels = c("Persistent", "Non-persistent"),
values = c("Persistent" = "black",
"Non-persistent" = "#D3D3D3"),
na.value = "grey90") +
new_scale_fill()
# Ring 2 - ST (middle ring)
t2 <- gheatmap(t2, meta.st, width = 0.06, offset = 0.015,
colnames = FALSE) +
scale_fill_viridis_d(option = "H", name = "ST",
limits = st_levels) +
new_scale_fill()
# Ring 3 - Facility Type (outer ring)
t2 <- gheatmap(t2, meta.facility, width = 0.06, offset = 0.025,
colnames = FALSE) +
scale_fill_manual(name = "Facility Type",
values = c("Meat" = "#404788",
"Dairy" = "#55C667",
"Fish" = "#FDE725"),
na.value = "grey90")
# Legend arrangement
t2 <- t2 +
theme(
legend.position = "right",
legend.box = "vertical",
legend.text = element_text(size = 14), # bigger text
legend.title = element_text(size = 18, face = "bold"), # bigger title
legend.key.size = unit(0.6, "cm"),
legend.spacing.y = unit(0.5, "cm")
)
ggsave(
"~/Desktop/listeria_tree.png",
plot = t2,
width = 16,
height = 14,
dpi = 300,
units = "in"
)