In [3]:

    

%pylab inline
from scipy import signal

def coherence(original_values):
    """Calculate coherence with an idea ribo-seq signal,
    this function is equivalent to phasescore implemented
    above, but we expect scipy has better handling of numerical
    precision.
    Parameters
    ----------
    values : array like
             List of value
    Returns
    -------
    coh : float
          Periodicity score calculated as
          coherence between input and idea 1-0-0 signal
    valid: int
           number of valid codons used for calculation
    """
    coh, valid = 0.0, -1
    for frame in [0, 1, 2]:
        values = original_values[frame:]
        normalized_values = []
        i = 0
        while i+2 < len(values):
            if values[i] == values[i+1] == values[i+2] == 0:
                i += 3
            else:
                real = values[i] + values[i+1] * cos(2*pi/3) + values[i+2] * cos(4*pi/3)
                image = values[i+1] * sin(2*pi/3)+values[i+2] * sin(4*pi/3)
                norm = sqrt(real**2 + image**2)
                if norm == 0:
                    norm = 1
                normalized_values += [values[i]/norm,
                                      values[i+1]/norm,
                                      values[i+2]/norm
                                      ]
                i += 3

        length = len(normalized_values)//3*3
        if length == 0:
            coh, valid = (0.0, 0)
        else:
            normalized_values = normalized_values[:length]
            uniform_signal = [1, 0, 0] * (len(normalized_values)//3)
            f, Cxy = signal.coherence(normalized_values,
                                      uniform_signal,
                                      window=[1.0, 1.0, 1.0],
                                      nperseg=3,
                                      noverlap=0)
            periodicity_score = Cxy[np.argwhere(np.isclose(f, 1/3.0))[0]][0]
            if periodicity_score > coh:
                coh = periodicity_score
                valid = length//3
            if valid == -1:
                valid = length//3
    return np.sqrt(coh), valid


    

    




Populating the interactive namespace from numpy and matplotlib

In [5]:

    

coherence([100, 2, 1, 100, 2, 1])


    

    
Out[5]:





(1.0, 2)

In [6]:

    

coherence([3, 2, 1, 3, 2, 1])


    

    
Out[6]:





(1.0, 1)

In [ ]: