<reponame>sauloal/cnidaria
#!/usr/bin/env python
"""MotifFinder that searches for over-represented k-words in an alignment."""

from __future__ import division
from cogent.motif.util import Motif, Module, ModuleFinder, ModuleConsolidator, \
    MotifFinder, Location, ModuleInstance, MotifResults
from cogent.core.bitvector import PackedBases
from cogent.maths.stats.test import combinations, multiple_comparisons
from cogent.maths.stats.distribution import poisson_high
from numpy import array, fromstring

__author__ = "<NAME>"
__copyright__ = "Copyright 2007-2012, The Cogent Project"
__credits__ = ["<NAME>", "<NAME>", "<NAME>"]
__license__ = "GPL"
__version__ = "1.5.3"
__maintainer__ = "<NAME>"
__email__ = "<EMAIL>"
__status__ = "Development"

class KWordModuleFinder(ModuleFinder):
    """ModuleFinder that finds all k-words in an alignment.
    """
    
    def __init__(self,Alignment,MolType):
        """Initializing KWordModuleFinder.
        """
        self.ModuleDict = {}
        self.ModuleOrder = []
        self.Alignment = Alignment
        self.MolType = MolType
        self._make_char_array_aln()
    
    def _make_char_array_aln(self):
        """Turns self.Alignment into a character array.
        """
        for k,v in self.Alignment.items():
            self.Alignment[k]=array(v,'c')
        
            
    def __call__(self,word_length):
        """Builds a dict of all Modules and a list of their order.
        
            - module_dict is {module pattern:Module object}
            - module_order is a list in descending order of their count.
        """
        #Dictionary keying k-word to Module
        self.ModuleDict = {}
        #For each sequence in the alignment
        for key,seq in self.Alignment.items():
            #For each position in seq till end - word_length
            for i in range(0,len(seq)-word_length+1):
                #Get the current k-word
                word = seq[i:i+word_length].tostring()
                #Create a location object
                location = Location(key,i,i+word_length)
                #Create a ModuleInstance
                curr_instance = ModuleInstance(word,location)
                #Check to see if pattern is already in dict
                if word in self.ModuleDict:
                    #Add instance to Module
                    self.ModuleDict[word][(key,i)]=curr_instance
                #Not in dict
                else:
                    #Create a new module and add to dict
                    self.ModuleDict[word]=Module({(key,i):curr_instance},\
                        MolType=self.MolType)
        #Get list of counts
        module_counts = \
            [(len(mod.Names),word) for word,mod in self.ModuleDict.items()]
        #Sort and put in descending order
        module_counts.sort()
        module_counts.reverse()
        #Get list of only the words in descending order
        self.ModuleOrder = [word for i,word in module_counts]

    
class KWordModuleConsolidatorNucleotide(ModuleConsolidator):
    """Consolidates module instances obtained from an alignment into modules.
    
        - Must be initialized with an instance of KWordModuleFinder
    """
    
    def __init__(self,KFinder): 
        """Initializing KWordModuleConsolidatorNucleotide.
        """
        self.Modules=[]
        self.KFinder=KFinder
    
    def __call__(self,mismatches):
        """Consolidates ModuleInstances in KFinder.ModuleInstances into Modules.
        
            - mismatches accounts for the difference between two ModuleInstances
                in terms of purine vs pyrimidine or strong vs weak bonding.
        """
        #New list to hold consolidated modules
        consolidated_list = []
        #Bind module_dict and module_order locally
        module_dict = self.KFinder.ModuleDict
        module_order = self.KFinder.ModuleOrder
        #Check that dict is not empty
        if module_dict:
            #Iterate through modules in order
            for pattern in module_order:
                #Create a Bitvector representation of pattern
                pat_vec = PackedBases(pattern)
                added=False
                #Iterate through consolidated_list
                for curr_module,curr_vec in consolidated_list:
                    #If pat_vec and curr_vec are in the allowed mismatch cutoff
                    if sum(pat_vec ^ curr_vec) <= mismatches:
                        #Add module information to curr_module
                        curr_module.update(module_dict[pattern])
                        added=True
                        break
                if not added:
                    #Add new module to consolidated list
                    consolidated_list.append([module_dict[pattern],pat_vec])
        #self.Modules should have only the list of modules
        self.Modules = [i[0] for i in consolidated_list]
        for mod in self.Modules:
            mod.Template = str(mod)

class KWordModuleConsolidatorProtein(ModuleConsolidator):
    """Consolidates module instances obtained from an alignment into modules.
    
        - Must be initialized with an instance of KWordModuleFinder
    """
    
    def __init__(self,KFinder): 
        """Initializing KWordModuleConsolidatorProtein.
        """
        self.Modules=[]
        self.KFinder=KFinder
    
    def __call__(self,mismatches):
        """Consolidates ModuleInstances in KFinder.ModuleInstances into Modules.
        
            - mismatches accounts for the difference between two ModuleInstances
                in terms of purine vs pyrimidine or strong vs weak bonding.
        """
        #New list to hold consolidated modules
        consolidated_list = []
        #Bind module_dict and module_order locally
        module_dict = self.KFinder.ModuleDict
        module_order = self.KFinder.ModuleOrder
        #Check that dict is not empty
        if module_dict:
            #Iterate through modules in order
            for pattern in module_order:
                #Create a Bitvector representation of pattern
                pat_array = fromstring(pattern,'c') 
                added=False
                #Iterate through consolidated_list
                for curr_module,curr_array in consolidated_list:
                    #If pat_vec and curr_vec are in the allowed mismatch cutoff
                    if sum(pat_array != curr_array) <= mismatches:
                        #Add module information to curr_module
                        curr_module.update(module_dict[pattern])
                        added=True
                        break
                if not added:
                    #Add new module to consolidated list
                    consolidated_list.append([module_dict[pattern],pat_array])
        #self.Modules should have only the list of modules
        self.Modules = [i[0] for i in consolidated_list]

class KWordModuleFilterProtein(ModuleConsolidator):
    """Filters list of Modules by number of Modules and seqs where Module is in.
    
        - This is a strict ModuleConsolidator.  Does not allow mismatches, but
        rather consolidates based on a minimum allowed modules and minimum
        number of sequences that a module must be present in.
    """
    
    def __init__(self,KFinder,Alignment): 
        """Initializing KWordModuleFilterProtein.
        """
        self.Modules=[]
        self.KFinder=KFinder
        self.Alignment=Alignment
    
    def __call__(self,min_modules,min_seqs):
        """Adds modules to self.Modules from self.KFinder.ModuleDict as follows:
        
            - modules that have at least min_modules and are present in at least
            min_seqs will be added to self.Modules.
        """
        #First filter by min_modules
        for k,v in self.KFinder.ModuleDict.items():
            #If number of modules is greater than the minimum
            if len(v) >= min_modules:
                #Check to see that module is present in at least min_seqs
                curr_seqs = {}
                for curr in v.keys():
                    curr_seqs[curr[0]]=True
                if len(curr_seqs) >= min_seqs:
                    self.Modules.append(self.fixModuleSequence(v))
    
    def fixModuleSequence(self,module):
        """Remaps original (non-reduced) sequence string for each ModuleInstance
        """
        module_len=len(str(module))
        module.Template=str(module)
        for k,v in module.items():
            seq_id, module_start = k
            module_end = module_start+module_len
            loc = Location(seq_id,module_start,module_end)
            curr_str = \
                self.Alignment[seq_id][module_start:module_end]
            curr_instance = ModuleInstance(curr_str,loc)
            module[k]=curr_instance
        return module


class KWordMotifFinder(MotifFinder):
    """Constructs a list of Motifs from a list of Modules.
    
        - Must be initialized with a list of Modules and an Alignment.
    """
    
    def __init__(self,Modules,Alignment,Mismatches,BaseFrequency):
        """Initializing KWordMotifFinder.
        """
        self.Modules=Modules        #List of modules from Alignment
        self.Alignment=Alignment    #Alignment used to find modules
        self.Mismatches=Mismatches  #Number of mismatches used for consolidation
        self.BaseFrequency=BaseFrequency
        

    def __call__(self,threshold,alphabet_len=None,use_loose=True):
        """Returns a MotifResults object containing a list of significant Motifs
        
            - When a Module is found that is less than the threshold, a new 
                Motif object is created and added to the list of significant
                motifs.
        """
        #Create new MotifResults object
        k_word_results = MotifResults()
        k_word_results.Modules=[]
        k_word_results.Motifs=[]
        #Add alignment to MotifResults
        k_word_results.Alignment = self.Alignment
        #Give each module a unique ID
        module_id = 0
        #For each module        
        for i,module in enumerate(self.Modules):
            
            if use_loose:
                p_curr = self.getPValueLoose(module,alphabet_len)
            else:
                p_curr = self.getPValueStrict(module,alphabet_len)
            #If the P value is less than or equal to the threshold
            if p_curr <= threshold:
                #Add the P value to the current module
                module.Pvalue = p_curr
                module.ID = module_id
                module_id+=1
                #Append the module to the k_word_results Modules list
                k_word_results.Modules.append(module)
                #Create a Motif from the module and append to Motifs list
                k_word_results.Motifs.append(Motif(module))
        return k_word_results
        
    
    def getPValueStrict(self,module, alphabet_len=None):
        """Returns the Pvalue of a module.
        
            - pass alphabet_len if alphabet other than module.MolType was used
            to find modules (i.e. using a protein reduced alphabet)
        """
        #Length of the module
        module_len = len(str(module))
        
        #if moltype length has not been passed, get it from module.MolType
        #
        if not alphabet_len:
            alphabet_len = len(module.MolType.Alphabet)
        
        #Total length of the alignment
        aln_len=sum(map(len,self.Alignment.values()))
        #Number of sequences in the alignment
        num_seqs = len(self.Alignment)
        #get module_p
        module_p = 1
        for char in module.Template:
            module_p *= self.BaseFrequency.get(char,1)
        
        #Mean for passing to poisson_high NOT using loose correction
        strict_mean = \
            (aln_len + num_seqs*(1-module_len))*(module_p)
            
        #Strict P value from poisson_high
        strict_p_value = poisson_high(len(module)-1,strict_mean)
        #Correct P value for multiple comparisons
        strict_p_corrected = \
            multiple_comparisons(strict_p_value,alphabet_len**module_len)

        return strict_p_corrected

    
    def getPValueLoose(self,module, alphabet_len=None):
        """Returns the Pvalue of a module.
        
            - pass alphabet_len if alphabet other than module.MolType was used
            to find modules (i.e. using a protein reduced alphabet)
        """
        #Length of the module
        module_len = len(module.Template)
        
        #if moltype length has not been passed, get it from module.MolType
        #
        if not alphabet_len:
            alphabet_len = len(module.MolType.Alphabet)
        
        #Total length of the alignment
        aln_len=sum(map(len,self.Alignment.values()))
        #Number of sequences in the alignment
        num_seqs = len(self.Alignment)
        

        #get module_p
        module_p = 1
        for char in module.Template:
            module_p *= self._degen_p(char,module.MolType)
        
        #Mean for passing to poisson_high using loose correction
        
        loose_mean = \
            (aln_len + num_seqs*(1-module_len))*(module_p)
        #Loose P value from poisson_high
        loose_p_value = poisson_high(len(module)-1,loose_mean)
        #Correct P value for multiple comparisons
        loose_p_corrected = \
            multiple_comparisons(loose_p_value,alphabet_len**module_len)
       
        return loose_p_corrected
    
    def _degen_p(self,char,alphabet):
        """Returns the sum of the probabilities of seeing each degenerate char.
        """
        all = alphabet.Degenerates.get(char,char)
        return sum([self.BaseFrequency.get(x,1) for x in all])
            
        
