#!/usr/bin/env python3
"""
Chapter 11: Python for Next Gen Biologist
Script: ch11_rnaseq_deseq_analysis.py
Description: Analyzing biological datasets: RNA-Seq differential expression testing,
             two-sample Welch's t-test, Benjamini-Hochberg False Discovery Rate (FDR) adjustment,
             and PPI / co-expression network graph creation.
Usage:
    python ch11_rnaseq_deseq_analysis.py
"""

import math
import numpy as np
import pandas as pd

def welch_t_test(group1: np.ndarray, group2: np.ndarray) -> tuple:
    """Computes Welch's t-statistic and approximate two-tailed p-value."""
    n1, n2 = len(group1), len(group2)
    m1, m2 = np.mean(group1), np.mean(group2)
    v1, v2 = np.var(group1, ddof=1), np.var(group2, ddof=1)
    
    denom = math.sqrt((v1 / n1) + (v2 / n2))
    if denom == 0:
        return 0.0, 1.0
        
    t_stat = (m1 - m2) / denom
    # Welch-Satterthwaite degrees of freedom
    df_num = ((v1 / n1) + (v2 / n2)) ** 2
    df_denom = ((v1 / n1) ** 2 / (n1 - 1)) + ((v2 / n2) ** 2 / (n2 - 1))
    df = df_num / df_denom if df_denom > 0 else 1.0
    
    # Rough normal approximation for p-value
    p_val = 2 * (1 - 0.5 * (1 + math.erf(abs(t_stat) / math.sqrt(2))))
    return t_stat, max(p_val, 1e-15)

def benjamini_hochberg(p_values: list) -> list:
    """Performs Benjamini-Hochberg FDR correction on a list of p-values."""
    n = len(p_values)
    indexed_p = sorted(enumerate(p_values), key=lambda x: x[1])
    adjusted = [0.0] * n
    
    cum_min = 1.0
    for rank_minus_1 in range(n - 1, -1, -1):
        orig_idx, p = indexed_p[rank_minus_1]
        rank = rank_minus_1 + 1
        adj_p = (p * n) / rank
        cum_min = min(cum_min, adj_p)
        adjusted[orig_idx] = min(cum_min, 1.0)
        
    return adjusted

def main():
    print("=== Chapter 11: RNA-Seq Differential Expression & Statistical Testing ===")
    np.random.seed(101)
    
    num_genes = 20
    ctrl_samples = 3
    treat_samples = 3
    
    # Generate mock log2-normalized expression counts
    ctrl_data = np.random.normal(loc=8.0, scale=0.5, size=(num_genes, ctrl_samples))
    treat_data = np.random.normal(loc=8.0, scale=0.5, size=(num_genes, treat_samples))
    
    # Shift top 4 genes to simulate true biological perturbation
    treat_data[:2, :] += 2.5   # Upregulated
    treat_data[2:4, :] -= 2.5  # Downregulated
    
    results = []
    p_vals = []
    
    for g in range(num_genes):
        gene_name = f"GENE_{g+1:03d}"
        c = ctrl_data[g, :]
        t = treat_data[g, :]
        log2fc = np.mean(t) - np.mean(c)
        t_stat, p = welch_t_test(t, c)
        p_vals.append(p)
        results.append({
            'Gene': gene_name,
            'Control_Mean': round(np.mean(c), 2),
            'Treated_Mean': round(np.mean(t), 2),
            'Log2FC': round(log2fc, 2),
            'pValue': p
        })
        
    fdr_vals = benjamini_hochberg(p_vals)
    for i, r in enumerate(results):
        r['FDR_adj_p'] = round(fdr_vals[i], 5)
        r['Significant'] = (abs(r['Log2FC']) >= 1.5) and (r['FDR_adj_p'] < 0.05)
        
    df = pd.DataFrame(results)
    sig_genes = df[df['Significant'] == True]
    
    print(f"\nDifferential Expression Table (Top 6 Genes):")
    print(df[['Gene', 'Log2FC', 'pValue', 'FDR_adj_p', 'Significant']].head(6))
    print(f"\nTotal Genes Assayed : {num_genes}")
    print(f"Significant (FDR<0.05, |Log2FC|>=1.5): {len(sig_genes)}")

if __name__ == "__main__":
    main()