diff --git a/Fig2_S1_S2/Circadian_amplitude_period_exp.py b/Fig2_S1_S2/Circadian_amplitude_period_exp.py deleted file mode 100644 index b3461a2..0000000 --- a/Fig2_S1_S2/Circadian_amplitude_period_exp.py +++ /dev/null @@ -1,164 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 12:54:22 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer -from pyboat import ensemble_measures as em - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -path = '...Raw_data/' -output = '...Fig2/' -dose = ['untreated', '5uM', '10uM'] -density = ['high'] -channel1 = 'circadian' -channel2 = 'cell_cycle' -colors = ['tab:blue', 'tab:orange', 'tab:green'] - -dt = 0.5 -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -#%%Changing dose -for i in density: - plt.figure(figsize = (12,10)) - periods_all = pd.DataFrame() - periods_rd = pd.DataFrame() - count = 0 - for j in dose: - - condition = str(i)+'_density_'+str(j)+'_' - print(condition) - - #Circadian data - data1 = pd.read_csv(path+condition+channel1+r'_filtered.csv', index_col=0) - #Cell-cycle data - data2 = pd.read_csv(path+condition+channel2+r'_filtered.csv', index_col=0) - - time = data1.index.values*0.5 - - ridge_results1 = {} - ridge_results2 = {} - - for column in data1: - norm_signal = data1[column].dropna() - signal = data2[column].dropna() - - if len(norm_signal)>100: - wAn.compute_spectrum(norm_signal,do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd.set_index(norm_signal.index,inplace=True) - ridge_results1[column] = rd - - wAn.compute_spectrum(signal,do_plot=False) - rd2 = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd2.set_index(signal.index,inplace=True) - ridge_results2[column] = rd2 - - powers_series = em.average_power_distribution(ridge_results1.values(), signal_ids = ridge_results1.keys()) - high_power_ids = powers_series[powers_series > 10].index - - osc = 0 - for col in high_power_ids: - osc += data1[col].fillna(0) - - signal = wAn.sinc_smooth(osc/len(high_power_ids), T_c=10) - plt.plot(time, signal, linewidth=5, label=condition) - - wAn.compute_spectrum(signal, do_plot=False) - rd = wAn.get_maxRidge(power_thresh = 10, smoothing_wsize=50) - rd_all = wAn.get_maxRidge(power_thresh = 0, smoothing_wsize=50) - periods_rd[count] = rd['periods'] - periods_all[count] = rd_all['periods'] - count+=1 - - plt.xlabel('Time [hours]') - plt.xticks([0,24,48,72,96,120]) - plt.ylabel('Mean circadian signal [a.u.]') - plt.legend(loc='best') - plt.savefig(output+'Mean_circ_osc_'+str(i)+'.svg') - plt.show() - - plt.figure(figsize=(10,8)) - for j in range(len(dose)): - plt.plot(periods_rd[j].index.values*0.5, periods_rd[j], linewidth=5, color=colors[j], label=str(dose[j])) - # plt.plot(time, periods_all[j], '--', linewidth=5, color=colors[j], alpha=0.5, label=str(dose[j])) - plt.xlabel('Time [hours]') - plt.xticks([0,24,48,72])#,96,120]) - plt.ylabel('Period') - plt.legend(loc='best') - plt.savefig(output+'Period_mean_circ_osc_'+str(i)+'.svg') - plt.show() - -#%% Changing density -density = ['high', 'medium', 'low'] - -fig = plt.figure(figsize = (12,10)) -count = 0 - -periods_rd = pd.DataFrame() - -for jj in density: - condition = str(jj)+'_density_untreated_' - print(condition) - - #Circadian data - data1 = pd.read_csv(path+condition+channel1+r'_filtered.csv', index_col=0) - - time=data1.index.values*0.5 - - ridge_results1 = {} - df1 = pd.DataFrame() - - for column in data1: - norm_signal = data1[column].dropna() - - if len(norm_signal)>100: - wAn.compute_spectrum(norm_signal,do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd.set_index(norm_signal.index,inplace=True) - ridge_results1[column] = rd - rd['traceId'] = column - - powers_series = em.average_power_distribution(ridge_results1.values(), signal_ids = ridge_results1.keys()) - high_power_ids = powers_series[powers_series > 10].index - - osc = 0 - for col in high_power_ids: - osc += data1[col].fillna(0) - - signal = wAn.sinc_smooth(osc/len(high_power_ids), T_c=10) - plt.plot(time, signal, linewidth=5, label=str(jj)+' density') - - wAn.compute_spectrum(signal, do_plot=False) - rd = wAn.get_maxRidge(power_thresh = 10, smoothing_wsize=50) - rd_all = wAn.get_maxRidge(power_thresh = 0, smoothing_wsize=50) - periods_rd[count] = rd['periods'] - count+=1 - -plt.xlabel('Time [hours]') -plt.xticks([0,24,48,72,96,120]) -plt.ylabel('Mean circadian signal [a.u.]') -plt.legend(loc='best') -plt.savefig(output+'Mean_circ_osc_densities.svg') -plt.show() - -plt.figure(figsize=(10,8)) -for j in range(len(density)): - plt.plot(periods_rd[j].index.values*0.5, periods_rd[j], linewidth=5, color=colors[j], label=str(density[j])+' density') -plt.xlabel('Time [hours]') -plt.xticks([0,24,48,72])#,96,120]) -plt.ylabel('Period') -plt.legend(loc='best') -plt.savefig(output+'Period_mean_circ_osc_densities.svg') -plt.show() diff --git a/Fig2_S1_S2/Circadian_cellcycle_entrainment_regions.py b/Fig2_S1_S2/Circadian_cellcycle_entrainment_regions.py deleted file mode 100644 index f002287..0000000 --- a/Fig2_S1_S2/Circadian_cellcycle_entrainment_regions.py +++ /dev/null @@ -1,117 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 12:49:54 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -from statistics import NormalDist -import pandas as pd -from mpl_toolkits.axes_grid1 import make_axes_locatable - -plt.rcParams.update({'font.size': 20}) -plt.rcParams['svg.fonttype'] = 'none' - -output = '.../' - -coupling = np.arange(0, 0.42, step=0.02) -detuning = np.arange(0, 8.1, step=0.1) - -arn_tong = np.zeros((len(coupling), len(detuning))) -det_xticks = [] -coupl_yticks = [] - -#Coupled oscillators parameters -N = 200 -dt = 0.1 -tf = 200 - -gg = 1 -ll = 1 -a0 = 1 -acc = 1 - -kappa = 0.1 #extracellular coupling - -for mm in range(len(coupling)): - for nn in range(len(detuning)): - delta = detuning[nn] - mu1, sigma1 = 20+delta, np.sqrt(4) - circ_per = np.random.normal(mu1, sigma1, N) - mu2, sigma2 = 28-delta, np.sqrt(4) - cell_per = np.random.normal(mu2, sigma2, N) - - overlap = NormalDist(mu=mu1, sigma=sigma1).overlap(NormalDist(mu=mu2, sigma=sigma2)) - print('overlap', overlap*100) - - eps = coupling[mm] - print('coupling, detuning:', eps, mu2-mu1) - - t = np.arange(0, tf, step=dt) - X = np.zeros((N, len(t))) - Y = np.zeros((N, len(t))) - XX = np.zeros((N, len(t))) - YY = np.zeros((N, len(t))) - - X[:,0] = np.random.uniform(-1, 1, size=(N)) - Y[:,0] = np.random.uniform(-1, 1, size=(N)) - XX[:,0] = np.random.uniform(-1, 1, size=(N)) - YY[:,0] = np.random.uniform(-1, 1, size=(N)) - - for i in range(len(t)-1): - sum_x = 0 - sum_y = 0 - for m in range(N): - middle = -gg*(np.sqrt(X[m,i]**2+Y[m,i]**2)-a0) - X[m,i+1] = X[m,i]+dt*(middle*X[m,i]-2*np.pi*Y[m,i]/circ_per[m]+sum_x*kappa/(2*N)) - Y[m,i+1] = Y[m,i]+dt*(middle*Y[m,i]+2*np.pi*X[m,i]/circ_per[m]+sum_y*kappa/(2*N)) - - mid = -ll*(np.sqrt(XX[m,i]**2+YY[m,i]**2)-acc) - XX[m,i+1] = XX[m,i]+dt*(mid*XX[m,i]-2*np.pi*YY[m,i]/cell_per[m]+eps*(X[m,i]+XX[m,i])/2) - YY[m,i+1] = YY[m,i]+dt*(mid*YY[m,i]+2*np.pi*XX[m,i]/cell_per[m]+eps*(Y[m,i]+YY[m,i])/2) - for k in range(N): - sum_x = sum_x+X[k,i] #global coupling - sum_y = sum_y+Y[k,i] #global coupling - - ph_diff = pd.DataFrame(columns=np.arange(0, N, step=1)) - std_all = [] - count_sync = 0 - for m in range(N): - difference = [] - for n in range(len(t)): - diff = np.arctan2(Y[m,n], X[m,n])-np.arctan2(YY[m,n], XX[m,n]) - phasediff = np.arctan2(np.sin(diff), np.cos(diff)) - difference.append(phasediff) - ph_diff[m] = difference - std_ph = np.std(difference[int(len(t)/2):len(t)]) - std_all.append(std_ph) - if std_ph < 0.5: - count_sync+=1 - arn_tong[mm,nn] = count_sync - - print('synch',count_sync) - print('---') - -arn_tong = np.flipud(arn_tong) -np.savetxt('Arn_tongue_intr_N200_k005.txt', arn_tong) - -file_path = '/Users/nicagutu/Nextcloud/Manuscripts/Circadian-CellCycle/Code/Poincare_model/Arn_tongue_intr_N200_k005.txt' -arn_tong = np.loadtxt(file_path) - -fig = plt.figure(figsize=(12, 10)) -color_map = plt.imshow(arn_tong, interpolation='bicubic', aspect='auto') -color_map.set_cmap('jet') -plt.yticks([0, 10, 20], [0.4, 0.2, 0]) -plt.xticks([0, 10, 20, 30, 40, 50, 60, 70, 80], [-8, -6, -4, -2, 0, 2, 4, 6, 8]) -plt.xlabel('Detuning') -plt.ylabel('Intracellular coupling') -ax = plt.gca() -divider = make_axes_locatable(ax) -cax = divider.append_axes("right", size="5%", pad=0.1) -cb = plt.colorbar(cax=cax) -cb.set_label('# locked oscillators') -# plt.savefig(output+'Arn_tongue_intra'+str(k)+'.svg') -plt.show() diff --git a/Fig2_S1_S2/Circadian_entrainment_regions.py b/Fig2_S1_S2/Circadian_entrainment_regions.py deleted file mode 100644 index 1c202c6..0000000 --- a/Fig2_S1_S2/Circadian_entrainment_regions.py +++ /dev/null @@ -1,107 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 12:30:00 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from mpl_toolkits.axes_grid1 import make_axes_locatable -from pyboat import WAnalyzer - -plt.rcParams.update({'font.size': 26}) -plt.rcParams['svg.fonttype'] = 'none' - -output = '.../' -path2 = '.../' - -dt = 0.1 -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -coupling = np.arange(0, 0.0505, step=0.0005) -dispersion = np.arange(0, 10.1, step=0.1) -arn_tong = np.zeros((len(coupling),len(dispersion))) - -#Simulation parameters -N = 20 -dt = 0.1 -tf = 150 - -#Coupled oscillators parameters -gg = 1 -ll = 1 -a0 = 1 -acc = 1 - -eps = 0 #Intracellular coupling - -for mm in range(len(coupling)): - for nn in range(len(dispersion)): - delta = dispersion[nn] - print('dispersion', delta) - - mu1, sigma1 = 24, np.sqrt(delta) - circ_per=np.random.normal(mu1, sigma1, N) - - kappa = coupling[mm] - print(kappa) - - t = np.arange(0, tf, step=dt) - X = np.zeros((N, len(t))) - Y = np.zeros((N, len(t))) - - X[:,0] = np.random.uniform(-1, 1, size=(N)) - Y[:,0] = np.random.uniform(-1, 1, size=(N)) - - for i in range(len(t)-1): - sum_x = 0 - sum_y = 0 - for k in range(N): - sum_x = sum_x+X[k,i] #global coupling - sum_y = sum_y+Y[k,i] #global coupling - - for m in range(N): - middle = -gg*(np.sqrt(X[m,i]**2+Y[m,i]**2)-a0) - X[m,i+1] = X[m,i]+dt*(middle*X[m,i]-2*np.pi*Y[m,i]/circ_per[m]+sum_x*kappa/(2*N)) - Y[m,i+1] = Y[m,i]+dt*(middle*Y[m,i]+2*np.pi*X[m,i]/circ_per[m]+sum_y*kappa/(2*N)) - - -#%%Synchronization index - mean_cell_sqr_time = np.mean((np.mean(X, axis=0))**2) - mean_cell_time_sqr = (np.mean(np.mean(X, axis=0)))**2 - avg = 0 - for i in range(N): - avg = avg+np.mean((X[i,:])**2)-(np.mean(X[i,:]))**2 - R = (mean_cell_sqr_time-mean_cell_time_sqr)/(avg/N) - arn_tong[mm,nn] = R #np.mean(R_all[int(len(t)*0.75):len(t)]) #to save - - -arn_tong = np.flipud(arn_tong) - -np.savetxt(path2+'Arn_tongue_extracell_coupling.txt', arn_tong) - -file_path = path2+'Arn_tongue_extracell_N20_k005.txt' -arn_tong = np.loadtxt(file_path) -arn_tong_dupl = np.hstack((np.fliplr(arn_tong), arn_tong)) - -fig = plt.figure(figsize=(14, 10)) -plt.xticks(np.linspace(0, 200, 9), [10, 7.5, 5, 2.5, 0, 2.5, 5, 7.5, 10]) -color_map = plt.imshow(arn_tong_dupl, interpolation='bicubic', aspect='auto') -color_map.set_cmap('jet') -plt.xlabel('Dispersion circadian periods [hours]') -plt.yticks(np.linspace(0,100,11), np.round(np.linspace(0.05, 0, 11), 5)) -plt.locator_params(axis='y', nbins=20) -plt.ylabel('Extracellular coupling') -ax = plt.gca() -divider = make_axes_locatable(ax) -cax = divider.append_axes("right", size="5%", pad=0.1) -cb = plt.colorbar(cax=cax) -cb.set_label('Synchronization index R') -# plt.savefig(output+'Extra_coupled_circadian_arntongue.svg') -plt.show() diff --git a/Fig2_S1_S2/Circadian_osc_cells.py b/Fig2_S1_S2/Circadian_osc_cells.py deleted file mode 100644 index 7301773..0000000 --- a/Fig2_S1_S2/Circadian_osc_cells.py +++ /dev/null @@ -1,120 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 13:00:17 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer -from scipy.signal import find_peaks - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -path = '.../Raw_data/' -output = '...' -dose = ['untreated', '5uM', '10uM'] -density = ['high'] -channel1 = 'circadian' -channel2 = 'cell_cycle' -colors = ['gold', 'tab:green', 'tab:blue'] - -dt = 0.5 -lowT = 16 -highT = 40 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -#%%Changing dose within same density - -num_peaks_cond = [] -for i in density: - - for j in dose: - count = 0 - - condition = str(i)+'_density_'+str(j)+'_' - - #Circadian data - data1 = pd.read_csv(path+condition+channel1+r'_filtered.csv', index_col=0) - - time = data1.index.values*0.5 - - for column in data1: - norm_signal = data1[column].dropna() - - if len(norm_signal) > 0: - wAn.compute_spectrum(norm_signal,do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd.set_index(norm_signal.index,inplace=True) - mean_power = (np.mean(rd['power'])) - - peaks = find_peaks(norm_signal, height=0, distance=32, prominence=0.1) - - if len(peaks[0]) > 1 and mean_power > 10: - count += 1 - - - num_peaks_cond.append((count/len(data1.columns))*100) - - - print(num_peaks_cond) - - fig = plt.figure(figsize=(10,10)) - barlist = plt.bar(dose, num_peaks_cond) - for ii in range(len(dose)): - barlist[ii].set_color(colors[ii]) - plt.ylabel('% oscillating cells') - plt.ylim([0,100]) - plt.savefig(output+'Number_osc_cells_high_density.svg') - plt.show() - - - -#%% Changing density -density = ['high', 'medium', 'low'] - -num_peaks_cond = [] -for jj in density: - count = 0 - - condition = str(jj)+'_density_'+str(dose[0])+'_' - print(condition) - - #Circadian data - data1 = pd.read_csv(path+condition+channel1+r'_filtered.csv', index_col=0) - - - for column in data1: - norm_signal = data1[column].dropna() - - if len(norm_signal) > 0: - wAn.compute_spectrum(norm_signal,do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd.set_index(norm_signal.index,inplace=True) - mean_power = (np.mean(rd['power'])) - - peaks = find_peaks(norm_signal, height=0, distance=32, prominence=0.1) - - if len(peaks[0]) > 1 and mean_power > 10: - count += 1 - - - num_peaks_cond.append((count/len(data1.columns))*100) - - -print(num_peaks_cond) - -fig = plt.figure(figsize=(10,10)) -barlist = plt.bar(density, num_peaks_cond) -for ii in range(len(density)): - barlist[ii].set_color(colors[ii]) -plt.ylabel('% oscillating cells') -plt.ylim([0,100]) -plt.savefig(output+'Number_osc_cells_condition_changing_density.svg') -plt.show() - diff --git a/Fig2_S1_S2/Circadian_properties_Poincare.py b/Fig2_S1_S2/Circadian_properties_Poincare.py deleted file mode 100644 index ea37b76..0000000 --- a/Fig2_S1_S2/Circadian_properties_Poincare.py +++ /dev/null @@ -1,160 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 12:36:30 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -output = '.../Fig2/' - -dt = 0.1 -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -#Parameters -Nosc = 200 -dt = 0.1 -tf = 120 - -gg = 1 -ll = 1 -a0 = 1 -acc = 1 - -kappa_vect = [0, 0.05, 0.075, 0.5] -eps1 = 0 -eps2 = 0.01 - -periods = pd.DataFrame() -phase_coh = pd.DataFrame() -mean_osc = pd.DataFrame() - -plt.figure(figsize=(12,10)) - -for kk in range(len(kappa_vect)): - kappa = kappa_vect[kk] - - #Vectors initialization - t = np.arange(0,tf,step=dt) - X = np.zeros((Nosc,len(t))) - XX = np.zeros((Nosc,len(t))) - Y = np.zeros((Nosc,len(t))) - YY = np.zeros((Nosc,len(t))) - - #Initial position - X[:,0] = np.random.uniform(-0.3,1,size=(Nosc)) - XX[:,0] = np.random.uniform(-1,1,size=(Nosc)) - Y[:,0] = np.random.uniform(-0.5,1,size=(Nosc)) - YY[:,0] = np.random.uniform(-1,1,size=(Nosc)) - - #Period distribution - mu1, sigma1 = 22, np.sqrt(6) - period_circ = np.random.normal(mu1, sigma1, Nosc) - - pert_resset = 0 - - # Integration over time and population - for i in range(len(t) - 1): - sum_x = 0 - sum_y = 0 - for k in range(Nosc): - sum_x = sum_x+X[k,i] - sum_y = sum_y+Y[k,i] - - - if int((tf*0.0)/dt)100: - wAn.compute_spectrum(norm_signal,do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd.set_index(norm_signal.index,inplace=True) - ridge_results1[column] = rd - rd['traceId'] = column - df1 = df1.append(rd) - - for ii in rd['periods'].index.values: - if 16 < rd['periods'][ii] < 50: - periods_mean.append(rd['periods'][ii]) - - wAn.compute_spectrum(signal,do_plot=False) - rd2 = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd2.set_index(signal.index,inplace=True) - ridge_results2[column] = rd2 - rd2['traceId'] = column - df2 = data2.append(rd2) - - periods_cond[j] = periods_mean - - powers_series = em.average_power_distribution(ridge_results1.values(), signal_ids = ridge_results1.keys()) - high_power_ids = powers_series[powers_series > 10].index - high_power_ridge_results = [ridge_results1[i] for i in high_power_ids] - res = em.get_ensemble_dynamics(high_power_ridge_results) - # pl.ensemble_dynamics(*res) - # plt.show() - - print('Sample size: ', len(ridge_results1)) - - print('Maximum ph coh: ', np.where(res[a][str(value)]==max(res[a][str(value)]))[0]) - print('Minimmum ph coh: ', np.where(res[a][str(value)]==min(res[a][str(value)]))[0]) - - plt.plot(time, res[a][str(value)], color=colors[count], label=str(j), linewidth=5) - # plt.fill_between(time, res[a]['Q1'], res[a]['Q3'], color=colors[count], alpha=.1) - count+=1 - plt.xlabel('Time [hours]') - plt.ylabel(str(magn)) - plt.legend(loc='best') - # plt.savefig(output+'ED_'+str(a)+'_'+str(i)+'0_10uM.svg') - plt.show() - -periods_cond_df = pd.DataFrame().from_dict(periods_cond, orient='index') -periods_cond_df = periods_cond_df.T - -print(periods_cond_df) - -plt.figure(figsize=(12,10)) -sns.boxplot(data=periods_cond_df, showfliers = False) -plt.ylabel('Periods [hours]') -plt.xlabel('Inhibitor dose') -# plt.savefig(output+'Boxplot_instantaneous_periods_high_density.svg') -plt.show() - - -#%% Changing density -density = ['high', 'medium', 'low'] - -fig = plt.figure(figsize = (12,10)) -count = 0 - -periods_cond = {} - -for jj in density: - condition = str(jj)+'_density_'+str(dose[0])+'_' - print(condition) - - #Circadian data - data1 = pd.read_csv(path+condition+channel1+r'_filtered.csv', index_col=0) - #Cell-cycle data - data2 = pd.read_csv(path+condition+channel2+r'_filtered.csv', index_col=0) - - time=data1.index.values*0.5 - - ridge_results1 = {} - df1 = pd.DataFrame() - ridge_results2 = {} - df2 = pd.DataFrame() - - periods_mean = [] - - for column in data1: - norm_signal = data1[column].dropna() - signal = data2[column].dropna() - - if len(norm_signal)>100: - wAn.compute_spectrum(norm_signal,do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd.set_index(norm_signal.index,inplace=True) - ridge_results1[column] = rd - rd['traceId'] = column - df1 = df1.append(rd) - - for ii in rd['periods'].index.values: - if 16 < rd['periods'][ii] < 50: - periods_mean.append(rd['periods'][ii]) - - wAn.compute_spectrum(signal,do_plot=False) - rd2 = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - rd2.set_index(signal.index,inplace=True) - ridge_results2[column] = rd2 - rd2['traceId'] = column - df2 = data2.append(rd2) - - periods_cond[jj] = periods_mean - - powers_series = em.average_power_distribution(ridge_results1.values(), signal_ids = ridge_results1.keys()) - high_power_ids = powers_series[powers_series > 10].index - high_power_ridge_results = [ridge_results1[i] for i in high_power_ids] - res = em.get_ensemble_dynamics(high_power_ridge_results) - # pl.ensemble_dynamics(*res) - # plt.show() - - print('Sample size: ', len(ridge_results1)) - - print('Maximum ph coh: ', np.where(res[a][str(value)]==max(res[a][str(value)]))[0]) - print('Minimmum ph coh: ', np.where(res[a][str(value)]==min(res[a][str(value)]))[0]) - - plt.plot(time, res[a][str(value)], color=colors[count], label=str(jj), linewidth=5) - # plt.fill_between(time, res[a]['Q1'], res[a]['Q3'], color=colors[count], alpha=.1) - count+=1 - -plt.xlabel('Time [hours]') -plt.ylabel(str(magn)) -plt.legend(loc='best') -# plt.savefig(output+'ED_'+str(a)+'_changing_density.svg') -plt.show() - -periods_cond_df = pd.DataFrame().from_dict(periods_cond, orient='index') -periods_cond_df = periods_cond_df.T - -print(periods_cond_df) - -colors = {'high': 'tab:orange', 'medium': 'tab:green', 'low': 'tab:blue'} - -plt.figure(figsize=(12,10)) -sns.boxplot(data=periods_cond_df, showfliers = False) -plt.ylabel('Periods [hours]') -plt.xlabel('Initial density') -# plt.savefig(output+'Boxplot_instantaneous_periods_changing_density.svg') -plt.show() diff --git a/Fig2_S1_S2/Heatmap_extra_vs_intra_Poincare.py b/Fig2_S1_S2/Heatmap_extra_vs_intra_Poincare.py deleted file mode 100644 index 0f67072..0000000 --- a/Fig2_S1_S2/Heatmap_extra_vs_intra_Poincare.py +++ /dev/null @@ -1,150 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 13:06:00 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from matplotlib.colors import LinearSegmentedColormap - -def create_custom_colormap(): #input: maximum number divisions - all_colors = [ - # (220, 20, 60), # Crimson - (199, 21, 133), # Medium Violet Red - (218, 112, 214), # Orchid - (240, 128, 128), # Light Coral - (221, 160, 221), # Plum - (219, 112, 147), # Pale Violet Red - (186, 85, 211), # Medium Orchid - (153, 50, 204), # Dark Orchid - (102, 51, 153), # Rebecca Purple - (147, 112, 219), # Medium Purple - (123, 104, 238), # Medium Slate Blue - (72, 61, 139) # Dark Slate Blue - ] - - normalized_colors = [(r / 255, g / 255, b / 255) for r, g, b in all_colors] - cmap_name = 'custom_colormap' - - return LinearSegmentedColormap.from_list(cmap_name, normalized_colors, N=256) - -plt.rcParams.update({'font.size': 20}) -plt.rcParams['svg.fonttype'] = 'none' - -output = '.../' - -extra = np.arange(0, 1.05, step=0.05) -intra = np.arange(0, 0.45, step=0.05) -arn_tong = np.zeros((len(extra), len(intra))) -det_xticks = [] -coupl_yticks = [] - -#Coupled oscillators -gg = 1 -ll = 1 -a0 = 1 -acc = 1 -N = 100 -dt = 0.2 -tf = 200 - -eps1 = 0 - -for mm in range(len(extra)): - coupl_yticks.append("%.3f"%extra[mm]) - for nn in range(len(intra)): - mu1, sigma1=20, np.sqrt(4) - circ_per = np.random.normal(mu1, sigma1, N) - mu2, sigma2 = 28, np.sqrt(4) - cell_per = np.random.normal(mu2, sigma2, N) - det_xticks.append("%.3f"%intra[nn]) - - kappa = extra[mm] - eps2 = intra[nn] - print(kappa, eps2) - - t = np.arange(0, tf, step=dt) - X = np.zeros((N, len(t))) - Y = np.zeros((N, len(t))) - XX = np.zeros((N, len(t))) - YY = np.zeros((N, len(t))) - - X[:,0] = np.random.uniform(-1, 1, size=(N)) - Y[:,0] = np.random.uniform(-1, 1, size=(N)) - XX[:,0] = np.random.uniform(-1, 1, size=(N)) - YY[:,0] = np.random.uniform(-1, 1, size=(N)) - - for i in range(len(t)-1): - sum_x = 0 - sum_y = 0 - for k in range(N): - sum_x = sum_x+X[k,i] #global coupling - sum_y = sum_y+Y[k,i] #global coupling - - for m in range(N): - middle = -gg*(np.sqrt(X[m,i]**2+Y[m,i]**2)-a0) - X[m,i+1] = X[m,i]+dt*(middle*X[m,i]-2*np.pi*Y[m,i]/circ_per[m]+sum_x*kappa/(2*N)+eps1*(XX[m,i])) - Y[m,i+1] = Y[m,i]+dt*(middle*Y[m,i]+2*np.pi*X[m,i]/circ_per[m]+sum_y*kappa/(2*N)+eps1*(YY[m,i])) - - mid = -ll*(np.sqrt(XX[m,i]**2+YY[m,i]**2)-acc) - XX[m,i+1] = XX[m,i]+dt*(mid*XX[m,i]-2*np.pi*YY[m,i]/cell_per[m]+eps2*(X[m,i])) - YY[m,i+1] = YY[m,i]+dt*(mid*YY[m,i]+2*np.pi*XX[m,i]/cell_per[m]+eps2*(Y[m,i])) - - ph_diff = pd.DataFrame(columns=np.arange(0, N, step=1)) - for m in range(N): - difference = [] - for n in range(len(t)): - diff = np.arctan2(Y[m,n], X[m,n])-np.arctan2(YY[m,n], XX[m,n]) - phasediff = np.arctan2(np.sin(diff), np.cos(diff)) - difference.append(phasediff) - ph_diff[m] = difference - - ph_diff = ph_diff.T - R_all = np.zeros(len(t)) - for col in ph_diff: - sum_phi_all = 0 - count = 0 - for ii in range(len(ph_diff[0])): - if str(ph_diff[col][ii]) != 'nan': - count+=1 - sum_phi_all = sum_phi_all+np.exp(ph_diff[col][ii]*1j) - - R_all[col] = np.abs(sum_phi_all/count) - R_all[len(t)-1] = R_all[len(t)-2] - print('mean ph coh', np.mean(R_all[int(len(t)*0.5):len(t)])) - arn_tong[mm,nn] = np.mean(R_all[int(len(t)*0.5):len(t)]) - - print('---') - -arn_tong = np.flipud(arn_tong) - -np.savetxt('Arn_tongue_extra_intr2.txt', arn_tong) - -#%%Plotting -# file_path = '...' -# arn_tong = np.loadtxt(file_path) -# # arn_tong_dupl = np.hstack((np.fliplr(arn_tong), arn_tong)) - -# own_colors = ["indigo", "rebeccapurple", "mediumorchid", "orchid", "violet", "lightpink", "pink", "lavenderblush"] -# cmap = mcolors.LinearSegmentedColormap.from_list("", own_colors) - -# fig = plt.figure(figsize=(12, 8)) -# plt.xticks([0, 20, 40],[0, 0.1, 0.2]) -# color_map = plt.imshow(arn_tong, interpolation='bicubic', aspect='auto')#, norm=PowerNorm(gamma=5)) -# color_map.set_cmap(cmap) -# plt.xlabel('Intracellular coupling') -# plt.yticks([0, 10, 20, 30, 40, 50],[0.05, 0.04, 0.03, 0.02, 0.01, 0]) -# plt.ylabel('Extracellular coupling') -# ax = plt.gca() -# divider = make_axes_locatable(ax) -# cax = divider.append_axes("right", size="5%", pad=0.1) -# cb = plt.colorbar(cax=cax) -# cb.set_label('Mean phase coherence of phase difference') -# plt.savefig(output+'Arn_tongue_extra_intra.svg') -# plt.show() - - diff --git a/Fig2_S1_S2/Ph_coh_ph_diff_Poincare.py b/Fig2_S1_S2/Ph_coh_ph_diff_Poincare.py deleted file mode 100644 index ae10c14..0000000 --- a/Fig2_S1_S2/Ph_coh_ph_diff_Poincare.py +++ /dev/null @@ -1,110 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 13:10:05 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -output = '.../' - -dt = 0.1 -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -#Parameters -Nosc = 200 -dt = 0.1 -tf = 120 - -gg = 1 -ll = 1 -a0 = 1 -acc = 1 - -kappa_vect = [0, 0.05] #extracellular circadian coupling -eps_vect = [0, 0.15] #intracellular circadian - cell cycle #0.1 or 0.15 -eps1 = 0 - -plt.figure(figsize=(12,10)) - -for kappa in kappa_vect: - for eps2 in eps_vect: - #Vectors initialization - t = np.arange(0,tf,step=dt) - X = np.zeros((Nosc,len(t))) - XX = np.zeros((Nosc,len(t))) - Y = np.zeros((Nosc,len(t))) - YY = np.zeros((Nosc,len(t))) - - #Initial position - # np.random.seed(0) - X[:,0] = np.random.uniform(-1,1,size=(Nosc)) - XX[:,0] = np.random.uniform(-1,1,size=(Nosc)) - Y[:,0] = np.random.uniform(-1,1,size=(Nosc)) - YY[:,0] = np.random.uniform(-1,1,size=(Nosc)) - - #Period distribution - mu1, sigma1 = 22, np.sqrt(6) - period_circ = np.random.normal(mu1, sigma1, Nosc) - mu2, sigma2 = 26, np.sqrt(6) - period_cell = np.random.normal(mu2, sigma2, Nosc) - - - #Integration over time and population - for i in range(len(t) - 1): - sum_x = 0 - sum_y = 0 - for k in range(Nosc): - sum_x = sum_x+X[k,i] - sum_y = sum_y+Y[k,i] - - for m in range(Nosc): - middle = -gg*(np.sqrt(X[m,i]**2+Y[m,i]**2)-a0) - X[m,i+1] = X[m,i]+dt*(middle*X[m,i]-2*np.pi*Y[m,i]/period_circ[m]+sum_x*kappa/(2*Nosc)+eps1*(XX[m,i])) - Y[m,i+1] = Y[m,i]+dt*(middle*Y[m,i]+2*np.pi*X[m,i]/period_circ[m]+sum_y*kappa/(2*Nosc)+eps1*(YY[m,i])) - - mid = -ll*(np.sqrt(XX[m,i]**2+YY[m,i]**2)-acc) - XX[m,i+1] = XX[m,i]+dt*(mid*XX[m,i]-2*np.pi*YY[m,i]/period_cell[m]+eps2*(X[m,i])) - YY[m,i+1] = YY[m,i]+dt*(mid*YY[m,i]+2*np.pi*XX[m,i]/period_cell[m]+eps2*(Y[m,i])) - - - #Phase coherence of phase differences - ph_diff = pd.DataFrame(columns=np.arange(0, Nosc, step=1)) - for m in range(Nosc): - difference=[] - for n in range(len(t)): - diff = np.arctan2(Y[m,n], X[m,n])-np.arctan2(YY[m,n], XX[m,n]) #which oscillator - phasediff = np.arctan2(np.sin(diff), np.cos(diff)) - difference.append(phasediff) - ph_diff[m] = difference - - ph_diff = ph_diff.T - R_all = np.zeros(len(t)) - for col in ph_diff: - sum_phi_all = 0 - count = 0 - for ii in range(len(ph_diff[0])): - if str(ph_diff[col][ii]) != 'nan': - count+=1 - sum_phi_all = sum_phi_all+np.exp(ph_diff[col][ii]*1j) - R_all[col] = np.abs(sum_phi_all/count) - R_all[len(R_all)-1] = R_all[len(R_all)-2] - - plt.plot(t, R_all, linewidth=5, label='k:'+str(kappa)+' eps:'+str(eps)) -plt.xlabel('Time [hours])') -plt.ylabel('Phase coherence of phase differences') -plt.legend(loc='best') -plt.savefig(output+'Ph_coh_ph_diff_4extremecases.svg') -plt.show() - diff --git a/Fig2_S1_S2/Ph_lock_Poincare.py b/Fig2_S1_S2/Ph_lock_Poincare.py deleted file mode 100644 index db804c9..0000000 --- a/Fig2_S1_S2/Ph_lock_Poincare.py +++ /dev/null @@ -1,137 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 12:42:44 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -output = '.../' - -def generate_green_scale(colormap, N): - cmap = plt.cm.get_cmap(colormap) - values = np.linspace(0, 1, N) - colors = cmap(values) - return colors - -dt = 0.1 -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -#Parameters -Nosc = 200 -dt = 0.1 -tf = 300 - -gg = 1 -ll = 1. #works for [0.5, 1.5] -a0 = 1 -acc = 1. # works for [0.5, 1.5] - -kappa = 0.05 #extracellular circadian coupling -eps2 = 0.01 #intracellular circadian - cell cycle #0.01 or 0.015 -eps1 = 0 - -#Vectors initialization -t = np.arange(0,tf,step=dt) -X = np.zeros((Nosc,len(t))) -XX = np.zeros((Nosc,len(t))) -Y = np.zeros((Nosc,len(t))) -YY = np.zeros((Nosc,len(t))) - -#Initial position -X[:,0] = np.random.uniform(-0.3,1,size=(Nosc)) -XX[:,0] = np.random.uniform(-1,1,size=(Nosc)) -Y[:,0] = np.random.uniform(-0.5,1,size=(Nosc)) -YY[:,0] = np.random.uniform(-1,1,size=(Nosc)) - -#Period distribution -mu1, sigma1 = 22, np.sqrt(6) -period_circ = np.random.normal(mu1, sigma1, Nosc) -mu2, sigma2 = 26, np.sqrt(6) -period_cell = np.random.normal(mu2, sigma2, Nosc) - -pert_resset = 0 - -# Integration over time and population -for i in range(len(t) - 1): - sum_x = 0 - sum_y = 0 - for k in range(Nosc): - sum_x = sum_x+X[k,i] - sum_y = sum_y+Y[k,i] - - if int((tf*0.0)/dt) 100: - wAn.compute_spectrum(norm_signal, do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0, smoothing_wsize=5) - rd.set_index(norm_signal.index, inplace=True) - ridge_results1[column] = rd - rd['traceId'] = column - df1 = df1.append(rd) - phases_circ[column] = rd['phase'] - - - wAn.compute_spectrum(signal, do_plot=False) - rd2 = wAn.get_maxRidge(power_thresh=0, smoothing_wsize=5) - rd2.set_index(signal.index, inplace=True) - ridge_results2[column] = rd2 - rd2['traceId'] = column - df2 = df2.append(rd2) - - diff_rel = (rd['phase'] - rd2['phase']) - diff = np.arctan2(np.sin(diff_rel), np.cos(diff_rel)) - - phases_diff[column] = diff - - powers_series = em.average_power_distribution(ridge_results1.values(), signal_ids = ridge_results1.keys()) - powers_series = powers_series.reset_index() - - - #%% Periods --------------------------- - df_period_diff = pd.DataFrame(index=data1.index.values) - for col in data1: - period_circ = ridge_results1[col]['periods'] - period_cc = ridge_results2[col]['periods'] - df_period_diff[col] = period_circ-period_cc - - - #%% Phase coherence ----------------------------- - phases_circ = phases_circ.T - phases_diff = phases_diff.T - R_circ = np.zeros(len(time)) - R_diff = np.zeros(len(time)) - - for col in phases_circ: - sum_phi_circ = 0 - sum_phi_all = 0 - count = 0 - for ii in range(len(phases_circ[col])): - if str(phases_circ[col][ii]) != 'nan': - count+=1 - sum_phi_circ = sum_phi_circ+np.exp(phases_circ[col][ii]*1j) - sum_phi_all = sum_phi_all+np.exp(phases_diff[col][ii]*1j) - R_circ[col] = np.abs(sum_phi_circ/count) - R_diff[col] = np.abs(sum_phi_all/count) - R_circ[len(R_circ)-1] = R_circ[len(R_circ)-2] - R_diff[len(R_diff)-1] = R_diff[len(R_diff)-2] - - - #%%Overlapped plot -------------------------------- - - max_value = max(R_circ) - - tp = 0.5*int(len(time)/2)#0.5*min(range(len(list(R_circ))), key=lambda i: abs(list(R_circ)[i]-max_value/2)) - print(tp) - - detuning = wAn.sinc_smooth(df_period_diff.median(axis=1), T_c=20) - - fig, ax1 = plt.subplots(figsize=(14, 12)) - ax1.plot(time, R_circ/max(R_circ), linewidth=5, color='tab:brown', label='Circadian clocks') - ax1.plot(time, R_diff/max(R_diff), linewidth=5, color='tab:purple', label='Phase differences') - ax1.axvline(x=tp, color='grey', linestyle='--') - ax1.set_ylabel('Relative phase coherence') - ax2 = ax1.twinx() - ax2.plot(time, detuning/(-min(detuning)), linewidth=5, color='tab:grey', label='Detuning circadian cell cycle') - ax2.set_ylabel('Relative detuning [hours]') - ax1.set_xlabel('Time [hours]') - ax1.set_xticks([0,24,48,72,96,120]) - lines, labels = ax1.get_legend_handles_labels() - lines2, labels2 = ax2.get_legend_handles_labels() - ax2.legend(lines + lines2, labels + labels2, loc='upper left') - plt.savefig(output+'Median_period_ph_coh_circ_ph_diff_weak_coupling_experimental.svg') - plt.show() - - - #%%Split the experiment in two parts and choose the first/second part of the data - - #1st part - select_index = powers_series.loc[powers_series[0]>10]['index'] - circadian = df1.loc[df1['traceId'].isin(select_index)] - cell_cycle = df2.loc[df2['traceId'].isin(select_index)] - - select_index_2=cell_cycle.loc[cell_cycle['time']<60 - cell_cycle=cell_cycle.loc[cell_cycle['time'].isin(select_index_2)] - circadian=circadian.loc[circadian['time'].isin(select_index_2)] - - ph_cell = (cell_cycle['phase']/(2*np.pi)).to_numpy() - ph_circ = (circadian['phase']/(2*np.pi)).to_numpy() - - fig = plt.figure(figsize = (12,10)) - hist2D = plt.hist2d(ph_circ, ph_cell, bins=np.arange(0,1.01,0.05), cmap='plasma', density=True, vmin=0.5, vmax=1.5) - plt.colorbar() - plt.xlabel(r'Circadian phase ($\theta$/2$\pi$)') - plt.ylabel(r'Cell cycle phase ($\theta$/2$\pi$)') - # plt.savefig(output+'Phase_lock'+str(condition)+'1st.svg') - plt.show() - - #2nd part - select_index = powers_series.loc[powers_series[0]>10]['index'] - circadian = df1.loc[df1['traceId'].isin(select_index)] - cell_cycle = df2.loc[df2['traceId'].isin(select_index)] - - select_index_2=cell_cycle.loc[cell_cycle['time']>tp]['time'] # change ><60 - cell_cycle=cell_cycle.loc[cell_cycle['time'].isin(select_index_2)] - circadian=circadian.loc[circadian['time'].isin(select_index_2)] - - ph_cell = (cell_cycle['phase']/(2*np.pi)).to_numpy() - ph_circ = (circadian['phase']/(2*np.pi)).to_numpy() - - fig = plt.figure(figsize = (12,10)) - hist2D = plt.hist2d(ph_circ, ph_cell, bins=np.arange(0,1.01,0.05), cmap='plasma', density=True, vmin=0.5, vmax=1.5) - plt.colorbar() - plt.xlabel(r'Circadian phase ($\theta$/2$\pi$)') - plt.ylabel(r'Cell cycle phase ($\theta$/2$\pi$)') - # plt.savefig(output+'Phase_lock'+str(condition)+'2nd.svg') - plt.show() - \ No newline at end of file diff --git a/Fig2_S1_S2/Ph_lock_exp.py b/Fig2_S1_S2/Ph_lock_exp.py deleted file mode 100644 index 5b7483a..0000000 --- a/Fig2_S1_S2/Ph_lock_exp.py +++ /dev/null @@ -1,90 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Tue Feb 27 12:58:03 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer -from pyboat import ensemble_measures as em - -plt.rcParams.update({'font.size': 20}) -plt.rcParams['svg.fonttype'] = 'none' - -path = '.../Raw_data/' -output = '.../Fig2/' -dose = ['untreated','5uM','10uM'] -density = ['high','medium','low'] -channel1 = 'circadian' -channel2 = 'cell_cycle' - -dt = 0.5 -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -for i in density: - for j in dose: - - condition = str(i)+'_density_'+str(j)+'_' - print(condition) - - #Circadian data - data1 = pd.read_csv(path+condition+channel1+r'_filtered.csv', index_col=0) - #Cell-cycle data - data2 = pd.read_csv(path+condition+channel2+r'_filtered.csv', index_col=0) - time = data1.index.values*0.5 - - ridge_results1 = {} - df1 = pd.DataFrame() - ridge_results2 = {} - df2 = pd.DataFrame() - - for column in data1: - norm_signal = data1[column].dropna() - signal = data2[column].dropna() - - if len(norm_signal) > 100: - wAn.compute_spectrum(norm_signal, do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0, smoothing_wsize=5) - rd.set_index(norm_signal.index, inplace=True) - ridge_results1[column] = rd - rd['traceId'] = column - df1 = df1.append(rd) - - wAn.compute_spectrum(signal, do_plot=False) - rd2 = wAn.get_maxRidge(power_thresh=0, smoothing_wsize=5) - rd2.set_index(signal.index, inplace=True) - ridge_results2[column] = rd2 - rd2['traceId'] = column - df2 = df2.append(rd2) - - powers_series = em.average_power_distribution(ridge_results1.values(), signal_ids = ridge_results1.keys()) - powers_series = powers_series.reset_index() - - select_index = powers_series.loc[powers_series[0]>10]['index'] - circadian = df1.loc[df1['traceId'].isin(select_index)] - cell_cycle = df2.loc[df2['traceId'].isin(select_index)] - - # #Split the experiment in two parts and choose the first/second part of the data - # select_index_2=cell_cycle.loc[cell_cycle['time']>60]['time'] # change ><60 - # cell_cycle=cell_cycle.loc[cell_cycle['time'].isin(select_index_2)] - # circadian=circadian.loc[circadian['time'].isin(select_index_2)] - - #Normalization - ph_cell = (cell_cycle['phase']/(2*np.pi)).to_numpy() - ph_circ = (circadian['phase']/(2*np.pi)).to_numpy() - - fig = plt.figure(figsize = (12,10)) - hist2D = plt.hist2d(ph_circ, ph_cell, bins=np.arange(0, 1.01, 0.05), cmap='plasma', density=True, vmin=0.5, vmax=1.5) - plt.colorbar() - plt.xlabel(r'Circadian phase ($\theta$/2$\pi$)') - plt.ylabel(r'Cell cycle phase ($\theta$/2$\pi$)') - plt.savefig(output+'Phase_lock'+str(condition)+'.svg') - plt.show() - diff --git a/Fig2_S1_S2/amplitude_conditions.py b/Fig2_S1_S2/amplitude_conditions.py deleted file mode 100644 index cd5eb4d..0000000 --- a/Fig2_S1_S2/amplitude_conditions.py +++ /dev/null @@ -1,129 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Mon Oct 7 13:21:50 2024 - -@author: nicagutu -""" - -import numpy as np -import matplotlib.pyplot as plt -import pandas as pd -from pyboat import WAnalyzer -from scipy.signal import savgol_filter -from scipy.signal import find_peaks -import seaborn as sns - -plt.rcParams.update({'font.size': 24}) -plt.rcParams['svg.fonttype'] = 'none' - -def remove_outliers(df, column): - Q1 = df[column].quantile(0.25) # First quartile (25th percentile) - Q3 = df[column].quantile(0.75) # Third quartile (75th percentile) - IQR = Q3 - Q1 # Interquartile range - lower_bound = Q1 - 1.5 * IQR - upper_bound = Q3 + 1.5 * IQR - return df[column][(df[column] >= lower_bound) & (df[column] <= upper_bound)] - -path = 'RawData/' -output = 'Figures/' -dose = ['untreated', '5uM', '10uM'] # '1.25uM','2.5uM', -density = ['high'] #,'medium','low' -channel1 = 'circadian' -channel2 = 'cell_cycle' - -dt = 0.5 -lowT = 16 -highT = 32 -periods = np.linspace(lowT, highT, 200) -wAn = WAnalyzer(periods, dt, time_unit_label='hours') - -#Changing the drug concentration -amplitudes_dict = {} -for j in dose: - condition = str(density[0])+'_density_'+str(j)+'_' - print(condition) - - ## Circadian data - # data1 = pd.read_csv(path+condition+channel1+r'_filtered.csv', index_col=0) - data1 = pd.read_csv(path+condition+channel1+r'.csv', index_col=0) - - amplitude_circadian = [] - for column in data1: - signal = data1[column].dropna() - detrended = wAn.sinc_detrend(signal, T_c = 50) - # norm_signal = wAn.normalize_amplitude(signal, window_size=50) - - wAn.compute_spectrum(detrended, do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - power = np.mean(rd['power']) - amplitude_circadian.extend(rd['amplitude']) - - amplitudes_dict[j] = amplitude_circadian - -df_ampl = pd.DataFrame({key: pd.Series(value) for key, value in amplitudes_dict.items()}) -print(df_ampl) - -cleaned_df_ampl = pd.DataFrame() -for col in df_ampl: - cleaned_df_ampl[col] = remove_outliers(df_ampl, col) - -plt.figure(figsize = (10, 8)) -sns.violinplot(data=cleaned_df_ampl, linewidth=2.5) -plt.ylabel('Amplitude [a.u.]') -plt.xlabel('Condition') -plt.savefig(output+'Amplitude_inhibitor.svg') -plt.show() - -plt.figure(figsize = (10, 8)) -sns.boxplot(data=cleaned_df_ampl, linewidth=2.5) -plt.ylabel('Amplitude [a.u.]') -plt.xlabel('Condition') -plt.show() - -#%% Changing the seeding density -dose = ['untreated']#, '5uM', '10uM'] # '1.25uM','2.5uM', -density = ['high' ,'medium','low'] - -amplitudes_dict = {} -for j in density: - condition = str(j)+'_density_'+str(dose[0])+'_' - print(condition) - - ## Circadian data - # data1 = pd.read_csv(path+condition+channel1+r'_filtered.csv', index_col=0) - data1 = pd.read_csv(path+condition+channel1+r'.csv', index_col=0) - - amplitude_circadian = [] - for column in data1: - signal = data1[column].dropna() - detrended = wAn.sinc_detrend(signal, T_c = 50) - # norm_signal = wAn.normalize_amplitude(signal, window_size=50) - - wAn.compute_spectrum(detrended, do_plot=False) - rd = wAn.get_maxRidge(power_thresh=0,smoothing_wsize=5) - power = np.mean(rd['power']) - amplitude_circadian.extend(rd['amplitude']) - - amplitudes_dict[j] = amplitude_circadian - -df_ampl = pd.DataFrame({key: pd.Series(value) for key, value in amplitudes_dict.items()}) - -cleaned_df_ampl = pd.DataFrame() -for col in df_ampl: - cleaned_df_ampl[col] = remove_outliers(df_ampl, col) - - -plt.figure(figsize = (10, 8)) -sns.violinplot(data=cleaned_df_ampl, linewidth=2.5) -plt.ylabel('Amplitude [a.u.]') -plt.xlabel('Condition') -plt.savefig(output+'Amplitude_density.svg') -plt.show() - -plt.figure(figsize = (10, 8)) -sns.boxplot(data=cleaned_df_ampl, linewidth=2.5) -plt.ylabel('Amplitude [a.u.]') -plt.xlabel('Condition') -plt.show() - diff --git a/Fig2_S1_S2/voronoi_diagrams.py b/Fig2_S1_S2/voronoi_diagrams.py deleted file mode 100644 index b839094..0000000 --- a/Fig2_S1_S2/voronoi_diagrams.py +++ /dev/null @@ -1,179 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Created on Fri Sep 20 12:00:10 2024 - -@author: nicagutu -""" - -import numpy as np -import pandas as pd -from scipy.spatial import Voronoi, voronoi_plot_2d, distance -from scipy.spatial.distance import pdist, squareform -import matplotlib.pyplot as plt -from scipy.spatial import ConvexHull -from scipy.signal import savgol_filter -import seaborn as sns - -plt.rcParams.update({'font.size': 22}) - -# Function to compute areas of Voronoi cells -def compute_voronoi_areas(vor): - def polygon_area(vertices): - if len(vertices) < 3: # Not a polygon - return 0 - hull = ConvexHull(vertices) - return hull.volume - - areas = [] - for region in vor.regions: - if not -1 in region and len(region) > 0: # Exclude infinite regions - polygon = [vor.vertices[i] for i in region] - areas.append(polygon_area(polygon)) - - return np.array(areas) - -# Mean, Variance, Standard Deviation, Coefficient of Variation (CV) -def compute_voronoi_statistics(voronoi_areas): - mean_area = np.mean(voronoi_areas) - median_area = np.median(voronoi_areas) - variance_area = np.var(voronoi_areas) - sd_area = np.std(voronoi_areas) - cv_area = (sd_area / mean_area) * 100 # Coefficient of Variation in percentage - - # Create a dictionary with all the statistics - stats = { - "Statistic": ["Mean", "Median", "Variance", "Standard Deviation", "Coefficient of Variation (CV)"], - "Value": [mean_area, median_area, variance_area, sd_area, cv_area] - } - - stats_table = pd.DataFrame(stats) - return stats_table - -# Function to remove outliers using the IQR method -def remove_outliers(data): - q1 = np.percentile(data, 25) # First quartile (25th percentile) - q3 = np.percentile(data, 75) # Third quartile (75th percentile) - iqr = q3 - q1 # Interquartile range (IQR) - lower_bound = q1 - 1.5 * iqr # Lower bound - upper_bound = q3 + 1.5 * iqr # Upper bound - - # Filter the data within the bounds - filtered_data = [x for x in data if lower_bound <= x <= upper_bound] - return filtered_data - -path = 'RawData/' -output = 'Figures/' -list_cond = ['100', '50', '25'] - -statistics_table_all = {} -areas_all = {} -for cond in list_cond: - centroid_col_df = pd.read_excel(path+'Absolute_centroid_col_'+str(cond)+'.xlsx', index_col=0) - centroid_row_df = pd.read_excel(path+'Absolute_centroid_row_'+str(cond)+'.xlsx', index_col=0) - - # Extract the x and y coordinates, removing the first index column - # x_coords = centroid_col_df.iloc[0, :].values#.flatten() - # y_coords = centroid_row_df.iloc[0, :].values#.flatten() - - x_coords = centroid_col_df.median(axis=0) - y_coords = centroid_row_df.median(axis=0) - - # plt.plot(y_coords, x_coords, 'o') - # plt.show() - - # Combine x and y coordinates into an array of points - points = np.column_stack((y_coords, x_coords)) - - # Compute the Voronoi diagram - vor = Voronoi(points) - - # Plot the Voronoi diagram - fig, ax = plt.subplots(figsize=(12, 12)) - voronoi_plot_2d(vor, ax=ax, show_vertices=False, line_colors='tab:purple', alpha=0.5, line_width=1.5, point_size=2) - - # Display the plot - plt.title('Voronoi Diagram for '+str(cond)+'% density') - plt.savefig(output+'Vornoi_diagrams_'+str(cond)+'.svg') - plt.show() - - - # Compute Voronoi cell areas - voronoi_areas = compute_voronoi_areas(vor) - - # Compute statistics of Voronoi cell areas - statistics_table = compute_voronoi_statistics(voronoi_areas) - statistics_table_all[cond] = statistics_table - - # Convert areas to numpy array for easier manipulation - areas = np.array(voronoi_areas) - - # Store the areas in the dictionary, keyed by condition - areas_all[cond] = list(remove_outliers(voronoi_areas)) - - # Plot histogram of areas for this condition - plt.figure(figsize=(8, 6)) - plt.hist(remove_outliers(voronoi_areas), bins=30, color='skyblue', edgecolor='black', density=True) - plt.title(f'Distribution of Voronoi Cell Area for {cond}% density') - # plt.xlim([0,30]) - plt.xlabel('Area of Voronoi Cells') - plt.ylabel('Density') - plt.savefig(output+'Distribution_Voronoi_cell_area_'+str(cond)+'.svg') - plt.show() - -print(statistics_table_all) - -# Find the maximum length of all the lists -max_length = max(len(areas) for areas in areas_all.values()) - -# Pad each list with None values to make them the same length -for cond in areas_all: - areas_all[cond] = areas_all[cond] + [None] * (max_length - len(areas_all[cond])) - -# Convert areas_all dictionary into a DataFrame -df = pd.DataFrame(areas_all) - -# Plot boxplot and violin plot using seaborn -fig, axs = plt.subplots(1, 2, figsize=(12, 6)) - -# Boxplot -sns.boxplot(data=df, ax=axs[0]) -axs[0].set_title('Boxplot') - -# Violin plot -sns.violinplot(data=df, ax=axs[1]) -axs[1].set_title('Violin Plot') - -# Show the plots -plt.tight_layout() -plt.show() - -# Calculate the mean and standard deviation for each condition -means = df.mean() -std_devs = df.std() - -# Create a bar plot with error bars showing standard deviation -plt.figure(figsize=(10, 8)) -plt.bar(means.index, means.values, yerr=std_devs.values, capsize=5, color='skyblue', edgecolor='black') -plt.xlabel('Seeding density') -plt.ylabel('Mean Voronoi Cell Area [um]') -plt.savefig(output+'Barplot_mean_std_seeding_density_Voronoi_cell_area.svg') -plt.show() - -# Convert areas_all dictionary into a long-format DataFrame -df_long = pd.DataFrame([ - {'Condition': cond, 'Area': area} - for cond, areas in areas_all.items() - for area in areas if area is not None -]) - -# Create a strip plot with error bars showing the mean and standard deviation -plt.figure(figsize=(10, 8)) -sns.stripplot(x='Condition', y='Area', data=df_long, jitter=True, color='skyblue', alpha=0.5) -sns.pointplot(x='Condition', y='Area', data=df_long, ci='sd', capsize=0.2, color='tab:blue', join=False) -plt.xlabel('Seeding density [%]') -plt.ylabel('Voronoi Cell Area [um^2]') -plt.savefig(output+'Striplot_mean_std_seeding_density_Voronoi_cell_area.svg') -plt.show() - -