## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.show = "hold",
  fig.pos = "h",
  out.width = "42%",
  fig.retina = 2
)

## ----setup, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE---------------
# Loading all packages required
# Data manipulation and reshaping
library(dplyr) # For data filtering, grouping, and summarising.
library(tidyr) # For reshaping data (e.g., pivot_longer, pivot_wider).

# Plotting and visualization
library(ggplot2) # For creating all the ggplot-based visualizations.
library(gridExtra) # For arranging multiple plots on a single page.
library(ggrepel) # For improved label placement in plots (e.g., volcano plots).
library(pheatmap) # For heatmap.2, which is used to generate heatmaps.
library(plot3D) # For creating 3D scatter plots in PCA and sPLS-DA analyses.
library(reshape2) # For data transformation (e.g., melt) in cross-validation plots.

# Statistical analysis
library(mixOmics) # For multivariate analyses (PCA, sPLS-DA, etc.).
library(e1071) # For computing skewness and kurtosis.
library(pROC) # For ROC curve generation in machine learning model evaluation.

# Machine learning
library(xgboost) # For building XGBoost classification models.
library(randomForest) # For building Random Forest classification models.
library(caret) # For cross-validation and other machine learning utilities.

# Package development and document rendering
library(knitr) # For knitting RMarkdown files and setting chunk options.
library(devtools) # For installing the development version of the package from GitHub.

# Load in the CytoProfile package
library(CytoProfile)

# Loading in data
data("ExampleData1")
data_df <- ExampleData1

## ----EDA1,  echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE---------------
# Generating boxplots to check for outliers for raw values
cyt_bp(data_df[, -c(1:3)], pdf_title = NULL)
# Removing the first 3 columns to retain only continuous variables.

# Generating boxplots to check for outliers for log2 values
cyt_bp(data_df[, -c(1:3)], pdf_title = NULL, scale = "log2")


## ----EDA2, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE----------------
# Raw values for group-specific boxplots
cyt_bp2(data_df[, -c(3, 5:28)], pdf_title = NULL, scale = NULL)

# Log2-transformed group-specific boxplots
cyt_bp2(data_df[, -c(3, 5:28)], pdf_title = NULL, scale = "log2")



## ----EDA3, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE----------------
# Histogram for overall raw data
cyt_skku(data_df[, -c(1:3)], pdf_title = NULL, group_cols = NULL)

# Histogram with grouping (e.g., "Group")
cyt_skku(ExampleData1[, -c(2:3)], pdf_title = NULL, group_cols = c("Group"))


## ----EDA4,  echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE---------------
# Generating basic error bar plots
data_df <- ExampleData1
cyt_errbp(
  data_df[, c("Group", "CCL.20.MIP.3A", "IL.10")],
  group_col = "Group",
  p_lab = FALSE,
  es_lab = FALSE,
  class_symbol = FALSE,
  x_lab = "Cytokines",
  y_lab = "Concentrations in log2 scale",
  log2 = TRUE
)

## ----EDA5,  echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE---------------
# Generating Error Bar Plot enriched with p-value and effect size
data_df <- ExampleData1
cyt_errbp(
  data_df[, c("Group", "CCL.20.MIP.3A", "IL.10")],
  group_col = "Group",
  p_lab = TRUE,
  es_lab = TRUE,
  class_symbol = TRUE,
  x_lab = "Cytokines",
  y_lab = "Concentrations in log2 scale",
  log2 = TRUE
)


## ----Univariate1,  echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE--------
# Performing Test
data_df <- ExampleData1[, -c(3)]
data_df <- dplyr::filter(data_df, Group != "ND", Treatment != "Unstimulated")
# Test example
cyt_ttest(
  data_df[, c(1:2, 5:6)],
  scale = "log2",
  verbose = TRUE,
  format_output = TRUE
)

## ----Univariate2,  echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE--------
# Perform ANOVA comparisons test (example with 2 cytokines)
data_df <- ExampleData1[, -c(3)]
cyt_anova(data_df[, c(1:2, 5:6)], format_output = TRUE)

## ----Multivariate1, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE, fig.show = "hold", out.width = "25%"----
# cyt_plsda function.
data <- ExampleData1[, -c(3)]
data_df <- dplyr::filter(data, Group != "ND" & Treatment == "CD3/CD28")
cyt_splsda(
  data_df,
  pdf_title = NULL,
  colors = c("black", "purple"),
  bg = FALSE,
  scale = "log2",
  ellipse = TRUE,
  conf_mat = FALSE,
  var_num = 25,
  cv_opt = "loocv",
  comp_num = 2,
  pch_values = c(16, 4),
  group_col = "Group",
  group_col2 = "Treatment",
  roc = TRUE
)

## ----Multivariate2, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE, fig.show = "hold", out.width = "25%", fig.ncol = 2, fig.width=7, fig.height=6----
# cyt_mint_plsda function.
data_df <- ExampleData5[, -c(2, 4)]
data_df <- dplyr::filter(data_df, Group != "ND")

cyt_mint_splsda(
  data_df,
  group_col = "Group",
  batch_col = "Batch",
  colors = c("black", "purple"),
  ellipse = TRUE,
  var_num = 25,
  comp_num = 2,
  scale = "log2",
  verbose = FALSE
)

## ----Multivariate3, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE,  fig.show = "hold", out.width = "25%"----
data <- ExampleData1[, -c(3, 23)]
data_df <- dplyr::filter(data, Group != "ND" & Treatment != "Unstimulated")
cyt_pca(
  data_df,
  pdf_title = NULL,
  colors = c("black", "red2"),
  scale = "log2",
  comp_num = 2,
  pch_values = c(16, 4),
  group_col = "Group"
)

## ----EDA6, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE----------------
# Generating Volcano Plot
data_df <- ExampleData1[, -c(2:3)]
cyt_volc(
  data_df,
  group_col = "Group",
  cond1 = "T2D",
  cond2 = "ND",
  fold_change_thresh = 2.0,
  top_labels = 15
)

## ----EDA7, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE, fig.width=5, fig.height=4----
# Generating Heat map
cyt_heatmap(
  data = data_df,
  scale = "log2", # Optional scaling
  annotation_col = "Group",
  title = NULL
)

## ----EDA8, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE----------------
# Generating dual flashlights plot
data_df <- ExampleData1[, -c(2:3)]
dfp <- cyt_dualflashplot(
  data_df,
  group_var = "Group",
  group1 = "T2D",
  group2 = "ND",
  ssmd_thresh = -0.2,
  log2fc_thresh = 1,
  top_labels = 10
)
# Print the plot
dfp
# Print the table data used for plotting
print(dfp$data, n = 25)

## ----ML1, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE-----------------
# Using XGBoost for classification
data_df0 <- ExampleData1
data_df <- data.frame(data_df0[, 1:3], log2(data_df0[, -c(1:3)]))
data_df <- data_df[, -c(2:3)]
data_df <- dplyr::filter(data_df, Group != "ND")

cyt_xgb(
  data = data_df,
  group_col = "Group",
  nrounds = 500,
  max_depth = 4,
  min_split_loss = 0,
  learning_rate = 0.05,
  nfold = 5,
  cv = TRUE,
  objective = "multi:softprob",
  eval_metric = "auc",
  early_stopping_rounds = NULL,
  top_n_features = 10,
  verbose = 0,
  plot_roc = TRUE,
  print_results = FALSE
)

## ----ML2, echo=TRUE, warning=FALSE, message=FALSE, cache=TRUE-----------------
# Using Random Forest for classification
cyt_rf(
  data = data_df,
  group_col = "Group",
  k_folds = 5,
  ntree = 1000,
  mtry = 4,
  run_rfcv = TRUE,
  plot_roc = TRUE,
  verbose = FALSE
)

