CIRCADIPY
This is the pipeline for basic execution of the CircadiPy library. Here you will find information needed to execute several features of the library.
1 - IMPORT LIBRARIES TO USE THE NOTEBOOK
[29]:
## IMPORT LIBRARIES ----------------------------------------------------------------------------------------------------
%matplotlib qt
%reload_ext autoreload
%autoreload 3
import matplotlib.pyplot as plt # Import matplotlib.pyplot to plot figures
import tkinter as tk # Import TK to open folder dialogs to select folders
from tkinter import filedialog # Import filedialog to open folder dialogs to select folders
import numpy # Import numpy to work with arrays and make calculations
import time # Import time to measure time
import warnings
from sys import platform # Import warnings to ignore warnings
warnings.filterwarnings('ignore') # Ignore warnings
## IMPORT CIRCADIPY ----------------------------------------------------------------------------------------------------
import chrono_reader as chr # Import chrono_reader to read data
import chrono_plotter as chp # Import chrono_plotter to plot data
import chrono_rithm as chrt # Import chrono_rithm to make calculations
import chrono_simulation as chs # Import chrono_simulation to simulate data
2 - SELECT THE FOLDER TO SAVE THE SIMULATION
In this first step we will generate a representative time series of the circadian activity of an experiment animal. The activity will be mimicked by sinusoids with noise that represent the number of movements made by an animal. To do this, we will first need to select a directory in which to save the simulated data
[10]:
## SET THE ENVIRONMENT -------------------------------------------------------------------------------------------------
if platform == "win32":
root = tk.Tk() # Create the root window
root.attributes('-topmost',True) # Keep the window on top of others
root.iconify() # Hide the root window
main_save_folder = filedialog.askdirectory(title="Select the folder to save the simulated data", parent=root) # Ask for the folder to save the simulated data
root.destroy() # Destroy the root window
else:
## If you are using a Mac, you will not be able to select the folder using the code above.
## In this case, please set the folder manually here:
main_save_folder = "Paste here the PATH/TO/SAVE/SIMULATED/DATA"
print(main_save_folder) # Print the folder to save the simulated data
E:/github/circadipy/src/simulated_data
3 - RUN AND SAVE A SIMPLE SIMULATION
Each simulation has a given number stages (num_stages), each stage can be configured with a single activity period. In this example we will simulate a protocol that has 3 stages and a sample rate of 5 minutes (this means that each 5 minutes, a measure is done):
A: a 5-day control period with a 12-12 hour light-dark (LD) cycle,
B: a 10-day free-running period with a 24-hour dark (DD) cycle,
C: a further 5 days of readaptation, with a 12-12 hour light-dark (LD) cycle.
Following a biological pattern, animals regulated to a 12-12 hour light-dark cycle have an activity period of 24 hours, however, when turning off the light track, leaving it in a 24-hour dark cycle, these animals tend to present a shorter period (in the example we will use 23 hours and 30 minutes). These animals tend to synchronize again when the light-dark cycle is reactivated, returning to a 24-hour period.
[70]:
## CONFIGURE THE SIMULATION --------------------------------------------------------------------------------------------
sampling_frequency = '5T' # Sampling frequency of the simulated data (30T = 30 minutes)
num_stages = 3 # Number of cycles in each simulation
days_per_stage = [5, 10, 5] # Number of days in each stage
activity_period = [24, 23.5, 24] # Activity period in each stage
cycles_per_stage = ['LD', 'DD', 'LD'] # Type of cycles in each stage
stage_labels = ['A', 'B', 'C'] # Labels of each stage
To perform the simulation, the simulation module (chrono_simulation imported as chs) of CircadiPy is used and the parameters discussed above are passed as arguments. The simulation can be generated in sine, square, square with low-pass filter, sawtooth and triangular format. In addition, Gaussian noise and its signal-to-noise ratio can be added. The user can also select whether negative values will be considered or not.
Note: the simulated data was saved in an .asc file, which is a simple text file. As long as the pattern of this file is followed, CircadiPy is able to read and process your data.
[104]:
## SIMULATE A EXPERIMENT AND SAVE THE DATA -----------------------------------------------------------------------------
file_name = 'simulated_data' # Create the file name for the simulated data with noise
raw_data = chs.simulate_protocol(file_name, main_save_folder, sampling_frequency, days_per_stage, activity_period,
signal_type = 'square', noise = True, snr_db = 20, only_positive = True,
remove_file = True) # Simulate the data with noise (chrono_simulation.py)
Saved successfully in E:/github/circadipy/src/simulated_data/simulated_data.asc
Run the next block to visualize the generated data.
[105]:
## PLOT THE RAW DATA ---------------------------------------------------------------------------------------------------
fig = plt.figure(figsize = (12, 5)) # Create a figure
ax = fig.add_subplot(111) # Add a subplot
ax.plot(raw_data, color = 'maroon', linewidth = 2) # Plot the data with noise
ax.set_xlabel('Time (Hours)') # Set the x label
ax.set_ylabel('Activity') # Set the y label
ax.set_title('Simulated Data With Noise') # Set the title
plt.show() # Show the plot
4 - IMPORT THE SIMULATION AND BUILD THE PROTOCOL TO BE STUDIED
To study a certain protocol, simply use the read_protocol function of the chrono_reader module (imported as chr). Just configure the following parameters: - Protocol name: a representative name (e.g. ‘simulation_example’) - File containing the data: the .asc file path (can be other formats). It was previously saved using the variables main_save_folder and file_name. It can also be accessed in one of the previously executed blocks. - ZT0 reference: ZT is the time the light was switched on (in this example we will use 00:00) - Labels dictionary: a dictionary containing the characteristics of each experimental stage. It is necessary to configure the type of cycle (e.g. light-dark , LD), a stage label (arbitrary name that will differentiate the stages) and the number of days of the stage. In our case, we will use the same variables used to create the simulation - Type of data: the type of file to be read (it can be a generic .asc file, or a predefined template such as intellicage) - Flag to consider the first day: if True, the first day of the experiment will be considered, if False, no
At the end, the protocol object (or variable) will be saved, it contains all the information needed to apply various features of CircadiPy. The data is stored in this object in the form of a table, which can be seen as a result of the block execution.
[106]:
## SET THE DATA FILE ---------------------------------------------------------------------------------------------------
# simulation_file = 'insert/your/file/directory.asc # If you want to use other dataset, please uncomment this line and set yout data path
simulation_file = main_save_folder + "/" + file_name + '.asc' # Set the path to the simulation description file
print('Data read from: ' + simulation_file)
## CONFIGURE THE EXPERIMENT --------------------------------------------------------------------------------------------
zt_0_time = 0 # Set the time for the first ZT (0)
type = 'generic' # Set the type of the data decoding (endocrino or intellicage)
labels_dict = {'cycle_types': ['LD', 'DD', 'LD'], 'test_labels': ['A', 'B', 'C'], 'cycle_days': [5, 10, 5]} # Create a dictionary to store to pass as argument to the read_protocol function
## CREATE THE PROTOCOL OBJECT WITH THE EXPERIMENT ----------------------------------------------------------------------
protocol = chr.read_protocol('simulation', simulation_file, zt_0_time = zt_0_time, labels_dict = labels_dict,
type = type, consider_first_day = True) # Read the protocol (chrono_reader.py)
print('\nTable generated after the data importing')
display(protocol.data) # Display the data
Data read from: E:/github/circadipy/src/simulated_data/simulated_data.asc
Table generated after the data importing
| values | is_night | cycle_types | test_labels | real_date | day | |
|---|---|---|---|---|---|---|
| 2022-01-01 00:00:00 | 101.0 | False | LD | A | 2022-01-01 00:00:00 | 2022-01-01 |
| 2022-01-01 00:05:00 | 0.0 | False | LD | A | 2022-01-01 00:05:00 | 2022-01-01 |
| 2022-01-01 00:10:00 | 10.0 | False | LD | A | 2022-01-01 00:10:00 | 2022-01-01 |
| 2022-01-01 00:15:00 | 2.0 | False | LD | A | 2022-01-01 00:15:00 | 2022-01-01 |
| 2022-01-01 00:20:00 | 0.0 | False | LD | A | 2022-01-01 00:20:00 | 2022-01-01 |
| ... | ... | ... | ... | ... | ... | ... |
| 2022-01-20 23:35:00 | 0.0 | True | LD | C | 2022-01-20 23:35:00 | 2022-01-20 |
| 2022-01-20 23:40:00 | 8.0 | True | LD | C | 2022-01-20 23:40:00 | 2022-01-20 |
| 2022-01-20 23:45:00 | 4.0 | True | LD | C | 2022-01-20 23:45:00 | 2022-01-20 |
| 2022-01-20 23:50:00 | 3.0 | True | LD | C | 2022-01-20 23:50:00 | 2022-01-20 |
| 2022-01-20 23:55:00 | 2.0 | True | LD | C | 2022-01-20 23:55:00 | 2022-01-20 |
5760 rows × 6 columns
CircadiPy allows the user to resample the signal (modify the number of samples per second), apply filters for noise removal and normalize the data. In the following example, the data will be resampled from 5 minutes to 10 minutes, in addition, a second order moving average filter will be applied using a 3-sample (in this case 30 minutes) window and, finally, a normalization will be made, adjusting the minimum and maximum values of each day in a range from 0 to 1. The result of the modifications is shown by the graph generated at the end of the block.
To show how the modifications that preprocessing causes to the signal, we will use the basic plot functions of the chrono_plotter module. The time_series function will show the whole time series, while time_series_sum_per_day will sum all the values corresponding to each day.
Note: If you want to skip this step and not apply any pre-processing, move on to the next block.
[107]:
## PLOT THE DATA BEFORE PREPROCESSING ----------------------------------------------------------------------------------
chp.time_serie(protocol, title = 'Time Series (Before)', x_label = 'Time (Days)', y_label = 'Amplitude',
color = 'midnightblue', save_folder = None, save_suffix = '') # Plot the data (using chrono_plotter.py)
## RESAMPLE, FILTER AND NORMALIZE --------------------------------------------------------------------------------------
protocol.resample('10T', method = 'sum') # Resample the data (chrono_reader.py)
protocol.apply_filter(window = 3, type = 'moving_average', order = 2, reverse = False) # Apply a filter to the data (chrono_reader.py)
protocol.normalize_data(type = 'minmax', per_day = True) # Normalize the data (chrono_reader.py)
## PLOT THE DATA AFTER PREPROCESSING -----------------------------------------------------------------------------------
chp.time_serie(protocol, title = 'Time Series (After)', x_label = 'Time (Days)', y_label = 'Amplitude',
color = 'midnightblue', save_folder = None, save_suffix = '') # Plot the data (using chrono_plotter.py)
chp.time_serie_sum_per_day(protocol, title = 'Sum of Time Series Per Day', x_label = 'Time (Days)',
y_label = 'Amplitude', color = 'midnightblue', save_folder = None, save_suffix = '') # Plot the sum of the data per day (using chrono_plotter.py)
5 - RUN THE COSINOR MODELING
In this block, the pre-processed (or not) data will be submitted to a parameterization process of several cosine models. The user will select those periods (or frequencies) that he wants to test and, for each of them, the best model will be calculated and its parameters defined, being them: the amplitude, acrophase, mesor. At the end of this step, the best model among all the periods tested will be selected and its parameters used to characterize data. this process will be carried out for each stage separately. Internally, CircadiPy will organize the necessary table so that the CosinorPy library can run and estimate the best model, in addition, it will use the best model to produce graphs and tables representing the parameterization.
To set up the modeling, a dictionary (dict) will be used as parameter of the fit_cosinor function. In this dictionary the following variables should be configured.
The time_shape can be ‘continuous’, ‘meadian’ or ‘mean’. If it is ‘continuous’, all samples will be cosider to calculate the model. If it is ‘median’ or ‘mean’, the samples of activity will be grouped by mean/median along one day.
The step is the time step to calculate the model (in hours).
The start/end_time defines the period that the model will try to fit. E.g. if the start time is 24, the end time is 26, and the step is 0.5, the model will try to fit the data to a cosine curve with period equal to 24, 24.5, 25, 25.5 and 26 hours.
The n_components is the number of components that the model will try to fit. If it is a list, the model will try to fit all the number of components in the list. E.g. if n_components = [1,2,3], the model will try to fit the data with 1, 2 and 3 components.
Then, the best model obtained for each stage can be used to generate daily models. In this functionality, the period selected as the best model for each stage will be fixed and the acrophase, mesor and amplitude parameters will be obtained for each day.
Note: CosinorPy provides as output all parameters and their statistical significance.
[108]:
## SET THE PARAMETERS FOR THE COSINOR FITTING --------------------------------------------------------------------------
dict = {'time_shape': 'continuous',
'step': 0.01,
'start_time': 22,
'end_time': 26,
'n_components': [1]} # Create a dictionary to pass as argument to the fit_cosinor function
## FIT THE COSINOR TO THE DATA -----------------------------------------------------------------------------------------
init = time.time() # Get the initial time
best_models, best_models_file = chrt.fit_cosinor(protocol, dict = dict, save_folder = main_save_folder) # Fit the cosinor to the data (chrono_rithm.py)
best_models_fixed, best_models_fixed_file = chrt.fit_cosinor_fixed_period(protocol, best_models,
save_folder = main_save_folder) # Fix the best period calculated and fit the cosinor for each day using this period (chrono_rithm.py)
end = time.time() - init # Get the time elapsed
print("Cosinor fitted to " + protocol.name + " and results saved!") # Print that the cosinor was fitted and the results saved
print("Time elapsed: " + "{:.2f}".format(end) + " seconds") # Print the time elapsed
display(best_models) # Display the best models
Cosinor fitted to simulation and results saved!
Time elapsed: 33.13 seconds
| significant | test | period | n_components | p | q | p_reject | q_reject | amplitude | CI(amplitude) | ... | CI(acrophase) | p(acrophase) | q(acrophase) | mesor | CI(mesor) | p(mesor) | q(mesor) | acrophase_zt | acrophase_zt_lower | acrophase_zt_upper | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | A | 24.14 | 1.0 | 1.110223e-16 | 1.158494e-16 | NaN | NaN | 0.580089 | [0.5596956018197057, 0.6004821916021363] | ... | [1.628973968528041, 1.7044476659115793] | 0.0 | 0.0 | 0.481137 | [0.4677809779739576, 0.49449351542240333] | 0.0 | 0.0 | 17.596496 | 17.741481 | 17.451512 |
| 1 | 1 | B | 23.51 | 1.0 | 1.110223e-16 | 1.158494e-16 | NaN | NaN | 0.586008 | [0.5680602716599603, 0.6039556246030819] | ... | [0.8993948938205021, 0.9748685912040402] | 0.0 | 0.0 | 0.482428 | [0.47275207399766567, 0.4921039752479259] | 0.0 | 0.0 | 20.493504 | 20.634705 | 20.352302 |
| 2 | 1 | C | 23.96 | 1.0 | 1.110223e-16 | 1.158494e-16 | 1.110223e-16 | NaN | 0.578468 | [0.5536516768067955, 0.6032844665068849] | ... | [2.660447832769734, 2.735921530153274] | 0.0 | 0.0 | 0.481172 | [0.4675651208797714, 0.4947779483121282] | 0.0 | 0.0 | 13.710871 | 13.854775 | 13.566967 |
3 rows × 23 columns
5a - COMPARE THE MODEL RESULTS FOR EACH CYCLE
In this section, we can compare the results of the cosinor fitting with the ground-truth period used to generate the simulated data. The mean and standard deviation (stage A, B and C) of the error between the period calculated by the cosinor fitting and the period used to generate the data is printed.
[109]:
## COMPARE THE COSINOR FITTING WITH THE GORUND-TRUTH PERIODS -----------------------------------------------------------
model_error = [] # Create an empty list to store the model error
activity_period_model = list(best_models['period']) # Get the activity period calculated by the model
animal_model_error = [abs(a - b) for a,b in zip(activity_period_model, activity_period)] # Calculate the error between the model and the real activity period used to generate the data
model_error.extend(animal_model_error) # Store the model error in the list
print("Mean error: " + "{:.2f}".format(numpy.mean(model_error)) + " +/- " +
"{:.2f}".format(numpy.std(model_error)) + " hours") # Print the mean and standard deviation of the model error
Mean error: 0.06 +/- 0.06 hours
6 - PLOT THE RESULTS
Now, we will plot actograms, periodograms and the results of the cosinor fitting. To do this, we will use the chrono_plotter module.
Note: If you want to save the images in the directory, please change the plot variable to False, if you want only to plot the images, please change the plot variable to True
[110]:
## CONFIGURE IF THE PLOTS WILL BE SAVED OR SHOWN -----------------------------------------------------------------------
plot = False # Set the plot parameter
if plot == False: # If the plot parameter is False, save the images in the directory
print("The plots will be saved in the folder: " + main_save_folder)
save_folder = main_save_folder
else:
print("The plots will be shown in the screen")
save_folder = None # If the plot parameter is True, save the images in the directory
The plots will be saved in the folder: E:/github/circadipy/src/simulated_data
[111]:
## PLOT ACTOGRAMS ------------------------------------------------------------------------------------------------------
chp.actogram_bar(protocol, first_hour = 18, save_folder = save_folder, save_suffix = 'bar') # Plot the actogram with bars (using chrono_plotter.py)
chp.actogram_colormap(protocol, first_hour = 18, save_folder = save_folder, save_suffix = 'colormap', norm_color = None)# Plot the actogram with colormap (using chrono_plotter.py)
## PLOT PERIODOGRAMS ---------------------------------------------------------------------------------------------------
chp.data_periodogram(protocol, time_shape = 'continuous', method = 'periodogram', max_period = 48,
unit_of_measurement = 'Amplitude', save_folder = save_folder, save_suffix = 'periodogram') # Plot the periodogram (using chrono_plotter.py)
chp.data_periodogram(protocol, time_shape = 'continuous', method = 'welch', max_period = 48,
unit_of_measurement = 'Amplitude', save_folder = save_folder, save_suffix = 'welch') # Plot the welch periodogram (using chrono_plotter.py)
[112]:
## PLOT THE COSINOR FITS -----------------------------------------------------------------------------------------------
chp.model_overview_detailed(protocol, best_models_fixed, save_folder = main_save_folder) # Plot representative parameters of the model fitted to the data for each day (using chrono_plotter.py)
chp.model_overview(protocol, best_models, save_folder = main_save_folder) # Plot representative parameters of the model fitted to the data for each stage (using chrono_plotter.py)
chp.model_over_signal(protocol, best_models, position = 'head', mv_avg_window = 1, save_folder = main_save_folder) # Plot representative figure with the data and the respective model (using chrono_plotter.py)
7 - RUN THE MODELING FOR SPECIFIED DAY WINDOWS
In this last section the user can parameterize models for each day separately, in addition, the user can also select a moving window of days to run the modeling (e.g. if the window is 3, the model will be generated using data from day 1 to 3, then from day 2 to 4 and so on).
[113]:
## SET THE PARAMETERS FOR THE COSINOR FITTING --------------------------------------------------------------------------
dict = {'day_window': 2,
'step': 0.01,
'start_time': 22,
'end_time': 26,
'n_components': [1]} # Create a dictionary to pass as argument to the fit_cosinor_per_day function
## FIT A COSINOR MODEL FOR EACH DAY ------------------------------------------------------------------------------------
best_models_per_day, best_models_fixed_file = chrt.fit_cosinor_per_day(protocol, dict = dict, plot = True,
save_folder = save_folder) # Fit the cosinor to the data for each day or window of days (chrono_rithm.py)