Preparing the dataset (TM R-User Manual sections 1.1 to 1.5)

Before conducting any analysis, the TissueMiner automated workflow reads three configuration files that contain (1) user-defined regions of interest (ROI’s), (2) time offsets for movie synchronization, and (3) the rotation angle used to align the tissue to a standard orientation (Figure 1—figure supplement 1). 

Defining regions of interest (howto Video 1)

As cellular behaviors may be spatially patterned, one should have the ability to quantify and compare cell dynamics within different ROI’s. TissueMiner provides a Fiji macro (draw_n_get_ROIcoord.ijm) to manually define a set of ROI’s directly on one given image of the movie. This program manages several lists of ROI’s, which the user can create, modify and delete. These lists help maintain the consistency of ROI labels, which is essential for subsequent analysis (Video 1).

In addition, defining ROI’s of different shapes and following them backwards and forwards in time (Figure 1A–E’) is a useful method to visualize tissue deformations (Figure 1D–D’, Video 2). These ROI’s can be defined at any frame within the movie. Thus, it is even possible to specify a region based on morphological features that only arise late in the morphogentic process under study, which is true of wing veins for example (see Figure 1E–E’). ROI definition allows the user to define morphologically relevant regions of interest and compare the behavior of cells in the different regions.

Video 1.

Click here to view.^{(14M, mp4)}

HOWTO: drawing ROI’s.

DOI: http://dx.doi.org/10.7554/eLife.14334.006

Video 2.

Click here to view.^{(12M, mp4)}

Visualizing tissue deformation by using vertical stripes.

DOI: http://dx.doi.org/10.7554/eLife.14334.007

Open in a separate window

Figure 1.

Regions of interest are followed in time by browsing the cell lineages.

(A) Largest population of cells (purple) that remains visible throughout the entire time-lapse. Two cell rows in contact to margin cells were discarded as margin cells are usually not well segmented. (B) Largest blade cell population (green) that remains visible throughout the entire time-lapse. The blade region of interest (yellow line) was defined on the last frame of the time-lapse using a custom Fiji macro (https://github.com/mpicbg-scicomp/tissue_miner/blob/master/fiji_macros/). The underlying cell population was then subset using our lineage browser algorithm. (C) One can define veins and inter-vein regions of interest and apply the same algorithm as in (B). (D–D') Regularly spaced regions of interest automatically selected and followed over time to visualize tissue deformation. (E–E') Here, we make use of the lineage browser routine to trace back the vein positions at 15 hAPF, as they aren't visible yet at 15 hAPF. Scale bar 50 microns.

DOI: http://dx.doi.org/10.7554/eLife.14334.004

Figure 1—figure supplement 1.

Open in a separate window

Flow chart of TissueMiner.

Solid lines depict the three main steps to analyze epithelial morphogenesis within TissueMiner. Dashed lines indicate additional inputs to the automated workflow: red boxes represent required inputs and black boxes indicate optional inputs. Arabic numbers indicate the order in which the tools are described in the main text. Cumulative time of the movie must be listed in a text file called cumultimesec.txt and located along with the movie images. The snakemake automated workflow is described in Figure 7.

DOI: http://dx.doi.org/10.7554/eLife.14334.005

By default, TissueMiner generates two regions of interest – raw and whole_tissue – in order to select cell populations by name. The raw ROI corresponds to all segmented and tracked cells. However cells located at the tissue margin may move in and out of the field of view of the microscope lens. TissueMiner identifies the population of cells (whole_tissue) whose entire lineage lies within the field of view throughout the movie. To identify this population, we developed a filtering tool to discard in each movie frame margin cells located at the edge of the segmentation mask and one additional row of cells that contact the margin cells. The choice of discarding two rows of cells is motivated by the fact that segmentation quality drops near the margin. We iterate over all time points to ensure that we discard all cells moving in and out the field of view (see Materials and methods). User-defined ROI’s are also subjected to this filtering.

Visualizing cell area, cell shape and cell packing on the entire tissue (TM R-User Manual sections 2.2 and 2.6, Py-tutorial sections 2.1 to 2.3)

An important step in analyzing tissue morphogenesis is to quantify cell state properties over time. These properties include cell area, shape anisotropy and packing geometry. In this section, we demonstrate the analysis and visualization tools of TissueMiner by comparing how these state properties evolve during wing morphogenesis in vein and inter-vein regions.

Cell area and elongation (TM R-User Manual sections 2.2 – 2.5, Py-tutorial sections 2.1 – 2.2)

Morphogenesis is often characterized by changes in cell area and elongation. In the TissueMiner workflow, these properties are calculated from the original segmentation masks and stored in the database (Materials and methods). To visualize the evolution of the cell area pattern at the scale of the whole tissue, we map the area values of each individual cell to a gradient color scale (see Figure 2A–A’, Video 4). Each cell contour is filled with a color that corresponds to its area. Figure 2A’ shows the pattern of cell areas in the wing at the end pupal wing blade elongation. This visualization scheme reveals that cells in the proximal hinge and in wing veins have a smaller cross-sectional area (blue) at this time.

Open in a separate window

Figure 2.

Patterned cell state properties in the developing pupal wing of Drosophila.

(A–D') Cell state patterns at 22 hr and 31 hr after puparium formation (hAPF). (A–A') Color-coded cell area. (B–B') Color-coded cell elongation. The magnitude of cell elongation corresponds to the norm of the cell elongation nematic tensor. (C–C'') Coarse-grained pattern of cell elongation nematics and (C'') cell elongation nematics represented as bars on each individual cell. The wing was divided into adjacent square-grid elements of 33x33 microns in which cell elongation nematics were averaged. (D–D') Color-coded representation of the cell neighbor number. (E) Time evolution of the average cell area in different regions of interest: wing blade (Figure 1B), veins (Figure 1E), and inter-vein regions. (F) Time evolution of the average cell elongation magnitude in the blade, veins and inter-vein regions. Scale bar: 50 microns.

DOI: http://dx.doi.org/10.7554/eLife.14334.009

Cell elongation is characterized by a nematic tensor describing the axis and magnitude of the elongation (Aigouy et al., 2010). As with cell area, we map the magnitude of cell elongation to a color scale (Figure 2B–B’, Video 5). This fine-grained quantification of cell elongation highlights striking differences between inter-vein and vein cells. Inter-vein cells are more elongated than vein cells at 22hr after puparium formation (hAPF), but this pattern is reversed by 31hAPF. 

The color scale above reveals only the magnitude of the tensor. To visualize both the magnitude and direction of cell elongation, we represent the elongation nematic as a line whose length and angle correspond to the magnitude and angle of cell elongation, respectively. Nematics can also be averaged across multiple cells in a region in order to coarse-grain the patterns and highlight the main features (Figure 2C–C’’, Video 6). For example, the coarse-grained elongation nematics shown in Figure 2C, highlight the global alignment of cell elongation in the proximal-distal direction at 22hAPF.

Packing geometry (TM R-User Manual section 2.6, Py-tutorial section 2.3)

Cells in the wing become progressively more hexagonal during pupal wing morphogenesis (Classen et al., 2005). To visualize packing geometry, we map the neighbor number of each cell to a discrete color code (Figure 2D–D’, Video 7). This makes changes in packing geometry during morphogenesis immediately obvious (22and 31hAPF).

Plotting temporal evolution of average cell properties (TM R-User manual sections 3.3 to 3.6, Py-tutorial section 3)

The visualization tools described above effectively reveal detailed spatial patterns of cell properties. To highlight how average cell properties change over time, and to facilitate comparison between movies and ROI’s, TissueMiner also provides tools to create plots of average quantities as a function of time. In Figure 2E and Figure 2F, we compare the time evolution of the average cell area and the average cell elongation in movies of the 3 WT wings (blue, green, red) used in (Etournay et al., 2015). The plots in Figure 2 compare the time evolution of average cell elongation and area values for vein and inter-vein cells. We previously showed that average cell area in the wing blade decreases during morphogenesis, but that cell area decrease is balanced by cell divisions to maintain wing blade area. Quantifying average area values in vein and inter-vein ROI’s reveals that vein cells contract over a longer period of time than inter-vein cells, and thus have a smaller cross-sectional area at the end of morphogenesis (Figure 2F). As previously described, cells in the wing blade elongate and then relax their shapes during pupal wing morphogenesis (Etournay et al., 2015) (Figure 2E, blade part). Plotting elongation in vein and inter-vein ROI’s reveals that vein cells elongate more slowly and also relax their elongation more slowly than inter-vein cells. These differences suggest that vein and inter-vein cells have different mechanical properties.

Visualizing patterns of cell divisions (TM R-User Manual sections 2.7 – 2.9, Py-tutorial section 2.4)

Oriented tissue morphogenesis may reflect the number, orientation and spatio-temporal pattern of cell divisions. TissueMiner provides several tools to visualize these events. Overlaying color-coded generation number on a pupal wing movie reveals patterns of cell divisions as they occur (Video 8), and examining the last frame of the movie (Figure 3A) reveals the cumulative pattern of cell divisions. This analysis is largely consistent with the cell division timing inferred from classical BrdU pulse-chase experiments (Schubiger and Palka, 1987; Garcia-Bellido et al., 1994; Milan et al., 1996), but also reveals unexpected additional features. The pattern of cell divisions correlates with veins: most cells in the wing blade divide only once during pupal morphogenesis, whereas in some parts of inter-vein regions they divide twice. These include the cells lying adjacent to veins L3, L4 and L5, and the region posterior to L5. We estimate the median cell-cycle length between the first and second rounds of cell divisions to be (5.25 ± 1.50) hr.

Open in a separate window

Figure 3.

Visualization of cell generations and cell divisions.

(A) Color-coded pattern of cell generations. The wing cartoon on the bottom right shows the names of subregions that we analyze in panel B. Scale bar 50 microns. (B) Cell division rate in different regions of interest. To smooth fluctuations, these rates were averaged in discrete time intervals of one hour (TM R-User Manual, section 3.7). We further averaged these rates amongst the three wild-type wings. Error bars depict the standard deviation between wings. Cells divide earlier in veins L2 and L4 than in L3 and L5. Two maxima corresponding to two rounds of divisions are visible in inter-vein regions: interL2-L3, distInterL3-L4 and postL5. (C–C') A dividing cell with its unit nematic depicting the division orientation. Scale bar 10 microns. (D) Coarse-grained pattern of cell division orientation (grid size of 33x33 microns). Scale bar 50 microns.

DOI: http://dx.doi.org/10.7554/eLife.14334.015

To further investigate how cell divisions are patterned in the blade, we quantified the time evolution of cell division rates in each vein and inter-vein region (Figure 3B). This analysis reveals differences in the timing and numbers of cell divisions in these different ROI’s. Cells in veins L2 and L4 divide before those in L3 and L5. These divisions are followed by a second peak of division in the inter-vein regions distInterL3-L4, interL2-L3 and postL5 (see cartoon in Figure 3A).

To more easily visualize the spatio-temporal pattern of divisions in veins only, the powerful tools available in TissueMiner allow us to assign vein cells a color corresponding to the time at which they divide: blue for 16–18hAPF and red for 18–20 hAPF (see Video 9). This analysis reveals more detailed patterning in division timing. Cell divisions in vein regions that protrude ventrally (L2 and proximal L4), peak at the same time and earlier than those that protrude dorsally (L3, distal L4 and L5). Precise correlation of cell divisions with specific vein and inter-vein regions suggests that they are autonomously controlled by signaling associated with veins.

To measure the orientation of cell divisions, we define a unit nematic tensor (see Materials and methods). For each cell division, the orientation of this unit nematic is defined by the line connecting the centers of mass of the two daughter cells when they first appear (see Figure 3C–C', and TM R-User Manual section 2.8). Each nematic is assigned a position on the tissue that corresponds to the center of combined mass of the two daughter cells. To visualize division orientation patterns, unit nematics can be added within different regions and averaged over different time intervals (Figure 3D, Video 10, TMR-User Manual section 2.9).

Visualizing cell junction dynamics (TM R-User Manual sections 2.10 – 2.12, 3.8–3.9)

Epithelial tissues can be reshaped by cell rearrangements, or T1 transitions (for review [Walck-Shannon and Hardin, 2014]). In the simplest case, a T1 transition involves two pairs of cells, that exchange neighbors by disassembling one cell-cell contact and replacing it by another – bringing together two previously separated cells (Figure 4A). In reality, cell contacts may undergo multiple rounds of shrinkage and regrowth before resolving (Figure 4A'). Furthermore some epithelia undergo the related process of rosette formation where multiple cell junctions are disassembled before new neighbors are brought into contact. By separately quantifying the orientation with which cell contacts are gained and lost, one can reveal whether there is a net directionality to cell junction assembly and disassembly. To identify gained and lost cell contacts, we compare cell neighbor relationships between 2 subsequent frames. We exclude changes in neighbor relationships resulting from cell division, extrusion or a cell moving in and out of the field of view. The remaining neighbor relationship changes are used to define cell contacts that have appeared or disappeared.

Open in a separate window

Figure 4.

Visualization and quantification of T1 transitions.

(A–A') Cartoon depicting an effective T1 transition (A) that corresponds to cell-contact loss and gain in different directions. Each contact loss or gain is assigned a unit nematic describing its orientation. (B–B') Pattern of cells losing contact (green), gaining contact (red) or both (blue). (C) Rate of neighbor change per cell and per hour in the blade, veins and inter-vein regions of interests. Rates were averaged within discrete time intervals of one hour and further averaged among the 3 WT wings (TM R-User Manual, section 3.8). Error bars depict the standard deviation amongst wings. (D) Coarse-grained pattern of neighbor exchange orientation at 17 hAPF. Cell neighbor change nematics were obtained by summing up unit nematics in each grid elements of 33x33 microns and further averaged in time using a 50 min time window. Scale bar 50 microns.

DOI: http://dx.doi.org/10.7554/eLife.14334.018

Figure 4—figure supplement 1.

Open in a separate window

T1 and cell elongation nematic orientation.

(A) Cell neighbor change nematics were averaged at each frame within each region of interest and are represented as bars in a circular diagram. The bar angle indicate the average T1 orientation, and its length (nematic norm) reflects how ordered cell neighbor change nematics are in a given region of interest. Their color depicts the developmental time in hours after puparium formation. (B) Cell elongtation nematics were also averaged at each frame within each region of interest. The average T1 nematic orientation starts to match the average cell elongation nematic orientation from about 22 hAPF (peak of cell stretch) on, when stress-induced PD-oriented T1 dominate over autonomous AP-oriented T1.

DOI: http://dx.doi.org/10.7554/eLife.14334.019

We characterize the orientation of contact gains and losses by assigning them a unit nematic tensor. For contact loss, the orientation of the nematic is defined by the axis intersecting the two cell centers. For contact gain, the orientation of the nematic is perpendicular to the axis intersecting the two cell centers (Figure 4A–A’). If there is a simple disappearance and reappearance of a single cell contact, corresponding nematics will cancel out. Therefore, the sum of contact gain and contact loss nematics over time and/or space will represent an effective T1 nematic describing net direction of contact assembly/disassembly.

The rate of contact gain and loss can be visualized in different ways. Cell contact dynamics can be viewed directly on movies of tissue morphogenesis by assigning colors to cells as they gain (red) or lose (green) contacts. Those cells that simultaneously gain and lose different cell contacts are colored blue (Figure 4B–B’).

The frequency of contact gain and loss, independent of orientation, can be plotted over time. Figure 4C compares the frequency of contact assembly/disassembly in vein and inter-vein regions. In both regions, this rate begins to decrease in the second half of morphogenesis.

To visualize the pattern of orientation of T1 transitions throughout the wing, we sum contact gain and loss nematics over square grid elements, and average over a chosen time window (about 50 min in Figure 4D, Video 11, see TM R-User Manual section 2.12). 

Finally, the average orientation of effective T1 nematics in sub-regions over time can be visualized using circular diagrams, where nematics are color-coded to indicate developmental time. Figure 4—figure supplement 1A reveals that the orientation of effective T1’s is along the anterior-posterior (AP) axis early (blue) and shifts to the proximal-distal (PD) axis in the second half of morphogenesis (red). A similar approach can be used to illustrate average cell elongation nematics over time (Figure 4—figure supplement 1B).

Quantification of tissue deformation and the contribution of different cellular events (TM R-user manual section 3.10)

While it is useful to quantify the number and orientation of cellular events like elongation, rearrangement, extrusion and division, this by itself does not provide quantitative information about the amount of tissue shape change contributed by each type of event. We therefore devised a method to measure deformation caused by these cellular processes such that they sum to the measured tissue deformation.

Tissue deformation can be decomposed into isotropic and anisotropic parts that distinguish changes in area (compression/expansion) from changes in aspect ratio (pure shear, for details see also Materials and methods). The quantities describing area changes are scalar, whereas the quantities describing shear rate in a 2D-network are nematic tensors harboring two distinct components that describe the orientation and magnitude of the shear.

Tissue area changes can be calculated based on cell area change and the number of cells gained and lost by divisions and extrusions – information that is all available in the TissueMiner database (Etournay et al., 2015).

To quantify the cellular contributions to anisotropic tissue deformation, TissueMiner uses the so-called Triangle Method, which is based on a triangular tiling of the junctional network (Etournay et al., 2015; Merkel et al., 2016). Triangle elongation is a proxy for cell elongation, and topological changes in the network result in redrawing of triangles (Figure 5A–C). The resulting change in average triangle elongation can be used to calculate the shear due to the topological changes (Etournay et al., 2015). In addition to contributions from divisions, cell rearrangements, extrusions and cell shape changes, the method also takes into account deformation caused by correlations between elongation and both area change and rotation.

Open in a separate window

Figure 5.

Visualization and quantification of anisotropic cell and tissue deformation.

(A) Triangulation of the cell network: each triangle vertex corresponds to a cell center. (B–B') Cartons depicting triangle pure shear and total tissue shear along the x axis. (C) Cartons depicting shear due to T1 transition, cell division and extrusion. (D) Pattern of local tissue shear rate obtained from the triangulation method. Scale bar 50 microns. (E) shows the average rate of tissue shear (blue) in the blade, interveins and veins, and the corresponding cellular shear contributions (other colors). Shaded regions indicate the standard deviation amongst wings. (F) shows the accumulated tissue shear over time and the accumulated contributions of each type of cellular event. The tissue shear (blue) in veins is orientated along the PD axis and it is higher than in inter-vein regions during most of pupal morphogenesis. It leads to an extension along the PD axis and to a narrowing along the anterior-posterior (AP) direction. By the end of the movie, accumulated tissue shear (blue) is almost twice as high in veins as in inter-vein regions. Shaded regions represent the standard deviation amongst wings.

DOI: http://dx.doi.org/10.7554/eLife.14334.021

Figure 5—figure supplement 1.

Open in a separate window

Measurements of cell and tissue deformation from two computer-generated sheets of hexagonal cells.

(A–D) One dataset corresponds to hexagonal cells undergoing a constant isotropic expansion rate of 3.50 10^–2 per frame, and the other corresponds to hexagonal cells undergoing constant pure shear rate of 1.75 10^–2 per frame. These datasets are termed iso.exp movie and shear movie respectively in graphs. There isn't any topological change. To keep consistent sets of cells in time, we filtered out cells that become in contact to the image border. We then performed our measurement on these tracked regions of about 50 cells in the shear movie and about 100 cells in the iso.exp movie. (A) Relative tissue area changes (blue) and its decomposition into cell area changes (green), cell number increase by divisions (orange) and cell number descrease by extrusions (cyan). Their corresponding cumulative sums are shown in (B). (C) shows the average tissue shear (blue) and its decomposition into cellular shear contributions (other colors). Their corresponding cumulative sums are shown in (D).

DOI: http://dx.doi.org/10.7554/eLife.14334.022

Figure 5—figure supplement 2.

Open in a separate window

Tissue isotropic deformation and cellular contributions in different regions.

(A) Relative rates of tissue area changes (blue) averaged over 3 WT wings for the blade, veins and interveins, and its decomposition into cell area changes (green), cell number increase by divisions (orange) and cell number descrease by extrusions (cyan). Their corresponding cumulative sums are shown in (B). (B) Cumulative tissue area changes and its cellular contributions. Shaded regions represent the standard deviation amongst wings.

DOI: http://dx.doi.org/10.7554/eLife.14334.023

Figure 5—figure supplement 3.

Open in a separate window

Comparison of patterns of cell event orientation with their correponding quantitative patterns of shear.

(A–A') Coarse-grained patterns of cell division orientation (A) and of shear contributed by cell division (A'). The pattern shown in (A) was obtained by summing up cell division nematics in each grid element and by further averaging in time. The pattern shown in (A') was obtained by averaging the shear nematics in each grid element and by further averaging in time. (B–B') Coarse-grained patterns of neighbor-change orientation (B) and of shear contributed by neighbor changes (B'). These patterns were obtained similarly as for cell divisions. Only the shear patterns (A' and B') obtained with the triangulation method provide a quantitative measurement of the local deformation induced by each type of cellular event. Square-grid size of 26x26 microns. Time averaging covering about 55 min (11 frames) in each grid element. Scale bar 50 microns.

DOI: http://dx.doi.org/10.7554/eLife.14334.024

Validation of tissue deformation measurements using computer-generated cells

To test the reliability of TissueMiner in calculating large cell and tissue deformations, we created two computer-generated movies of hexagonal cells sheets (Videos 12, 13). In one movie, we imposed a constant isotropic expansion rate of 3.50 10⁻² per frame, without any anisotropic deformation. In the second movie, we imposed a constant pure shear along the x-axis with a rate of 1.75 10⁻² per frame, and without any isotropic expansion. The amounts of isotropic expansion and pure shear have been chosen to be at least 10 times higher than what we measure between subsequent frames of pupal wing movies.

We then asked if TissueMiner could quantitatively recapitulate the respectively imposed deformation rates. In each dataset, TissueMiner automatically defines a 'whole_tissue' region of interest that corresponds to a consistent set of cells that are always visible (about 100 cells in the isotropic expansion movie and about 50 cells in the pure shear movie, green labels in Videos 12 and 13). All measurements are done in this ROI to avoid measuring deformation due to inward and outward cell flows. Figure 5—figure supplement 1 shows the time evolution of the measured tissue expansion rate (panel A) and tissue shear rate (panel C) that were averaged over the 'whole_tissue' ROI, and their respective cellular contributions. Panels B and D show the corresponding cumulated curves. As expected, in the isotropic expansion movie we observe a nearly constant isotropic expansion rate, which is accounted for by the cell area change contribution. We measure an average expansion rate of (3.53 ± 0.04) 10^–2 per frame, which is consistent with the value imposed when creating the movie. The measured uncertainty is the 95% confidence interval of the standard error of the mean. The pure shear rate and its cellular contributions nearly vanish in this movie (Figure 5—figure supplement 1C,D).

For the pure shear movie, we measure an approximately constant horizontal component of the pure shear rate of (1.74 ± 0.02) 10^–2 per frame, which is consistent with the value imposed when creating the movie. This pure shear rate is entirely accounted for by cell elongation change. The isotropic expansion rate and its cellular contributions nearly vanish (Figure 5—figure supplement 1A, B). Other contributions to expansion and shear rates are negligible in both movies.

The pixelated nature of individual cell contours contributes to fluctuations of our measured values. Moreover, we find that these fluctuations cancel out when cumulating the deformation (Figure 5—figure supplement 1B and D). Thus, the current implementation of TissueMiner captures the tissue isotropic expansion and pure shear rates as well as the corresponding cellular contributions with a good precision in these computer-generated movies.

A relational database to store the history of cells, their lineage and their constituent bonds and vertices

The visualization and quantification of cell dynamics underlying tissue morphogenesis rely on the ability to extract information about cell geometry, cell neighbor topology and cell histories from time-lapse movies (Aigouy et al., 2010; Etournay et al., 2015). We use TissueAnalyzer to segment and track the cell network over time. This results in a series of digital images that contain this information (Figure 6—figure supplement 1). To facilitate its access and use, we developed tools in the TissueMiner framework to extract and convert this information initially stored in images into a specific database format (see details in appendix 1), which we call 'TM-DB' (schematically outlined in Figure 6A).

Open in a separate window

Figure 6.

Construction of the relational database of TissueMiner.

(A) Conceptual scheme of the database. Entities (square boxes) are related to other entities by associations (rounded boxes). Each entity contains an identifier (underlined) that uniquely defines each record. The database can be implemented by converting entities into tables (see appendix 1 and Figure 6—figure supplement 2). (B) Cell lineage trees are stored in the database: upon division a mother cell identifier a gives rise to two new daughter cell identifiers b and c. {a,b,c,d,e,f,g} defines one lineage group. (C) A pixelated cell contour in the 2D cell network: green=bond pixels, red=vertex pixels, white=other cell network pixels. (D) Vectorized representation of the cell shown in (C). To preserve the topology of the cell network, directed bonds (cyan) are defined from within a given cell alpha and ordered anticlockwisely along the cell contour. Each directed bond is complemented by a conjugated bond (magenta) and is linked to it next counter-clockwise follower (dashed).

DOI: http://dx.doi.org/10.7554/eLife.14334.027

Figure 6—figure supplement 1.

Open in a separate window

Tracked cells identified by unique colors in TissueAnalyzer.

(A) shows two consecutive frames depicting colored-tracked cells from a time-lapse movie processed with TissueAnalyzer. Each cell is assigned a color identifier that uniquely defines it in the course of the time-lapse. One pixel wide cell-cell interfaces are visible in white on the raster image.

DOI: http://dx.doi.org/10.7554/eLife.14334.028

Figure 6—figure supplement 2.

Open in a separate window

Logical scheme of the relational database.

(A) The conceptual scheme shown in Figure 6A can be automatically converted to a logical scheme shown here by using softwares such as IntelliJDEA or MySQL workbench. The rules of conversion are briefly evoked in appendix 1. The entities defined in the conceptual scheme are converted into tables containing one primary key (upper part of the table) that uniquely defines each record in the table, the properties of each record, and the foreign keys (arrows). Foreign keys are properties of one table pointing to the primary key of a related table (ex: conj_dbond_id:dbond_id means that the conj_dbond_id column is a foreign key whose values must be defined in the dbond_id column of the dbonds table). As a consequence of logical contraints by foreign keys, tables harbor more columns that one expected from looking at Figure 6A. This logical scheme now shows all tables and columns of the database. This scheme is implemented in physical SQLite tables can are indexed for the sake of performance (see CreateDbFromParser.R on https://github.com/mpicbg-scicomp/tissue_miner).

DOI: http://dx.doi.org/10.7554/eLife.14334.029

First, the history of each tracked cell in the movie is stored as a separate row in the cell_histories table of the TM-DB (Figure 6A). This includes the movie frames in which it first appears and disappears and why, along with its lineage relationship to other cells (see appendix 1). The reason for cell appearance and disappearance is inferred by the parser. A primary reason could be a cell division, which results in the disappearance of the mother cell and in the appearance of two daughter cells. It could be a cell extrusion that results in its disappearance. It could also be that cells move in and out of the field of view of the microscope lens, resulting in gain and loss of cells. Furthermore, we use this information to establish the lineage relationship that corresponds to each group of cells related by ancestry (Figure 6B). Each cell within the lineage group is assigned a generation number. The lineage group and generation number for each cell are listed in the cell_histories table.

We store the time points at which the movie images were recorded into a frames table that links each movie frame to its corresponding time point. For each movie frame, we need to store geometrical and topological information about cells within the cellular network. Geometrical information includes position and shape descriptors, whereas topological information indicates the arrangement of neighboring cells around each cell. We use cell histories, geometry and topology to understand how individual cells contribute to the whole tissue deformation during morphogenesis (Etournay et al., 2015).

The geometrical information is stored in three tables of the TM-DB: cells, bonds and vertices. They correspond to the 3 generic entities - cells, cell-cell contacts and intersections between cell-cell contacts, respectively illustrated in Figure 6C. These entities are commonly used in vertex model simulations (for review [Fletcher et al., 2014]). The cells table contains cell geometrical data (center of mass, area, shape anisotropy) and the polarized distribution of proteins along the cell circumference, as represented by a polarity nematic tensor (Aigouy et al., 2010). The bonds table informs about bond length, and the vertices table about vertex position in each movie frame.

The directed_bonds table exclusively stores the cell neighbor topological information at each frame, i.e. how bonds are organized around each cell along with the cell neighbor relationship information. To store the cell neighbor topology in an unambiguous manner, we define for each cell a directed path of consecutive bond vectors oriented counterclockwise, which forms the oriented circumference of the cell (Figure 6D, see also [Kachalo et al., 2015]). We link each directed bond to its counter clockwise follower (left directed bond) in the same cell. To store the cell neighbor relationship, we link each directed bond to its corresponding directed bond (conjugated bond) of the neighboring cell (Figure 6D, and appendix 1).

The TM-DB is relational, which means that it establishes contextual relationships between items stored in one ore more tables (see appendix 1). These relationships are outlined in rounded boxes in the conceptual scheme of the TM-DB (Figure 6A). Technically, each item in a table is stored in a separate row and is given a unique number as identifier. For a relationship between two tables, one of the tables contains an additional column, which refers to items in the other table by holding their identifier number. Such additional columns for the TM-DB format are shown in blue in Figure 6—figure supplement 2. When extracting information from a database using so-called queries, these columns serve as bridges connecting the information stored about related items.

In essence, this structure creates a generic relational model to represent complex cell tracking data in 2D. In practice, the data for each movie is stored in a separate SQLite database file. Since all movie files are stored using the same database structure, automated data mining and visualization are greatly facilitated. For the same reason, usage of the TissueMiner database format encourages exchange of both movie data and analysis tools.

An automated workflow compliant with high performance computing platforms

To help the user to perform complex tissue morphogenesis analysis, we developed an automated pipeline that uses the tracked data from TissueAnalyzer as an input to build the database and perform all downstream analyses described above. To do so, we use the snakemake workflow engine developed by Koster and Rahmann (Koster and Rahmann, 2012). This engine channels the different processing steps into a well-formed workflow graph. Snakemake automatically determines the execution order, provides means for error recovery and job control, and supports High Performance Computing (HPC) environments. By using snakemake we enable the user to easily run and monitor TissueMiner, while maintaining a proper decoupling of tools as independent executables.

Practically, the user defines a workflow definition file in which processing steps are defined as a set of execution rules, namely a list of scripts to be run along with required input(s) and expected output(s). Snakemake automatically builds a directed graph from which the execution order of processing steps is inferred. If only one branch of the graph needs to be run, the engine will ensure that all input data are present and will automatically run upstream steps if necessary. This engine also provides the possibility to visualize a directed acyclic execution dependency and execution state graph (DAG) for a given workflow (see Figure 7).

Open in a separate window

Figure 7.

Automated workflow using snakemake.

(A) The snakemake engine can generate a directed acyclic graph (DAG) where we show an example here. This graph represents both the execution dependency (grey arrows) and the execution state of the workflow (solid or dashed line). Each box corresponds to an execution rule, namely a program to be run along with required input(s) and expected output(s). This DAG can be generated at any time when running the workflow (see documentation). Solid lines indicate the rules that have not been executed yet, whereas dashed lines depict completed jobs. The first rule to be executed is called 'prepare_movie': it prepares the tracked images from TissueAnalyzer to be converted by the parser into tables of values containing all necessary entities along with their properties ('parse_tables' rule). Then the 'make_db' rule is executed for building the database. Following the grey arrows can one navigate into the next steps of the workflow. The 'roi_tracking' rule filters out cells in contact to margin cells including user-defined regions of interest, and the 'roi_movie' rule allows us to visualize regions of interest over time. The 'deformation_movies' and 'db_elongation_movies' rules generate annotated movies to visualize the deformation of the tissue and the cell state properties (area, elongation). The 'four_way' rule detects four-way vertices and performs basic statistics on vertices. The 'tri_create' rule performs the triangulation of the network for further shear calculation and visualization ('shear_calculate' and 'shear_movies'). It also enables triangle tracking and mapping to each type of cell event ('tri_categorize'). The 'topo_countT1' rule detects neighbor changes that are not due to division or extrusion, and categorizes them into gained or lost neighbors. The 'topo_movies' rule allows one to visualize the coarse-grained rates of division and neighbor changes on the tissue. The 'topo_unbalance' rule is a quality check to verify that the number of gained neighbors is similar to the number of lost neighbors. The 'polygon_class' rule performs the cell-neighbor number count. The 'lineage_colors' rule allows us to optimize the color of each lineage group such that adjacent lineage groups always have different colors. Finally, the 'lineage_movies' rule allows one to visualize lineage groups and cell generations on the tissue. The rule 'all' checks that all upstream jobs have been completed.

DOI: http://dx.doi.org/10.7554/eLife.14334.030

One major advantage of a workflow engine such as snakemake is that it can run the workflow on various architectures - from single-core workstations to multi-core servers and clusters - without the need to modify the rules, thereby facilitating reproducible research. To simplify the TissueMiner installation procedure, we provide a pre-configured system to be loaded in the docker software available at http://docker.com. The TissueMiner docker image can be run without any setup using provided example data or custom user data as detailed out on the TissueMiner GitHub project page. More advanced users can use TissueMiner directly from the command-line with or without snakemake and can thus perform simultaneous analyses of multiple movies.

A user-friendly data-mining library to easily collect information for comparing multiple datasets

After applying our automated workflow to different movies, the results can be easily compared using a collection of command-line tools written in R and Python. These tools aggregate different experiments for plotting and performing comparative analysis. Here we describe the tools written in R. Python tools are described in the corresponding tutorial. The R tools are designed to be used in an integrated development environment such as RStudio, which provides a user-friendly environment to assist the user in writing and executing command lines. These command line tools are organized in the spirit of a grammar of data manipulation and they can be combined with the existing R tools like dplyr (Francois and Wickham, 2015) or ggplot2 (Wickham, 2009) for manipulating and visualizing data (https://mpicbg-scicomp.github.io/tissue_miner/user_manual/Learning_the_R_basics_for_TissueMiner.html).

We developed generic 'multi-query functions' (mqf) to collect specific information for individual movies. These mqf tools are organized into fine-grained and coarse-grained categories according to the type of analysis to be carried out. The fine-grained tools aggregate data at cellular level, namely individual cell properties inside regions of interest. These tools are prefixed with 'mqf_fg_'. The coarse-grained mqf tools are further separated into 'roi' and 'grid' categories to distinguish between regions moving with the tissue and static square regions tiled into a grid. They allow one to visualize and quantify average cell properties at different tissue locations and various spatial scales, and are prefixed with 'mqf_cg_roi_' and 'mqf_cg_grid_' respectively.

To compare fine-grained and coarse-grained cell properties amongst movies we developed a 'multi-db-query' tool, which streamlines the application of the mqf tools to a set of movies. To use this tool, the user should first align the movies in time, using convenient morphological or cellular landmarks. As for the Drosophila wing, we align movies such that the peaks of cell elongation coincide in the different movies. The user can then apply a chosen mqf tool to multiple movies and multiple ROI’s. All mqf tools, alone or in combination with the 'multi-db-query' tool, generate a table that contains individual or averaged measurements to be visualized on the tissue (Figure 1 A–E’, Figure 2A–D’, Figure 3A,D, Figure 4B,D, Figure 5D) or in graphs (Figure 2E–F, Figure 3B, Figure 4C, Figure 5E–F). This library of tools is described in detail in the TM R-User Manual, which also provides many examples. These tools can be easily extended to address project specific questions.

Detecting gain and loss of cell contacts

To detect cell neighbor changes, we developed a routine in R that queries the DB and establishes the cell-neighbor relationship at each frame. By comparing the list of neighboring cell identifiers for a given cell between two consecutive frames [f, f + 1], can one identify and count the changes in neighbor relationships. These can be subdivided into those caused by cell divisions, cell extrusions or the simple gain or loss of a cell contact (not due to division or extrusion). We call these half-T1’s because they resemble the gain and loss of cell contacts that occurs during a T1 transition – although they may also be generated by other events such as rosette formation. To assign a neighbor change to the half-T1 category, the corresponding cell identifiers must be present in both frames, ruling out extrusions and cells moving in and out of the field of view. To detect half-T1’s that occur simultaneously with divisions, we mask neighbor changes due to divisions by propagating the mother cell identifier (frame f) to the two daughter cells (frame f+1) that we fuse into one fake cell having the mother cell identifier. We iterate over each pair of consecutive frames and store the half-T1 events due to a gain and a loss of cell neighbors.

Cell lineages and lineage browsing to follow ROI’s forward and backward in time

We pool all lineage information (as contained in the cell_id, left_daughter_cell_id and right_daughter_cell_id columns from the cell_histories table) into a directed lineage graph (Nepusz, 2006) from which we infer a lineage group identifier and a generation number. By definition the root of each lineage tree is considered as the F₀ generation and is thus given a generation value of 0. We follow ROI’s backward and forward in time by browsing lineage graphs that were selected based on the regions drawn by using the draw_n_get_ROIcoord.ijm Fiji macro. However cells may be lost or not detected when browsing the lineages. One primary reason is that extruding cells are not detected when browsing the lineage backward in time. Cells could also be lost due to possible tracking mistakes. To improve spatial continuity of ROI’s we have implemented a method to reassign lost cells to ROI’s when located within ROI’s. To identify lost cells for a frame within a given ROI, we first distinguish bonds that connect two cells within the ROI, only one cell within the ROI or none. All corresponding cell-pairs define an undirected graph on which a connected component analysis (Nepusz, 2006) allows to identify the ROI and non-ROI regions. All cells of non-ROI regions, except for the largest one, are reassigned to become part of the ROI. By doing so, we make the assumption that the largest non-ROI component is defined by the tissue surrounding the ROI.

Nematic tensors to describe cell elongation and the orientation of cellular processes

When analyzing and visualizing single cell properties, we use the same cell elongation definition as in Aigouy et al. (2010). For a given Cartesian xy coordinate system, the elongation of a given cell is defined by the nematic tensor

with

Here, A_c is the area of the given cell, and the integrals are carried out over all points r within the cell. The angle is the angle between the vector r − r_c and the x axis, where r_c is the cell center defined as

Here, the integral is again carried out over all points r within the cell. The magnitude of the elongation is given by ϵ = (ϵ_xx²+ϵ_xy²)^½ and the elongation angle φ is given by the following two equations

Note that this definition of cell elongation is different from the triangle-based definition that is also discussed in this article. However for the fruit fly wing, both cell elongation definitions yield very similar results.

To characterize the axes of cell divisions and T1 transition, we introduce the unit nematic tensors ${\displaystyle {\tilde{n}}_{\mathrm{C}\mathrm{D}}}$, ${\displaystyle {\tilde{n}}_{\mathrm{T}1+}}$, and ${\displaystyle {\tilde{n}}_{\mathrm{T}1-}}$. The orientation of a single cell division is quantified by the unit nematic ${\displaystyle {\tilde{n}}_{\mathrm{C}\mathrm{D}}}$ defined by:

Here, the angle ϕ is the angle of the line connecting both cell centers with respect to the x axis, measured in counter-clockwise sense. The orientation for a half-T1 transition during which two cell lose neighborship is characterized by:

where ϕ_T1+ is the angle of the line connecting the centers of the cells losing neighborship. The orientation for a half-T1 transition during which two cell gain neighborship is characterized by:

where ϕ_T1− is the angle of the line connecting the centers of the cells that gain neighborship. The axes of the nematics ${\displaystyle {\tilde{n}}_{\mathrm{C}\mathrm{D}}}$, ${\displaystyle {\tilde{n}}_{\mathrm{T}1+}}$ and ${\displaystyle {\tilde{n}}_{\mathrm{T}1-}}$ roughly correspond to the axis along which the tissue extends due to the respective cell division or half-T1 transition. In particular, because of the minus sign in the definition of ${\displaystyle {\tilde{n}}_{\mathrm{T}1-}}$, when the same two cells gain neighborship and lose it again along the same axis, the total effect adding ${\displaystyle {\tilde{n}}_{\mathrm{T}1+}}$ and ${\displaystyle {\tilde{n}}_{\mathrm{T}1-}}$ is zero.

This work was a truly collaborative effort and the authors jointly wrote the manuscript. We are grateful to Christian Dahmann, Marcus Michel and Jacques Boutet de Monvel for critical reading of the manuscript, Benoit Lombardo for his help in Fiji macroing, and Peter Steinbach for useful discussion about workflow engines. We thank Franz Gruber, Vincent Michel and Nathalie Gourreau for testing the quickstart tutorials. RE acknowledges a Marie Curie fellowship from the EU 7th Framework Programme (FP7). This work was supported by the Max Planck Gesellschaft, and by the BMBF. SE acknowledges funding from the ERC.