Abstract
Analysis of coupled variables is a core concept of cell biological inference, with colocalization of two molecules as a proxy for protein interaction being a ubiquitous example. However, external effectors may influence the observed colocalization independently from the local interaction of two proteins. Such global bias, although biologically meaningful, is often neglected when interpreting colocalization. Here, we describe DeBias, a computational method to quantify and decouple global bias from local interactions between variables by modeling the observed colocalization as the cumulative contribution of a global and a local component. We showcase four applications of DeBias in different areas of cell biology, and demonstrate that the global bias encapsulates fundamental mechanistic insight into cellular behavior. The DeBias software package is freely accessible online via a webserver at https://debias.biohpc.swmed.edu.
https://doi.org/10.7554/eLife.22323.001eLife digest
Cell biologists often use microscopes to look closely at cells and see what is happening during an experiment. Cell biology experiments typically involve measuring more than one aspect of the cells, for example, the forces a cell is experiencing and the direction it is moving, or the locations of two different components in the cell. The task is then to decipher the interactions between these independent variables to better understand the inner workings of a living cell.
This task, however, can be challenging because other variables can mask the interactions between the pairs of variables being studied. For example, it is difficult to know if two components of cells overlap in microscopy images because they directly interact or simply because the overall organization of the cell makes it more likely they will end up in the same place. Several statistical methods have been developed to estimate and eliminate such confounding effects to reveal specific interactions. However, these confounding effects – termed “global biases” – may themselves contain valuable information about how cells work. Ignoring the possible roles of global biases may limit our understanding of biological processes.
Zaritsky et al. have now addressed this issue by developing an algorithm called DeBias that takes global bias into account by decoupling it from true interactions between two variables. DeBias uses global bias as a second measurement to understand how two biological variables interact via confounding variables. Zaritsky et al. then used DeBias with data from four different cell biology experiments to show that the algorithm works. For example, animal cells contain several proteins that form filaments, which give them their shape and help them move. Socalled vimentin filaments and microtubules had previously been seen to occur in the same place within cells but it was not clear whether their alignment was due to local interactions or determined by the overall shape of the cell. DeBias was used to analyze cells that were either still or moving in a particular direction. The algorithm could tease apart the effect of this movement and showed that coalignment of vimentin filaments and microtubules was caused more by the movement and shape of the cell and less by specific interactions.
Overall, DeBias is a new mathematical tool for cell biologists that is freely accessible online. Other researchers can now use this tool in future studies to identify local interactions and global biases in a wide range of cell biology experiments and interpret the data in a meaningful way.
https://doi.org/10.7554/eLife.22323.002Introduction
Interpretation of the relations among coupled variables is a classic problem in many areas of cell biology. One example is the spatiotemporal colocalization of molecules – a critical clue to interactions between molecular components; another example is alignment of molecular structures, such as filamentous networks. However, colocalization or alignment may also occur because the observed components are associated with external effectors. For example, the internal components of a polarized cell are organized along the polarization axis, making it difficult to quantify how much of the observed alignment between two filamentous networks is related to common organizational constraints, and how much of it is indeed caused by direct interaction between filaments. Another example is introduced with protein colocalization, where their intensity distributions may be heavily biased to specific levels regulated by the cell state. The combined effects of global bias with local interactions are manifest in the joint distribution of the spatially coupled variables. The contribution of global bias to this joint distribution can be recognized from the deviation of the marginal distributions of each of the two variables from an (unbiased) uniform distribution.
Although global bias can significantly mislead the interpretation of colocalization and coorientation measurements, most studies do not account for this effect (Adler and Parmryd, 2010; Bolte and Cordelières, 2006; Costes et al., 2004; Das et al., 2015; Dunn et al., 2011; Kalaidzidis et al., 2015; Rizk et al., 2014; SerraPicamal et al., 2012; Tambe et al., 2011). Previous approaches indirectly assessed spatial correlations (e.g., [Drew et al., 2015; Karlon et al., 1999]), variants of mutual information (e.g., [Krishnaswamy et al., 2014; Reshef et al., 2011]) or spatial biases (Helmuth et al., 2010) but did not explicitly quantify the contribution of the global bias to the observed joint distribution. These methods approach the global bias as a confounding factor (VanderWeele and Shpitser, 2013) that must be eliminated for more accurate assessment of the true local interaction, but ignore the possibility that the global bias contains byitself valuable mechanistic information to cell behavior.
Here, we present DeBias as an algorithm to decouple the global bias (represented by a global index) from the bona fide local interaction (represented by a local index) in colocalization and coorientation of two independentlymeasured spatial variables. The decoupling enables simultaneous investigation of processes that drive global bias and local interactions between spatiallymatched variables. Our method is dubbed DeBias because it Decouples the global Bias from local interactions between two variables.
To highlight its capabilities, DeBias was applied to data from four different areas in cell biology, ranging in scale from macromolecular to multicellular: (1) alignment of vimentin fibers and microtubules in the context of polarized cells; (2) alignment of cell velocity and traction stress during collective migration; (3) ﬂuorescence resonance energy transfer of Protein Kinase C; and (4) recruitment of transmembrane receptors to clathrincoated pits during endocytosis. These examples demonstrate the generalization of the method and underline the potential of extracting global bias as an independent functional measurement in the analysis of multiplex biological variables.
Results
Similarity of observed coorientation originating from different mechanisms
The issue of separating contributions from global bias and local interactions is best illustrated with the alignment of two sets of variables that carry orientational information. Examples of coorientation include the alignment of two filament networks (Drew et al., 2015; Gan et al., 2016; Nieuwenhuizen et al., 2015), or the alignment of cell velocity and traction stress, a phenomenon referred to as plithotaxis (Das et al., 2015; Tambe et al., 2011; Trepat and Fredberg, 2011). In these systems, global bias imposes a preferred axis of orientation on the two variables, which is independent of the local interactions between the two variables (Figure 1A).
Similar observed alignments may arise from different levels of global bias and local interactions. This is demonstrated by simulation of two independent random variables X and Y, representing orientations (Figure 1B, left), from which pairs of samples x_{i} and y_{i} are drawn to form an alignment angle θ_{i} (Figure 1B, middle). Then, a local interaction between the two variables is modeled by coaligning θ_{i} by ζ_{i} degrees, resulting in two variables x_{i}’ and y_{i}’ with an observed alignment θ_{i}  ζ_{i} (Figure 1B, right).
We show the joint distribution of X, Y for four simulations (Figure 1C) where X and Y are normally distributed with identical means but different standard deviations (σ), truncated to $\left[{90}^{\circ},{90}^{\circ}\right]$, and different magnitudes of local interactions (ζ). The latter is defined as ζ = αθ (Figure 1B, α = 1 for perfect alignment). Throughout the simulations both σ and α are gradually increased (Figure 1C, lefttoright), implying that the global bias in the orientational variables is reduced while their local interactions increase. As a result, all simulations display similar observed alignments (mean values, 18.9°−19.5°). Figure 1D visualizes 100 samples from each of the two most distinct scenarios: low σ and no local interaction (σ = 17°, α = 0) leads to tendency of X and Y to align independently to one direction (left); higher variance together with increased interaction (σ = 40°, α = 0.5) leads to more diverse orientations of X and Y (right), while maintaining similar mean alignment. This simple example highlights the possibility of observing similar alignments arising from different mechanisms. While the described properties are well known and many others have used statistical postprocessing to eliminate confounding factors for accurate assessment of local interactions (Drew et al., 2015; Helmuth et al., 2010; Karlon et al., 1999; Krishnaswamy et al., 2014; Reshef et al., 2011), we aim at directly quantifying the global bias, with the goal of extracting encapsulated information that is fundamental to the biological question.
DeBias: a method to assess the global and local contribution to observed coalignment
DeBias models the observed marginal distributions X’ and Y’ as the sum of contributions by a common effector, i.e., the global bias, and by local interactions that effect the coalignment of the two variables in every data point (Figure 2A).
In a scenario without any global bias or local interaction between X’ and Y’, the observed alignment would be uniformly distributed (denoted uniform). Hence any deviation from the uniform distribution would reflect contributions from both the global bias and the local interactions. To extract the contribution of the global bias we constructed a resampled alignment distribution (denoted resampled) from independent samples of the marginal distributions X’ and Y’, which decouples matched pairs (x_{i}’, y_{i}’), and thus excludes their local interactions. The global bias is defined as the dissimilarity between the uniform and resampled distributions and accordingly describes to what extent elements of X’ and Y’ are aligned without local interaction (Figure 2B). If a local interaction exists then the distribution of the observed alignment angles will differ from the independently resampled alignment distribution. Hence, the uniform distribution will be less similar to the experimentally observed alignment distribution (denoted observed) than to the resampled distribution. Accordingly, the local interaction is defined by the difference of dissimilarity between the observed and uniform distributions and dissimilarity between the resampled and uniform distributions (Figure 2B).
The Earth Mover's Distance (EMD) (Peleg et al., 1989; Rubner et al., 2000) was used to calculate dissimilarities between distributions. The EMD of 1dimensional distributions is defined as the minimal 'cost' to transform one distribution into the other (Kantorovich and Rubinstein, 1958). This cost is proportional to the minimal accumulated number of moving observations to adjacent histogram bins needed to complete the transformation. Formally, we calculate $EMD\left(A,B\right)={\sum}_{i=1,\dots ,K}{\sum}_{j=1,\dots ,i}{a}_{j}{\sum}_{j=1,\dots ,i}{b}_{j}$, with K  number of histogram bins, a_{j} and b_{j}  fraction of observations in bin j for distributions A and B correspondingly. Introducing the EMD defines scalar values for the dissimilarities and allows us to define the EMD between resampled and uniform alignment distributions as the global index (GI) and the local index (LI) as the difference of EMD between observed and uniform and the GI (Figure 2B).Figure 2C, demonstrates how the GI and LI recognize the global bias and local interactions between the matched variable pairs (x_{i}’, y_{i}’) established in Figure 1C. For a scenario with no local interaction (α = 0) DeBias correctly reports LI~0 and GI~3. For a scenario with gradually wider distributions X,Y, i.e., less global bias, and gradually stronger local interactions (α > 0), the LI increases while the GI decreases. Additional simulations showed that similar properties apply for negative local interactions ζ = αθ (Figure 1B) were α < 0 (Figure 2—figure supplement 1).
In the previous illustrations, changes in spread of the distributions X and Y were compensated by changes in the local interactions. When leaving the interaction parameter α constant while changing the spread of X and Y, a weak, but intrinsically negative correlation between LI and GI becomes apparent (Figure 2D–E). Thus, while DeBias can correctly distinguish scenarios with substantial shifts from global bias to local interactions, the precise numerical values estimating the contribution of LI varies between scenarios with a low versus high global bias. To address this issue we propose to exploit the variation between experiments for modeling the relation between LI and GI. This is demonstrated by comparing two distinct values of the interaction parameter, α, emulating different experimental settings (Figure 2F). Within experiments variation was obtained by drawing multiple values of σ from a normal distribution. Due to the negative correlation between LI and GI the experimental patterns can only be discriminated by combining LI and GI into a twodimensional descriptor (Figure 2F). This point will be further demonstrated in one of the following case studies and in the Discussion.
Theory and limiting cases of DeBias
To characterize the properties of DeBias we used theoretical statistical reasoning. The first limiting case is set by the scenario in which observations from X and Y are independent. The expected values of the observed and resampled alignments are identical; accordingly, LI converges to 0 for large N (Appendix 1, Theorem 1). The second limiting case is set by the scenario in which X and Y are both uniform. The corresponding resampled alignment is also uniform; accordingly, GI converges to 0 for a large N (Appendix 1, Theorem 2). The third limiting case occurs with perfect alignment, i.e., x_{i} = y_{i} for all i. In this case the observed alignment distribution is concentrated in the bin containing θ = 0. We examine two scenarios of perfect alignment: (1) When all the locally matched measurements are identical (x_{i} = y_{j} for all i, j), the resampled distribution is also concentrated in the bin θ = 0 implying that LI = 0 and GI assumes the maximal possible value: $GI=\frac{1}{K}{\displaystyle \sum}_{i=1,\mathrm{..},K}\left(i1\right)=\frac{K1}{2}$, where K is the number of quantization bins (Appendix 1, Theorem 3.I). (2) When X, Y are uniform (and x_{i} = y_{i} for all i), the resampled distribution is uniform, thus GI = 0 and LI reaches its maximum value: $LI=\frac{1}{K}{\displaystyle \sum}_{i=1,\mathrm{..},K}\left(i1\right)=\frac{K1}{2}$, (Appendix 1, Theorem 3.II). Generalizing this case, we prove that LI is a lower bound for the actual contribution of the local interaction to the observed alignment (Appendix 1, Theorem 4). Complementarily, GI is an upper bound for the contribution of the global bias to the observed alignment.
Last, we show that when X and Y are truncated normal distributions, or when the alignment distribution is truncated normal, GI reduces to a limit of 0 as σ → ∞, when σ is the standard deviation of the normal distribution before truncation (Appendix 1, Theorem 5). Simulations complement this result demonstrating that σ and GI are negatively associated, i.e., GI decreases with increasing σ (Figure 2E). This final property is intuitive, because resampling from more biased distributions (smaller σ) tends to generate high agreement between (x_{i}, y_{i}) leading to reduced alignment angles and increased GI.
The modeling of the observed alignment as the sum of GI and LI allowed us to assess the performance of DeBias from synthetic data. By using a constant local interaction parameter ζ (ζ = c), we were able to retrieve the portion of the observed alignment that is attributed to the local interaction and to compare it with the true predefined ζ (Appendix 2, Appendix 2—figure 1). These simulations demonstrated again the need for a GIdependent interpretation of LI (first shown in Figure 2E–F). Simulations were also performed to assess how the choice of the quantization parameter K (i.e., number of histogram bins) and number of observations N affect GI and LI (Appendix 2, Appendix 2—figures 2–3), and this was also verified in our experimental data (Figure 3—figure supplement 1). In summary, by combining theoretical considerations and simulations we demonstrated the properties and limiting cases of DeBias in decoupling paired matching variables from orientation data.
Local alignment of vimentin and microtubule filaments
We applied DeBias to investigate the degree of alignment between vimentin intermediate filaments and microtubules in polarized cells. Recent work using genomeedited Retinal Pigment Epithelial (RPE) cells with endogenous levels of fluorescently tagged vimentin and αtubulin showed that vimentin provides a structural template for microtubule growth, and through this maintains cell polarity (Gan et al., 2016). The effect was strongest in cells at the wound front where both vimentin and microtubule networks collaboratively align with the direction of migration (Figure 3A–C). An open question remains as to how much of this alignment is caused by the extrinsic directional bias associated with the collective migration of cells into the wound as opposed to a local interaction between the two cytoskeleton systems.
Analysis of the GI and LI revealed that most of the discrepancy in vimentinmicrotubule alignment originated from a shift in the global bias (Figure 3D), suggesting that the local interaction between the two cytoskeletons is unaffected by the cell position or knockdown of vimentin. Instead, the reduced alignment between the two cytoskeletons is caused by a loss of cell polarity in cells away from the wound edge, probably associated with the reduced geometric constraints imposed by the wound edge. In a similar fashion, reduction of vimentin expression relaxes global cell polarity cues that tend to impose alignment.
To corroborate our conclusion that the global state of cell polarity is encoded by the GI, we performed a live cell imaging experiment, in which single cells at the edge of a freshly inflicted wound in a RPE monolayer were monitored for 80 min after scratching. DeBias was applied to calculate a time sequence of LI and GI. Cells at the wound edge tended to gradually increase their polarity and started migrating during the imaging time frame (Figure 3E, Video 1). Accordingly, the GI increased over time (Figure 3F). We also used this data set to verify that the reported shifts in GI are independent of the number of data points in and the binning of the distribution (Figure 3—figure supplement 1). This demonstrates the capacity of DeBias to distinguish fundamentally different effectors of cytoskeleton alignment.
Identifying molecular factors in alignment of cell velocity and mechanical forces during collective cell migration
Collective cell migration requires intercellular coordination, achieved by mechanical and chemical information transfer between cells. One mechanism for cellcell communication is plithotaxis, the tendency of individual cells to align their velocity with the maximum principal stress orientation (He et al., 2015; Tambe et al., 2011; Zaritsky et al., 2015). As in the previous example of vimentin and microtubule interaction, much of this alignment is associated with a general directionality of velocity and stress field parallel to the axis of collective migration (Zaritsky et al., 2015).
Using a wound healing assay, Das et al. (Das et al., 2015) screened 11 tightjunction proteins to identify pathways that promote motionstress alignment (Figure 4A). Knockdown of Merlin, Claudin1, Patj and Angiomotin (Amot) reduced the alignment of velocity direction and stress orientation (Das et al., 2015). Further inspection of these hits showed that the stress orientation remained stable upon depletion of these proteins, but the velocity direction distribution was much less biased towards the wound edge (Zaritsky et al., 2015). Here, we further analyze this data to demonstrate the capacity of DeBias to pinpoint tightjunction proteins that alter specifically the global or local components that induce velocitystress alignment.
By distinguishing GI and LI we generated a refined annotation of the functional alteration that depletion of these tightjunction components caused in mechanical coordination of collectively migrating cells (Figure 4B–C). First, we confirmed that the four hits reported by (Das et al., 2015) massively reduced the GI, consistent with the notion that absence of these proteins diminished the general alignment of velocity to the direction induced by the migrating sheet (Figure 4B, red dashed rectangle). Merlin, Patj and Angiomotin reduced the LI to values close to 0, suggesting that the local dependency between stress orientation and velocity direction was lost. Depletion of Claudin1, or of its paralog Claudin2, which was not reported as a hit in the Das et al. screen, reduced the LI to a lesser extent, similarly for both proteins, but had very different effects on the GI (Figure 4B, purple dashed rectangle). This suggested that the analysis by (Das et al., 2015) missed effects that do not alter the general alignment of stress or motion, and implied the existence of a local velocitystress alignment mechanism that does not immediately change the collective aspect of cell migration.
When assessing the marginal distributions of stress orientation and velocity direction we observed that depletion of Claudin1 reduced the organization of stress orientations and of velocity direction, while Claudin2 reduced only the latter. The LI values of depletion conditions were similar and lower than control (Figure 4D). Merlin depletion is characterized by an even lower LI and marginal distributions with aligned stress orientation and almost uniform alignment distribution (Figure 4D). Since we think that aligned stress is transformed to aligned motion (He et al., 2015; Zaritsky et al., 2015), we propose that in this data the LI quantifies the effect of local mechanical communication on parallelizing the velocity among neighboring cells. Accordingly, stressmotion transmission mechanism is impaired to a similar extent by reduction of Claudin1 and Claudin2, albeit less than by reduction of Merlin.
Using LI as a discriminative measure also allowed us to identify a group of new hits (Figure 4C). ZO1, ZO2, ZO3, Occludin and ZONAB are all characterized by small reductions in GI but a substantial increase in LI relative to control (Figure 4B, orange dashed rectangle). A quantitative comparison of control and ZO1 depleted cells provides a good example for the type of information DeBias can extract: both conditions yield similar observed alignment distributions with nearly identical means, yet ZO1 depletion has an 83% increase in LI and 8% reduction in GI, i.e., the mild loss in the marginal alignment of velocity or stress is compensated by enhanced local alignment in ZO1 depleted cells (Figure 4D). This might point to a mechanism, in which stress orientation is reduced by tightjunction depletion, but enhanced by transmission of stress orientation into motion orientation, leading to comparable alignment. Notably, all paralogs, ZO1, ZO2 and ZO3 fall into the same cluster of elevated LI and slightly reduced GI relative to control experiments. This phenotype is in agreement with the outcomes of a screen that found ZO1 depletion to increase both motility and celljunctional forces (Bazellières et al., 2015).
Using DeBias to assess proteinprotein colocalization
Proteinprotein colocalization is another ubiquitous example of correlating spatially matched variables in cell biology. To quantify GI and LI for proteinprotein colocalization, we normalized each channel to intensity values between 0 and 1. The 'alignment' θ_{i} of matched observations (x_{i},y_{i}) was replaced by the difference in normalized fluorescent intensities x_{i}  y_{i} (Materials and methods). Simulations demonstrated that stronger interactions in colocalization are translated to larger LI values and validated that the choice of K (number of histogram bins) and N (number of observations) marginally affect GI and LI (Appendix 3, Appendix 3—figure 1). While LI could serve as a measure to assess colocalization, the interpretation of GI is less intuitive. In the following, we present two examples of applying DeBias for proteinprotein colocalization, and demonstrate the type of information that can be extracted from the combined GI and LI analysis.
PKC FRET: a simple example of pixelbased proteinprotein colocalization
To test the potential of DeBias in quantification of pixelbased colocalization, we analyzed the effect of ﬂuorescence resonance energy transfer (FRET) in the C kinase activity reporter (CKAR), which reversibly responds to PKC activation and deactivation (Violin et al., 2003). Reduced PKC activity leads to energy transfer from CFP to YFP_{CFP}, resulting in reduced FRET ratio ($\frac{CFP}{YF{P}_{CFP}}$) (Figure 5A). Assuming that the CFP signal is dominant (CFP > YFP_{CFP}), this alteration should reduce the difference between the CFP and YFP_{CFP} channels, which would in DeBias yield an increased LI (Figure 5A, Materials and methods).
To test this we labeled hTERTRPE1 cells with CKAR and imaged CFP and YFP_{CFP} channels before and after specific inhibition of PKC with HA100 dihydrochloride (Figure 5B, Materials and methods), leading to reduced pixel differences in their normalized fluorescent intensities (Figure 5C). As expected, the $\frac{CFP}{YF{P}_{CFP}}$ ratio decreased (Figure 5D), LI values increased (Figure 5E) and seemed more sensitive to the FRET. Surprisingly, DeBias indicated a shift in the GI values (Figure 5F), reflected in a more homogeneous marginal distribution of both channels before inhibition (Figure 5G). Control experiments with cytoplasmic GFP and mCherry expression did not show the shifts observed in LI or GI (Figure 5H). Thus, we conclude that PKC inhibition changes the localization of PKC towards a more random spatial distribution. One possible mechanism for this behavior is that deactivation releases the kinase from the substrate. This example illustrates DeBias’ capabilities to simultaneously quantify changes in local interaction and global bias in pixelbased colocalization.
Inferring colocalization of molecular cargo and clathrincoated pits during endocytosis
Clathrinmediated endocytosis (CME) is the major pathway for entry of cargo receptors into eukaryotic cells. Cargo receptor composition plays an important role in regulating clathrincoated pit (CCP) initiation and maturation (Liu et al., 2010; Loerke et al., 2009). The clustering of transferrin receptors (TfnR), the classic cargo receptor used to study CME, promotes CCP initiation, in concert with clathrin and adaptor proteins (Liu et al., 2010). Recent evidence suggests a diversity of mechanisms regulating endocytic trafficking, including crosstalks between signaling receptors and components of the endocytic machinery (Di Fiore and von Zastrow, 2014). For example, the oncogenic protein kinase Akt has been shown to play an important role in mediating CME in cancer cells (Liberali et al., 2014; Reis et al., 2015), but not in normal epithelial cells (Reis et al., 2015). Here we tested how the decoupling by DeBias of global and local contributions to the overall intensity alignment of clathrin and TfnR, can be used to simultaneously investigate colocalization and predict CCP dynamics, using fixed cell fluorescence imaging.
We used fluorescence images of fixed non‐small lung cancer cells (H1299) or untransformed human retinal pigment epithelial cells (ARPE19) expressing clathrin light chain A fused to eGFP (eGFPCLCa) as a CCP marker (Figure 6A–B). Cells were either treated with DMSO or with an AKT inhibitor (Akt inhibitor X, ‘ten’), and imaged by Total Internal Reﬂection Fluorescence Microscopy (TIRFM). CCPs were reported in the eGFPCLCa channel and TfnR was visualized by immunofluorescence in a second channel (Materials and methods). For single cells, the location of fluorescent signals of CLCa and TfnR were recorded and the data were pooled and processed by DeBias (Materials and methods).
LI values, indicative of the colocalization between TfnR and CLCa, were significantly lower in Aktinhibited H1299 compared to control cells (Figure 6C). In contrast, Akt inhibition increased the LI values in ARPE19 cells but this effect was less prominent (Figure 6D). Akt inhibition resulted in increased GI values for both cell lines, to a much greater degree in H1299 cells (Figure 6C–D). To test whether GI enhances the ability to distinguish between control and Aktinhibited cells, we applied Linear Discriminative Analysis (LDA) classification to calculate the true positive rate versus the false positive rate for LI alone (black lines, Figure 6E–F) or the pair (GI, LI), (orange lines, Figure 6E–F). The area under these curves (AUC) provided a direct measure of the ability of each method to accurately classify the experimental condition of single cells. AUC for the (GI, LI) representation was superior to using LI alone for both cell lines (H1299: 0.96 versus 0.88, Figure 6E; ARPE19: 0.83 versus 0.72, Figure 6F). Similar benefit in discrimination was achieved when using the GI to complement Pearson’s correlation as an alternative to LI for measuring local interaction (Figure 6—figure supplement 1A–F). Such improved discrimination is indicative of distinct molecular processes that were altered upon Akt inhibition. We also used this data set to experimentally validate the independence of GI and LI of the number of observations (N, Figure 6—figure supplement 1G,H) and the choice of the number of histogram bins (K, Figure 6—figure supplement 1I,J).
To interpret the increased GI values for Aktinhibited cells, we examined the joint and marginal distributions of CLCa and TfnR. Upon Aktinhibition, the joint distributions were more biased toward regions of low TfnR intensities (Figure 6G–H). This was clearly observed in the marginal distributions (Figure 6I–J). Hence, although the CLCa distribution appeared not to change upon AKT inhibition, the frequency of CCPs with fewer TfnRs increased. Given the positive relation between TfnR cargo quantities and CCPs maturation (Loerke et al., 2009), we wondered whether the increased frequencies of CCPs containing less TfnR might alter CCPs dynamics. Indeed, liveimaging of H1299 cells showed a higher frequency of CCPs with shorter lifetimes upon Aktinhibition, which was not seen in normal ARPE19 cells (Figure 6K–L, Materials and methods). It has previously been shown that Akt inhibition reduces the rate of TfnR CME uptake in H1299 cells, but not in ARPE19 cells ([Reis et al., 2015], see also Figure 6M–N); therefore, these findings indicate that the reduced levels of TfnR in CCPs upon Akt inhibition results in an increase in shortlived, most likely abortive events, and hence a decrease in CME efficiency.
Altogether, DeBias could distinguish alterations in the regulation of CME between two cell types. The decoupling to GI and LI indicated that upon Akt inhibition, both untransformed and cancer cells showed a global bias towards CCPs with lowered TfnR intensities. This conclusion could not have been reached by considering only the LI, which increased for normal and decreased for transformed cell lines.
Discussion
We introduce DeBias as a new method to assess global bias and local interactions between coupled cellular variables. Although the method is generic, we show here specific examples of DeBias analysis in coorientation and colocalization studies. The source code is available, https://github.com/DanuserLab/DeBias, as well as via a webbased platform, https://debias.biohpc.swmed.edu. The website also provides detailed instructions for the operation of the user interface.
DeBias defines a generalizable framework for eliminating confounding factors in the analysis of interacting variables. Our examples demonstrate that the distinction of global and local contributions to the level of variable coupling eliminates much of the global confounder bias in the analysis of more direct interactions and can unearth in the form of global bias mechanisms that are missed by a single parameter analysis (Figures 1–2). In the example of vimentinmicrotubule alignment (Figure 3), the significant decrease in GI as opposed to the LI upon partial vimentin knockdown indicated that the reduction in alignment between the two cytoskeleton systems is associated with a reduction of cell polarity as the global cue. In the example of stressvelocity alignment (Figure 4), depletion of some tight junction proteins increased LI, suggestive of enhanced local stressmotion transmission; knockdown of others decreased GI indicating an overall impaired alignment of velocity in the direction of wound closure. In the example of FRET experiments (Figure 5), PKC inhibition lead to increased LI, validating the FRET response, while a reduced GI was indicative of weaker interactions of the inactivated kinase with its substrates. In the example of Tfn receptor (TfnR) colocalization with CCPs during CME (Figure 6), the increased GI in response to Akt inhibition related to a higher fraction of CCPs containing less TfnR. Moreover, Akt inhibition induced opposite shifts in LI for normal and cancer cells, reflecting differential alterations in colocalization between cell types. Thus, DeBias provided insight into the regulation of cargopit association by kinase activity that depended on a proper deconvolution of local and global effects on the interaction of the clathrin and receptor signal. We then validated our conclusions by further analyses of the marginal distributions, liveimaging and uptake assays (Figure 6). Overall, the four applications shown in this work first emphasize the general need for a confounder analysis when dealing with coupled biological variables and second indicate that the global bias may be linked to mechanistically meaningful properties of the studied system. These properties were either ignored or eliminated by previous methods, and now can be assessed directly by DeBias.
Other approaches have been used to address global confounders for assessment of local interactions between biological variables. For the specific example of objectbased colocalization, Helmuth et al. simulated the spatial distribution of objects in the absence of local interactions to calibrate colocalization measurement (Helmuth et al., 2010). Other methods mostly used secondorder spatial statistics on distances between neighbor points to exclude confounders for better colocalization sensitivity (reviewed in Lagache et al., 2015). Importantly, we show applications of DeBias on colocalization that do not require initial object detection (Figure 5). While the phenomenon of confounder bias is independent of object versus pixelbased colocalization, we distinguish the peculiarities of the two scenarios in Appendix 4.
An important and more general approach to revealing local interactions masked by global biases was recently proposed by (Krishnaswamy et al., 2014), using applications to single cell mass cytometry data as examples. The authors developed a measure referred to as conditionalDensity Resampled Estimate of Mutual Information (DREMI) to quantify the influence of a protein X on protein Y based on the conditional probability P(YX). DREMI takes advantage of the abundant mass cytometry data to equally weigh data at different intervals along the range of X values using > 10,000 cells per experimental condition. This approach is less reliable when limited data is available, because of the low confidence in the conditional probability of observations with low data abundance. Thus, DREMI is not well suited for image data, which typically has fewer observations.
DeBias estimates GI and LI assuming a constant global bias and local interaction for all observations. Moreover, its quantification power is relative. For example, a twofold increase in the direct interaction of two variables would not necessarily result in a twofold increase in LI. Another limitation is the absence of complete orthogonality of GI and LI values (Figure 2E–F, Appendix 3—figure 1), which complicates the interpretation of GI and LI in certain scenarios. These three limitations apply to all current approaches for quantifying interactions between coupled variables. The main conceptual advance DeBias seeks to make relates to the explicit integration of confounding factors in the analysis of coupled variables, which implies an expansion of the coupling metric from a scalar to a twodimensional score. A forth limitation in the current implementation of DeBias is the linear normalization of multiple intensity variables in colocalization applications. Future versions may include nonlinear normalization methods, although such normalization is usually highly specific to a particular data set. Last, the mechanism encoded by the GI is not always obvious. Sometimes it requires additional experiments to unveil the information contained by the GI. For example, we combined fixed cell dualcolor imaging with liveimaging and uptake assays to show that shifts in the GI encode a shift in the relative populations of short and longlived CCPs between conditions (Figure 6). Despite some of the discussed complexities, DeBias offers a simple means to quantify and interpret mechanisms that alter confounders in the coupling of two variables and to largely exclude such global biases from the quantification of direct interactions.
Materials and methods
DeBias procedure
Request a detailed protocolThe DeBias procedure is depicted in Figure 2A. The marginal distributions X and Y are estimated from the experimental data, $\forall i,{x}_{i},{y}_{i}\in \left[0,90\xb0\right]$. The experimentally observed alignment distribution (denoted observed) is calculated from the alignment angles θ_{i} of matched (x_{i},y_{i}) paired variables, for all i.
The resampled alignment distribution (denoted resampled) is constructed by independent sampling from X and Y. N random observations (where N = X is the original sample size) from X and Y are independently sampled with replacement, arbitrarily matched and their alignment angles calculated to define the resampled alignment. This type of resampling precludes the local dependencies between the originally matched (x_{i},y_{i}) paired variables.
The uniform alignment distribution (denoted uniform) is used as a baseline for comparison between distributions. This is the expected alignment distribution when neither global bias (reflected by uniform X, Y distributions) nor local interactions exist. The Earth Mover's Distance (EMD) (Peleg et al., 1989; Rubner et al., 2000) was used as a distance metric between alignment distributions. The EMD for two distributions, A and B, is defined as follows:
$EMD\left(A,B\right)=\sum _{j=1,\dots ,K}\sum _{j=1,\dots ,i}{a}_{j}\sum _{j=1,\dots ,i}{b}_{j}$, where ${a}_{j}$ and ${b}_{j}$ are the frequencies of observations in bins $j$ of the histograms of distributions A and B, respectively, each containing K bins.
The global index (GI) is defined as the EMD between the uniform distribution and the resampled alignment:
GI = EMD(uniform,resampled)
The local index is determined by subtraction of the global index from the EMD between the uniform distribution and the experimentally observed alignment distribution:
LI = EMD(uniform,observed)  global index
DeBias for proteinprotein colocalization
Request a detailed protocolThe following adjustments to this procedure are implemented to allow DeBias to quantify proteinprotein colocalization:
Levels of fluorescence are not comparable between different channels due to different expression levels and imaging parameters. Thus, each channel is normalized to [0,1] by the fifth and 95th percentiles of the corresponding fluorescence intensities to achieve a stable and robust distribution.
The alignment angle θ_{i} of the matched observation (x_{i},y_{i}) is calculated as the difference in normalized fluorescence intensities x_{i}  y_{i} and the alignment distribution is thus defined on the interval [−1,1].
The number of histogram bins for the alignment distributions (observed, resampled and uniform) was K = 15 for orientational data, 19 for PKC and 40 for CME colocalization data.
Automated selection of number of histogram bins, K
Request a detailed protocolThe FreedmanDiaconis rule (Freedman and Diaconis, 1981) was used to automate the selection of histogram bin width: $bin\text{}size=\frac{{Q}_{3}\left(x\right){Q}_{1}\left(x\right)}{\sqrt[3]{n}}$, where Q_{i} is the i^{th} quartile of the empirical distribution x and n is the number of data points contained. A function to calculate K is included in our publicly available source code and this functionality was also integrated to the webserver implementation. Importantly, GI and LI across experiments can be compared only when evaluated with the same K value and this is enforced by the webserver. It is the responsibility of the sourcecode user to validate using the same K values when comparing different experimental conditions.
Simulating synthetic data
Simulating coalignment data
Request a detailed protocolLet us define $X$ and $Y$ as the angular probability distribution functions, with angle instances denoted ${x}_{i}$ and ${y}_{i}$, respectively. When simulating local relations, for each pair of angles, one of the angles will be shifted towards the other by $\zeta $ degrees (Figure 1B), unless $\left{x}_{i}{y}_{i}\right<\zeta $, in which case it will be shifted by ${x}_{i}{y}_{i}$ degrees. The angle to be shifted (either ${x}_{i}$ or ${y}_{i}$) is chosen by a Bernoulli random variable, $p$, with probability 0.5. The observed angles for pixel $i$ will therefore be
and
The alignment of angles at pixel $i$ will be:
For example, for our simulations we choose $X,\text{}Y$ to be truncated normal distributions on $(90,90)$ with $\mu =0$ and varying values of $\sigma $.
$\zeta $ is modeled in two ways: either as a constant value, e.g. $\zeta =5\xb0$ (Appendix 2—figures 1–3), or as a varying value dependent on $\left{x}_{i}{y}_{i}\right$ (Figures 1–2). For the latter, $\zeta $ is defined as a fraction $0<\alpha <1$ from ${x}_{i}{y}_{i}$ for each pair of observations; namely, ${\zeta}_{i}=\alpha {x}_{i}{y}_{i}$ (see Figure 1B). Note, that the observed marginal distributions X’, Y’ may be slightly different from X, Y.
Simulating colocalization data
Request a detailed protocolLet us define $X$ as a probability distribution function, with instances denoted ${x}_{i}$. Local interactions were simulated as ${y}_{i}={x}_{i}{\zeta}_{i}$, where ${\zeta}_{i}$ is an instance of a probability distribution $Z$. For our simulations we chose $X$ to be a truncated normal distribution on $[0,1]$ with μ_{ x} = 0.5 and $Z$ to be a normal distribution with μ_{ζ} = 1. This model of interaction assumes on average a onetoone interaction between X and Y, deviation of ${\zeta}_{i}$ from one implies reduced interaction. ${y}_{i}$ samples were also truncated to [0,1]. This ensures that $\forall i,0\le {y}_{i},{x}_{i}\le 1$ and accordingly, $\forall i,1\le {x}_{i}{y}_{i},\le 1$, making unnecessary the normalization step in DeBias colocalization calculation.
When simulating scenarios where only subgroups of the observations undergo interactions, we sampled the noneinteracting observations ${y}_{i}$ from Y the truncated normal distribution on $(0,1)$, with μ_{y} = 0.5 and σ_{y} = σ_{x} (same as X).
Vimentin and microtubule filaments experiments and analysis
Cell model
Request a detailed protocolhTERTRPE1 cells (ATCC, RRID: CVCL_4388) were TALENgenome edited to endogenously label vimentin with mEmerald and αtubulin with mTagRFPt, and validated for protein expression levels (Gan et al., 2016). Cells were stably transfected with shRNA against vimentin to knock down vimentin and the knockdown efficiency was validated as ~75% (Gan et al., 2016). The cell line has been tested negative for mycoplasma contamination.
Fixed cell imaging
Request a detailed protocolhTERTRPE1 mEmeraldvimentin/mTagRFPtαtubulin cells expressing shRNAVIM or scrambled control shRNA Scr were plated into MatTek (Ashland, MA) 35 mm glassbottom dishes (P35G0–20 C) coated with 5 µg/mL fibronectin. Cells were incubated overnight to allow them to adhere and form monolayers. Monolayers were scratched with a pipette tip to form a wound. Cells were incubated for 90 min, washed briefly and fixed with methanol at −20°C for 15 min. Cells were imaged at the wound edge (denoted ‘front’ cells), and at 2–3 cell rows from the wound edge (denoted ‘back’ cells, only for control condition). Images were acquired using a Nikon Eclipse Ti microscope, equipped with a Nikon Plan Apo Lambda 100x/1.45 N.A. objective. Images were recorded with a Hamamatsu ORCA Flash 4.0 with 6.45 μm pixel size (physical pixel size: 0.0645 × 0.0645 μm). All microscope components were controlled by Micromanager software.
Live cell imaging
Request a detailed protocolhTERTRPE1 mEmeraldvimentin/mTagRFPtαtubulin cells expressing scrambled control shRNA Scr were plated into MatTek (Ashland, MA) 35 mm glassbottom dishes (P35G0–20 C) coated with 5 µg/mL fibronectin. Cells were incubated overnight to allow them to adhere and form monolayers. Monolayers were scratched with a pipette tip to form a wound. Imaging started 30 min after scratching with an Andor Revolution XD spinning disk microscope mounted on a Nikon Eclipse Ti stand equipped with Perfect Focus, a Nikon Apo 60 × 1.49 N.A. oil objective and a 1.5x optovar for further magnification. Images were recorded with an Andor IXON Ultra EMCCD camera with 16 μm pixel size (physical pixel size: 0.16 × 0.16 μm). Lasers with 488 nm and 561 nm light emission were used for exciting mEmerald and mTagRFPt, respectively. The output powers of the 488 nm and 561 nm lasers were set to 10% and 20% of the maximal output (37 mW and 23 mW, respectively). The exposure time was 300 ms per frame for both channels and images were collected at a frame rate of 1 frame per minute. During acquisition, cells were kept in an onboard environmental control chamber. All microscope components were controlled by Metamorph software.
Filament extraction and spatial matching
Request a detailed protocolWe applied the filament reconstruction algorithm reported in (Gan et al., 2016). Briefly, multiscale steerable filtering is used to enhance curvilinear image structures, centerlines of candidate filament fragments are detected, clustered to high and low confidence sets and iterative graph matching is applied to connect fragments into complete filaments. Each filament is represented by an ordered chain of pixels and the local filament orientation derived from the steerable filter response. Spatial matching was performed as follows: each pixel belonging to a filament detected in the MT channel is recorded to the closest pixel that belongs to a filament in the VIM channel. If the distance between the two pixels is less than 20 pixels, then the pair of VIM and MT orientations at this pixel is recorded for analysis. The same process is repeated to record matched pixels from VIM to MT filaments.
Collective cell migration experiments and analysis
Coupled measurements of velocity direction and stress orientation were taken from the data originally published by Tamal Das et al. (Das et al., 2015). Particle image velocimetry (PIV) was applied to calculate velocity vectors, and monolayer stress microscopy (Tambe et al., 2011) was used to extract stress orientations. Velocity and stress measurements were recorded 3 hr after collective migration was induced by lifting off the cultureinsert in which the cells have grown to confluence. Validated siRNAs were used for gene screening. Detailed experimental settings can be found in Das et al., 2015.
Statistical test
Request a detailed protocolWe devised a permutation test to determine statistical significance of differences in LI values between experiments (conditions) (Figure 4C). For every pair of conditions (i,j), the following procedure was repeated for 100 iterations. 50% of the velocitystress observations were randomly selected for each condition and the LI (and GI) were calculated from this subsampling. Without loss of generality, for LI_{i} < LI_{j} (based on Figure 4B) the pvalue was recorded as the fraction of iterations in which the subsampled LI value for condition i was greater than the LI value for condition j. A fraction of 0 thus implies pvalue < 0.01.
PKC FRET experiments
Request a detailed protocolhTERTRPE1 cells (ATCC, RRID: CVCL_4388) expressing GFP and mCherry (for the control experiment) or C kinase activity reporter (CKAR, for the PKC activation experiment) (Violin et al., 2003) were plated with DMEM/F12 medium containing 10% fetal bovine serum and 1% penicillinstreptomycin in fibronectincoated 35 mm MatTek plates (P35G0–10 c). The cell line has been tested negative for mycoplasma contamination. Cells were incubated overnight and imaged with a custombuilt DeltaVision OMX widefield microscope (GE healthcare life sciences, Chicago, IL) equipped with an Olympus PLAN 60 × 1.42 N.A. oil objective and CoolSNAP HQ2 interline CCD cameras. FRET experiments were performed with a 445 nm laser, and control experiments were performed with 488 and 561 nm lasers. 478/35, 541/22, 528/48 and 609/37 emission bandpass filters were used for the CFP, YFP, GFP and mCherry channels, respectively. The output powers of the lasers were set to 10% the maximal output (100 mW). The exposure time was 200 ms per frame for both channels. 10 µM of the PKC inhibitor HA100 dihydrochloride (Santa Cruz Biotechnology, Dallas, TX) was added to the media after the first image was recorded and a second image was recorded 20 min later.
Single cells were manually selected for analysis. In particular, cells with higher intensities in the CFP channel were found to provide reproducible changes in their FRET intensity. Each cell was manually annotated and analyzed with the Biosensor Processing Software 2.1 to produce the ratio images (Hodgson et al., 2010). Briefly, the field of view was corrected for uneven illumination, background was subtracted, the image was masked with the single cell annotation, and the ratio image was calculated as CFP/YFP_{CFP}. Statistics was determined using the nonparametric Wilcoxon signedrank test.
Clathrin mediated endocytosis experiments
Cells, cell culture and chemicals
Request a detailed protocolARPE19 (retinal pigment epithelial cells) (ATCC, RRID: CVCL_0145) stably expressing eGFPCLCa were grown under 5% CO2 at 37°C in DMEM high glucose medium (Life Technologies, ), supplemented with 20 mM HEPES, 10 mg/ml streptomycin, 66 ug/ml penicillin and 10% (v/v) fetal calf serum (FCS, HyClone). H1299 (nonsmall cell lung carcinoma) (RRID: CVCL_0060, a generous gift from Dr. J. Minna at the UT Southwestern Medical Center) stably expressing eGFPCLCa were grown under 5% CO_{2} at 37°C in RPMI, supplemented with 20 mM HEPES, 10 mg/ml streptomycin, 66 ug/ml penicillin and 5% (v/v) fetal calf serum (FCS, HyClone). STR profiling was performed to ensure cell identity. No mycoplasma contamination was found. The AKT inhibitor (Akt inhibitor X, ‘ten’) was purchased from Calbiochem.
Transferrin receptor internalization
Request a detailed protocolTfnR uptake was measured by an ‘incell’ ELISA assay using the antiTfnR monoclonal antibody, HTRD65 (Schmid and Smythe, 1991), as ligand, exactly as previously described (Reis et al., 2015). Internalized D65 was expressed as the percentage of the total surfacebound D65 at 4°C (i.e., without acid wash step), measured in parallel.
Fixed cell imaging
Request a detailed protocolTransferrin receptor (TfnR) surface levels were measured using the antiTfnR mAb (HTRD65). ARPE19 and H1299 cells (1.5 × 105 cells per well in a 6well plate) were grown overnight on glass cover slips and further preincubated with 4 ug/ml of D65 in TfnR assay buffer (PBS4+: PBS supplemented with 1 mM MgCl2, 1 mM CaCl2, 5 mM glucose and 0.2% bovine serum albumin) at 4°C for 30 min. After being washed with PBS^{4+}, cells were fixed in 4% PFA for 30 min at 37°C, permeabilized with 0.1% Triton X100 for 5 min and further blocked with QPBS (2% BSA, 0.1% lysine, pH 7.4) for 30 min. After three washes with PBS, cells were incubated with a 1:500 dilution of goat antimouse Alexa568 labelled secondary antibody (Life Technologies) for 30 min, washed an additional three times with PBS before TIRFM imaging using a 100 × 1.49 NA Apo TIRF objective (Nikon, Japan) mounted on a TiEclipse inverted microscope equipped with the Perfect Focus System (Nikon). Images were acquired with an exposure time of 150 ms for both channels using a pcoedge 5.5 sCMOS camera with 6.5 um pixel size. For inhibition studies, cells were initially preincubated in the presence of Akt inhibitor X (10 uM) for 30 min at 37°C, followed by preincubation with 4 ug/ml of D65 at 4°C for 30 min, in continued presence of the inhibitor.
Livecell imaging and analysis
Request a detailed protocolDuring TIRFM imaging, cells were maintained in DMEM lacking phenol red and supplemented with 2.5% fetal calf serum. Timelapse image sequences were acquired at a frame rate of 1 frame/s with exposure time of 150 ms using a pcoedge 5.5 sCMOS camera with 6.5 um pixel size. CCP detection, tracking and construction of lifetime distributions were performed with the custom CME analysis software described in (Aguet et al., 2013). Lifetime distribution was defined at the resolution of 1 s and limited to 160 s. Longer CCP trajectories were excluded from the analysis. To compare lifetime distributions for single field of views (FOVs) between WT and AKTinhibited cells we measured the heterogeneity as the EMD distance to the uniform distribution. FOV’s score were compared between the different experimental conditions using the nonparametric Wilcoxon ranksum test.
Image analysis for fixed cell experiments
Request a detailed protocolSingle cell masks were manually annotated in each FOV. We applied the approach described in (Aguet et al., 2013) to automatically detect CCPs from the CLC channel. Briefly, CLC ﬂuorescence was modeled as a twodimensional Gaussian approximation of the microscope PSF above a spatially varying local background. CCP candidates were first detected via filtering, followed by a modelfitting for subpixel localization. The fluorescent intensity of the CLC and any other acquired channel were recorded in the detection coordinates to define the matched observations for DeBias. GI and LI were calculated independently for each single cell. Linear Discriminant Analysis (LDA) classification (Fisher, 1936) was applied to assess single cell classification accuracy. Every cell constituted an observation, a label was assigned based on the experimental condition and the representation was either the LI or the pair (GI,LI). The LDA classifier was trained on a labeled dataset consisting of WT and AKTinhibition for H1299 or ARPE19 cells. The area under the Receiver Operating Characteristic (ROC) curve was recorded to assess and compare the discriminative accuracy of different measures. The truepositive rate (TPR) is the percentage of control cells classified correctly. The falsepositive rate (FPR) is the percent of AKT inhibited cells classified as control. When comparing the potential accuracy of several classification algorithms, a measure that has higher truepositive rate for any fixed falsepositive rate values is proved to be the better one. Thus, higher curves (larger areas under the ROC curve, or AUC) correspond to more discriminative measures. Statistical significance for comparing classification performance of LDA classifiers that were trained for scalar measures with or without the GI was calculated by bootstrapping (Figure 6E–F). The following process was repeated 1000 times and the frequency for which the scalarbased classifier outperformed the classifier trained on pairs of measures was reported as the pvalue. Random resampling with replacement was performed to obtain a sample size identical to that of the observed dataset. The area under the curve (AUC) of the competing pretrained LDA classifiers was assessed for this resampled dataset and recorded when the model that was trained without the GI predicted better.
Webserver
Request a detailed protocolThe DeBias code was implemented in Matlab, compiled with Matlab complier SDK and transferred to a webbased platform to allow public access for all users at https://debias.biohpc.swmed.edu. The graphical user interface (GUI) was designed to be simple and easy to use. The user uploads one or more datasets to the DeBias webserver and selects the mode of operation (colocalization/orientation). GI and LI values are calculated and the results are displayed and emailed to the user. ‘DeBias Analyst’ enables to group experiments into two experimental conditions (usually control versus treatment), visualizes and outputs statistics on the alterations of GI and LI. The software’s flow chart and a detailed user manual are available in the online user manual. Source code is publicly available, https://github.com/DanuserLab/DeBias Danuser, 2017 (with a copy archived at https://github.com/elifesciencespublications/DeBias).
Appendix 1
Theoretical results for coorientation data
Terms and definitions
Let X, Y be the distribution functions of two random variables representing angles on $[90\xb0,90\xb0]$. Spatially matched random variables from these distributions are denoted ${x}_{i}$ and ${y}_{i}$, $i=1,2\dots N$, where $N$ is the number of observations. ${x}_{i}^{*}$ and ${y}_{i}^{*}$ are random variables sampled from X and Y independently (without considering the spatial matching). The observed and resampled alignment distributions are denoted $A$ and ${A}^{*}$, respectively. The alignment distributions represent angles on $[0\xb0,90\xb0]$, and are functions of $X,Y$, the interaction between $X$ and $Y$, and $N$. Random variables from $A$ are denoted ${\theta}_{i}$, $i=1,2\dots N$. Random variables from ${A}^{*}$are denoted ${\theta}_{i}^{*}$. Let $K$ be the number of bins in the alignment histogram. Histogram bins are denoted $bi{n}_{i},i=0,\dots ,K1$ where $bi{n}_{0}$ contains the lowest values (including 0°) and $bi{n}_{k1}$ the highest values (including 90°). $U$denotes the uniform distribution on $[0\xb0,90\xb0]$ with the same K bins as $A,{A}^{*}$.
Theorem 1: Local index of independent variables
If $X,Y$ are independent, then $E\left(LI(X,Y)\right)=0$
Proof:
* Resampling does not change the expectation of the difference of two independent variables.
Theorem 2: Global index of uniform distributions
Let $X,Y$ be uniform distributions, then ${\mathrm{l}\mathrm{i}\mathrm{m}}_{N\to \mathrm{\infty}}GI=0$.
Proof:
Since $X$ and $Y$ are uniform distributions, the density of the resampled distributions are ${f}_{{x}_{i}^{\ast}}\left(\omega \right)={f}_{{y}_{i}^{\ast}}\left(\omega \right)=\frac{1}{180},90\text{}\text{}\omega \text{}\text{}90$. We can compute the distribution of the difference by $\begin{array}{r}{f}_{{x}_{i}^{\ast}{y}_{i}^{\ast}}\left(\omega \right)=\underset{\mathrm{\infty}}{\overset{\mathrm{\infty}}{\int}}{f}_{{x}_{i}^{\ast}}\left(x\right){f}_{{y}_{i}^{\ast}}\left(x\omega \right)dx=\underset{\mathrm{\infty}}{\overset{\mathrm{\infty}}{\int}}\frac{1}{{180}^{2}}{I}_{x\in \left(90,90\right)}{I}_{x\omega \in \left(90,90\right)}dx=\\ {\displaystyle \{\begin{array}{cc}\frac{1}{{180}^{2}}\underset{90}{\overset{90+\omega}{\int}}1dx& 0<\omega <180\\ \frac{1}{{180}^{2}}\underset{90+\omega}{\overset{90}{\int}}1dx& 0<\omega <180\end{array}=\{\begin{array}{cc}\frac{1}{{180}^{2}}(180+\theta )& 180<\omega <0\phantom{\underset{90}{\overset{90+\omega}{\int}}1dx}\\ \frac{1}{{180}^{2}}(180\theta )& 0<\omega <180\phantom{\underset{90}{\overset{90+\omega}{\int}}1dx}\end{array}}\end{array}$
We conclude that the distribution of the difference is
Therefore, when taking the absolute value of the angle difference we get:
Finally, since the alignment is limited to $[{0}^{\circ},{90}^{\circ}]$ we apply the function $\displaystyle g\left(\omega \right)=\{\begin{array}{lc}180\omega & 90\text{}\text{}\omega \le 180\\ \omega & 0\text{}\text{}\omega \le 90\end{array}$ so that
Therefore ${f}_{g\left(\left{x}_{i}^{*}{y}_{j}^{*}\right\right)}={A}^{*}$ is a uniform distribution.
For $N$ sufficiently large, the histogram has approximately $\frac{1}{K}$ of the observations in each of the $K$ equally spaced intervals between 0 and 90 and thus ${lim}_{N\to \mathrm{\infty}}\text{}GI={lim}_{N\to \mathrm{\infty}}\text{}EMD\left({A}^{\ast},U\right)=0$.
Theorem 3: Perfect alignment
If $\forall i,j,{x}_{i}={y}_{j}$ then $GI=\frac{\left(K1\right)}{2},\text{}\underset{N\to \mathrm{\infty}}{lim}LI=0$
If $X$ and $Y$ are uniform distributions, and $\forall i,{x}_{i}={y}_{i}$ then $\underset{N\to \infty}{\text{lim}}GI=0,\underset{N\to \infty}{\text{lim}}LI=\frac{\left(K1\right)}{2}$
Proof:
(I)
$A={A}^{*}$ because $\forall i,j{a}_{i}={a}_{j}^{*}=0$. For large $N$ the random variable drawn from the alignment distribution $A$ will be approximately:
The EMD of the alignment from the uniform distribution is therefore simply 'moving' $\frac{1}{K}$ observations from $bi{n}_{0}$ to every other bin, which sums up to
Therefore, for large $N$, $LI=EMD\left(A,U\right)EMD\left({A}^{*},U\right)=0$, $GI=EMD\left({A}^{*},U\right)=\frac{K1}{2}$.
(II)
Since $\forall i,{x}_{i}={y}_{i}$ we get that, similarly to part (I), $EMD\left(A,U\right)=\frac{K1}{2}$.
On the other hand, since $X$ and $Y$ are uniform distributions, we get from theorem 2 that ${lim}_{N\to \mathrm{\infty}}\text{}EMD\left({A}^{\ast},U\right)=0$.
Therefore, for infinite observations, $LI=\frac{K1}{2}$, $GI=0$.
Theorem 4: $LI$ is a lower bound for the local contribution to the observed alignment
Assuming that the observed alignment distribution $A$ is cumulatively explained by a global bias and a local interaction, we construct a new alignment distribution ${A}_{\zeta}$ encoding the true cumulative local contribution to the observed alignment and demonstrate that $LI\le EMD\left(U,{A}_{\zeta}\right)GI$ to conclude that LI is a lower bound for the local contribution to the observed alignment.
Proof:
We first define ${A}^{}$, the alignment distribution corresponding to $A$ that does not include any local interaction. Thus, ${A}^{}$, can be interpreted as an alignment distribution constructed from ${X}^{}$and ${Y}^{}$, denoting $X$ and $Y$ after elimination of the (unknown) alignment correction due to local interactions between the observations $({x}_{i},{y}_{i})$. The construction of $A}_{\zeta$ is based on the corresponding matching pairs $({x}_{i}^{}\in {X}^{},{y}_{i}^{}\in {Y}^{})$ with alignment correction by the local interaction ${\zeta}_{i}$ (see Figure 1B for as a schematic depiction). Such local interaction exists in our model (although it might not be explicitly known) and can be represented as a vector $\zeta \in {\mathbb{R}}^{N},{\zeta}_{i}\ge 0\forall i$. Note, that this construction supports different ${\zeta}_{i}$ values for every observation $i$ and thus can provide a more detailed platform than the single measure LI that DeBias outputs (which assumes ${\zeta}_{i}={\zeta}_{j}\forall i,j$). Also note, that when ${\zeta}_{i}>{\theta}_{i}^{}$ (${\theta}_{i}^{}$ is the alignment angle between $({x}_{i}^{},{y}_{i}^{}$)), then the observed alignment ${\theta}_{i}^{}{\zeta}_{i}<0$.
Accordingly, ${A}_{\zeta}$ is defined as the alignment distribution of ${\theta}_{i}^{}{\zeta}_{i}$. As described above, ${A}_{\zeta}$ can contain negative values for ${\zeta}_{i}>{\theta}_{i}^{}$. A, the experimentally observed alignment, thus can be generated from ${A}_{\zeta}$ as well, by truncating the ‘saturated’ observations (where ${\zeta}_{i}>{\theta}_{i}^{}$) to the value 0. More formally, the elements in A are defined by
We can get an upper bound for $EMD\left(A,{A}^{}\right)$ in the form of:
$\begin{array}{c}EMD\left(A,{A}^{}\right)\le EMD\left({A}_{\zeta},\text{}{A}^{}\right)\le \sum _{i=1}^{N}\frac{1}{N}\lceil \frac{{\zeta}_{i}}{bin}\rceil \end{array}$.
Where $\leftbin\right$ defines the size of the angular interval of a bin in the alignment histogram.
This equation is intuitively interpreted as every observation i is locally aligned by ${\zeta}_{i}$, and therefore is translocated $\lceil \frac{{\zeta}_{i}}{bin}\rceil$ bins, at most.
Note that a decreased bin size reduces this bound as close as needed to the value of $EMD\left({A}_{\zeta},{A}^{}\right)$.
Finally,
Thus, the LI is a lower bound on the contribution of the direct interaction between $X$ and $Y$ on the alignment distribution.
Additionally, we get that
Implying that the GI is an upper bound of the contribution of the global bias.
* by corollary 2
Corollary 2: For any alignment distribution A, $LI\le EMD\left(A,{A}^{*}\right)$
Proof:
Let ${A}_{i},{A}_{i}^{*},{U}_{i}$ denote the relative frequency of observations in $bi{n}_{i},$$0\le i\le k1$ for $A,{A}^{*},U$, respectively.
* triangle inequality
Theorem 5: GI limits for highly variant truncated normal distributions
${lim}_{\genfrac{}{}{0ex}{}{\sigma \to \mathrm{\infty}}{N\to \mathrm{\infty}}}GI=0$ for the following scenarios:
The resampled alignment is a truncated normal distribution with variance parameter ${\sigma}^{2}$.
$X$ and $Y$ are truncated normal distributions, each with variance parameter ${\sigma}^{2}$.
Proof:
(I)
Let ${A}_{\sigma}^{*}$ be the truncated normal resampled alignment distribution, defined by the parameters $\mu =0$, $\sigma $, with the support interval $(a,b)$, such that $a\le \mu \le b$. Let ${\varphi}^{\sigma}\left(x\right)$be the probability density function (PDF) of the truncated normal distribution and $u\left(x\right)$, $U$ respectively, the PDF and CDF (cumulative distribution function)of the uniform distribution function on $(a,b)$. The PDF and CDF of the normal distribution function is denoted in the standard notation of $\varphi $ and $\text{\Phi}$ respectively.
First we prove that $\underset{\sigma \to \mathrm{\infty}}{\mathrm{l}\mathrm{i}\mathrm{m}}{\varphi}^{\sigma}\left(x\right)=u\left(x\right)$ and use this to conclude that ${lim}_{\genfrac{}{}{0ex}{}{\sigma \to \mathrm{\infty}}{N\to \mathrm{\infty}}}\text{}\mathrm{E}\mathrm{M}\mathrm{D}\left({\text{A}}_{\sigma}^{\text{*}},\text{U}\right)={lim}_{\genfrac{}{}{0ex}{}{\sigma \to \mathrm{\infty}}{N\to \mathrm{\infty}}}\text{GI}=0$.
Therefore, ${lim}_{\sigma \to \mathrm{\infty}}\text{}{\varphi}_{\text{t}}^{\sigma}\left(x\right)=Constant$. Since the support of $\varphi}_{\mathrm{t}}^{\sigma$ is $(\text{a},\text{b})$, the only constant satisfying that ${lim}_{\sigma \to \mathrm{\infty}}\text{}{\varphi}_{t}^{\sigma}\left(\text{x}\right)$ is the probability distribution, $\frac{1}{\text{b}\text{a}}=u\left(\text{x}\right).$ Therefore, ${lim}_{\genfrac{}{}{0ex}{}{\sigma \to \mathrm{\infty}}{N\to \mathrm{\infty}}}\text{}\mathrm{E}\mathrm{M}\mathrm{D}\left({\text{A}}_{\sigma}^{\text{*}},\text{U}\right)={lim}_{\genfrac{}{}{0ex}{}{\sigma \to \mathrm{\infty}}{N\to \mathrm{\infty}}}\text{GI}=0$.
(II)
Let $\text{X}$, $\text{Y}$ be truncated normal distributions. In part (I) we prove that ${lim}_{\sigma \to \mathrm{\infty}}\mathrm{X}={lim}_{\sigma \to \mathrm{\infty}}\mathrm{Y}=\mathrm{u}(\mathrm{x})$. Theorem 2 implies that when $\text{X}$ and $\text{Y}$ are uniform distributions ${lim}_{\genfrac{}{}{0ex}{}{\sigma \to \mathrm{\infty}}{N\to \mathrm{\infty}}}\mathrm{G}\mathrm{I}=0$.
Appendix 2
Simulations with constant ζ
To assess the performance of DeBias we tested its ability to retrieve a predetermined local interaction parameter $\zeta$ (see Figure 1B) from simulated synthetic data. $\mathrm{X}$ and $\mathrm{Y}$ were modeled as truncated normal distributions on (−90, 90), with $\mu $=0 and changing $\sigma}_{\mathrm{x}},{\sigma}_{\mathrm{y}$. Pairs of ($\mathrm{X}}_{\mathrm{i}},{\mathrm{Y}}_{\mathrm{i}$) were sampled from $\mathrm{X},\mathrm{Y}$ and shifted towards each other by $\zeta$ degrees (similar to Figure 1B, but with a constant cumulative $\zeta$) to construct the observed alignment angles. To avoid confusion we denote $\mathrm{X},\mathrm{Y},{\sigma}_{\mathrm{x}},{\sigma}_{\mathrm{y}}$ as the observed values postsimulation. For a given constant $\zeta$, we exhaustively explored the $\sigma}_{\mathrm{x}},{\sigma}_{\mathrm{y}$ space. For each $\sigma}_{\mathrm{x}},{\sigma}_{\mathrm{y}$, we performed 20 independent simulations with N = 1600 observations ($\mathrm{X}}_{\mathrm{i}},{\mathrm{Y}}_{\mathrm{i}$). For each simulation we constructed the resampled distribution 10 times based on 400 observations drawn from the marginal $\mathrm{X},\mathrm{Y}$ distributions, and used the mean GI, LI. The final recorded GI, LI were averaged over the independent simulations.
The expected mean alignment when neither global bias nor local interactions exist is 45°. We begin by examining the deviation of the mean observed alignment from this value (45°  θ_{mean}). Better alignment is reflected by higher 45°  θ_{mean} values implying a larger deviation from the unbiased and nointeractions scenario. Low standard deviations $\sigma$, correspond to better alignment, improving with growing $\zeta$, as expected (Appendix 2—figure 1A). The GI follows a similar pattern and remains relatively stable for small changes in $\zeta$ (Appendix 2—figure 1B). The similar patterns between Appendix 2—figure 1A and B indicate that the global bias has a prominent role in determining the observed alignment.
The LI grows with $\zeta$ (Appendix 2—figure 1C) and its relative contribution to the observed alignment grow with increasing $\zeta$ (Appendix 2—figure 1D, quantified by LI/(LI+GI)), as expected. This relative contribution can be harnessed to restore an estimated $\zeta$ as the corresponding fraction from (45°  θ_{mean}) (Appendix 2—figure 1E). The estimated $\zeta$ is a lower bound for the actual value (Appendix 1, Theorem 4). Estimation is more accurate for larger $\zeta$ and for large $\sigma$ (Appendix 2—figure 1F). These results again highlight the importance of exploiting the GI for better interpretation of the LI (first introduced in Figure 2D–E).
We also investigated the effect of the choice of the number of bins K, used for sampling and computation of the EMD between distributions. Increased K induces linear growth in LI and GI values, as expected (Appendix 1 Theorem 5, Appendix 2—figure 2A–B) and stabilized its accuracy in predicting $\zeta$ for K ≥ 11 (Appendix 2—figure 2C–E). Large K will require more observations to estimate the true distribution. Using a constant K for a specific application assures fair comparison between different cases. Varying N, the number of observations, did not have a major effect on these measurements (Appendix 2—figure 3A–D), but increasing N reduced the noise which increased the accuracy in predicting ζ (Appendix 2—figure 3E).
Appendix 3
Simulating proteinprotein colocalization
Colocalization of molecules is commonly used to predict potential local interactions under the assumption that the stoichiometry of interacting molecules remains constant across all observations. We demonstrated this by simulating a random variable $\mathrm{X}$, representing molecular counts, and a random variable $\mathrm{Z}$, representing the local interaction. Pairs of samples $\mathrm{x}}_{\mathrm{i}$ from $\mathrm{X}$ and $\zeta}_{i$ from $\mathrm{Z}$were drawn, and $y}_{i}={x}_{i}{\zeta}_{i$ represented the molecular counts of the colocalized molecule (Appendix 3—figure 1A). The joint distributions of $\mathrm{X},\mathrm{Y}$ for four simulations are shown in Appendix 3—figure 1B. X is normally distributed with mean $\mu$ = 0.5 and standard deviations $\sigma$ = 0.2, truncated to [0,1]. Z is normally distributed with mean $\mu$ = 1 and different standard deviations, simulating gradually increasing local interactions (Appendix 3—figure 1B, lefttoright). Applying DeBias on these data demonstrated that increased local interactions (reduced $\sigma}_{\zeta$) are translated to increased LI (Appendix 3—figure 1C).
We next simulated a scenario, in which a partial subset of observations $\mathrm{y}}_{i$ interacted with $\mathrm{x}}_{i$; the remaining $\mathrm{y}}_{i$ were drawn independently from the distribution $\mathrm{Y}=\mathrm{X}$. We found that LI values increased with the fraction of interacting observations (Appendix 3—figure 1D).
Finally, we assessed how variation in the quantization parameter (K) and number of observations (N) alter GI and LI (Appendix 3—figure 1EF) and this was also demonstrated in our experimental data (Figure 6—figure supplement 1G–J).
Appendix 4
Pixel versus objectbased colocalization
Many methods have been developed to quantify proteinprotein colocalization. Pixelbased methods measure pixelwise correlation coefficients (Adler and Parmryd, 2010; Bolte and Cordelières, 2006; Costes et al., 2004; Manders et al., 1993; Pearson, 1901), exploiting the notion that fluorescence levels of colocalized proteins are correlated. They suffer from background signal where no colocalization exists. Objectbased methods first detect objects of interest and then assess colocalization based on secondorder statistics of the spatial distributions of the detections (Helmuth et al., 2010; Kalaidzidis et al., 2015; Lagache et al., 2015; Rizk et al., 2014; Semrau et al., 2011). Objectbased methods remove the background pixels but lose the information contained by the fluorescence levels at the detected objects. Thus, objectbased methods are best applicable for the colocalization of binary signals (Lagache et al., 2015), but less suited for applications in which colocalization accounts for coupling of molecular counts on a continuous spectrum. Moreover, objectbased methods require detection of objects in both channels, which often limits their applicability. Pros and cons for using either of these methods are presented in Appendix 4—Table 1. In the examples of receptorCCP colocalization we implemented a hybrid of the two approaches: colocalization analysis by DeBias focused on the intensity of fluorescent readouts within detected CCPs to decompose the coupling of the two intensity variables into LI and GI (Figure 6). The same decomposition was demonstrated for a diffuse signal without objects in the PKCFRET example (Figure 5).
References

Control of cellcell forces and collective cell dynamics by the intercellular adhesomeNature Cell Biology 17:409–420.https://doi.org/10.1038/ncb3135

A guided tour into subcellular colocalization analysis in light microscopyJournal of Microscopy 224:213–232.https://doi.org/10.1111/j.13652818.2006.01706.x

Endocytosis, signaling, and beyondCold Spring Harbor Perspectives in Biology 6:a016865.https://doi.org/10.1101/cshperspect.a016865

Metrics for assessing cytoskeletal orientational correlations and consistencyPLoS Computational Biology 11:e1004190.https://doi.org/10.1371/journal.pcbi.1004190

A practical guide to evaluating colocalization in biological microscopyAJP: Cell Physiology 300:C723–C742.https://doi.org/10.1152/ajpcell.00462.2010

The use of multiple measurements in taxonomic problemsAnnals of Eugenics 7:179–188.https://doi.org/10.1111/j.14691809.1936.tb02137.x

On the histogram as a density estimator: l 2 theoryProbability Theory and Related Fields 57:453–476.

Biosensors for characterizing the dynamics of rho family GTPases in living cellsCurrent Protocols in Cell Biology Chapter 14:14.11. 11–14.11.14.https://doi.org/10.1002/0471143030.cb1411s46

Measurement of orientation and distribution of cellular alignment and cytoskeletal organizationAnnals of Biomedical Engineering 27:712–720.https://doi.org/10.1114/1.226

Statistical analysis of molecule colocalization in bioimagingCytometry Part A 87:568–579.https://doi.org/10.1002/cyto.a.22629

Local clustering of transferrin receptors promotes clathrincoated pit initiationThe Journal of Cell Biology 191:1381–1393.https://doi.org/10.1083/jcb.201008117

Measurement of colocalization of objects in dualcolour confocal imagesJournal of Microscopy 169:375–382.https://doi.org/10.1111/j.13652818.1993.tb03313.x

A unified approach to the change of resolution: space and graylevelIEEE Transactions on Pattern Analysis and Machine Intelligence 11:739–742.https://doi.org/10.1109/34.192468

The earth mover's distance as a metric for image retrievalInternational Journal of Computer Vision 40:99–121.https://doi.org/10.1023/A:1026543900054

Stagespecific assays for coated pit formation and coated vesicle budding in vitroThe Journal of Cell Biology 114:869–880.https://doi.org/10.1083/jcb.114.5.869

Mechanical waves during tissue expansionNature Physics 8:628–634.https://doi.org/10.1038/nphys2355

Collective cell guidance by cooperative intercellular forcesNature Materials 10:469–475.https://doi.org/10.1038/nmat3025

Plithotaxis and emergent dynamics in collective cellular migrationTrends in Cell Biology 21:638–646.https://doi.org/10.1016/j.tcb.2011.06.006

On the definition of a confounderThe Annals of Statistics 41:196–220.https://doi.org/10.1214/12AOS1058

A genetically encoded fluorescent reporter reveals oscillatory phosphorylation by protein kinase CThe Journal of Cell Biology 161:899–909.https://doi.org/10.1083/jcb.200302125

Seeds of locally aligned motion and stress coordinate a collective cell migrationBiophysical Journal 109:2492–2500.https://doi.org/10.1016/j.bpj.2015.11.001
Decision letter

Fabrice CordelièresReviewing Editor; CNRS, France
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Decoupling global biases and local interactions between cell biological variables" for consideration by eLife. Your article has been favorably evaluated by Naama Barkai (Senior Editor) and three reviewers, one of whom, Fabrice Cordelières (Reviewer #1) served as Reviewing Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Chong Zhang (Reviewer #2); Perrine PaulGilloteaux (Reviewer #3).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Summary:
Zaritsky et al. propose a new approach to characterise interactions between two variables as a result of two influences: a global effect, independent of the two variables and influencing both, and a local effect that depends on their "direct" dependency. This paper raises the point of differentiating situations where correlations are observed, not due to a direct link between variables, but rather being the result of an overall behaviour. The authors built a framework in which both effects are reflected by two indicators: global and local indexes (GI and LI). The paper is divided into a theoretical part, where simulations are used to setup the framework, and a series of 4 biological examples where the "DeBias" method is applied.
Essential revisions:
1) The authors have performed simulation in which the bin widths of distributions or the number of observations have been progressively increased. Conclusions seem to be that a 15 quantization bins and a number of samples around 100 is acceptable. However, how should these parameters be set depending on the nature of the data to analyse? How should these parameters be set, depending on the distributions to analyse? Is there a generic rule? If so, this should be explained. If not, why did the author make those choices? The Methods section includes information about the number of quantization bins: why is it 89 for angular data, where as the value is 10 for PKC experiments and 39 for colocalization? From theoretical results, it seems LI or GI tends to (K1)/2 when the other indicator is zero, K being the number of quantization bins. How would the authors recommend proceeding when the K parameter is not the same from one experiment to another? I fear we are missing proper guidelines.
2) For example of lack of usability: I did test it also with my own data of spot colocalization (simulation) and interestingly the GIs were inversely related to the spot density (but which was also related to the image size or N), and LIs was less obvious (my N was 100^2 and 256^2 with 100 spots each time with no colocalisation or with 50% of them colocalizing). I’ve found
GI no coloc 50% coloc
(squareroot size) 100 4.47 4.32
(squareroot size)256 3.67 3.62
LI no coloc 50% coloc
(squareroot size)100 1.59 2.24
(squareroot size)256 1.69 1.57
Before any conclusion, interpretation would require more characterisation work, or maybe provided with a simulation script in order to better apprehend the significance of it. I would have had to repeat the experiment, too laborious from the server, and in addition I do not know which number of bins was used in that case and if it was constant.
3) In terms of local interaction ζ, it is simulated as a ratio of θ by α, is this a generic assumption that local interaction is a variable linearly related to the alignment?
4) Subsection “Simulating synthetic data”, the authors demonstrated in the simulation the effectiveness of GI and LI. While it seems true by simulating local interaction of alignment by shifting one of the angles towards the other, how would GI and LI behave when the local interaction is not alignment, i.e., by shifting one away from the other? Would GI and LI show the opposite effect so as to be able to discriminate opposite local interactions?
5) Subsection “DeBias procedure”, for the adjustment #1 for the colocalization quantification, such normalization implicitly assumes the channel signals/intensities have linear relationship. This could be a strong assumption. What would be the adjustment for those signals that do not have linear relationship?
6) Probably some additional validation to highlight the method's advantage over other conventional intensity based colocalization methods would be useful. Such additional evaluation is probably not necessarily within the scope of a "Tools" submission, but that if such data were available, it would be a valuable scientific addition to the paper. Or at least some kind of generic guidelines about under which situations when DeBias may or may not be applied would be instructive.
7) The quantification power is actually relative: in the same conditions (meaning here same K, and potentially same N see point 2), the GIs and LIs would actually give relative information, but a twofold interaction factor would not show a twofold LI. In addition each experiment shows very different range of LIs and GIs. These values are theoretically in the range of [0; (k1)/2], k been the number of bins, and does not denotes directly to the strength of the interaction, and have to be compared. Normalizing these descriptors by their upper bound may help to interpret.
8) In each of these experiments, the conclusions were drawn with a different path of reasoning and different constructions of statistical tests to compare and assess the significance of observed differences between these 2 descriptors for different conditions. I do believe that the manuscript would gain in impact by providing a general workflow (maybe a scheme with different possibilities), including the assessment of significance of the descriptors differences. Otherwise, the potential user of this powerful technique may struggle with the analysis and be led to erroneous conclusions.
9) The Discussion of this article is quite disappointing. The first part is a summary of previous conclusions, the second part being a short comparison to already published, approaching methods. The manuscript should benefit from this comparison being included into the proper section of the manuscript, namely the colocalisation one.
https://doi.org/10.7554/eLife.22323.019Author response
Essential revisions:
1) The authors have performed simulation in which the bin widths of distributions or the number of observations have been progressively increased. Conclusions seem to be that a 15 quantization bins and a number of samples around 100 is acceptable. However, how should these parameters be set depending on the nature of the data to analyse? How should these parameters be set, depending on the distributions to analyse? Is there a generic rule? If so, this should be explained. If not, why did the author make those choices? The Methods section includes information about the number of quantization bins: why is it 89 for angular data, where as the value is 10 for PKC experiments and 39 for colocalization? From theoretical results, it seems LI or GI tends to (K1)/2 when the other indicator is zero, K being the number of quantization bins. How would the authors recommend proceeding when the K parameter is not the same from one experiment to another? I fear we are missing proper guidelines.
Setting K is a tradeoff between the number of observations and the range of observed values. In the manuscript, the narrow dynamic range in the PKC imaging and wider range in the CME colocalization experiments implied different selection of K.
In response to these valid questions by the reviewers, we have further validated that changing K does not alter the interpretation of simulated or experimental data analysis for both coorientation and colocalization scenarios (Figure 3—figure supplement 1DE, Figure 6—figure supplement 1IJ, Appendix 3—figure 1F).
The FreedmanDiaconis rule is designed to calculate the histogram bin width so to minimize the difference between the area under the observed and the theoretical distribution (D Freedman, P Diaconis et al., 1981). Following the reviewers’ comments we implemented and integrated this rule to our source code and webserver implementation as a means to determine a default value of K and the marginal distributions’ bin width.
K must remain constant for a given application, for testing the alterations of LI/GI as response to different perturbations / experimental conditions. However, different applications may use different K values. This is emphasized in the Methods section (“Importantly, GI and LI across experiments can be compared only when evaluated with the same K value and this is enforced by the webserver.”).
2) For example of lack of usability: I did test it also with my own data of spot colocalization (simulation) and interestingly the Gis was inversely related to the spot density (but which was also related to the image size or N), and Lis was less obvious ((my N was 100^2 and 256^2 with 100 spots each time with no colocalisation or with 50% of them colocalizing. I’ve found
GI no coloc 50%coloc
(squareroot size) 100 4.47 4.32
(squareroot size)256 3.67 3.62
LI no coloc 50%coloc
(squareroot size)100 1.59 2.24
(squareroot size)256 1.69 1.57before any conclusion, interpretation would require more characterisation work, or maybe provided with a simulation script in order to better apprehend the significance of it. I would have had to repeat the experiment, too laborious from the server, and in addition I do not know which number of bins was used in that case and if it was constant.
It is not clear to us how these simulations were performed. Below we list possible explanation to the reviewer’s observations:
GI is inversely associated with the density of colocalized pixels. If the background pixels are drawn from a uniform distribution then more “background” pixels (in this case 6.5 fold) imply reduced GI, as observed.
LI is not increased when N = 256^2 (assuming 100 colocalized pixels). A fraction of ~0.0015 (100/255^2) colocalized pixels could be too subtle for DeBias to identify. In the case N = 100^2 and corresponding fraction of 1% colocalized pixels, this is easily identified (LI: 1.59 → 2.24).
In our new colocalization simulations we have mixed interacting and noninteracting observations and demonstrated that LI associates with an increased fraction of locally interacting observations (Appendix 3—figure 1D). In these simulations, the non interacting observations from the 2^{nd} channel (Y) were drawn from the distribution of the first channel (X). Following the reviewer concern, we have also simulated very low fractions of interacting observations of 0%, 0.1%, 0.5% and 1% to highlight the sensitivity of DeBias. See Author response image 1 (parameters were: σ_{ζ} = 0.1, σ_{x} = 0.2). We did not include these results in the manuscript because Appendix 3—figure 1D already captures a wide range of fractions of interacting observations. Its purpose was to demonstrate DeBias’ ability to identify scenarios of partial indications and this was achieved with the new figure panel.
Following the reviewers’ suggestion, we have made our colocalization and coorientation simulation publicly available.
3) In terms of local interaction ζ, it is simulated as a ratio of θ by α, is this a generic assumption that local interaction is a variable linearly related to the alignment?
The calculation of LI and GI is based on the assumption that the local interaction is constant, implemented by modeling the LI as an additive constant by complementing GI (calculated from the marginal distributions) to the observed coalignment / colocalization distribution. This assumption might be unrealistic in many biological systems. As DeBias is the first direct measure to quantify both local interactions and global bias, we simplified the estimation of LI and used it as a proxy to complex and more biologicalrealistic interactions. Our coalignment simulations modeled the local interactions as a fraction of the alignment θ. We were able to show, even with this simplifying assumption, that we can discriminate and estimate ζ accurately (Appendix 2—figure 1EF). We have now also included a simulation for colocalization, where the local interactions are modeled as associations between the channels (y = αx). These results are shown in Appendix 3. Adjusting DeBias to nonconstant bias is outofscope of this manuscript and will be pursued as future work.
4) Subsection “Simulating synthetic data”, the authors demonstrated in the simulation the effectiveness of GI and LI. While it seems true by simulating local interaction of alignment by shifting one of the angles towards the other, how would GI and LI behave when the local interaction is not alignment, i.e. by shifting one away from the other? Would GI and LI show the opposite effect so as to be able to discriminate opposite local interactions?
We have implemented simulations of two locally repulsive vector variables and demonstrated that yields negative LI values, indeed (Figure 2—figure supplement 1).
5) Subsection “DeBias procedure”, for the adjustment #1 for the colocalization quantification, such normalization implicitly assumes the channel signals/intensities have linear relationship. This could be a strong assumption. What would be the adjustment for those signals that do not have linear relationship?
This assumption was made to adjust the channels’ intensity to the [0,1] range. Indeed, this is a simplifying assumption that is now discussed in the new section on DeBias’ limitations in the Discussion. Adjusting the normalization to possibly nonlinear interactions is outofscope of this manuscript.
6) Probably some additional validation to highlight the method's advantage over other conventional intensity based colocalization methods would be useful. Such additional evaluation is probably not necessarily within the scope of a "Tools" submission, but that if such data were available, it would be a valuable scientific addition to the paper. Or at least some kind of generic guidelines about under which situations when DeBias may or may not be applied would be instructive.
We already demonstrated that including GI increases discriminative abilities for the experimental CME colocalization data (Figure 6EF) and for simulated coalignment data (Figure 2F – only subjectively). We now complement this by demonstrating enhanced discrimination performance of Pearson correlation by including the GI as a coupled measure (Figure 6—figure supplement 1CF). We also demonstrated the high association of LI and Pearson’s correlations (Pearson’s coefficient > 0.95), validating LI as a measure for local interactions.
7) The quantification power is actually relative: in the same conditions (meaning here same K, and potentially same N see point 2), the GIs and Lis would actually give relative information, but a twofold interaction factor would not show a twofold Li. In addition each experiment shows very different range of Lis and Gis. These values are theoretically in the range of [0; (k1)/2], k been the number of bins, and does not denotes directly to the strength of the interaction, and have to be compared. Normalizing these descriptors by their upper bound may help to interpret.
It is true that the quantification power is relative. K can be set independently for every application to discriminate between different experimental conditions. It is also true that a twofold increase in the direct interaction would not necessarily imply a twofold increase in LI – this requires a priori knowledge of the nature of the interaction and adjustment of the calculation of LI. Tackling this issue for all of the four application areas seems outside the scope of this manuscript. When absolute values are needed, we will in the future refer to the DeBias framework and calibrate the LI change on an applicationby application basis. It should be noted that alternative measures for pixelbased co localization are relative in nature as well. Normalizing by the upper bound will not solve the second concern (twofold increased interaction mapped to a twofold increase in LI). Since dividing by the constant K(K1) will induce very low numbers for LI and GI (e.g., for K = 15 and LI = 1 → < 0.005) we have decided to avoid it in our implementation.
We excluded the possibility that N has a major effect on the calculated GI and LI values (point #2).
8) In each of these experiments, the conclusions were drawn with a different path of reasoning and different constructions of statistical tests to compare and assess the significance of observed differences between these 2 descriptors for different conditions. I do believe that the manuscript would gain in impact by providing a general workflow (maybe a scheme with different possibilities), including the assessment of significance of the descriptors differences. Otherwise, the potential user of this powerful technique may struggle with the analysis and be led to erroneous conclusions.
We included a general assessment of discrimination and integrated it into our source code and webserver implementation: the pvalue of a twosided Wilcoxon rank sum test for LI and GI of two distinct experimental conditions. In the webserver this part is now referred as ‘DeBias Analyst’.
In the manuscript, the only data assessed differently were the velocitystress alignment measures (Figure 4) – because we had a single observation per condition and thus could not perform statistical analysis between independent experimental instances. This data was not produced in our lab (previously published data) so we could not retrieve more replicates and thus had to come up with a different analysis. Other methods were also proposed in the manuscript (e.g., classification accuracy, Figure 6EF) but are not provided in our code distribution.
9) The Discussion of this article is quite disappointing. The first part is a summary of previous conclusions, the second part being a short comparison to already published, approaching methods. The manuscript should benefit from this comparison being included into the proper section of the manuscript, namely the colocalisation one.
We feel that the summary part starting the Discussion is necessary because of the nature of the work. Showcasing DeBias on four quite distinct examples emphasizes the generalization of the method, especially, in extracting insightful mechanistic information that was not available in earlier colocalization quantification methods. This point is now emphasized in the ‘summary’ of the Discussion. Another reason for a summary is that many people read only the Discussion before they decide whether to invest time in reading the full text.
We did not follow the suggestion to move the comparison between pixel and object based colocalization to the colocalization part of the manuscript. We feel that it will be confusing to follow a colocalization discussion in the Results section immediately before going into a general discussion of the features of DeBias analysis in the Discussion section. Having the comparison at the beginning of the colocalization results would interfere with the flow of the manuscript. Instead, we have rewritten the Discussion part addressing this issue and inserted a new Appendix 4 that discusses colocalization approaches in greater depth.
We also added a section on limitations of DeBias to the Discussion. It includes the issues raised by the reviewers: linearity assumption for colocalization normalization, adjusting DeBias to nonconstant bias, relative quantification power, LI and GI are slightly coupled and that in some cases it takes additional means to interpret what biological properties are encoded in the GI.
https://doi.org/10.7554/eLife.22323.020Article and author information
Author details
Funding
EMBO (postdoctoral fellowship)
 Uri Obolski
National Institute of General Medical Sciences (R01 GM073165)
 Gaudenz Danuser
 Sandra L Schmid
Cancer Prevention and Research Institute of Texas (R1225)
 Gaudenz Danuser
National Institute of General Medical Sciences (P01 GM103723)
 Gaudenz Danuser
National Institute of General Medical Sciences (P01 GM096971)
 Gaudenz Danuser
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We thank the BioHPC team at UTSW and especially Liqiang Wang for the help in implementing the DeBias webserver and making it freely available for public use. We are grateful to Tamal Das and Joachim Spatz for providing us with the motion and stress data from their tightjunction screen. We thank Andrew Jamieson for helping to package the source code and creating the repository. We thank Dana Reed for assistance with cell characterization and mycoplasma testing. We thank Sangyoon Han, Claudia Schaefer, Meghan Driscoll, Erik Welf, Phillippe Roudot and Marcel Mettlen for critically reading the manuscript and Marcel Mettlen for fruitful discussions and advice. This work was supported by the Cancer Prevention and Research Institute of Texas (CPRIT R1225 to GD), by NIH P01 GM103723, P01 GM096971 (to GD) and by GM713165 (to SLS and GD). UO was supported by an EMBO LongTerm fellowship. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.
Reviewing Editor
 Fabrice Cordelières, CNRS, France
Publication history
 Received: October 13, 2016
 Accepted: March 10, 2017
 Accepted Manuscript published: March 13, 2017 (version 1)
 Version of Record published: May 2, 2017 (version 2)
Copyright
© 2017, Zaritsky et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,956
 Page views

 417
 Downloads

 10
 Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.