PiQSARS: A pipeline for quantitative and statistical analyses of ratiometric fluorescent biosensors.

Genetically encoded ratiometric fluorescent probes are cutting-edge tools in biology. They allow precise and dynamic measurement of various physiological parameters within cell compartments. Because data extraction and analysis are time consuming and may lead to inconsistencies between results, we describe here a standardized pipeline for•Semi-automated treatment of time-lapse fluorescence microscopy images.•Quantification of individual cell signal.•Statistical analysis of the data.First, a dedicated macro was developed using the FIJI software to reproducibly quantify the fluorescence ratio as a function of time. Raw data are then exported and analyzed using R and MATLAB softwares. Calculation and statistical analysis of selected graphic parameters are performed. In addition, a functional principal component analysis allows summarizing the dataset. Finally, a principal component analysis is performed to check consistency and final analysis is presented as a visual diagram. The method is adapted to any ratiometric fluorescent probe. As an example, the analysis of the cytoplasmic HyPer probe in response to an acute cell treatment with increasing amounts of hydrogen peroxide is shown. In conclusion, the pipeline allows to save time and analyze a larger amount of samples while reducing manual interventions and consequently increasing the robustness of the analysis.

Specifications table

Method details

Introduction and general background information

Since their discovery, fluorescent probes have been engineered to accurately monitor and quantify in real-time an increasing variety of physiological parameters in living cells such as ion concentration variations [1], cysteines [2], pH [3], enzymatic activity [4] or oxidative stress [5], [6], [7], [8], but also various cellular compounds [9], [10], [11], [12]. Genetically-encoded fluorescent proteins are of particular interest as they can be addressed to specific organelles, such as the nucleus or the mitochondria, allowing a subcellular resolution [13,14]. For some probes, the readout is a simple change in their fluorescence intensity, whereas others show a shift in their excitation or emission spectrum. These latter are the most reliable since their spectral shift is independent of the illumination intensity, the probe concentration or photobleaching [15], [16], [17], [18]. To record a cell parameter in real time with these probes, either simultaneous excitation at two different wavelengths or simultaneous signal recovery at two emission wavelengths is performed depending if the shift deals with the excitation or the emission spectrum. The ratio of both signals is then calculated, which is independent of the probe local concentration. These probes are thus called ratiometric.

Image acquisition with ratiometric probes generates a tremendous amount of data. Data processing and analysis are laborious, yet crucial for a reproducible and efficient exploitation of these biological tools. Macros, like FluoQ [19] and softwares, like Friker's one [18] are available for analysis of multidimensional images. Fricker's software, for example, enables to quantify ratio data from 4D-images (x,y,z,t) with up to 5 fluorescence channels in a single analysis round [18]. However, this software is rather dedicated to optical sections of entire plants or fungi and does not perform individual cell segmentation. In the same way, the FluoQ macro allows the rapid processing of large quantities of images, but the cells are segmented on a t-projection of the timelapse, which is efficient only for cells with no or little mobility [19]. Consequently, both tools do not allow easy tracking of individual moving mammalian cells in culture. To circumvent this drawback, we developed a complementary approach that considers morphological changes and movements of individual cells over time during image acquisition, in response to an exogenous stimulus. Derived from Belousov's methods related to a redox ratiometric probe [20,21], we developed an automated FIJI-encoded, user-friendly macro that preprocesses and quantifies data from any ratiometric probe. Secondly, we developed a set of tools for statistical analyses of the generated data using R software. This method paper describes this pipeline step-by-step, as summarized in Fig. 1.

Fig. 1

Schematic diagram of the analysis pipeline. Starting from raw images, a two-step analysis procedure is shown within the gray rectangle. The FIJI macro extracts quantitative data from images, while R and MATLAB scripts perform and display data statistical analysis.

As a proof-of-concept, we used time-lapse data collected from cells treated with varying H2O2 concentrations and expressing the ratiometric H2O2-specific HyPer probe, which is routinely used in our laboratory [22,23]. HyPer is a widely used, genetically encoded fluorescent redox biosensor based on the Yellow Fluorescent Protein (YFP). It senses H2O2 concentration variations within living cells allowing oxidative stress monitoring [24]. Of note, other YFP-derived or GFP-derived probes may also be used, such as roGFP2 [17]. The reliability and robustness of our method is discussed.

Experimental setup

Cell line and performed experiments

We stably expressed nucleus-excluded cytoplasmic HyPer under the control of a doxycycline-inducible promoter in rabbit kidney epithelial cells (RK13 cell line [25]). HyPer fluorescence was monitored for 65 min and 0, 10, 40 or 75 µM H2O2 were added to the cells 5 min after beginning the acquisition. Each field contained 3 to 6 cells, and each experiment was repeated 3 times. Between 10 and 20 curves for each H2O2 concentration tested were thus obtained.

Cell culture

Before time-lapse acquisitions, cells were plated in glass-bottom dishes using a weakly fluorescent red phenol-free culture medium. Cell density must be low (300 – 500 cells.cm-2) for easier single cell analysis.

Cell imaging

Perform acquisition at each excitation peak of your probe of interest. Table 1 provides the theoretical values of the excitation wavelengths for HyPer, as well as possible experimental acquisition setups.

Table 1

Theoretical fluorescence peaks of HyPer and possible experimental setups for fluorescence microscopy acquisition. Theoretical maxima are compiled from our observations and from Meyer and Dick [17].

		Theoretical maxima	Possible setups
Excitation	Oxidized probe	500 nm	488 nm laser469 nm LED511 nm LED
	Reduced probe	420 nm	405 nm laser423 nm LED
Emission	Oxidized and reduced probe	516 nm	> 500 nm

Set acquisition parameters so that cells are bright enough to ensure a satisfactory dynamic range without any noticeable phototoxicity. The image dynamic is satisfactory when it reaches at least half of the maximum dynamic range (i.e. 2048 gray levels for a 12-bit image). Dynamic can be improved by post acquisition additive binning using the provided macro (see Step 2: image processing).

With HyPer, perform excitation near 420 nm first for each frame, as this reverses photo-conversion that occurs when exciting the probe at 500 nm [20].

Acquisition parameters optimization was based on previous work [20]. Initial fluorescence intensities in both channels must be similar and the ratio curve of cells exposed to 0 µM H2O2 has to remain flat.

Save or format the acquired image as a 2-channel stack. Channel 1 must correspond to the near-420-nm acquisition and channel 2 to the near-500-nm acquisition. (Fig. 2A).

Fig. 2

Image formating, macro opening and editing. A: The image is a 2-channel stack containing as many frames as time points (here: 131). B: Screenshot of the menu to open the macro. C: Screenshot of FIJI macro editor. The language “IJ1 Macro” has to be selected to edit or run the macro.

For quantification, images must not be compressed and a high dynamic range (at least 12 bits) is preferable. Any image format supported by FIJI is acceptable, such as tiff or czi.

Required softwares for image quantification and data analysis

FIJI or ImageJ: Use either ImageJ (imagej.nih.gov/) or FIJI (https://fiji.sc/) to run the macro.

R and MATLAB: Most of statistical analyses are performed using R (cran.r-project.org) through the Rstudio interface (https://rstudio.com/). Parameters calculation is performed using MATLAB (https://uk.mathworks.com/products/matlab.html), but free analysis softwares like Scilab (https://www.scilab.org/) or spreadsheet programs like LibreOffice Calc can also be used.

Further details regarding softwares are described in the Materials and Methods section.

Data extraction from raw acquired images

We developed a macro, based on the image analysis protocol established by the creators of the HyPer probe [20,21]. Briefly, the raw 2-channel image is split into two images (STEP 1), which are then processed to improve their signal-to-noise ratio as well as dynamic range (STEP 2). The ratio of the two denoised images is then calculated. For each single cell, a rectangle containing the cell on each slice is duplicated on each image (STEP 3). The cell is then segmented, and both fluorescence intensities are quantified within the cell for each time step, as well as ratio values (STEP 4). The macro eventually exports (i) a summary table containing the time values in minutes and the fluorescence intensities for each analyzed cell and (ii) a report including all the choices made by the user (STEP 5). At every step, the macro saves regions of interest (ROIs) and/or intermediates images. These files are referred as “control data” in the following text.

STEP 0: Macro installation

Download a recent version of FIJI or ImageJ software.

Save the provided macro in the file of your choice.

Open the macro using the Macros → Edit menu and select the IJ1 macro Language (Fig. 2B). Alternatively, the macro can simply be run using the Macros → Run… menu.

If desired, edit the macro parameters. Then launch it by clicking run (Fig. 2C).

STEP 1: Image import and working folder selection

Click “OK” on the 1st dialog box and choose the folder to save the control data (ROIs positions, intermediates images…), later referred as “working folder” (Fig. 3A).

Fig. 3

Importing images and choosing the export folder using the macro. Dialog boxes enable the user to choose the folder for saving the control data (A) and to open the raw image (B). If the Bio-Format importer is used, there is no need for further action (B). The image is automatically split into two channels, which are renamed “C1.tif” and “C2.tif”. They are displayed using the HiLo LUT (C).

Click “OK” on the 2nd dialog box and open the image to analyze. If the bioformat importer is used, choose “view stack with Hyperstack” and do not tick any box (Fig. 3B). The image is automatically split by the macro into two channels, which are renamed “C1.tif” and “C2.tif” respectively. Both images are automatically displayed using the HiLo LUT to easily visualize the zero and saturated pixels (colored in blue and red respectively) (Fig. 3C).

STEP 2: Image processing

Image processing is a crucial step as it improves image quality and signal-to-noise ratio, and makes cell segmentation easier.

Observe the images and decide whether or not to use additive binning before clicking “OK” on the dialog box (Fig. 4A). The analyse → Histogram menu enables to check the dynamic range. Additive binning increases the signal dynamic range but decreases resolution (Figure S1).

Fig. 4

Image processing options. A: Before clicking OK on the 1st dialog box, the user is free to check image dynamics and resolution to decide whether using binning or not. B: A dialog box enables to choose among several options of image processing. C: To subtract background, the user defines a background-containing rectangle on the channel of his choice. D: The quantification will be performed on denoised images. Potentially blurred images are used for cell segmentation. E: The ratio image is calculated as a 32 bit stack from both denoised images.

Select the desired image processing options (Fig. 4B):

Ratio to be calculated: for instance, choose for a HyPer-based probe.

Post-acquisition additive 2 × 2 binning: useful to improve signal dynamic range for low-fluorescence cells with sufficient resolution (Figure S1). If this box is checked, both binned images are saved in the working folder as C1_binning.tif and C2_binning.tif.

Background subtraction (recommended): subtraction of the mean gray value of a cell-devoid rectangle from the gray value of each pixel (on the binned images, if so) on the corresponding slice. Background mean gray values as a function of slice number are saved in the working folder as C1_background.csv and C2_background.csv, as well as the images from which the background noise has been subtracted (as C1_background_substraction.tif and C2_background_substraction.tif). The position of the background rectangle is saved as “background_rectangle.zip”.

Median filter: useful to decrease the electronic noise of images acquired with a confocal microscope (Figure S2). If the box is checked, please adjust the filter radius (default: 2 pixels). The denoised images are saved as C1_med.tif and C2_med.tif respectively.

Gaussian filter: useful to smooth punctiform signal of mitochondrial probes to facilitate cell segmentation (Figure S2). If this box is checked, please adjust the filter radius (default: 4 pixels). The smoothed images are saved in the working folder as C1_gauss.tif and C2_gauss.tif respectively.

After this step, only denoised images (for quantification) and potentially blurred images (for segmentation) are opened (Fig. 4D). The fluorescence ratio of both denoised images is then calculated as a 32-bit image in order to maximize its dynamic range (Fig. 4E). This image is exported in the working folder as ratio.tif.

STEP 3: Creation of regions of interest with a single cell inside

Count the cells to be analyzed and fill the dialog box accordingly (Fig. 5A and B).

Fig. 5

Selection of individual cell-containing regions of interest. A: The user can adjust the contrast and scroll through the frames to count the cells to be analyzed before clicking “OK” on this dialog box. B: The user indicates how many cells he wants to analyze (default: 1). C: The user can choose on which channel to draw the ROIs. D: For each cell, the user draws a rectangle before clicking “OK” on this dialog box. E: The numbered ROIs are displayed on the active image as overlays. The ROI manager (right) displays the ROI position.

Choose the channel on which you want to draw the ROI (Fig. 5C).

For each cell, draw a rectangular ROI that contains the cell on all frames of the stack (Fig. 5D) before clicking “OK” on the dialog box (Fig. 5E).

The ROIs are displayed on the image as numbered overlays. Their positions and numbering are saved in the working folder as “analyzed_cells.zip”, it enables to check afterwards the numbering of each cell Fig. 5E.

Each ROI is automatically duplicated on both denoised acquisitions, on both blurred acquisitions and on the ratio image.

STEP 4: Cell segmentation

Identifying cells as separate objects different from the image background (cell segmentation) is a determinant step for successful further analysis, especially since changes in cells shapes or cells movements may occur during the acquisition.

Next steps are performed sequentially on each duplicated ROI.

Threshold the 1st channel (i.e. choose a gray value above which the signal is considered as cell and below which it is considered as background) (Fig. 6A).

Fig. 6

Cell segmentation using the macro. Both channels are thresholded (A, B, D, E) and binarized (C, F). Then, according to the user's choice, cells are detected either on one binary image or on the sum of both (G).

Click “OK” on the dialog box (Fig. 6B) to binarize the 1st channel (Fig. 6C).

Repeat these two steps on the 2nd channel (Fig. 6D, E, F)

Choose the binary image that will be used for segmentation: it can be either one of the two binary images or the sum of both. Choose whether you want to analyze edge-touching cells (Fig. 6G).

Summing both binary images usually optimizes cell segmentation (Figure S 3A). However, better results may be obtained with one of the two binary images only, e.g. if the fluorescence of the probe on a whole channel is insufficient to achieve a clear segmentation. It is recommended to exclude edge-touching cells, to ensure that the whole cell is considered for all measurements.

The default parameters for cell detection are:

Size (µm²): > 100 µm².

Circularity: 0.00 – 1.00 (1 corresponding to a perfectly circular particle).

No hole inclusion (to exclude the nucleus region).

Depending on cell size, cell morphology and probe localization, these parameters can be modified in the script of the macro before running it (Figure S 3B). We advise to perform segmentation tests by varying size and circularity parameters to check whether these settings are accurate for a proper segmentation of the user's cells.

At the end of this step, for each cell, both binary images are automatically saved in the working folder as i_C1_binary.tif and i_C2_binary.tif, with i being the number of the considered cell. The mask, i.e. the binary image corresponding to the segmented cell after size and circularity filtering, is saved as “i_mask.tif”, to check cell segmentation if desired.

STEP 5: Fluorescence intensity quantification

Choose the folder where the quantification results will be saved (Fig. 7A).

Fig. 7

Quantification results. A: A dialog box enables the user to choose the file to save the quantification results. B: The quantification results are a 5-column matrix containing the fluorescence intensities of both channels and on ratio image for each frame of each cell (3rd to 5th columns). Time values in minutes are extracted from the image metadata (2nd column). C: Report of all image processing and segmentation parameters that have been chosen for each cell.

The macro automatically performs the quantification.

The images are automatically closed.

These results are stored in a 5-column matrix (Fig. 7B):

Column 1 (entitled “cell_number”): number of the cell

Column 2 (entitled “time”): acquisition time extracted from the image metadata

Column 3 (entitled “i_C1”): mean fluorescence intensities on the 1st channel

Column 4 (entitled “i_C2”): mean fluorescence intensities on the 2nd channel

Column 5 (entitled “ratio”): mean fluorescence intensities on the ratio image

This stacked data format allows to gather data from time-lapses with different time bases and is adequate for analysis using R software.

A report of all image processing and segmentation parameters chosen for each cell is automatically exported in the result folder (Fig. 7C).

The mean fluorescence intensity quantification is performed as follows (Figure S 3B): each image of the mask stack is divided by 255 so that the pixel value within the cell equals 1 (=255/255) and background pixels equal 0 (=0/255). Each image of the three stacks to be analyzed is divided by the corresponding image of the mask stack. This operation does not alter the signal inside the segmented cell but assigns “not a number” (NaN) to each background pixel. The mean fluorescence intensity over the ROI is calculated for each image of the three stacks. This value corresponds to the mean gray value of each image.

Data analysis

STEP 0: R configuration

Download a recent version of R (www.r-project.org) and Rstudio (https://rstudio.com/).

Save the provided R script in the folder of your choice.

Open the script using the File → open menu (Fig. 8A to C).

Fig. 8

Opening the script and installing the packages with Rstudio. A: Rstudio interface overview. The left window enables to script directly into the R console. The “environment” tab (top right window) gives an overview of all the objects (here, the a, b and c variables). Plots are displayed in the bottom right window. B: The File → Open file menu enables to open any R script in a new window above the console. C: Overview of the script window. A single line can be run using the Ctrl + enter shortcut. The script can be run either by clicking on “run” or by using the Ctrl + Enter shortcut. D: The “Packages” tab enables to install and load the required packages.

Install the required packages: {tidyr} [26], {xlsx} [27], {reshape2} [28], {ggplot2} [29], {ggpubr} [30], {gridExtra} [31], {FactoMineR} [32] and {factoextra} [33]. This step can be performed either by using the Packages tab (Fig. 8D) or the install.packages function. The installed packages will stay installed for your next working session using R.

Load the required packages by running the 1st lines of the script (Fig. 8C), for example by using the Ctrl + Enter shortcut. This step needs to be performed at each new working session.

STEP 1: Qualitative variable implementation and data compilation

Qualitative variables are required for data analysis. They must be implemented in each raw data file exported by the macro. Here, we added the following variables: H2O2 concentration, date, replicate number, cell passage number, image number, and a unique identifier for each cell (Table 2). Some variables like H2O2 concentration are expected to be important parameters for the analysis. Others like the replicate number or the cell passage number can be added to examine their potential implication in the final results. Based on these variables, data can easily be filtered and analyzed accordingly. The cell identifier is useful to retrieve the initial image. Then, these completed datasets must be compiled together.

Table 2

Meaning of the dataset columns.

	Column name	Variable's signification
Variables already present in the raw tables exported by FIJI	cell_number	Number of the cell in the image
	time	Time from the start of acquisition in minutes (extracted from the image metadata)
	i_C1	Mean cell fluorescence intensity on the 1st channel
	i_C2	Mean cell fluorescence intensity on the 2nd channel
	ratio	Mean cell fluorescence intensity on the ratio image
Additional qualitative variables	C_H2O2	H2O2 concentration applied on cells (levels: 0, 10, 40 and 75 µM)
	date	Date of the experiment (yyyymmdd format)
	replicate	Replicate number: e.g. 2 if the experiment is done for the 2nd time on the same day
	passage	Number of passages after stable transfection (see M&M)
	induction_time	Probe expression time (from doxycycline addition)
	image_num	Image number for a given date
	id_cell	Unique cell identifier, obtained by concatenating the date, the image number for a given date and the number of the cell in this image, separated by "_"

We provide a method to perform these actions using R (Fig. 9).

Fig. 9

Data completing and compiling prior to statistical analysis. A: Example of R script for raw data importing, completing with qualitative variables and renaming. B: Visualization of the 1st raw data file. C: Visualization of the 1st data file after qualitative variables implementation. D: Example of R script for data compiling and exporting. E: Visualization of the complete dataset.

1. Open a new R script using the Ctrl + Shift + n shortcut.

2. Set as working directory the folder where the raw data are saved using the setwd function.

setwd(“folder_where_your_data_are_saved”) # Sets the proper working directory

3. Import the 1st raw data file using the read.table function (Fig. 9A and B).

rawdata<-read.table(“data1.csv”, sep=”,”) # Imports the dataset and names it rawdata

4. Add the 1st required qualitative variable as shown on the following example; it adds a column entitled C_H2O2 which is filled of “10”. Quotation marks around “10” allow the software to consider the variable “C_H2O2” as a multi-level factor rather than a quantitative variable.

rawdata [“C_H2O2”] <- “10” # Adds a full column of “10” entitled C_H2O2

5. Repeat step 4 to add any other desired qualitative variable (Fig. 9A).

6. One of these qualitative variables must be a unique identifier for each cell (Id_cell column, Fig. 9C and E). This allows (i) to represent the results as curves using the {ggplot2} package, and (ii) to easily retrieve the cell and experiment associated to a particular curve.

7. Rename this 1st completed data frame using the following script line, which creates a new data frame called data1 that equals the rawdata data frame (Fig. 9A):

data1 <- rawdata

8. Perform steps 3 to 7 with the other raw data files, renaming each completed data file with a new name (for example data2, data3, …, datak).

9. Compile these completed data files using the rbind function as follows (Fig. 9D and E):

dataset <- rbind (data1, data2, …,datak)

10. Export the complete dataset in the chosen folder using the write.table function (Fig. 9D):

setwd (“the pathway of the chosen working directory”)

write.table (dataset, file = ”TheNameYouWant.txt”)

STEP 2: Data visualization and processing

Even though the provided script is customized to our dataset, with minor adjustments it enables the complete analysis of any dataset.

Go back to the R script.

Load the dataset by running the script lines entitled “raw data loading”.

To visualize data, the provided script (Fig. 10A) uses the {ggplot2} package to display the ratio curves colored according to the experiment date and replicate. H2O2 addition time is indicated by vertical bars. The facet_grid function is used to split the response curves as a function of the 1st qualitative variable: H2O2 concentration (Fig. 10D). These settings can directly be modified in the script (see https://ggplot2.tidyverse.org/ for more details).

Fig. 10

Data visualization and processing. A: Script to plot the ratio curves facetted by H2O2 concentrations applied on cells (0, 10, 40 or 75 µM H2O2). B: Script enabling to get the same time base for all curves. C: Spread data are exported into an excel file. D: Plot of the raw data ratio curves facetted by H2O2 concentrations. E: Plot of the interpolated data using the same settings.

If the time base is not the same for all cells, run the part of the script enabling to interpolate your data (Fig. 10B) using the spline function from the {stats} package (R Core Team 2019). The number of knots can be set to the number of data points (=131 in the provided example).

Check that the interpolation does not noticeably alter the shape of the curves by plotting the interpolated ratio curves (Fig. 10E).

Spread this interpolated dataset using the spread function of the {tidyr} package. This step changes the dataset from a long table, in which all the cells are one above the other (Fig. 9E), to a wide table with one cell per column (Fig. 10C).

Export the spread dataset as an excel file in the working directory, using the write.xlsx ({xlsx}) function (Fig. 10C).

STEP 3: Choosing kinetic parameters

Save the provided MATLAB script in the same folder as the excel file generated at previous step.

Open the script in MATLAB.

Adjust the following variables at the beginning of the script (Fig. 11A):

Fig. 11

Kinetic parameters calculation using the MATLAB script. A: Overview of the beginning of the MATLAB script. Strings preceded by “%” (in green) are comments. The string underlined in green has to be the name of the excel file containing your spread data. Numbers framed in red must be adjusted to your own number of time points. Numbers framed in purple enable to adjust the smoothing level. B: Parameters calculated using the MATLAB script are exported as an excel spreadsheet with one line per cell. C: Overview of the calculated parameters for one curve.

The imported file name if required

The number of time points before H2O2 addition

The smoothing parameters if required

Indeed, the curves are smoothed before the kinetic parameters are calculated. In our case, the slope of the curve increased sharply upon addition of H2O2. We therefore chose to smooth separately the two parts of each curve (before and after H2O2 addition) to take this into account. For this reason, the number of time points before H2O2 addition is required.

Run the script by clicking “run” (Fig. 11A).

The following parameters, summarized in Fig. 11C, are calculated on smoothed curves:

The maximum ratio, rmax

The time of the maximum ratio in minutes, tmax

The slope at the maximal inflexion point in minutes−1, noted infl

The time of the maximal inflexion point in minutes, tinfl

The ratio value at the maximal inflexion point, noted rinfl

The lag time, in minutes, noted tlag, is estimated by extending the tangent at tinfl to the initial baseline (y=rinit, where rinit is the median ratio before H2O2 addition).

Except the maximal ratio, these parameters cannot be estimated for flat curves, as observed in absence of treatment. The tmax, infl, tinfl, rinfl and tlag parameters are calculated for curves whose maximum ratio is arbitrarily at least two-fold the initial ratio and reached within 30 min. This filter can directly be customized in the script.

The calculated parameters are automatically saved as an excel spreadsheet in the working folder (Fig. 11B).

From the R script, run the script lines entitled “Analysis of the parameters calculated by MATLAB”. This part of the script merges the kinetic parameters with the complete dataset and plots each parameter as a function of qualitative variable 1 (H2O2 concentration in our example).

As our samples are small (less than 30 individuals per H2O2 concentration tested), we chose to superimpose boxplots and dotplots to visualize the parameters as a function of H2O2 concentration. The boxplots give a simple visualization of the variables by showing median and extreme values. The overlapping dotplots precisely show the data distribution since each cell is shown. The p-values associated with the Kruskal-Wallis multiple comparison test are calculated and displayed on each graph using the stat_compare_means ({ggpubr}) function. The null hypothesis tested by this nonparametric test is H0 = ”All H2O2 concentrations give identical results”. We rejected this hypothesis (i.e. we considered that there was at least one significant difference) as soon as the associated p-value was below 0.05. Mann Whitney pairwise comparison tests were also performed to compare H2O2 concentrations two by two. The plots are presented in Fig. 12.

Fig. 12

Kinetic parameters of ratio curve analysis after cell exposure to 0, 10, 40 and 75 µM H2O2 (Kruskal-Wallis and Mann-Whitney tests).

We observed that:

The maximal ratio is reached later for 10 µM H2O2 than for 40 and 75 µM H2O2.

The maximal ratio is minimal for cells treated with 0 µM H2O2. It is lower for 10 µM than for 40 and 75 µM H2O2.

The inflexion point is reached later for 10 than for 40 and 75 µM H2O2. Surprisingly, it is reached later for 75 than for 40 µM H2O2 (but with a higher p-value).

The slope at the inflexion point is lower for 10 µM than for 40 and 75 µM H2O2.

The lag time is longer for 10 than for 40 and 75 µM H2O2.

All these observations strongly suggest that the lower the H2O2 concentration, the slower the ratio increases and the lower the maximal ratio. These conclusions are valid for low H2O2 concentrations, but we do not observe any significant difference between the kinetic parameters of the curves obtained with 40 and 75 µM H2O2.

Overall, these results highlight that, in our cell model, HyPer ratio increase is lower and slower for cells treated with 10 µM H2O2 than for cells treated with 40 and 75 µM. Consistency of these results with those presented in the following sections is shown in Table III (see Supplementary).

STEP 4: End-point ratio analysis

Run the R script section entitled “Analysis of the final median ratios”. The median ratio between 60 and 65 min is calculated for each curve and this new variable, named r_end, is added to the dataset (Fig. 11C). Final ratios are then plotted as a function of H2O2 concentration (Fig. 13).

Fig. 13

Final ratios analysis after cells exposure to 0, 10, 40 and 75 µM H2O2 during  ≈  1 h (Kruskal-Wallis and Mann-Whitney tests).

We observed that the final ratio increased with H2O2 concentration, reaching a plateau at 40 µM H2O2. Therefore, in our cell model, HyPer ratio returned to its basal state more efficiently for cells treated with 10 µM H2O2 than for cells treated with 40 and 75 µM H2O2. These results are consistent with results of step 3 (Fig. 12).

STEP 5: Assigning a global score to each curve using functional principal components analysis

In this section, each curve is replaced by a single score using functional Principal Components Analysis (fPCA). This approach enables to preserve almost the entire variability of the initial dataset, allowing a fast and accurate comparison of the responses obtained under the different experimental conditions. Moreover, in contrast to the previous kinetic parameters, this approach can be used for all curves including flat shaped curves. According to fPCA principle, each curve is decomposed into a sum of components called harmonics representing their main variation modes.

Go to the R script section entitled “functional PCA”.

Load both the stacked and spread datasets by running the appropriate script lines.

Convert the spread dataset into a functional data object by executing the corresponding section of the script (Fig. 14A). A spline basis is created, whose parameters need be adjusted: rangeval must be set to the time range and nbasis to the number of spline functions of order norder that are summed to fit the curves. The higher the values of these two parameters, the lower the smoothing will be (Fig. 14A).

Fig. 14

Functional principal components analysis of the dataset. A: Script for the functional data object generation from the spread dataset. Adjust the numbers framed in red according to your dataset and desired smoothing level. B: Visualization of the dataset 1st two principal variation modes. The 1st harmonic (left) is a pharmacokinetic-type function whose score quantifies the amplitude of the variation and the overall height. This is the overwhelming mode of variation in our dataset as it accounts for 99% of its variability. The score on the second harmonic (right) indicates whether the curve tends to have a high maximum ratio and a low final ratio. This mode of variation only accounts for 0.4% of the variability in our dataset, so we chose to neglect it. C: Curves colorization according to their score on harmonic 1. The lower the score of a curve, the closer its color is to light orange. The higher the score of a curve, the closer its color is to dark blue. The higher the score of a curve, the higher it is and the greater its amplitude. D: Curves score on the 1st harmonic of fPCA after cells exposure to 0, 10, 40 and 75 µM H2O2 (Kruskal-Wallis and Mann-Whitney tests).

Perform the functional PCA and plot the types of variation represented by the 1st two harmonics (Fig. 14B).

The 1st harmonic, namely the function type that is able to best describe the set of curves, is a pharmacokinetics-like function with a lag time followed by a growing phase and by a slow decrease (Fig. 14B, left). This harmonic accounts for 99% of our dataset variability. The 2nd harmonic only accounts for 0.4% of the variability and can thus be neglected (Fig. 14B).

Plot all ratio curves colored according to this score (Fig. 14C).

The plot unambiguously shows that the higher the score of a curve on the first harmonic, the higher the overall ratio: cells with an overall high cytoplasmic H2O2 concentration have a high score and cells with an overall low cytoplasmic H2O2 concentration, e.g. cells treated with 0 µM H2O2, have a low score.

Plot the score on the 1st harmonic of fPCA as a function of the qualitative variable of your choice (the H2O2 concentration in our example) (Fig. 14D).

We noticed that the score increases with H2O2 concentration, except for 40 and 75 µM H2O2 where no significant difference is observed. The score on the 1st harmonic of the fPCA is thus an excellent way to summarize a complex dataset.

STEP 6: Principal components analysis of all parameters

To get an overview of the above kinetic parameters and scores, we performed a Principal Components Analysis (PCA) on all these variables. PCA is a procedure enabling to simplify a dataset by replacing an important number of variables by a limited number of new variables. These new variables are called axes and are computed to account for the variability of the original dataset [34].

Go to the R script section entitled “PCA with every qualitative parameter”.

The individuals factor maps (i.e. the projection of each cell in the plane formed by the 1st two axes) and the correlation circle (i.e. the projection of the variables in the same plane) are automatically displayed.

The {factoextra} package enables to plot these data in a more esthetic and customizable way. Individuals and variables were colored according to the H2O2 concentration (Fig. 15A), and to the quality of their representation (cos2), respectively (Fig. 15B). The closer to 1 the cos2 value is (green vectors), the better the corresponding variable is represented in the plane formed by the 1st two axes of the PCA. Variables with a cos2 close to 0 (red curves) are poorly represented.

Fig. 15

PCA on all quantitative variables. A: Individuals projection in the plane formed by the first 2 axes of the PCA with coloration according to the H2O2 concentration. B: Correlation circle with vectors coloration according to the quality of representation of variables. The shorter the vector representing a variable and the closer to 0 its cos2, the less well it is represented. Here, all variables are satisfactorily represented except for rinfl. C: Projection of the entire dataset (variables and individuals) in the plane formed by the first two axes of the PCA. Individuals are colored according to the H2O2 concentration. The center of gravity of each point cloud is presented as a larger point. Two independent (orthogonal) groups of variables collinear to each other were identified. Variables in the set 1 (blue vectors) describe the ratio increase kinetics, while variables in the set 2 (purple vectors) describe the overall ratio height.

Individuals and variables factor maps can also be superimposed (Fig. 15C).

From these plots, we note that the 1st two axes account for 78.39% of the total variability. The projection of the variables onto the plane formed by these two axes is thus a good image of the complete dataset Fig. 15. Excepting the ratio at the inflexion point, all the variables are well represented on this projection (Fig. 15C). The ratio at the inflexion point is thus set aside in the following interpretations.

The variable projection highlights two independent (orthogonal) groups of variables being collinear to each other (Fig. 15D). These sets of variables are complementary:

Set 1: kinetics of the ratio increase (blue vectors). The time the inflexion point is reached, the lag time, and the time the maximum is reached, are highly correlated. Expectedly, this set of variables is highly anticorrelated with the slope at the inflexion point, corresponding to the speed of ratio increase.

Set 2: overall ratio height (purple vectors). The maximal ratio, the final ratio and the score on the 1st harmonic of the fPCA are highly correlated. This set of variables is not correlated with the one described above. Both sets of variables are thus complementary with each other.

The joint observation of the circle of correlations and the projection of individuals leads to the following conclusions (Fig. 15D):

The cells treated by 0 µM H2O2 are grouped together along the set 1 of variables, because tinfl, tlag, tmax and infl were not calculated for these cells. Consistently, they exhibit low overall ratio heights.

The cells treated by 10 µM H2O2 on the one hand and by 40 or 75 µM H2O2 on the other hand are well discriminated along both sets of variables. Consistently, 10 µM H2O2-treated cells have both a slower (or more delayed) ratio increase and a lower overall ratio than 40 or 75 µM H2O2-treated cells.

Overall, these analyses show that the response is faster and higher as the H2O2 concentration increases. Individual cells treated with 0, 10, 40 or 75 µM H2O2 can be efficiently discriminated, using various statistical parameters as shown in Fig. 15D.

Method validation

Related to the widespread utilization of GFP- and YFP-based biosensors in the scientific community, methods have been proposed by several authors with the common objective of rendering microscopy data analysis less fastidious and time-consuming [18,19]. In this study, we describe a procedure, which is simple to use, reproducible and robust, to analyze microscopy data obtained from ratiometric biosensors. Considering cell morphology changes during the acquisition time is a major advantage of our provided macro. Following image acquisition and quantification, we further developed a complete set of tools for analysis, including the automatic determination of key “diagnostic” kinetics parameters (representative time, slope and ratio values, see STEP3) for appropriate curve mathematical definition. Functional PCA was introduced to simplify each curve by a single score yet retaining variability. This dataset can be analyzed using parametric or non-parametric statistical tests, depending on data size and distribution. Final PCA on all quantitative variables allows a global overview of the conclusions, through a clear and visual presentation of the results.

As previously mentioned, our free method mainly uses open source softwares (FIJI and R) and associated packages, allowing easy access to all. Although some parameters were calculated using the MATLAB software, calculation can be performed using any curve analysis free software like Scilab or spreadsheet program like LibreOffice Calc.. Moreover, scripts or macros are easily editable and intelligible thanks to many embedded commentaries. Modifications for particular needs are thus straightforward.

Through this presentation, the user is guided step by step. Moreover, since the process is semi-automated, the user can check and adjust the parameters at key steps, avoiding errors in data processing. In contrast, this could not be easily performed when using fully automated methods. Further, checking intermediary results or retrieving analysis parameters is made possible thanks to results saving at each step of the process.

To evaluate the performance of our procedure in analyzing data and reaching biological conclusions, we used an experimental scenario involving cultured epithelial rabbit kidney RK13 cells expressing the ratiometric HyPer probe in the cytoplasm. Cells were treated with increasing H2O2 concentrations (0, 10, 40 and 75 µM) and time-lapse fluorescence microscopy data were collected. Substantial changes occurred in the measured parameters when increasing the exogenous H2O2 concentration. Statistical significance was established for all parameters between 0 and 10 µM, and between 10 and 40 µM H2O2 (data are summarized in Table III, see Supplementary), showing an efficient discrimination between cells treated with various H2O2 concentrations ranging from 0 to 40 µM by our method. This demonstration in a standard experiment involving dose-dependent oxidation of the sensors validates that the set of tools is capable of reflecting biological changes, accurately and quantitatively.

Noticeably, no statistically significant difference was found between 40 and 75 µM H2O2 for most of the parameters. It is possible that HyPer fluorescence ratio evolution within this range of H2O2 concentrations is not linear anymore, consistently with previous studies that described that probe signal is linear within the 5–25 µM range in E. coli and eukaryotic cells [24] or within the 0–10 µM range in human cells [35]. It is also possible that the actual intracellular H2O2 concentration might not linearly correlate with increasing exogenous H2O2 concentration, as intracellular antioxidant proteins such as peroxiredoxins and catalases may act to metabolize H2O2 when it is above a certain threshold [36]. It is very likely that the absence of statistically significant difference between some of the parameters within the 40–75 µM range originates in cell physiology or in sensors oxidation dynamics in vivo.

One potential discrepancy also appeared. We identified that tinfl was higher for 75 than for 40 µM H2O2 (otherwise expected to remain equal or to decrease as H2O2 concentration increases), in a statistically significant manner (see Fig. 12). This discrepancy may very well be due to the parameters calculation imprecision, as the curves have a step of 0.5 min and the inflection point generally appears a few minutes after H2O2 addition. Decreasing time between time points during acquisition, if technically feasible, should likely resolve this point.

As a conclusion, we provide here a step-by-step semi-automated method to analyze microscopy data obtained from ratiometric biosensors. This method proved to be effective and accurate, and provides robust statistical analysis. Results are shown as a visual diagram. Its applications to various experimental contexts and within different organelles should be helpful to many researchers interested in the utilization of ratiometric sensors.

Construction of the HyPer-expressing cell line

Construction of HyPer cyto NES expression vector

The original pHyPer-cyto vector encoding HyPer was purchased from Evrogen [24]. Nucleus-excluded HyPer was constructed by amplifying the HyPer sequence using primers 719 and 720 (Table V, Supplementary) introducing a Nuclear Export Signal (NES). The PCR product, fusing NES in frame to the 3′ end of HyPer, was then cloned in pCMV/myc/cyto vector (Life Technologies) at PstI/NotI sites, leading to pCMV HyPer cyto NES plasmid. HyPer cyto NES encoding sequence was then amplified by PCR using primers HyPer cyto NES for and HyPer cyto NES rev (Table V, Supplementary), and cloned at the SalI/NdeI sites in the pTRE3G plasmid (Clontech laboratories) for inducible gene expression.

Construction and culture of the RK13 HyPer cyto cell line

We used the rabbit kidney epithelial lineage RK13 [37]. We first established a stable Tet3G expressing cell line by transfection with the pTetOn3G vector (Clontech) using Lipofectamine 2000 (ThermoFisher) and selection with 500 µg/mL G418 (Sigma-Aldrich). About three million RK13 pTetOn3G cells were then stably cotransfected with 2 µg pTRE3G HyPer cyto NES (not encoding any antibiotic resistance) and 100 ng of a selection plasmid encoding puromycin resistance by electroporation using the Nucleofector system (Lonza). The transfection was performed in Opti-MEM medium with fetal bovine serum (FBS) (Eurobio) but without antibiotics using the pre-recorded program optimized for MDCK cells. After selection of transfected cells with 1 µg.mL−1 puromycin (Sigma-Aldrich), a clonal cell line was isolated from the population by limiting dilution to maximize the level of probe expression as well as the proportion of probe-expressing cells. Aliquots of cells were frozen on the first passage after checking for the absence of mycoplasma using the MycoAlert mycoplasma detection kit (Lonza).

The cells were passaged to the fifth once a week in Opti-MEM medium (Gibco) supplemented with 10% FBS, 100 I.U.mL−1 penicillin, 100 µg.mL−1 streptomycin (PS) (Eurobio), 1 µg.mL−1 puromycin and 500 µg.mL−1 G418. Cells were kept at 37 °C under a water-saturated atmosphere containing 5% CO2.

Fluorescence microscopy

Two days before acquisition, 1000 cells were seeded on the coverslip of a 35 mm glass bottom dish (MatTek) in 500 µL OptiMEM supplemented with 10% FBS and PS. Cells were incubated until they adhere to the glass (90 min at 37 °C). Then, the medium was replaced by 2 mL FluoroBrite DMEM (Gibco) with 10% FBS, 1 I.U..mL−1 PS, 2 mM l-glutamine, 25 mM HEPES pH 7.3 and 1 µg.mL−1 doxycycline and the cells were incubated for 48 h ± 2 h at 37 °C before microscopy acquisitions.

Timelapse images were acquired using a Zeiss AxioVert.A1 wide field fluorescence microscope (LD Plan-Neofluar 40x magnification/0.6 Korr M27 dry objective) using LED illumination at 423/44 and 469/38 nm (colibri-7). Images were acquired in 12 bit using an AxioCam MR R3 imaging device (Zeiss). Two-channel time-lapse images were acquired by exciting cells sequentially at 423 and at 469 nm. Acquisition at 423 nm was performed first for each frame. For both channels, light emitted between 500 and 550 nm was collected.

The parameters used for the acquisitions must be the same for all the experiments compared with each other. To set the optimal conditions, the objective is to maximize the signal on both channels while maintaining a low phototoxicity due to the excitation itself. In particular, it is important to check that the ratio remains constant under the chosen acquisition conditions in the absence of exogenous redox treatment. The authors also recommend having an initially similar fluorescence intensity on both channels (see [20] for further details). The excitation parameters we used are summarized in Table IV (see Supplementary) for information purposes. We performed 65-minute acquisitions with 2 images per minute. After 5-min acquisition, 1 mL culture medium with a given H2O2 concentration was added to the cells. Culture medium alone was added to the control cells. The .czi stacks we obtained were directly processed and analyzed using the macro presented in the present paper.

Softwares used for macro development and subsequent data analysis

- Macro development

We used the FIJI distribution of ImageJ to develop the macro, that can be executed either with the simple ImageJ or with FIJI. FIJI is an open source project based on the ImageJ software [38], [39], [40].

- R, RStudio and associated packages

R is a powerful language for statistical computing and graphics based on an open source environment highly compatible with every operating system [41]. The Rstudio interface was used for programming [42].

The following packages were used:

{tidyr} for data spreading [26]

{xlsx} for excel spreadsheet exporting [27]

{reshape2} for data melting [28]

[43] for functional PCA performing [44]

{ggplot2} for esthetic plotting [29]

{ggpubr} for p-value associated with nonparametric tests displaying on plots [30]

{gridExtra} for graphs faceting [31]

{FactoMineR} for PCA performing [32]

{factoextra} for PCA results visualization [33]

- MATLAB

MATLAB is a particularly suitable language for matrices manipulation. To write the script for parameters calculation, the MATLAB version 9.5.0.1033004 [45] was used.

Subject Area	Bioinformatics
More specific subject area	Quantitative and statistical analyses of data from ratiometric fluorescent biosensors
Method name	PiQSARS
Name and reference of original method	Not applicable
Resource availability	https://fiji.sc/https://rstudio.com/https://uk.mathworks.com/products/matlab.html Macro and scripts are provided as separate documents as supplementary materials