# Install the full version with data from GitHub
if (!require("devtools")) install.packages("devtools")
devtools::install_github("RightSZ/BioThermR")

# Install dependencies
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("EBImage")

library(ggplot2)
library(ggpubr)
library(BioThermR)

# 1. Get the path to the 30 images included with the package
data_folder <- system.file("extdata", package = "BioThermR")

# 2. Read using the batch reading function
obj_list <- read_thermal_batch(data_folder)
length(obj_list)
#> [1] 29

# 3. Use the lapply function to batch call roi_segment_ebimage for segmentation
# Otsu's method is used by default
obj_list <- lapply(obj_list, roi_segment_ebimage)

# Display the thermal heatmap of one BioThermR object
p1 <- plot_thermal_heatmap(obj_list[[1]], use_processed = FALSE)
p1

# Display the segmented image of one BioThermR object
p2 <- plot_thermal_heatmap(obj_list[[1]], use_processed = TRUE)
p2

# Display the 3d image of one segmented BioThermR object
plot_thermal_3d(obj_list[[1]])

# Display the ROI density plot of one BioThermR object
p3 <- plot_thermal_density(obj_list[[1]])
p3

# Display the overall display montage plot
p4 <- plot_thermal_montage(obj_list, ncol = 5, padding = 2, text_size = 3)
p4

# Display the overall display cloud plot
p5 <- plot_thermal_cloud(obj_list) 
p5

ggpubr::ggarrange(p1,p2,p3,p4,p5,labels = c("A","B","C","D","E"))

# Batch calculate statistics for each BioThermR object
obj_list <- lapply(obj_list, analyze_thermal_stats)

# Compile statistics into a data frame
df <- compile_batch_stats(obj_list)

# Read grouping data and merge
pd_path <- system.file("extdata", "group.csv", package = "BioThermR")
pd <- read.csv(pd_path)
df <- merge_clinical_data(df,pd)

# Aggregate technical replicates (e.g., 3 photos per mouse)
df_new <- aggregate_replicates(df,
                               method = "mean", # use mean method
                               id_col = "Sample", # merge by sample
                               keep_cols = c("Group")
                               )
head(df_new)
#>   Sample Group n_replicates   Pixels      Min      Max     Mean   Median
#> 1  ND4_1   ND4            3 355.3333 29.27191 35.91896 32.84147 32.21257
#> 2  ND4_2   ND4            3 396.0000 29.21382 36.22407 33.27716 32.30919
#> 3  ND4_3   ND4            3 336.3333 30.35210 35.94674 34.28877 35.60987
#> 4  ND4_4   ND4            3 392.0000 30.04441 35.92882 33.42511 32.41518
#> 5  ND4_5   ND4            2 388.0000 31.11169 36.26280 34.78027 35.70538
#> 6   ND_1    ND            3 402.6667 29.36788 35.87173 31.70707 31.55465
#>         SD      Q25      Q75       IQR         CV Peak_Density Weight Glucose
#> 1 1.994384 31.38211 35.57433 4.1922170 0.06072339     32.02256   30.0    3.13
#> 2 2.169935 31.45619 35.67073 4.2145389 0.06520944     35.73361   27.4    2.61
#> 3 1.874250 32.23967 35.70550 3.4658298 0.05466743     35.69518   25.8    3.35
#> 4 2.023300 31.69513 35.66528 3.9701481 0.06053138     34.42208   28.9    3.55
#> 5 1.695506 32.56883 35.85285 3.2840147 0.04874905     35.80615   27.2    3.95
#> 6 1.311257 31.07301 31.86817 0.7951624 0.04135084     31.62128   22.1    5.69
#>   Rep
#> 1   2
#> 2   2
#> 3   2
#> 4   2
#> 5   2
#> 6   2

# Display grouping data plots, grouped by Group
my_comparisons <- list( c("ND4", "ND"))
# Mean
s1 <- viz_thermal_boxplot(data = df_new,
                          y_var = "Mean",
                          x_var = "Group") +
  ggpubr::stat_compare_means(comparisons = my_comparisons,
                             method = "t.test",
                             label = "p.signif",
                             hide.ns = FALSE
                             )

# Median
s2 <- viz_thermal_boxplot(data = df_new,
                          y_var = "Median",
                          x_var = "Group")+
  ggpubr::stat_compare_means(comparisons = my_comparisons,
                             method = "t.test",
                             label = "p.signif",
                             hide.ns = FALSE
  )

# IQR
s3 <- viz_thermal_boxplot(data = df_new,
                          y_var = "IQR",
                          x_var = "Group")+
  ggpubr::stat_compare_means(comparisons = my_comparisons,
                             method = "t.test",
                             label = "p.signif",
                             hide.ns = FALSE
  )
# Peak Density
s4 <- viz_thermal_boxplot(data = df_new,
                          y_var = "Peak_Density",
                          x_var = "Group")+
  ggpubr::stat_compare_means(comparisons = my_comparisons,
                             method = "t.test",
                             label = "p.signif",
                             hide.ns = FALSE
  )

# In this example, selected thermal metrics are tested against body weight and blood glucose using Spearman correlation. The results are visualized as a correlation heatmap and a representative scatter plot.
res <- correlate_thermal_traits(
  df_new,
  thermal_vars =  c("Max","Min","Mean","Median","IQR","Peak_Density"),
  external_vars = c("Weight","Glucose"),
  method = "spearman")

# Correlation heatmap summarizing the associations between thermal metrics and external traits.
s5 <- viz_cor_heatmap(res)

# Representative scatter plot illustrating the relationship between one selected thermal feature (IQR) and an external trait (Weight).
s6 <- viz_cor_scatter(df_new,
                x_col = "IQR",
                x_label = "Thermal Metric: IQR",
                y_col = "Weight"
                )

ggarrange(
  s1, s2, s3,
  s4, s5, s6,
  ncol = 3, nrow = 2,
  labels = c("A", "B", "C", "D", "E", "F"),
  widths = c(1, 1, 1),
  heights = c(1, 1.1)
)