%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
coherence([100, 2, 1, 100, 2, 1])
(1.0, 2)
coherence([3, 2, 1, 3, 2, 1])
(1.0, 1)