from plastid import BAMGenomeArray, GTF2_TranscriptAssembler
import pandas as pd
import numpy as np

def calculate_te(riboseq_bam, rnaseq_bam, gtf_path):
    '''Calculate translation efficiency per gene'''
    # Load transcripts
    transcripts = list(GTF2_TranscriptAssembler(gtf_path))

    # Load alignments
    ribo = BAMGenomeArray(riboseq_bam)
    rna = BAMGenomeArray(rnaseq_bam)

    results = []
    for tx in transcripts:
        if tx.cds_start is None:
            continue

        # Get CDS region
        cds = tx.get_cds()

        # Count reads
        ribo_counts = ribo.count_in_region(cds)
        rna_counts = rna.count_in_region(tx)  # Full transcript for RNA-seq

        # Normalize by length
        cds_length = sum(len(seg) for seg in cds)
        tx_length = tx.length

        ribo_rpk = ribo_counts / (cds_length / 1000)
        rna_rpk = rna_counts / (tx_length / 1000)

        if rna_rpk > 0:
            te = ribo_rpk / rna_rpk
        else:
            te = np.nan

        results.append({
            'gene': tx.get_gene(),
            'transcript': tx.get_name(),
            'ribo_counts': ribo_counts,
            'rna_counts': rna_counts,
            'te': te
        })

    return pd.DataFrame(results)

library(riborex)

# Load count matrices
# Rows = genes, columns = samples
ribo_counts <- read.csv('ribo_counts.csv', row.names = 1)
rna_counts <- read.csv('rna_counts.csv', row.names = 1)

# Sample information
sample_info <- data.frame(
    sample = colnames(ribo_counts),
    condition = factor(c('control', 'control', 'treated', 'treated'))
)

# Run riborex
results <- riborex(
    rnaCntTable = rna_counts,
    riboCntTable = ribo_counts,
    rnaCond = sample_info$condition,
    riboCond = sample_info$condition
)

# Significant differential TE
sig_te <- results[results$padj < 0.05, ]

library(DESeq2)

# Combine Ribo-seq and RNA-seq counts
counts <- cbind(ribo_counts, rna_counts)

# Design matrix with interaction term
coldata <- data.frame(
    condition = factor(rep(c('ctrl', 'ctrl', 'treat', 'treat'), 2)),
    assay = factor(rep(c('ribo', 'rna'), each = 4)),
    row.names = colnames(counts)
)

dds <- DESeqDataSetFromMatrix(
    countData = counts,
    colData = coldata,
    design = ~ condition + assay + condition:assay
)

dds <- DESeq(dds)

# The interaction term tests for differential TE
res_te <- results(dds, name = 'conditiontreat.assayribo')

def normalize_counts(counts_df, method='tpm'):
    '''Normalize count matrix'''
    if method == 'tpm':
        # TPM normalization
        rpk = counts_df.div(counts_df['length'] / 1000, axis=0)
        scale = rpk.sum(axis=0) / 1e6
        tpm = rpk.div(scale, axis=1)
        return tpm

    elif method == 'rpkm':
        # RPKM normalization
        total = counts_df.sum(axis=0)
        rpm = counts_df / total * 1e6
        rpkm = rpm.div(counts_df['length'] / 1000, axis=0)
        return rpkm

def calculate_te_matrix(ribo_tpm, rna_tpm):
    '''Calculate TE from normalized matrices'''
    # Add pseudocount to avoid division by zero
    te = (ribo_tpm + 0.1) / (rna_tpm + 0.1)
    return np.log2(te)  # Log2 TE