params <-
list(report_date = structure(20584, class = "Date"))

## ----simulate-data, message=FALSE, warning=FALSE------------------------------
library(XYomics)
library(dplyr)

# Load precomputed bulk dataset
expression_data <- readRDS(
  system.file("extdata", "bulk_expression.rds", package = "XYomics")
)

sex <- readRDS(
  system.file("extdata", "bulk_sex.rds", package = "XYomics")
)

phenotype <- readRDS(
  system.file("extdata", "bulk_phenotype.rds", package = "XYomics")
)

## ----stratified-analysis, eval=FALSE------------------------------------------
# res <- sex_stratified_analysis_bulk(expression_data, sex, phenotype)
# # The `min_samples` parameter defines the minimum number of cells required
# # to perform the interaction analysis. The default value is 3.
# 
# 

## ----load-stratified-results--------------------------------------------------
res <- readRDS(
  system.file("extdata", "bulk_results_degs.rds", package = "XYomics")
)
# Internal QC results from XYomics can be accessed through:
# res$validation
#
# This includes validation of sex and phenotype balance, group sample sizes,
# detection of design imbalances, and identification of groups below the
# minimum sample threshold required for reliable downstream analysis.
res$validation

## ----sex-specific-analysis, eval=FALSE----------------------------------------
# res_cat <- categorize_sex(res$male_DEGs, res$female_DEGs)

## ----load-categorized-results-------------------------------------------------
res_cat <- readRDS(
  system.file("extdata", "bulk_results_cat.rds", package = "XYomics")
)

cat("Top sex-stratified differentially expressed genes:\n")
head(res_cat)

table(res_cat$DEG_Type)

## ----interaction-term-analysis-bulk, eval=FALSE-------------------------------
# interaction_term_results_bulk <- sex_interaction_analysis_bulk(expression_data, phenotype, sex)
# # The `min_samples` parameter defines the minimum number of cells required
# # to perform the interaction analysis. The default value ofr interaction analysis is 20.
# 

## ----load-interaction-results-------------------------------------------------
interaction_term_results_bulk <- readRDS(
  system.file("extdata", "bulk_interaction_results.rds", package = "XYomics")
)

head(interaction_term_results_bulk$summary_stats)# signficant sex-modulated DEGs: interaction_term_results_bulk$sig_results

## ----visualization------------------------------------------------------------
top_male_gene <- res_cat %>%
  filter(DEG_Type == "male-specific") %>%
  arrange(Male_FDR) %>%
  pull(Gene_Symbols) %>%
  head(1)

top_female_gene <- res_cat %>%
  filter(DEG_Type == "female-specific") %>%
  arrange(Female_FDR) %>%
  pull(Gene_Symbols) %>%
  head(1)
top_dimorphic_gene <- res_cat %>%
  filter(DEG_Type == "sex-dimorphic") %>%
  arrange(Female_FDR) %>%
  pull(Gene_Symbols) %>%
  head(1)
top_neutral_gene <- res_cat %>%
  filter(DEG_Type == "sex-neutral") %>%
  arrange(Female_FDR) %>%
  pull(Gene_Symbols) %>%
  head(1)

top_genes = c(top_male_gene, top_female_gene, top_dimorphic_gene, top_neutral_gene )

generate_violinplot_bulk(expression_data,sex,phenotype , top_genes)

## ----pathway-analysis, eval=FALSE---------------------------------------------
# pathway_results <- categorized_enrich(res_cat, enrichment_db = "GO")
# 
# # The `enrichment_db` parameter can be set to "GO", "KEGG", "REACTOME".
# # If the user wants to use a custom database, a TERM2GENE data frame must be
# # provided via the `custom_db` parameter.
# #
# # Gene identifiers can be either any Gene identifier type supported by OrgDb
# #' (e.g., "SYMBOL", "ENTREZID", "ENSEMBL", "UNIPROT")
# # using the `gene_type` parameter.
# #
# # The `return_df` parameter controls the output format:
# # if set to TRUE, results are returned as data frames;
# # by default (FALSE), enrichResult objects are returned

## ----load-pathway-results-----------------------------------------------------
pathway_results <- readRDS(
  system.file("extdata", "bulk_pathway_results.rds", package = "XYomics")
)

head(pathway_results)

## ----celltype1_plot , results='asis', fig.keep='all', fig.width=15, fig.height=7, message=FALSE----
# visualize via dotplot
plot <- plot_enrichment_dotplots(pathway_results)
print(plot)

## ----get-string-network, message=FALSE----------------------------------------
# Fetch STRING network (can be replaced with a custom network)
# g <- get_string_network(organism = "9606", score_threshold = 900)

# Load a pre-existing network from a file
g <- readRDS(system.file("extdata", "string_example_network.rds", package = "XYomics")) # You can also download the hormonal network "hormonal_network_metacore.rds" instead of "string_example_network.rds".

## ----construct-pcsf-network, eval=FALSE---------------------------------------
# # Use dimorphic DE results to define prizes
# dimorphic_specific <- res_cat[res_cat$DEG_Type == "sex-dimorphic", ]
# dimorphic_prizes <- -log10((dimorphic_specific$Male_FDR + dimorphic_specific$Female_FDR) / 2)
# names(dimorphic_prizes) <-dimorphic_specific$Gene_Symbols
# 
# # Construct the PCSF subnetwork
# 
# dimorphic_network <- construct_ppi_pcsf(g = g, prizes = dimorphic_prizes)
# network <- ppi_pipeline(res_cat, g)

## ----load-network-results-----------------------------------------------------
network <- readRDS(
  system.file("extdata", "bulk_network_results.rds", package = "XYomics")
)

neutral_network <- readRDS(
  system.file("extdata", "bulk_neutral_network.rds", package = "XYomics")
)

neutral_specific <- res_cat[res_cat$DEG_Type == "sex-neutral", ]

## ----visualize-network--------------------------------------------------------


plot_network(neutral_network, "sex-neutral", neutral_specific, show_barplot = T)


#Generate plots for all categories

all_plots <- plot_network_pipeline(network, res_cat)


## ----report, eval=FALSE-------------------------------------------------------
# # This command generates a comprehensive HTML report
# 
# template_path = system.file("extdata", "Template_report_bulk.Rmd", package = "XYomics") # fetch the template for bulk from the package
# 
# generate_cat_report(res_cat, pathway_results, network, template_path = template_path)
# 
#