######################
#   FISH Paper       #
#   R code           #
#   2026             #
######################

# This R code was used to generate the figures and supplementary figures shown in the manuscript 'Autofluorescence and fluorescence-based detection of anaerobic gut fungi (phylum Neocallimastigomycota)'
# Description of the data files:
	# Autofluorescence_total_Data_for_R.txt --> Total autofluorescence images taken on CLSM, Data analysed with ImageJ following the protocol provided in supplementary
	# Autofluorescence_Data_lasX_for_R_including_stains.txt --> Lamda scans produced on CLSM, Data analysed with LasX
	# Mono-labelled FISH, Data analysed with ImageJ --> Mono-labelled FISH with DY480XL-labelled probes, images taken on CLSM, Data analysed with ImageJ following the protocol provided in supplementary
	# HCR-FISH_data_for_R.txt --> HCR FISH with Cy5-labelled probes, images taken on CLSM, Data analysed with ImageJ following the protocol provided in supplementary
	# HCR-FISH_data_Cy3.txt --> HCR FISH with Cy3-labelled probes, images taken on an epifluorescence microscope, Data analysed with ImageJ, Data analysed with ImageJ following the protocol provided in supplementary
	
library(tidyverse)
library(ggplot2)
library(rstatix)
library(ggpubr)
library(RColorBrewer) 

packageVersion("tidyverse")
# ‘1.3.1’
packageVersion("ggplot2")
# ‘3.3.5’
packageVersion("ggpubr")
# ‘0.4.0’
packageVersion("rstatix")
# ‘0.7.0’

set.seed(2023)

################################
### Data import and Clean-up ###
################################

### Total autofluorescence, Data analysed with ImageJ ###
AF_data <- read_tsv("Autofluorescence_total_Data_for_R.txt")
AF_data$Replicate <- factor(AF_data$Replicate)
AF_data$Medium <- factor(AF_data$Medium)
AF_data$Strain <- factor(AF_data$Strain, levels = c("Anaeromyces", "Caecomyces", "Pecoramyces", "none"))
AF_data$Fixation <- factor(AF_data$Fixation, levels = c("A", "4F", "20F", "37F"), labels = c("acetone", "4% formaldehyde", "20% formaldehyde", "37% formaldehyde"))
AF_data$Stain<- factor(AF_data$Stain)
AF_data$Laser <- factor(AF_data$Laser, levels = c("405", "488", "561", "633"), labels = c("405 nm", "488 nm", "561 nm", "633 nm"))
AF_data$nm <- factor(AF_data$nm)
AF_data$LPGain<- factor(AF_data$LPGain)

### Lamda scans, Data analysed with LasX ###
LasX_data <- read_tsv("Autofluorescence_Data_lasX_for_R_including_stains.txt")
LasX_data$Strain    <- factor(LasX_data$Strain,    levels = c("Anaeromyces", "Caecomyces", "Pecoramyces", "none", "stains"))
LasX_data$Replicate <- factor(LasX_data$Replicate)
LasX_data$Medium    <- factor(LasX_data$Medium)
LasX_data$Fixation  <- factor(LasX_data$Fixation,  levels = c("acetone", "4% formaldehyde", "20% formaldehyde", "37% formaldehyde"))
LasX_data$Stain     <- factor(LasX_data$Stain)
LasX_data$Laser     <- factor(LasX_data$Laser,
                               levels = c("405", "488", "561", "633"),
                               labels = c("405 nm", "488 nm", "561 nm", "633 nm"))


### Mono-labelled FISH, Data analysed with ImageJ ###
Mono_data <- read_tsv("Monolabelled-FISH_results.txt")
Mono_data$Strain <- factor(Mono_data$Strain, levels = c("ViSuPo", "Pecora"), labels = c("Caecomyces", "Pecoramyces"))
Mono_data$Fixation <- factor(Mono_data$Fixation)
Mono_data$Formamide <- factor(Mono_data$Formamide, levels = c("20", "30", "40", "50"), labels = c("20% Formamide", "30% Formamide", "40% Formamide", "50% Formamide"))
Mono_data$Probes <- factor(Mono_data$Probes, levels = c("nonsense-DY480XL", "GGNL4R-DY480XL"))
Mono_data$Stain <- factor(Mono_data$Stain)
Mono_data$Laser <- factor(Mono_data$Laser)
Mono_data$Detection <- factor(Mono_data$Detection)
Mono_data$Replicate <- factor(Mono_data$Replicate)
Mono_data$LPGain <- factor(Mono_data$LPGain)


### HCR FISH with Cy5 (confocal microscope), Data analysed with ImageJ ###
HCR5_data <- read_tsv("HCR-FISH_data_for_R.txt")
HCR5_data$Strain <- factor(HCR5_data$Strain, levels = c("V", "P", "Vws"), labels = c("Caecomyces", "Pecoramyces","Caecomyces on straw"))
HCR5_data$Fixation <- factor(HCR5_data$Fixation)
HCR5_data$Formamide <- factor(HCR5_data$Formamide, levels = c("00%F", "10%F", "20%F", "30%F", "35%F", "40%F", "45%F", "50%F"), labels = c("0% Formamide", "10% Formamide", "20% Formamide", "30% Formamide", "35% Formamide", "40% Formamide", "45% Formamide", "50% Formamide"))
HCR5_data$Additives <- factor(HCR5_data$Additives, levels = c("none", "5µM-initH", "10µM-initH", "1µM-AP", "2.5µM-AP", "10µMinitH-2.5µMAP", "140mM-EDTA", "250mM-EDTA", "250mMEDTA-5µMinitH", "250mMEDTA-10µMinitH", "2mgBSA", "10mgBSA"))
HCR5_data$InitiatorH <- factor(HCR5_data$InitiatorH)
HCR5_data$Amplifierprobes <- factor(HCR5_data$Amplifierprobes)
HCR5_data$Probes <- factor(HCR5_data$Probes, levels = c("none", "initH", "H1H2", "initHH1H2"))
HCR5_data$NuclearStain <- factor(HCR5_data$NuclearStain)
HCR5_data$Laser <- factor(HCR5_data$Laser)
HCR5_data$Detected <- factor(HCR5_data$Detected)
HCR5_data$Replicate <- factor(HCR5_data$Replicate)
HCR5_data$LPGain <- factor(HCR5_data$LPGain, levels = c("15%900V", "15%800V", "10%900V", "10%800V", "5%900V", "5%650V", "3%800V", "1%900V", "1%700V"))

### HCR FISH with Cy3 (epifluorescence microscope), Data analysed with ImageJ ###
HCR3_data <- read_tsv("HCR-FISH_data_Cy3.txt")
HCR3_data$Date <- factor(HCR3_data$Date)
HCR3_data$ROIs <- factor(HCR3_data$ROIs)
HCR3_data$Sample <- factor(HCR3_data$Sample)
HCR3_data$Probes <- factor(HCR3_data$Probes, levels = c("none", "initH", "H1", "H2", "H1H2", "initHH1", "initHH2", "initHH1H2", "initHH1-H1H2"))
HCR3_data$Stain <- factor(HCR3_data$Stain, levels = c("none", "NucBlue", "CalcofluorWhite"))
HCR3_data$Objective <- factor(HCR3_data$Objective)
HCR3_data$Replicate <- factor(HCR3_data$Replicate)
HCR3_data$Filter_Exposure <- factor(HCR3_data$Filter_Exposure, levels = c("BF", "D_50ms1x", "D_100ms1x", "D_500ms1x", "D_1s1x", "D_2s1x", "D_2s3.4x", "B_500ms1x", "B_1s1x", "B_1s1.8x", "B_2s1x", "B_2s3.4x", "B_5s1x", "B_5s3.4x", "B_5s9.3x", "5_1s1x", "5_2s1x", "5_2s3.4x", "5_5s1x", "5_5s3.4x", "5_5s9.3x"))

######################################
##		  Visualize data			##
######################################

##################################################################################################################
### Total autofluorescence, Data analysed with ImageJ ###
#########################################################

#filter data 
LPG <- AF_data %>% filter (LPGain == "15%800V" & Medium == "Minimal")
LPG_nostain <- LPG %>% filter (Stain == "none")

#check for outliers:
LPG_outliers <- LPG %>%  group_by(Laser, Strain, Stain) %>%  identify_outliers(Mean)
LPG_outliers <- LPG_outliers[,-c(18,19)]
LPG_or <- setdiff(LPG, LPG_outliers) 

LPG_nostain_outliers <- LPG_nostain %>%  group_by(Laser, Strain) %>%  identify_outliers(Mean)
LPG_nostain_outliers <- LPG_nostain_outliers[,-c(18,19)] 
LPG_nostain_or <- setdiff(LPG_nostain, LPG_nostain_outliers) 

#filter for each laser line to make separate plots (so the y axis can be adapted)
LPG_nostain_or_405 <- LPG_nostain_or %>% filter (Laser == "405 nm")
LPG_nostain_or_488 <- LPG_nostain_or %>% filter (Laser == "488 nm")
LPG_nostain_or_561 <- LPG_nostain_or %>% filter (Laser == "561 nm")
LPG_nostain_or_633 <- LPG_nostain_or %>% filter (Laser == "633 nm")

## Boxplot to check average fluorescence intensity vs. Strain, Laser, Fixation, Stain: ##
boxplot <- ggplot(LPG_or, aes(x = Stain, y = Mean, fill = Stain)) +
  geom_boxplot() +
  stat_summary(fun.data = function(x) c(y = 195, label = length(x)),
               geom = "text", size = 4) +
  facet_wrap(Strain ~ Laser, scales = "fixed", nrow = 3) +
  coord_cartesian(ylim = c(0, 200)) +
  theme_bw() +
  labs(title = "Mean Fluorescence Intensity by Stain and Laser Line",
       y = "Mean total fluorescence intensity (a.u.)",
       x = "Fixation") +
  theme(axis.text.x = element_blank(),
        legend.position = "right",
        text = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14))
		
boxplot <- boxplot + scale_fill_manual(values = c("grey", "lightblue", "salmon"))
ggsave("Boxplot_AFtotal_Fiji.png", plot = boxplot, width = 10, height = 7, dpi = 300)

## Boxplot to check average fluorescence intensity vs. Strain, Laser, Fixation: ##

# getting the right colors
colors <- brewer.pal(4, "PRGn")

boxplot <- ggplot(LPG_nostain_or, aes(x = Fixation, y = Mean, fill = Fixation)) +
  geom_boxplot() +
  scale_fill_manual(values = colors) +
  stat_summary(fun.data = function(x) c(y = 48, label = length(x)),
               geom = "text", size = 4) +
  facet_wrap(~ Laser, scales = "fixed", nrow = 1) +
  coord_cartesian(ylim = c(0, 50)) +  
  theme_bw() +
  labs(title = "Mean Autofluorescence Intensity by Fixation and Laser Line",
       y = "Mean total fluorescence intensity (a.u.)",
       x = "Fixation") +
  theme(axis.text.x = element_blank(), 
        legend.position = "right",
        text = element_text(size = 12),  
        axis.title = element_text(size = 12, face = "bold"), 
        strip.text = element_text(size = 12, face = "bold"), 
		plot.title = element_text(size = 16, face = "bold", hjust = 0.5)) 

ggsave("Total_autofluorescence_Fixations.png", plot = boxplot, width = 10, height = 6, dpi = 300)
ggsave("Total_autofluorescence_Fixations.jpeg", plot = boxplot, width = 10, height = 6, dpi = 300)
##################################################################################################################


#####################################################
###        Autofluorescence LasX data	    	  ###
#####################################################

# AvSpectra_AF_perStrain.png:

#filter data
none <- LasX_data %>% filter(Stain == "none" & Strain != "none")

bin_width <- 11
average_spectrum <- none %>%
  mutate(binned_wavelength = cut(Wavelength, breaks = seq(min(Wavelength), max(Wavelength), by = bin_width))) %>%
  group_by(Strain, Fixation, Laser, binned_wavelength) %>%
  summarise(mid_wavelength = mean(Wavelength),
            mean_intensity  = mean(nfi, na.rm = TRUE),
            .groups = "drop")

filtered_spectrum <- average_spectrum %>%
  filter(
    (Laser == "405 nm" & mid_wavelength <= 700) |
    (Laser == "488 nm" & mid_wavelength <= 700) |
    (Laser == "561 nm" & mid_wavelength >= 500 & mid_wavelength <= 700) |
    (Laser == "633 nm" & mid_wavelength >= 600 & mid_wavelength <= 700)
  )

##### Stain spectra #####
# In the stains rows, the stain name is stored in the Replicate column
stain_spectra <- LasX_data %>%
  filter(Strain == "stains") %>%
  rename(StainName = Replicate) %>%          # make it explicit
  group_by(StainName, Laser, Wavelength) %>%
  summarise(mean_intensity = mean(nfi, na.rm = TRUE), .groups = "drop") %>%
  # apply the same wavelength windows as the autofluorescence data
  filter(
    (Laser == "405 nm" & Wavelength <= 700) |
    (Laser == "488 nm" & Wavelength <= 700) |
    (Laser == "561 nm" & Wavelength >= 500 & Wavelength <= 700) |
    (Laser == "633 nm" & Wavelength >= 600 & Wavelength <= 700)
  )

##### Vertical lines to indicate laser lines #####
vertical_lines <- data.frame(
  Laser      = c("405 nm", "488 nm", "561 nm", "633 nm"),
  xintercept = c(405, 488, 561, 633)
)
  
##### Plot #####
AF_plot <- ggplot(filtered_spectrum,
                  aes(x = mid_wavelength, y = mean_intensity,
                      color = Fixation, linetype = Fixation)) +
  geom_line(size = 0.8) +
  # Stain spectra overlaid in each facet column (all Strain rows share the same Laser column)
  geom_line(data      = stain_spectra,
            aes(x     = Wavelength,
                y     = mean_intensity,
                group = StainName,
                color = StainName),   # separate colour scale entry
            linetype  = "solid",
            linewidth = 1,
            inherit.aes = FALSE) +
  facet_grid(Strain ~ Laser, scales = "free_x") +
  geom_vline(data = vertical_lines,
             aes(xintercept = xintercept),
             linetype = "dashed", color = "black") +
  scale_color_manual(
	  values = c(
		# Fixation colours (PRGn palette values preserved)
		"acetone"            = "#762a83",
		"4% formaldehyde"    = "#af7ab3",
		"20% formaldehyde"   = "#7fbf7b",
		"37% formaldehyde"   = "#1b7837",
		# Stain colours
		"DAPI"    = "lightblue",
		"NucBlue" = "darkblue",
		"Cy3"     = "darkcyan",
		"DY480XL" = "darkorange",
		"Cy5"     = "red",
		"NucRed"  = "darkred"
	  )
	) +
  labs(x     = "Wavelength (nm)",
       y     = "Mean normalized fluorescence intensity (a.u.)",
       title = "Mean Normalized Autofluorescence Emission Spectra",
       color = "Fixation / Stain") +
  theme_minimal() +
  theme(
    strip.text      = element_text(size = 12, face = "bold"),
    axis.title      = element_text(size = 12, face = "bold"),
    plot.title      = element_text(size = 16, face = "bold", hjust = 0.5),
    legend.title    = element_text(size = 12),
    legend.position = "right",
    panel.spacing   = grid::unit(2, "lines")
  )

ggsave("AvSpectra_AF_perStrain.png", plot = AF_plot, width = 10, height = 6, dpi = 300)
ggsave("AvSpectra_AF_perStrain.jpeg", plot = AF_plot, width = 10, height = 6, dpi = 300)
#########################################################################################################################

##################################################################################################################
### Mono-labelled FISH, Data analysed with ImageJ ###
#####################################################

#filter data
DY480XL <- Mono_data %>% filter (Laser == "488" | Laser == "498")

#check for outliers:
Mono_outliers <- DY480XL %>%  group_by(Strain, Stain, Formamide, Probes) %>%  identify_outliers(Mean)
Mono_outliers <- Mono_outliers[,-c(17,18)] 
Mono_or <- setdiff(DY480XL, Mono_outliers) 

# Calculate n for each group
n_labels <- Mono_or %>%
  group_by(Strain, Formamide, Probes) %>%
  summarise(
    n = sum(!is.na(Mean)),
    y = max(Mean, na.rm = TRUE) * 1.05,  
    .groups = "drop"
  )

  
#visualize
boxplot <- ggboxplot(Mono_or,
          x = "Probes",
          y = "Mean",
          add = "jitter",
          color = "Probes",
          ylab = "Mean fluorescence intensity (a.u.)",
          xlab = "Probes") +
  facet_grid(Strain ~ Formamide) +  
  scale_color_manual(values = c("#696969", "#FF7F50")) +
  coord_cartesian(ylim = c(0, 28)) +
  geom_text(data = n_labels,
            aes(x = Probes,
                y = Inf,
                label = paste0("n = ", n)),
            vjust = 1.5,
            inherit.aes = FALSE) +
  ggtitle(
	expression(
		paste("Mono-labelled FISH for ", italic("Neocallimastigomycota"),
		  " using DY480XL-labelled probes")
	)
  ) +
       theme(axis.text.x = element_blank(),  
           strip.text.y = element_text(face = "italic", size = 11),
		   strip.text.x = element_text(size = 11), 
           axis.title = element_text(size = 12, face = "bold"), 
           plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
           legend.title = element_text(size = 12),
           legend.text = element_text(size = 14),  
           panel.spacing = unit(2, "lines"))  

							
ggsave("Boxplot_Formamide_Mono.png", plot = boxplot, width = 14, height = 5, dpi = 300)
ggsave("Boxplot_Formamide_Mono.jpeg", plot = boxplot, width = 14, height = 5, dpi = 300)
##################################################################################################################


##################################################################################################################
### HCR FISH with Cy5 (confocal microscope), Data analysed with ImageJ ###
##########################################################################

#check for outliers
HCR5_outliers <- HCR5_data %>%  group_by(Strain, Formamide, Additives, Probes, NuclearStain, Laser, LPGain) %>%  identify_outliers(Mean)
HCR5_outliers <- HCR5_outliers[,-c(28,29)] 
HCR5_or <- setdiff(HCR5_data, HCR5_outliers) 

#filter data
HCR5_none <- HCR5_or %>% filter (Additives == "none")
HCR5_none_Cy5 <- HCR5_none %>% filter (Formamide != "45% Formamide" & Formamide != "35% Formamide" & Laser == "633er")
HCR5_opt <- HCR5_or %>% filter (Strain == "Caecomyces" & Laser == "633er" & Additives != "none")
HCR5_straw <- HCR5_or %>% filter (Strain == "Caecomyces on straw" & Laser == "633er")

filtered_data <- HCR5_none_Cy5 %>%
  filter(!(Strain == "Caecomyces" & LPGain == "10%800V" & Formamide %in% c("50% Formamide", "40% Formamide"))) %>%
  filter(!(Strain == "Caecomyces" & LPGain == "5%650V" & Formamide %in% c("0% Formamide", "10% Formamide"))) %>%
  filter(!(Strain == "Caecomyces" & LPGain == "10%800V" & Formamide == "30% Formamide" & Probes %in% c("H1H2", "initHH1H2"))) %>%
  filter(!(Strain == "Pecoramyces" & LPGain == "5%650V")) 

# Calculate n for each group
n_labels <- filtered_data %>%
  group_by(Strain, Formamide, Probes) %>%
  summarise(
    n = sum(!is.na(Mean)),
    y = max(Mean, na.rm = TRUE) * 1.05,  
    .groups = "drop"
  )

#visualize Formamide series
boxplot <- ggboxplot(filtered_data, x = "Probes", y = "Mean", 
           add = "jitter", color = "Probes", shape = "LPGain",
           palette = c("#696969", "#CD0000", "darkred"),
           ylab = "Mean fluorescence intensity (a.u.)",
           xlab = "Probes") +
     facet_grid(Strain ~ Formamide, scales = "free_y") + scale_shape_manual(values = c(21, 25)) + 
     coord_cartesian(ylim = c(0, 250)) +
     geom_text(data = n_labels,
               aes(x = Probes,
                   y = Inf,
                   label = paste0("n = ", n)),
               vjust = 1.5,
               inherit.aes = FALSE) +
     ggtitle(
	  expression(
		paste("HCR-FISH for ", italic("Neocallimastigomycota"),
			  " using Cy5-labelled probes")
	  )
	) +
     theme(axis.text.x = element_blank(),  
           strip.text.y = element_text(face = "italic", size = 11),
		   strip.text.x = element_text(size = 11), 
           axis.title = element_text(size = 12, face = "bold"), 
           plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
           legend.title = element_text(size = 12),
           legend.text = element_text(size = 14),  
           panel.spacing = unit(2, "lines"))  

ggsave("Boxplot_Formamide_HCR-Cy5.png", plot = boxplot, width = 14, height = 5, dpi = 300)
ggsave("Boxplot_Formamide_HCR-Cy5.jpeg", plot = boxplot, width = 14, height = 5, dpi = 300)


####### visualize HCR FISH for Caecomyces on straw ######

filtered_data2 <- HCR5_straw %>%
  filter(!(Additives == "10µMinitH-2.5µMAP" & LPGain %in% c("5%650V", "3%800V"))) %>%
  filter(!(Additives == "250mMEDTA-5µMinitH" & LPGain == "3%800V")) %>%
  filter(!(Additives == "2mgBSA" & LPGain %in% c("5%650V", "3%800V", "1%700V"))) %>%
  filter(!(Additives == "10mgBSA" & LPGain %in% c("5%650V", "3%800V", "1%700V")))
  
# Calculate n for each group
n_labels <- filtered_data2 %>%
  group_by(Strain, Additives, Probes) %>%
  summarise(
    n = sum(!is.na(Mean)),
    y = max(Mean, na.rm = TRUE) * 1.05,  # position above the box
    .groups = "drop"
  )

# visualize
boxplot <- ggboxplot(filtered_data2, x = "Probes", y = "Mean", 
		add = "jitter", color = "Probes",
		palette = c("#696969", "#CD0000", "darkred"),
		ylab = "Mean fluorescence intensity (a.u.)",
		xlab = "Probes") +
	facet_grid(Strain ~ Additives, scales = "free_y") + 
     coord_cartesian(ylim = c(0, 200)) +
     geom_text(data = n_labels,
               aes(x = Probes,
                   y = Inf,
                   label = paste0("n = ", n)),
               vjust = 1.5,
               inherit.aes = FALSE) +
     ggtitle(
	  expression(
		paste("HCR-FISH for ", italic("Caecomyces"),
			  " grown on wheat straw using Cy5-labelled probes")
	  )
	) +
     theme(axis.text.x = element_blank(), 
           strip.text.y = element_text(face = "italic", size = 11),
		   strip.text.x = element_text(size = 11), 
           axis.title = element_text(size = 12, face = "bold"), 
           plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
           legend.title = element_text(size = 12),
           legend.text = element_text(size = 14), 
           panel.spacing = unit(2, "lines")) 

ggsave("Boxplot_HCR_CaecoStraw.png", plot = boxplot, width = 14, height = 5, dpi = 300)
ggsave("Boxplot_HCR_CaecoStraw.jpeg", plot = boxplot, width = 14, height = 5, dpi = 300)

####### visualize HCR FISH Optimization #######

#filter data
filtered_data3 <- HCR5_opt %>%
  filter(!(Additives == "5µM-initH" & LPGain %in% c("5%650V", "3%800V"))) %>%
  filter(!(Additives == "10µM-initH" & LPGain == "5%650V")) %>%
  filter(!(Additives == "1µM-AP" & LPGain == "5%650V")) %>%
  filter(!(Additives == "2.5µM-AP" & LPGain == "3%800V")) %>%
  filter(!(Additives == "250mM-EDTA" & LPGain %in% c("5%650V", "3%800V", "15%800V"))) %>%
  filter(!(Additives == "250mMEDTA-5µMinitH" & LPGain == "1%700V")) %>%
  filter(LPGain == "10%800V")
  
# Calculate n for each group
n_labels <- filtered_data3 %>%
  group_by(Strain, Additives, Probes) %>%
  summarise(
    n = sum(!is.na(Mean)),
    y = max(Mean, na.rm = TRUE) * 1.05, 
    .groups = "drop"
  )

# visualize
boxplot <- ggboxplot(filtered_data3, x = "Probes", y = "Mean", 
		add = "jitter", color = "Probes",
		palette = c("#696969", "#CD0000", "darkred"),
		ylab = "Mean fluorescence intensity (a.u.)",
		xlab = "Probes") +
	facet_grid(Strain ~ Additives) +  
     coord_cartesian(ylim = c(0, 270)) +
     geom_text(data = n_labels,
               aes(x = Probes,
                   y = Inf,
                   label = paste0("n = ", n)),
               vjust = 1.5,
               inherit.aes = FALSE) +
     ggtitle(
	  expression(
		paste("Optimization of HCR-FISH for ", italic("Neocallimastigomycota"),
			  " using Cy5-labelled probes")
	  )
	) +
     theme(axis.text.x = element_blank(),  
           strip.text.y = element_text(face = "italic", size = 11),
		   strip.text.x = element_text(size = 10), 
           axis.title = element_text(size = 12, face = "bold"), 
           plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
           legend.title = element_text(size = 12),
           legend.text = element_text(size = 14),  
           panel.spacing = unit(2, "lines")) 
		
ggsave("Boxplot_Formamide_HCR_optimization.png", plot = boxplot, width = 14, height = 5, dpi = 300)
ggsave("Boxplot_Formamide_HCR_optimization.jpeg", plot = boxplot, width = 14, height = 5, dpi = 300)
##################################################################################################################

##################################################################################################################
### HCR FISH with Cy3 (epifluorescence microscope), Data analysed with ImageJ ###
#################################################################################

#filter data
Cy3_AF <- HCR3_data %>% filter (Sample == "FL1" | Sample == "ViSuPo" & ROIs == "withROIs")
Env <- HCR3_data %>% filter (Sample != "FL1" & Sample != "ViSuPo" & ROIs == "noROIs" )
Env_filt <- Env %>% filter (Filter_Exposure == "D_500ms1x" | Filter_Exposure == "B_2s1x" | Filter_Exposure =="5_2s3.4x")
Cy3_filt <- Cy3_AF %>% filter (Filter_Exposure == "D_500ms1x" | Filter_Exposure == "B_2s1x" | Filter_Exposure =="5_2s3.4x")
Cy3_filt_a <- Cy3_filt %>% filter (Date == "02.05.2023")
Cy3_filt_b <- Cy3_filt %>% filter (Date != "02.05.2023" & Date != "17.04.2023")

#rename for graphs
Cy3_filt_a$Sample <- factor(Cy3_filt_a$Sample, levels = c("ViSuPo", "FL1"), labels = c("Caecomyces", "Neocallimastix"))
Cy3_filt_b$Sample <- factor(Cy3_filt_b$Sample, levels = c("ViSuPo", "FL1"), labels = c("Caecomyces", "Neocallimastix"))
Env_filt$Sample <-factor(Env_filt$Sample, levels = c("Hundekacke", "Gamskacke"), labels = c("fecal sample dog", "fecal sample chamois"))

Cy3_filt_a$Filter_Exposure <- factor(Cy3_filt_a$Filter_Exposure, levels = c("D_500ms1x", "B_2s1x", "5_2s3.4x"), labels = c("NucBlue", "Autofluorescence (B-2A)", "Cy3"))
Cy3_filt_b$Filter_Exposure <- factor(Cy3_filt_b$Filter_Exposure, levels = c("D_500ms1x", "B_2s1x", "5_2s3.4x"), labels = c("NucBlue", "Autofluorescence (B-2A)", "Cy3"))
Cy3_filt_b$Date <- factor(Cy3_filt_b$Date, levels = c("19.04.2023", "09.05.2023"), labels = c("repeated washing", "washing by 'desalting'"))
Env_filt$Filter_Exposure <- factor(Env_filt$Filter_Exposure, levels = c("D_500ms1x", "B_2s1x", "5_2s3.4x"), labels = c("NucBlue", "Autofluorescence (B-2A)", "Cy3"))

####### Cy3 washing optimization experiments ######

#filter for graph
Cy3_filt_a_cy3 <- Cy3_filt_a %>% filter (Filter_Exposure == "Cy3")
Cy3_filt_b_cy3 <- Cy3_filt_b %>% filter (Filter_Exposure == "Cy3")
Env_filt_cy3 <- Env_filt %>% filter (Filter_Exposure == "Cy3")

# Calculate n for each group
n_labels <- Cy3_filt_b_cy3 %>%
  group_by(Sample, Date, Probes) %>%
  summarise(
    n = sum(!is.na(Mean)),
    y = max(Mean, na.rm = TRUE) * 1.05,  # position above the box
    .groups = "drop"
  )

#visualize data
boxplot <- ggboxplot(Cy3_filt_b_cy3, x = "Probes", y = "Mean", 
		add = "jitter", color = "Probes", 
		palette = c("#696969", "#8FBC8F", "#008000"),
		ylab = "Mean fluorescence intentsity (a.u.)",
		xlab = "Probes") +
	facet_grid(Sample ~ Date, scales = "free_y") +  
     coord_cartesian(ylim = c(0, 35)) +
     geom_text(data = n_labels,
               aes(x = Probes,
                   y = Inf,
                   label = paste0("n = ", n)),
               vjust = 1.5,
               inherit.aes = FALSE) +
     ggtitle(
	  expression(
		paste("HCR-FISH for ", italic("Neocallimastigomycota"),
			  " using Cy3-labelled probes: Optimizing washing steps")
	  )
	) +
     theme(axis.text.x = element_blank(),  
           strip.text.y = element_text(face = "italic", size = 11),
		   strip.text.x = element_text(size = 10), 
           axis.title = element_text(size = 12, face = "bold"), 
           plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
           legend.title = element_text(size = 12),
           legend.text = element_text(size = 14), 
           panel.spacing = unit(2, "lines"))  
ggsave("Boxplot_Cy3_washing-optimization.png", plot = boxplot, width = 14, height = 5, dpi = 300)
ggsave("Boxplot_Cy3_washing-optimization.jpeg", plot = boxplot, width = 14, height = 5, dpi = 300)


####### Cy3 environmental samples ######

# Calculate n for each group
n_labels <- Env_filt %>%
  group_by(Sample, Filter_Exposure, Probes) %>%
  summarise(
    n = sum(!is.na(Mean)),
    y = max(Mean, na.rm = TRUE) * 1.05,  # position above the box
    .groups = "drop"
  )
  
#visualize data
boxplot <- ggboxplot(Env_filt, x = "Probes", y = "Mean", 
		add = "jitter", color = "Probes", 
		palette = c("#696969", "#8FBC8F", "#008000"),
		ylab = "Mean fluorescence intensity (a.u.)",
		xlab = "Probes") +
	facet_grid(Sample ~ Filter_Exposure) +
	coord_cartesian(ylim = c(0, 115)) +
     geom_text(data = n_labels,
               aes(x = Probes,
                   y = Inf,
                   label = paste0("n = ", n)),
               vjust = 1.5,
               inherit.aes = FALSE) +
     ggtitle("HCR-FISH with Cy3-labelled probes on environmental samples") +
	theme(axis.text.x = element_blank(),  # Remove x-axis labels
        strip.text.y = element_text(size = 10),
		strip.text.x = element_text(size = 10), , 
        axis.title = element_text(size = 12, face = "bold"), 
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        legend.title = element_text(size = 12),
		legend.text = element_text(size = 14),  # Adjust legend text size
		panel.spacing = unit(2, "lines"))  # Increase spacing between facet panels

ggsave("Boxplot_Envsample_noROIs.jpeg", plot = boxplot, width = 14, height = 5, dpi = 300)
ggsave("Boxplot_Envsample_noROIs.png", plot = boxplot, width = 14, height = 5, dpi = 300)

####### Cy3 a closer look experiment ######

# Calculate n for each group
n_labels <- Cy3_filt_a_cy3 %>%
  group_by(Sample, Probes) %>%
  summarise(
    n = sum(!is.na(Mean)),
    y = max(Mean, na.rm = TRUE) * 1.05,  # position above the box
    .groups = "drop"
  )
boxplot <- ggboxplot(Cy3_filt_a_cy3, x = "Probes", y = "Mean", 
		add = "jitter", color = "Probes", shape = "Probes",
		palette = c("#808080", "#696969", "#8FBC8F", "#3CB371", "#228B22", "#008000", "#556B2F", "#006400" , "black"),
		ylab = "Mean fluorescence intensity (a.u.)",
		xlab = "Probes") +
	facet_wrap(~ Sample) + 
	scale_shape_manual(values = c(21, 22, 23, 24, 25, 21, 22, 23, 24)) +
	coord_cartesian(ylim = c(0, 27)) +
     geom_text(data = n_labels,
               aes(x = Probes,
                   y = Inf,
                   label = paste0("n = ", n)),
               vjust = 1.5,
               inherit.aes = FALSE) +
     ggtitle(
	  expression(
		paste("HCR-FISH for ", italic("Neocallimastigomycota"),
			  " using Cy3-labelled probes: A closer look")
	  )
	) +
	theme(axis.text.x = element_blank(), 
		strip.text.x = element_text(face = "italic", size = 12), , 
        axis.title = element_text(size = 12, face = "bold"), 
        plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
        legend.title = element_text(size = 12),
		legend.text = element_text(size = 14), 
		panel.spacing = unit(2, "lines")) 
		
ggsave("Boxplot_Cy3_probes_a-closer-look.png", plot = boxplot, width = 14, height = 5, dpi = 300)
ggsave("Boxplot_Cy3_probes_a-closer-look.jpeg", plot = boxplot, width = 14, height = 5, dpi = 300)