## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(warning = FALSE, message = FALSE, fig.width = 8, fig.height = 6) 

## ----install, eval=FALSE------------------------------------------------------
# # Install from CRAN
# install.packages("ZetaSuite")
# 
# # Load the package
# library(ZetaSuite)

## ----shiny_app, eval=FALSE----------------------------------------------------
# # Launch the Shiny app
# ZetaSuiteApp()
# 
# # Launch without opening browser automatically
# ZetaSuiteApp(launch.browser = FALSE)
# 
# # Launch on a specific port
# ZetaSuiteApp(port = 3838)

## ----load_data----------------------------------------------------------------
library(ZetaSuite)

# Load example data
data(countMat)
data(negGene)
data(posGene)
data(nonExpGene)
data(ZseqList)
data(SVMcurve)

# Display data dimensions
cat("Count matrix dimensions:", dim(countMat), "\n")
cat("Negative controls:", nrow(negGene), "genes\n")
cat("Positive controls:", nrow(posGene), "genes\n")
cat("Non-expressed genes:", nrow(nonExpGene), "genes\n")

## ----qc_analysis--------------------------------------------------------------
# Perform quality control analysis
qc_results <- QC(countMat, negGene, posGene)

# Display QC plots
cat("QC analysis completed. Generated", length(qc_results), "diagnostic plots.\n")

## ----qc_plot1, fig.height = 5, fig.width = 8----------------------------------
qc_results$score_qc

## ----qc_plot2, fig.height = 5, fig.width = 6----------------------------------
qc_results$tSNE_QC

## ----qc_plot3, fig.height = 4, fig.width = 8----------------------------------
grid::grid.draw(qc_results$QC_box)

## ----qc_plot4, fig.height = 5, fig.width = 6----------------------------------
qc_results$QC_SSMD

## ----zscore_analysis----------------------------------------------------------
# Calculate Z-scores
zscore_matrix <- Zscore(countMat, negGene)

# Display first few rows and columns
cat("Z-score matrix dimensions:", dim(zscore_matrix), "\n")
cat("First 5 rows and columns:\n")
print(zscore_matrix[1:5, 1:5])

## ----event_coverage-----------------------------------------------------------
# Calculate event coverage
ec_results <- EventCoverage(zscore_matrix, negGene, posGene, binNum = 100, combine = TRUE)

# Display event coverage plots
cat("Event coverage analysis completed.\n")

## ----ec_plots, fig.height = 5, fig.width = 10---------------------------------
# Decrease direction (exon skipping)
ec_results[[2]]$EC_jitter_D

## ----ec_plots2, fig.height = 5, fig.width = 10--------------------------------
# Increase direction (exon inclusion)
ec_results[[2]]$EC_jitter_I

## ----zeta_calculation---------------------------------------------------------
# Calculate zeta scores without SVM correction
zeta_scores <- Zeta(zscore_matrix, ZseqList, SVM = FALSE)

# Display summary statistics
cat("Zeta score summary:\n")
cat("Number of genes:", nrow(zeta_scores), "\n")
cat("Zeta_D range:", range(zeta_scores$Zeta_D), "\n")
cat("Zeta_I range:", range(zeta_scores$Zeta_I), "\n")

# Show top hits
cat("\nTop 10 genes by Zeta_D (decrease direction):\n")
top_decrease <- head(zeta_scores[order(zeta_scores$Zeta_D, decreasing = TRUE), ], 10)
print(top_decrease)

## ----svm_analysis, eval=FALSE-------------------------------------------------
# # Run SVM analysis (can be computationally intensive)
# svm_results <- SVM(ec_results)
# 
# # Calculate zeta scores with SVM correction
# zeta_scores_svm <- Zeta(zscore_matrix, ZseqList, SVMcurve = svm_results, SVM = TRUE)

## ----screen_strength----------------------------------------------------------
# Calculate FDR cutoffs and Screen Strength
fdr_results <- FDRcutoff(zeta_scores, negGene, posGene, nonExpGene, combine = TRUE)

# Display Screen Strength plots
cat("Screen Strength analysis completed.\n")

## ----ss_plot1, fig.height = 5, fig.width = 8----------------------------------
fdr_results[[2]]$Zeta_type

## ----ss_plot2, fig.height = 5, fig.width = 6----------------------------------
fdr_results[[2]]$SS_cutOff

## ----fdr_table----------------------------------------------------------------
# Display FDR cutoff results
fdr_table <- fdr_results[[1]]
cat("FDR cutoff results summary:\n")
cat("Number of thresholds tested:", nrow(fdr_table), "\n")
cat("Screen Strength range:", range(fdr_table$SS), "\n")

# Show optimal thresholds (SS > 0.8)
optimal_thresholds <- fdr_table[fdr_table$SS > 0.8, ]
cat("\nOptimal thresholds (SS > 0.8):\n")
print(head(optimal_thresholds, 10))

## ----hit_selection------------------------------------------------------------
# Example: Select threshold with SS > 0.8 and reasonable number of hits
selected_threshold <- fdr_table[fdr_table$SS > 0.8 & fdr_table$TotalHits > 50, ]
if(nrow(selected_threshold) > 0) {
  best_threshold <- selected_threshold[which.max(selected_threshold$SS), ]
  cat("Recommended threshold:", best_threshold$Cut_Off, "\n")
  cat("Screen Strength:", best_threshold$SS, "\n")
  cat("Total hits:", best_threshold$TotalHits, "\n")
  
  # Identify hits
  combined_zeta <- zeta_scores$Zeta_D + zeta_scores$Zeta_I
  hits <- names(combined_zeta[combined_zeta >= best_threshold$Cut_Off])
  cat("Number of hits identified:", length(hits), "\n")
}

## ----single_cell_example, eval=FALSE------------------------------------------
# # Example single cell analysis (requires single cell count matrix)
# # single_cell_results <- ZetaSuitSC(count_matrix_sc, binNum = 10, filter = TRUE)

## ----custom_params, eval=FALSE------------------------------------------------
# # Custom event coverage analysis
# ec_custom <- EventCoverage(zscore_matrix, negGene, posGene,
#                           binNum = 200,    # More bins for finer resolution
#                           combine = FALSE) # Separate decrease/increase directions
# 
# # Custom zeta score calculation with SVM
# zeta_custom <- Zeta(zscore_matrix, ZseqList,
#                    SVMcurve = svm_results,
#                    SVM = TRUE)  # Use SVM correction
# 
# # Custom FDR analysis
# fdr_custom <- FDRcutoff(zeta_scores, negGene, posGene, nonExpGene,
#                        combine = FALSE)  # Analyze directions separately

## ----batch_processing, eval=FALSE---------------------------------------------
# # Example: Process multiple datasets
# datasets <- list(dataset1 = list(countMat = countMat1, negGene = negGene1),
#                 dataset2 = list(countMat = countMat2, negGene = negGene2))
# 
# results <- lapply(datasets, function(ds) {
#   zscore <- Zscore(ds$countMat, ds$negGene)
#   ec <- EventCoverage(zscore, ds$negGene, ds$posGene, binNum = 100)
#   zeta <- Zeta(zscore, ZseqList, SVM = FALSE)
#   return(list(zscore = zscore, ec = ec, zeta = zeta))
# })

## ----sessionInfo--------------------------------------------------------------
sessionInfo()