Additional file 3 All R code used to generate results in this study

title: "S.ruminantium_analysis"
author: "benV"
date: "29/07/2023"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r}
##### Build PopPUNK tree

### Set up R env
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

### Load libraries
library(readxl)
library(tidyverse)
library(ggtree)
library(colorspace)
library(viridis)
library(ggnewscale)
library(ape)
library(aplot)

### Load data

# NJ PopPUNK tree
NJ_core_tree <- read.tree("PopPUNK/vis_refined_core_NJ.nwk")

# Metadata
mdata <- read_excel("../0.Paper/Data_S1.xlsx", sheet = "Sheet1")

### Manip dfs

# Root tree to mid-root point (due to no outgroup)  
NJ_core_tree_mid <- phytools::midpoint.root(NJ_core_tree)

# Prepare table for heatmap
location_table <- mdata %>%
  select(1, 4, 5) %>%
  mutate(Year = as.character(Year)) %>%
  pivot_longer(c(2:3), names_to = "grouping", values_to = "grouping_value")

### Plot

# Construct metadata tree
NJ_core_tree_mid_mdata <- ggtree(NJ_core_tree_mid, layout="rectangular", size = 0.5) %<+% mdata +
  geom_tiplab(align=TRUE, offset=.00009) +
  geom_tippoint(aes(shape = Source, fill = as.factor(PopPUNK_cluster)), size=4.5, alpha = 1) + # ST colour
  scale_shape_manual(values = c("Cow" = 24, "Sheep" = 22)) +
  scale_fill_viridis("PopPUNK cluster", option="turbo", discrete=TRUE) +
  geom_treescale(x = 0.001, y = 1, offset = 0.25, fontsize = 4) +
  theme_tree() +
  xlim(0, 0.00598) + # this one for labelled heatmap
  guides(fill = guide_legend(override.aes = list(shape = 21)))
NJ_core_tree_mid_mdata

# Make heatmap of PopPUNK clusters
tree_heatmap <- ggplot(location_table, (aes(x = grouping, y = Isolate))) +
  geom_tile(fill = NA, colour = "grey") +
  geom_text(aes(label = grouping_value), color = "black", size = 4) +
  # scale_fill_brewer("Location", palette = "RdYlBu") + 
  theme_tree2(axis.text.x = element_text(angle = 45, hjust=1)) +
  coord_fixed(ratio = 1/6) # this sets the total width of the heatmap, nice. # ratio = 0.06 xlim = c(1, 2), ylim = c(0, 41)
tree_heatmap

## Combine all with aplot

# Join heatmaps and tree together
ruminantium_tree_hmap <- tree_heatmap %>%
  insert_left(NJ_core_tree_mid_mdata)
ruminantium_tree_hmap
ggsave(file="phylogeny/ruminantium_tree_hmap.svg", plot=ruminantium_tree_hmap, width=14, height=9)

ggsave(file="phylogeny/ruminantium_tree_hmap.png", plot=ruminantium_tree_hmap, width=9, height=10)
ggsave(file="phylogeny/ruminantium_tree_hmap.pdf", plot=ruminantium_tree_hmap, width=9, height=10)


ggsave(file="03.02.22_phylogeny/reroot_NJ_core_tree_7882_mdata_2.png", plot=reroot_NJ_core_tree_7882_mdata_2, width=10, height=8)


```


```{r}
##### Examine distribution of SNVs from SKA by country

### Set up R env
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

### Load libraries
library(tidyverse)

### Load data

# Create vertice dataframe https://igraph.org/r/doc/graph_from_data_frame.html
isolates <- read.delim("SKA/SNV_network_vertices.txt")

# Create edge dataframe
relationships <- read.delim("SKA/SNV_network_edges.txt")

### Manip dfs

# Prep df for SNV distribution by country
SNV_distributions_df <- relationships %>%
  left_join(isolates, by = c("Isolate_1" = "Isolate")) %>%
  rename(Location_Isolate_1 = Location) %>%
  left_join(isolates, by = c("Isolate_2" = "Isolate")) %>%
  rename(Location_Isolate_2 = Location) %>%
  select(1:8, 12, 13, 17) %>%
  mutate(Same_location = Location_Isolate_1 == Location_Isolate_2) %>%
  filter(Same_location == TRUE) %>%
  filter(Isolate_1 != Isolate_2)

SNV_distributions_df_sard <- SNV_distributions_df %>%
  filter(Location_Isolate_1 == "Sardinia, Italy")

### Plot

# Histogram of SNVs
SNV_country_histo <- ggplot(SNV_distributions_df) +
  geom_histogram(aes(x = SNVs), binwidth = 15, colour = "dark red") +
  facet_grid( ~ Location_Isolate_1) +
  scale_x_continuous("Pairwise SNVs") +
  scale_y_continuous("Number of pairwise genomes") +
  theme_bw()
SNV_country_histo
ggsave("SKA/SNV_country_histo.png", plot = SNV_country_histo, width=12, height=6)
ggsave("SKA/SNV_country_histo.pdf", plot = SNV_country_histo, width=12, height=6)


```


```{r}
##### Build network of pairwise SNPs

### Set up R env
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

### Load libraries
library(igraph)
library(ggraph)
library(RColorBrewer)
library(ggnewscale)
library(ggforce)

### Load data

# Create vertice dataframe https://igraph.org/r/doc/graph_from_data_frame.html
isolates <- read.delim("SKA/SNV_network_vertices.txt")

# Create edge dataframe
relationships <- read.delim("SKA/SNV_network_edges.txt")

### Manip dfs
relationships <- relationships %>%
  filter(Rank <= 53)

### Plot

# Create igraph network
g <- graph_from_data_frame(relationships, directed=FALSE, vertices = isolates)

# Graph/visualise and annotate using ggraph
set.seed(9998) # set seed for consistent layout, otherwise randomness occurs
set.seed(NULL)

SNV_ruminantium_network <- ggraph(g,
                                  layout = 'kk',
                                  # weights = (relationships$SNVs * 0.00005)
) + #this is the layout, there are a bunch to choose from including: 'fr', 'kk', 'stress', 'lgl', 'graphopt', 'circle', 'mds', 'gem', 'star', 'dh', 'grid', 'randomly', 'drl', 'linear'. You can also use layout = 'kk', maxiter = n) for extra variability. For directed graphs, you can use 'partition'.
  geom_edge_link(aes(linetype = as.factor(Rank),
                     alpha = as.factor(Rank))) +
  scale_edge_linetype_manual("Transmission (SNVs)",
                             values = c("15" = 1, "53" = 2, "1000" = 0), 
                             labels = c("\u2264 15", "\u2264 53", "\u2265 1000")) +
  scale_edge_alpha_manual("Transmission (SNVs)",
                          values = c("15" = 1, "53" = 0.5, "1000" = 0),
                          labels = c("\u2264 15", "\u2264 53", "\u2265 1000")) +
  geom_node_point(aes(fill = isolates$Year), shape = 21, size = 3) +
  scale_fill_distiller("Year of isolation", palette = "Blues", direction = 1) +
  facet_nodes(~ Location) +
  theme_bw() +
  coord_equal() +
  theme_no_axes()
SNV_ruminantium_network
ggsave(file="SKA/SNV_ruminantium_network.png", plot= SNV_ruminantium_network, width=12, height=10)


## COULD ALSO USE EDGES AS WEIGHTS:
# ggraph(g, layout = 'fr', weights = (relationships$Pairwise.SNPs*0.05))



```


```{r}
##### Just do network for Italy

### Set up R env
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

### Load libraries
library(igraph)
library(ggraph)
library(ggplot2)
library(RColorBrewer)
library(ggnewscale)
library(ggforce)
library(GGally)

### Load data

# Create vertice dataframe https://igraph.org/r/doc/graph_from_data_frame.html
Italy_nodes <- read.delim("SKA/SNV_network_vertices.txt")


# Create edge dataframe
Italy_edges <- read.delim("SKA/SNV_network_edges_Italy.txt")

### Manip dfs
Italy_edges <- Italy_edges %>%
  filter(Rank <= 53)
# Sardinia only nodes
Italy_nodes <- Italy_nodes %>%
  filter(Location == "Sardinia, Italy")

### Plot

# Create igraph network
g <- graph_from_data_frame(Italy_edges, directed=FALSE, vertices = Italy_nodes)

# Graph/visualise and annotate using ggraph
set.seed(9998) # set seed for consistent layout, otherwise randomness occurs
set.seed(NULL)

Italy_SNV_ruminantium_network <- ggraph(g,
                                  layout = 'stress',
                                  # layout = "fr", weights = (Italy_edges$SNVs*0.00005)
                                  # weights = (1/relationships$Rank)
                                  ) + #this is the layout, there are a bunch to choose from including: 'fr', 'kk', 'stress', 'lgl', 'graphopt', 'circle', 'mds', 'gem', 'star', 'dh', 'grid', 'randomly', 'drl', 'linear'. You can also use layout = 'kk', maxiter = n) for extra variability. For directed graphs, you can use 'partition'.
  geom_edge_link(aes(# colour = factor(Location),
                     # width = edges$snp_dist,
                     linetype = as.factor(Italy_edges$Rank), alpha = as.factor(Italy_edges$Rank))) + # this is the line. Can use geom_edge_fan() for bent line instead of straight
  # scale_edge_width(name = "SNP distances between isolates", range = c(2,1)) +
  # scale_edge_linetype_manual("Transmission events", values = c("10" = "Yes", "\u003E 10" = "No", "1000" = 0)) +
  scale_edge_linetype_manual("Transmission (SNVs)",
                             values = c("15" = 1, "53" = 2, "1000" = 0), 
                             labels = c("\u2264 15", "\u2264 53", "\u2265 1000")) +
  scale_edge_alpha_manual("Transmission (SNVs)",
                          values = c("15" = 1, "53" = 0.5, "1000" = 0),
                          labels = c("\u2264 15", "\u2264 53", "\u2265 1000")) +
  geom_node_point(aes(fill = as.factor(Italy_nodes$Year)), shape = 21, size = 5) + # this is where you can use a node point (also node size) instead of the labels I used here!
  scale_fill_brewer("Year of isolation", palette="Blues") +
  theme_bw() +
  coord_equal() +
  theme_no_axes()
Italy_SNV_ruminantium_network
ggsave(file="SKA/Italy_SNV_ruminantium_network.png", plot= Italy_SNV_ruminantium_network, width=9, height=7)
ggsave(file="SKA/Italy_SNV_ruminantium_network.pdf", plot= Italy_SNV_ruminantium_network, width=9, height=7)

```

```{r}
##### Try network on map using GGally - https://ggobi.github.io/ggally/articles/ggnetworkmap.html: Overall, I am sorta not a huge fan. It works, but not pretty.

library(GGally)
library(network)
library(sna)
library(tidyverse)
library(geosphere)

# Create vertice dataframe https://igraph.org/r/doc/graph_from_data_frame.html
Italy_nodes <- read.delim("SKA/SNV_network_vertices.txt")

# Create edge dataframe
Italy_edges <- read.delim("SKA/SNV_network_edges_Italy.txt")

### Manip dfs

# Sardinia only nodes
Italy_nodes <- Italy_nodes %>%
  filter(Location == "Sardinia, Italy")
rownames(Italy_nodes) <- Italy_nodes$Isolate

## Generate a list of lat/long for each node

# Function to add random variation to a value
add_random_variation <- function(value, variation_range) {
  value + runif(1, -variation_range, variation_range)
}

# Specify the latitude and longitude variation range
latitude_variation <- 0.5
longitude_variation <- 0.5

# Add random variation to latitude and longitude columns
Italy_nodes_lat_long <- Italy_nodes %>%
  # mutate(Latitude = sapply(rep(39.75, nrow(.)), add_random_variation, latitude_variation),
  #        Longitude = sapply(rep(9.45, nrow(.)), add_random_variation, longitude_variation))
  mutate(Latitude = sapply(rep(40, nrow(.)), add_random_variation, latitude_variation),
         Longitude = sapply(rep(9, nrow(.)), add_random_variation, longitude_variation))

# Construct network
Italy_network <- network(Italy_edges, directed = FALSE)

# add geographic coordinates
Italy_network %v% "lat" <- Italy_nodes_lat_long[network.vertex.names(Italy_network), "Latitude"]
Italy_network %v% "lon" <- Italy_nodes_lat_long[network.vertex.names(Italy_network), "Longitude"]

# compute degree centrality
Italy_network %v% "degree" <- degree(Italy_network, gmode = "digraph")

# Add year values
year_values <- Italy_nodes$Year

# Assign the "Year" values as a vertex attribute in the Italy_network graph object
Italy_network %v% "Year" <- year_values

## Subset the Italian map via latitude and longitude!

# Load the map data for Italy
italy_map_data <- map_data("italy")

# Set min and max coords
latitude_min <- 38.5
latitude_max <- 41.5
longitude_min <- 8.0
longitude_max <- 10

# Subset the map data to include only the desired area
sardinia_map_data <- subset(italy_map_data,
                            lat >= latitude_min & lat <= latitude_max &
                            long >= longitude_min & long <= longitude_max)

# # Filter the data to include only polygons for Sardinia - alternative way by region.
# sardinia_map_data <- subset(italy_map_data, region == c("Sassari", "Nuoro", "Oristano", "Cagliari"))

# Plot Sardinia
sardinia_plot <- ggplot(sardinia_map_data, aes(x = long, y = lat)) +
  geom_polygon(aes(group = group), color = "grey65", fill = "#f9f9f9", size = 0.2) +
  coord_fixed(1.3) +  # Optional: Adjust the aspect ratio for better display
  scale_x_continuous("Longitude") +
  scale_y_continuous("Latitude")
sardinia_plot

# overlay network data to map
Italy_network_map <- ggnetworkmap(
  sardinia_plot, Italy_network, size = 4, great.circles = FALSE,
  node.group = Year, segment.color = "black"
  # ring.group = degree, weight = degree
)
Italy_network_map
ggsave(file="SKA/Italy_network_map.png", plot= Italy_network_map, width=12, height=10)

```


```{r}

### Load libraries
library(ggmap)
library(maps)
library(mapdata)
library(ggplot2)
library(viridis)

### Load data
isolate_long_lat <- read_delim("isolate_long_lat.txt") %>%
  arrange(PopPUNK_cluster_bgmm) # sort for nexts step of colour stuff

# Produce data for Sardinia
sardinia <- map_data("italy", region = c("Sassari", "Oristano", "Nuoro", "Cagliari"))

# Produce labels for region for map - central
cnames <- aggregate(cbind(long, lat) ~ region, data=sardinia, 
                    FUN=function(x)mean(range(x)))
cnames[2,2] <- cnames[2,2]+0.3 # move the Nuoro heading further East
### Produce map

# Visualise map
sardinia_map <- ggplot(sardinia) +
  geom_polygon(aes(x = long, y = lat, group = group), colour = "black", fill = "light grey") +
  geom_text(data = cnames, aes(x = long, y = lat, label = region), fontface ="bold") +
  scale_x_continuous("Longitude") +
  scale_y_continuous("Latitude")
sardinia_map


```

```{r}
##### Putative virulence gene heatmap

### Set up R env
setwd(dirname(rstudioapi::getSourceEditorContext()$path))

### Load libraries
library(tidyverse)
library(RColorBrewer)
library(cluster) 

### Load data

# virulence data
virul_df <- read.delim("virulence/virulence_df_2.txt")

### Manip df

# Convert to binary for clustering
virul_df <- virul_df %>%
  mutate_all(~ ifelse(. == "Present", 1, .))

# prep df for plotting
virul_df_aes <- virul_df %>%
  select(1, 7:8) %>%
  pivot_wider(names_from = Gene, values_from = Presence, id_cols = filename) %>%
  mutate_all(~ ifelse(is.na(.), 0, .)) %>%
  pivot_longer(c(2:ncol(.)), names_to = "Gene", values_to = "Presence")


# Perform clustering of plot..
virul_df_complete <- virul_df %>%
  select(1, 7:8) %>%
  pivot_wider(names_from = Gene, values_from = Presence, id_cols = filename) %>%
  mutate_all(~ ifelse(is.na(.), 0, .))

virul_df_complete <- virul_df_complete %>%
    t() %>%
    data.frame()

# Move row 1 to colnames
colnames(virul_df_complete) <- virul_df_complete[1,]
# Remove row 1
virul_df_complete <- virul_df_complete[-1, ] 
# Convert df to factor
virul_df_complete[] <- lapply(virul_df_complete, factor)

# Create df for KOs as row for clustering analysis
isolate_row <- virul_df_complete %>%
 # mutate_at(vars(1:212), factor) %>%
 t() %>%
 as.data.frame()
isolate_row[] <- lapply(isolate_row, factor) # must all be factor for daisy to analyse

# Create df for KOs as row for clustering analysis
gene_row <- virul_df_complete %>%
 # mutate_at(vars(1:212), factor) %>%
 as.data.frame()
gene_row[] <- lapply(gene_row, factor) # must all be factor for daisy to analyse

# Perform gower dissimilarity analysis for input to hclust
isolate.gower.dist <- daisy(isolate_row, metric = c("gower"))
gene.gower.dist <- daisy(gene_row, metric = c("gower"))

# Run hclust on gower output
isolate.aggl.clust.c <- hclust(isolate.gower.dist, method = "complete")
gene.aggl.clust.c <- hclust(gene.gower.dist, method = "complete")

### Plot

vir_hmap <- ggplot(virul_df_aes, aes(x = Gene, y = filename, fill = as.factor(Presence))) +
  geom_tile(colour = "light grey", linewidth = 0.1) +
  scale_fill_manual("Presence of putative\nvirulence factor", values = c("0" = "white", "1" = "black")) +
  scale_x_discrete("Putative virulence gene", limits = rownames(gene_row)[gene.aggl.clust.c$order], expand=c(0,0)) +
  scale_y_discrete("Isolate", limits = rownames(isolate_row)[isolate.aggl.clust.c$order], expand=c(0,0)) + # expand allows cells   to go to borders
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
vir_hmap
ggsave(file="virulence/vir_heatmap.png", plot = vir_hmap, width=15, height=7)
ggsave(file="virulence/vir_heatmap.pdf", plot = vir_hmap, width=15, height=7)

```