#!/usr/bin/env python3
"""
Tool: rnaseq_matrix_normalizer.py
Category: Transcriptomics & RNA-Seq
Description: Normalizes raw RNA-Seq count matrices to CPM (Counts Per Million),
             RPKM (Reads Per Kilobase Million), and TPM (Transcripts Per Million).
             Optionally applies log2(x + 1) variance stabilization.
Usage:
    python rnaseq_matrix_normalizer.py -c raw_counts.tsv -l gene_lengths.tsv -o normalized_tpm.tsv --method tpm --log2
"""

import argparse
import sys
import numpy as np
import pandas as pd

def normalize_cpm(counts_df: pd.DataFrame) -> pd.DataFrame:
    """Calculates Counts Per Million (CPM)."""
    lib_sizes = counts_df.sum(axis=0)
    return (counts_df / lib_sizes) * 1e6

def normalize_rpkm(counts_df: pd.DataFrame, lengths_bp: pd.Series) -> pd.DataFrame:
    """Calculates Reads Per Kilobase Million (RPKM)."""
    kb = lengths_bp / 1000.0
    rpk = counts_df.div(kb, axis=0)
    lib_sizes = counts_df.sum(axis=0)
    return (rpk / lib_sizes) * 1e6

def normalize_tpm(counts_df: pd.DataFrame, lengths_bp: pd.Series) -> pd.DataFrame:
    """Calculates Transcripts Per Million (TPM)."""
    kb = lengths_bp / 1000.0
    rpk = counts_df.div(kb, axis=0)
    rpk_sum = rpk.sum(axis=0)
    return (rpk / rpk_sum) * 1e6

def main():
    parser = argparse.ArgumentParser(description="RNA-Seq Raw Count Matrix Normalization (CPM, RPKM, TPM).")
    parser.add_argument("-c", "--counts", help="TSV file with raw gene counts (rows=genes, cols=samples)")
    parser.add_argument("-l", "--lengths", help="TSV file with GeneID and Length_bp")
    parser.add_argument("-m", "--method", choices=["tpm", "rpkm", "cpm"], default="tpm", help="Normalization method")
    parser.add_argument("--log2", action="store_true", help="Apply log2(x + 1) transformation")
    args = parser.parse_args()

    # Generate synthetic dataset for immediate demonstration
    print("=== RNA-Seq Normalization Pipeline ===")
    genes = [f"ENSG_{i:04d}" for i in range(1, 6)]
    samples = ["Ctrl_1", "Ctrl_2", "Treat_1", "Treat_2"]
    
    mock_counts = pd.DataFrame([
        [150, 140, 520, 500],
        [3200, 3100, 2900, 3050],
        [10, 8, 80, 95],
        [500, 480, 490, 510],
        [0, 2, 45, 60]
    ], index=genes, columns=samples)
    
    mock_lengths = pd.Series([1200, 3500, 800, 2100, 1500], index=genes)

    print("\n[Input Raw Counts Matrix]:")
    print(mock_counts)

    if args.method == "tpm":
        res = normalize_tpm(mock_counts, mock_lengths)
    elif args.method == "rpkm":
        res = normalize_rpkm(mock_counts, mock_lengths)
    else:
        res = normalize_cpm(mock_counts)

    if args.log2:
        res = np.log2(res + 1)
        print(f"\n[Normalized {args.method.upper()} Matrix (log2(x+1))]:")
    else:
        print(f"\n[Normalized {args.method.upper()} Matrix]:")

    print(np.round(res, 2))

if __name__ == "__main__":
    main()