classifiers/scoring_functions.py

import numpy as np
import torch
from classifiers.hl.hl import HalfLife
from classifiers.pa.translation_rate import TranslationRate
from classifiers.te.riboseq import RiboSeq
from src.mrna.tasks.lm.mdlm import Diffusion

class ScoringFunctions:
    def __init__(self, device, five_utr_ind, cds_length, model=None, ckpt_path=None, conf_path=None, score_func_names=None):
        """
        Class for generating score vectors given generated sequence

        Args:
            score_func_names: list of scoring function names to be evaluated
            score_weights: weights to scale scores (default: 1)
            target_protein: sequence of target protein binder
        """
        if score_func_names is None:
            self.score_func_names = []
        else:
            self.score_func_names = score_func_names

        self.device = device

        halflife = HalfLife()
        translation_rate = TranslationRate()
        ribo_seq = RiboSeq()
        self.five_utr_ind = five_utr_ind
        self.three_utr_ind = five_utr_ind + cds_length

        self.all_funcs = {'halflife': halflife, 
                          'translation_rate': translation_rate,
                          'ribo_seq': ribo_seq
                          }

        self.all_funcs = {k: v for k, v in self.all_funcs.items() if k in self.score_func_names}

        if model is not None:
            self.model = model
        else:
            self.model = Diffusion.load_from_pretrained(ckpt_path, conf_path)

        self.model.eval().to(self.device)  # Set to evaluation mode
        
    def generate_embeddings(self, seqs):
        """
        Converts a list of RNA sequences into hidden state embeddings using a pretrained DiffusionRNALanguageModel.
        
        Args:
            seqs (List[str]): List of RNA sequences (e.g., ['ACGU...', 'UGCA...']).
            model_ckpt_path (str): Path to the model checkpoint.
            
        Returns:
            np.ndarray: Array of shape (N, L, d_model), where
                        N = number of sequences,
                        L = sequence length,
                        d_model = embedding dimension.
        """
        # Load model
        five_utrs = [seq[:self.five_utr_ind] for seq in seqs]
        cds = [seq[self.five_utr_ind:self.three_utr_ind] for seq in seqs]
        three_utrs = [seq[self.three_utr_ind:] for seq in seqs] 

        # Tokenize sequences
        tokens = torch.tensor(self.model.tokenizer.batch_tokenize(five_utrs, cds, three_utrs), dtype=torch.int64, device=self.device)

        with torch.no_grad():
            with torch.autocast(device_type='cuda', dtype=torch.float16):
                hidden = self.model.net(tokens)['last_hidden_state']

        return hidden.mean(dim=1).cpu().numpy()
    
    def forward(self, input_seqs):
        scores = []
        
        input_embs = self.generate_embeddings(input_seqs)

        for i, score_func in enumerate(self.score_func_names): 
            if score_func == 'gc_content':
                score = np.array([get_gc_content(seq) for seq in input_seqs])
            elif score_func == 'sequence_complexity':
                score = np.array([idt_sequence_complexity(seq) for seq in input_seqs])
            else:
                score = self.all_funcs[score_func](input_embs = input_embs)

            scores.append(score)
            
        # convert to numpy arrays with shape (num_sequences, num_functions)
        scores = np.float32(scores).T
        
        return scores
    
    def __call__(self, input_seqs: list):
        return self.forward(input_seqs)

def idt_sequence_complexity(dna: str, window: int = 100, step: int = 50, homopolymer_threshold: int = 6) -> float:
    """
    Calculate a simple sequence complexity score inspired by IDT's approach.
    It penalizes low complexity regions with homopolymers and repeats.

    Args:
        dna (str): DNA sequence (A,T,C,G).
        window (int): Sliding window size.
        step (int): Step size for sliding.
        homopolymer_threshold (int): Minimum homopolymer length to count as low complexity.

    Returns:
        float: Complexity score (0 to 1), higher is more complex.
    """
    dna = dna.upper()
    
    def homopolymer_fraction(seq: str) -> float:
        """Fraction of sequence covered by homopolymers longer than threshold."""
        count = 0
        current_char = ''
        current_run = 0
        for base in seq:
            if base == current_char:
                current_run += 1
            else:
                if current_run >= homopolymer_threshold:
                    count += current_run
                current_char = base
                current_run = 1
        # Check last run
        if current_run >= homopolymer_threshold:
            count += current_run
        return count / len(seq)
    
    def simple_repeat_fraction(seq: str, k: int = 4) -> float:
        """Fraction of seq that consists of repeated k-mers."""
        if len(seq) < k:
            return 0.0
        kmers = [seq[i:i+k] for i in range(len(seq) - k + 1)]
        unique_kmers = set(kmers)
        repeat_count = sum(kmers.count(kmer) - 1 for kmer in unique_kmers if kmers.count(kmer) > 1)
        return repeat_count / len(seq)
    
    scores = []
    for i in range(0, len(dna) - window + 1, step):
        sub_seq = dna[i:i+window]
        hp_frac = homopolymer_fraction(sub_seq)
        rep_frac = simple_repeat_fraction(sub_seq)
        
        # Combine penalties (you can tune weights if desired)
        penalty = (hp_frac + rep_frac) / 2
        
        # Complexity = 1 - penalty
        complexity = max(0.0, 1.0 - penalty)
        scores.append(complexity)
    
    if not scores:
        # Fallback for short sequences
        hp_frac = homopolymer_fraction(dna)
        rep_frac = simple_repeat_fraction(dna)
        penalty = (hp_frac + rep_frac) / 2
        return max(0.0, 1.0 - penalty)
    
    return sum(scores) / len(scores)

def get_gc_content(dna:str) -> float:
    """
    Calculate the GC content of a DNA sequence.

    Args:
        dna (str): The DNA sequence.

    Returns:
        float: The GC content as a percentage.
    """
    if len(dna) == 0:
        return 0.0

    gc_count = sum(1 for base in dna if base in 'GCgc')
    gc_content = (gc_count / len(dna)) * 100

    return gc_content <= 0.62

# Write a function that finds the length of the longest consecutive sequence of the same character in a string
def longest_consecutive_sequence(dna: str) -> int:
    """
    Find the length of the longest consecutive sequence of the same character in a DNA sequence.

    Args:
        s (str): The input DNA sequence.

    Returns:
        int: The length of the longest consecutive sequence.
    """
    max_length = 0
    current_length = 1

    for i in range(1, len(dna)):
        if dna[i] == dna[i - 1]:
            current_length += 1
        else:
            max_length = max(max_length, current_length)
            current_length = 1

    max_length = max(max_length, current_length)

    return 1 if max_length < 6 else 0

def unittest():
    device = 'cuda'
    
    scoring = ScoringFunctions(device, score_func_names=['halflife', 
                          'translation_rate',
                          'ribosome'])
    
    seq = ['ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGGACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTGTTGCAGCCGATCATCAGCACATCTAGGTTTCGTCCGGGTGTGACCGAAAGGTAAGATGAGTAAAGGAGAAGAACTTTTCACTGGAGTTGTCCCAATTCTTGTTGAATTAGATGGCGATGTTAATGGGCAAAAATTCTCTGTCAGTGGAGAGGGTGAAGGTGATGCAACATACGGAAAACTTACCCTTAAATTTATTTGCACTACTGGGAAGCTACCTGTTCCATGGCCAACACTTGTCACTACTTTCTCTTATGGTGTTCAATGCTTTTCAAGATACCCAGATCATATGAAACAGCATGACTTTTTCAAGAGTGCCATGCCCGAAGGTTATGTACAGGAAAGAACTATATTTTACAAAGATGACGGGAACTACAAGACACGTGCTGAAGTCAAGTTTGAAGGTGATACCCTTGTTAATAGAATCGAGTTAAAAGGTATTGATTTTAAAGAAGATGGAAACATTCTTGGACACAAAATGGAATACAACTATAACTCACATAATGTATACATCATGGCAGACAAACCAAAGAATGGAATCAAAGTTAACTTCAAAATTAGACACAACATTAAAGATGGAAGCGTTCAATTAGCAGACCATTATCAACAAAATACTCCAATTGGCGATGGCCCTGTCCTTTTACCAGACAACCATTACCTGTCCACACAATCTGCCCTTTCCAAAGATCCCAACGAAAAGAGAGATCACATGATCCTTCTTGAGTTTGTAACAGCTGCTGGGATTACACATGGCATGGATGAACTATACAAATAACAATCTTTAATCAGTGTGTAACATTAGGGAGGACTTGAAAGAGCCACCACATTTTCACCGAGGCCACGCGGAGTACGATCGAGTGTACAGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAATTTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGGAGAATGACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA',
           'TACACACGAATAAAAGATAACAAAGATGAGTAAAGGAGAAGAACTTTTCACTGGAGTTGTCCCAATTCTTGTTGAATTAGATGGCGATGTTAATGGGCAAAAATTCTCTGTCAGTGGAGAGGGTGAAGGTGATGCAACATACGGAAAACTTACCCTTAAATTTATTTGCACTACTGGGAAGCTACCTGTTCCATGGCCAACACTTGTCACTACTTTCTCTTATGGTGTTCAATGCTTTTCAAGATACCCAGATCATATGAAACAGCATGACTTTTTCAAGAGTGCCATGCCCGAAGGTTATGTACAGGAAAGAACTATATTTTACAAAGATGACGGGAACTACAAGACACGTGCTGAAGTCAAGTTTGAAGGTGATACCCTTGTTAATAGAATCGAGTTAAAAGGTATTGATTTTAAAGAAGATGGAAACATTCTTGGACACAAAATGGAATACAACTATAACTCACATAATGTATACATCATGGCAGACAAACCAAAGAATGGAATCAAAGTTAACTTCAAAATTAGACACAACATTAAAGATGGAAGCGTTCAATTAGCAGACCATTATCAACAAAATACTCCAATTGGCGATGGCCCTGTCCTTTTACCAGACAACCATTACCTGTCCACACAATCTGCCCTTTCCAAAGATCCCAACGAAAAGAGAGATCACATGATCCTTCTTGAGTTTGTAACAGCTGCTGGGATTACACATGGCATGGATGAACTATACAAATAAATGTCCAGACTTCCAATTGACACTAAAGTGTCCGAACAATTACTAAATTCTCAGGGTTCCTGGTTAAATTCAGGCTGAGACTTTATTTATATATTTATAGATTCATTAAAATTTTATGAATAATTTATTGATGTTATTAATAGGGGCTATTTTCTTATTAAATAGGCTACTGGAGTGTAT',
           'CGCCTGCCTGAATCTGTTCTGCCCCCTCCCCACCCATTTCACCACCACCATGACACCGGGCACCCAGTCTCCTTTCTTCCTGCTGCTGCTCCTCACAGTGCTTACAGCTACCACAGCCCCTAAACCCGCAACAGTTGTTACGGGTTCTGGTCATGCAAGCTCTACCCCAGGTGGAGAAAAGGAGACTTCGGCTACCCAGAGAAGTTCAGTGCCCAGCTCTACTGAGAAGAATGCTTTTAATTCCTCTCTGGAAGATCCCAGCACCGACTACTACCAAGAGCTGCAGAGAGACATTTCTGAAATGTTTTTGCAGATTTATAAACAAGGGGGTTTTCTGGGCCTCTCCAATATTAAGTTCAGGCCAGGATCTGTGGTGGTACAATTGACTCTGGCCTTCCGAGAAGGTACCATCAATGTCCACGACGTGGAGACACAGTTCAATCAGTATAAAACGGAAGCAGCCTCTCGATATAACCTGACGATCTCAGACGTCAGCGTGAGTGATGTGCCATTTCCTTTCTCTGCCCAGTCTGGGGCTGGGGTGCCAGGCTGGGGCATCGCGCTGCTGGTGCTGGTCTGTGTTCTGGTTGCGCTGGCCATTGTCTATCTCATTGCCTTGGCTGTCTGTCAGTGCCGCCGAAAGAACTACGGGCAGCTGGACATCTTTCCAGCCCGGGATACCTACCATCCTATGAGCGAGTACCCCACCTACCACACCCATGGGCGCTATGTGCCCCCTAGCAGTACCGATCGTAGCCCCTATGAGAAGGTTTCTGCAGGTAATGGTGGCAGCAGCCTCTCTTACACAAACCCAGCAGTGGCAGCCACTTCTGCCAACTTGTAGGGGCACGTCGCCCGCTGAGCTGAGTGGCCAGCCAGTGCCATTCCACTCCACTCAGGTTCTTCAGGGCCAGAGCCCCTGCACCCTGTTTGGGCTGGTGAGCTGGGAGTTCAGGTGGGCTGCTCACAGCCTCCTTCAGAGGCCCCACCAATTTCTCGGACA']

    scores = scoring(input_seqs=seq)
    print(scores)
    print(len(scores))

if __name__ == '__main__':
    unittest()