#!/usr/bin/env Rscript
# ==============================================================================
# Chapter 8: R for Next Gen Biologist
# Script: r08_statistical_hypothesis_testing.R
# Description: Statistical Hypothesis Testing for Biologists in R:
#              Parametric vs non-parametric tests (t-test, Wilcoxon), ANOVA,
#              Fisher's exact test for variant enrichment, and FDR multiple testing.
# Usage: Rscript r08_statistical_hypothesis_testing.R
# ==============================================================================

cat("\n=======================================================\n")
cat(" OMICS.CO.IN - R FOR NEXT GEN BIOLOGIST: CHAPTER 8\n")
cat(" Biostatistics: Hypothesis Testing & FDR Correction\n")
cat("=======================================================\n\n")

set.seed(42)

# 1. Two-Sample Testing: Student's t-test vs Wilcoxon Rank-Sum Test
cat("[Step 1] Comparing Gene Expression Between Control & Treated Cohorts...\n")
# Expression values (TPM)
control_group <- rnorm(20, mean = 25.4, sd = 4.2)
treated_group <- rnorm(20, mean = 34.8, sd = 5.1)

# Check normality using Shapiro-Wilk test
shapiro_ctrl <- shapiro.test(control_group)
shapiro_treat <- shapiro.test(treated_group)
cat(sprintf("Normality Test (Shapiro-Wilk) -> Control p: %.3f | Treated p: %.3f\n",
            shapiro_ctrl$p.value, shapiro_treat$p.value))

# Parametric: Welch's Two-Sample t-test
t_res <- t.test(treated_group, control_group)
cat(sprintf("Welch's Two-Sample t-test: t = %.2f | p-value = %.4e (Difference: %+.2f)\n",
            t_res$statistic, t_res$p.value, diff(t_res$estimate)))

# Non-Parametric: Wilcoxon Rank-Sum / Mann-Whitney U test
w_res <- wilcox.test(treated_group, control_group)
cat(sprintf("Wilcoxon Rank-Sum test   : W = %.1f | p-value = %.4e\n",
            w_res$statistic, w_res$p.value))

# 2. Multi-Group Testing: One-Way ANOVA across Tumor Stages
cat("\n[Step 2] Multi-Group Testing: ANOVA & Post-hoc Tukey across Stages...\n")
n_per_stage <- 15
stage_data <- data.frame(
  Expression = c(rnorm(n_per_stage, 20, 3),  # Stage I
                 rnorm(n_per_stage, 28, 4),  # Stage II
                 rnorm(n_per_stage, 39, 5),  # Stage III
                 rnorm(n_per_stage, 48, 6)), # Stage IV
  Stage = factor(rep(c("Stage_I", "Stage_II", "Stage_III", "Stage_IV"), each = n_per_stage))
)

aov_fit <- aov(Expression ~ Stage, data = stage_data)
aov_summary <- summary(aov_fit)
print(aov_summary)

cat("\nPost-hoc Tukey HSD Pairwise Comparisons:\n")
tukey_res <- TukeyHSD(aov_fit)
print(head(tukey_res$Stage, 3))

# 3. Categorical Contingency: Fisher's Exact Test for Mutation Enrichment
cat("\n[Step 3] Fisher's Exact Test for Somatic Mutation Enrichment...\n")
# Contingency table: TP53 Mutation vs Drug Response
#                Responders  Non-Responders
# TP53 Mutated       18            4
# TP53 Wild-Type      6           22
contingency_mat <- matrix(c(18, 6, 4, 22), nrow = 2,
                          dimnames = list(
                            TP53_Status = c("Mutant", "WildType"),
                            Drug_Response = c("Responder", "NonResponder")
                          ))
print(contingency_mat)

fisher_res <- fisher.test(contingency_mat)
cat(sprintf("Fisher's Exact Test p-value: %.4e\n", fisher_res$p.value))
cat(sprintf("Odds Ratio                 : %.2f (95%% CI: %.2f - %.2f)\n",
            fisher_res$estimate, fisher_res$conf.int[1], fisher_res$conf.int[2]))

# 4. Multiple Testing Correction: FDR (Benjamini-Hochberg)
cat("\n[Step 4] Multiple Testing Correction across 10,000 Transcripts...\n")
n_hypotheses <- 10000
raw_pvals <- c(rbeta(500, 0.1, 10), runif(9500)) # 500 true positives + 9500 null uniform

bonferroni_p <- p.adjust(raw_pvals, method = "bonferroni")
bh_fdr_p     <- p.adjust(raw_pvals, method = "BH")

cat(sprintf("Raw p < 0.05 (Uncorrected False Positives): %d genes\n", sum(raw_pvals < 0.05)))
cat(sprintf("Bonferroni Corrected (alpha = 0.05)       : %d genes (Extremely Conservative)\n", sum(bonferroni_p < 0.05)))
cat(sprintf("Benjamini-Hochberg FDR (q < 0.05)         : %d genes (Standard for RNA-Seq)\n", sum(bh_fdr_p < 0.05)))

# 5. Gene Co-expression Correlation (Pearson vs Spearman)
cat("\n[Step 5] Gene Co-expression Correlation Testing...\n")
geneA <- rnorm(30, 50, 10)
geneB <- geneA * 0.8 + rnorm(30, 0, 5) # Co-regulated partner

cor_pearson  <- cor.test(geneA, geneB, method = "pearson")
cor_spearman <- cor.test(geneA, geneB, method = "spearman")

cat(sprintf("Pearson Correlation (r) : %+.3f | p-value: %.3e\n",
            cor_pearson$estimate, cor_pearson$p.value))
cat(sprintf("Spearman Rank Corr (rho): %+.3f | p-value: %.3e\n",
            cor_spearman$estimate, cor_spearman$p.value))

cat("\n[SUCCESS] Chapter 8 statistical testing script executed successfully.\n")