2.2 Detecting global directionality: hotellingsTest
Textor et al (2011) proposed using Hotelling’s T-square test to detect directionality in a dataset. This test computes displacement vectors of all steps in the data, and tests if the mean displacement vector is significantly different from the null vector.
In using this test, it is important to realise that steps from within a track are not independent: cells usually exhibit some kind of persistent motion at least on a short time scale (that is: the same cell does not tend to go in completely opposite directions at to subsequent timepoints). This means that blindly extracting all steps violates the “independent observations” assumption of Hotelling’s test.
For example, applying this to the TCellsRaw dataset:
hotellingsTest( TCellsRaw, col = "gray" )
##
## Hotelling's one sample T2-test
##
## data:
## T2 = 1.3213, df1 = 2, df2 = 872, p-value = 0.5172
## alternative hypothesis: true location is not equal to c(0,0)
we find (as expected) that there is no no evidence for global directionality among T cells in a healthy lymph node. However, applying the same to the BCellsRaw dataset:
hotellingsTest( BCellsRaw, col = "gray" )
##
## Hotelling's one sample T2-test
##
## data:
## T2 = 18.978, df1 = 2, df2 = 1276, p-value = 8.175e-05
## alternative hypothesis: true location is not equal to c(0,0)
there seems to be a significant directionality. We know, however, that this bias shouldn’t be there–so what is going on?
This is an example of the problem of steps being correlated. If, by chance, a few more cells happen to go in one direction, on the step-level this shows up as a “significant” effect because the individual steps are treated as independent (inflating the power of the test). If we compare the B and T cells in their persistence time in an autocovariance plot (see the vignette on analysis methods for details), we get:
# Compute autocovariance, normalize so the first point lies at 1
Tacov <- aggregate( TCellsRaw, overallDot )
Tacov$value <- Tacov$value / Tacov$value[1]
# The same for B cells:
Bacov <- aggregate( BCellsRaw, overallDot )
Bacov$value <- Bacov$value / Bacov$value[1]
# Compare autocovariances:
plot( Tacov )
points( Bacov, col = "red" )
We see here that the B cells are more persistent, and the correlation between step directions (on the y-axis) only disappears after 10 steps or more.
When we only consider steps that are some distance apart using the step.spacing argument:
hotellingsTest( BCellsRaw, col = "gray", step.spacing = 10 )
##
## Hotelling's one sample T2-test
##
## data:
## T2 = 0.962, df1 = 2, df2 = 135, p-value = 0.6214
## alternative hypothesis: true location is not equal to c(0,0)
there is no longer evidence for directionality. (Note, however, that the power of the test is also reduced because the larger step spacing reduces the total number of steps).
By contrast, directionality remains obvious in the dataset with drift, even at a large step spacing:
hotellingsTest( TCells.drift, plot = TRUE, col = "gray", step.spacing = 10 )
##
## Hotelling's one sample T2-test
##
## data:
## T2 = 36.856, df1 = 2, df2 = 99, p-value = 1.799e-07
## alternative hypothesis: true location is not equal to c(0,0)
Here we also visualize the effect by setting plot=TRUE; the blue ellipse indicates the global direction (which is indeed the direction we simulated drift in).
! Note : The appropriate value of step.spacing depends on the level of persistence expected in the cells of interest, as well as on the time resolution of the experiment. An autocorrelation or autocovariance plot may provide insights in the persistence time for the cells of interest (see the vignette on analysis methods for details). If tracks of static structures are used to detect and correct for drift – as is recommended to distinguish between drift and truly directed movement – a step.spacing of zero is safe to use because there should be no inherent persistence for something that does not move.
Of note, by default hotellingsTest() is performed on a projection of the tracks on only the \(x\) and \(y\) dimensions. For 3D tracking data, we can specify that all three dimensions should be taken into account (although plotting is not supported in 3D):
hotellingsTest( TCellsRaw, dim = c("x","y","z"), step.spacing = 10 )
##
## Hotelling's one sample T2-test
##
## data:
## T2 = 1.5623, df1 = 3, df2 = 98, p-value = 0.6761
## alternative hypothesis: true location is not equal to c(0,0,0)
hotellingsTest( TCells.drift, dim = c("x","y","z"), step.spacing = 10 )
##
## Hotelling's one sample T2-test
##
## data:
## T2 = 47.594, df1 = 3, df2 = 98, p-value = 2.396e-08
## alternative hypothesis: true location is not equal to c(0,0,0)
However, this is not always recommended; the z-dimension in many two-photon experiments is much smaller than the x and y dimensions, and the resolution in this dimension is smaller as well. The z-coordinates of these data are therefore not always reliable and it can be better to leave them out for some track analyses.
2.3 Detecting global directionality: angle analysis
Beltman et al (2009) proposed an analysis of angles versus distance between cell pairs to detect global directionality in a dataset. Use the function analyzeCellPairs to get a dataframe with for each pair of tracks in the dataset the angle (between their displacement vectors) and the distance (min distance between the tracks at any timepoint). Note that the distance is not defined for pairs of tracks that have no overlap in timepoints – these will get an NA value.
# compute for both original as drift data
df.drift <- analyzeCellPairs( TCells.drift )
df.norm <- analyzeCellPairs( TCellsRaw )
# Plot
p.norm <- ggplot( df.norm, aes( x = dist, y = angle ) ) +
geom_point( color = "gray70", size = 0.5 ) +
stat_smooth( span = 1, fill = "blue" )+
labs( x = "distance between cell pairs",
y = "angle between cell pairs",
title = "original data") +
geom_hline( yintercept = 90, color = "red" ) +
theme_classic()
p.drift <- ggplot( df.drift, aes( x = dist, y = angle ) ) +
geom_point( color = "gray70", size = 0.5 ) +
stat_smooth( span = 1, fill = "blue" )+
labs( x = "distance between cell pairs",
y = "angle between cell pairs",
title = "data with drift") +
geom_hline( yintercept = 90, color = "red" ) +
theme_classic()
gridExtra::grid.arrange( p.norm, p.drift, ncol = 2 )
In the original data, the mean angle is roughly 90 degrees on average for cell pairs at any distance from each other, which is what we would expect if there is no global directionality in the data. In the data with drift, the average angle is lower than we would expect even at large distances between cell pairs – which indicates that there is global directionality (even cells far apart move roughly in the same direction).
We can perform a similar analysis on the level of single steps to boost the power. The function analyzeStepPairs() finds all pairs of steps from different cells that occur at the same timepoint, and then computes distances (between the step starting points) and angles between them. Note that we use a function in the argument filter.steps to take into account only steps with a minimum displacement, to avoid noise from steps where the cell is pausing.
# compute for both original as drift data.
df.drift <- analyzeStepPairs( TCells.drift, filter.steps = function(x) displacement(x)>2 )
df.norm <- analyzeStepPairs( TCellsRaw, filter.steps = function(x) displacement(x)>2 )
# Plot
p.norm <- ggplot( df.norm, aes( x = dist, y = angle ) ) +
geom_point( color = "gray70", size = 0.5 ) +
stat_smooth( span = 0.75, fill="blue" )+
labs( x = "distance between step pairs",
y = "angle between step pairs",
title = "original data") +
geom_hline( yintercept = 90, color = "red" ) +
theme_classic()
p.drift <- ggplot( df.drift, aes( x = dist, y = angle ) ) +
geom_point( color = "gray70", size = 0.5 ) +
stat_smooth( span = 0.75, fill="blue")+
labs( x = "distance between step pairs",
y = "angle between step pairs",
title = "data with drift") +
geom_hline( yintercept = 90, color = "red" ) +
theme_classic()
gridExtra::grid.arrange( p.norm, p.drift, ncol = 2 )
Again we see a mean step pair angle that is below 90 degrees for the data with drift, which is not the case for the original data.
2.4 Correcting drift
If Hotelling’s test and/or angle analyses provide evidence for drift in the dataset, this can be corrected by computing the mean step displacement:
# Get steps and find their displacement vectors
steps.drift <- subtracks( TCells.drift, 1 )
step.disp <- t( sapply( steps.drift, displacementVector ) )
# Get the mean
mean.displacement <- colMeans( step.disp )
# Divide this by the mean timestep to get a drift speed
drift.speed <- mean.displacement/timeStep( TCells.drift )
drift.speed
## [1] 0.05411294 0.04793798 0.04538616
This is indeed roughly the drift speed we introduced in the data. Now remove this:
correct.drift <- function( x, drift.vector )
{
# separate timepoints and coordinates
tvec <- x[,1]
coords <- x[,-1]
# compute movement due to drift.
drift.matrix <- matrix( rep( drift.vector, nrow(coords) ),
ncol = ncol(coords), byrow= TRUE )
drift.matrix <- drift.matrix * tvec
# Add drift to coordinates
x[,-1] <- coords - drift.matrix
return(x)
}
# Create data with drift
TCells.corrected <- as.tracks( lapply( TCells.drift, correct.drift, drift.speed ) )
# Compare (zoom in on part of the field to see better)
par( mfrow = c(1,2) )
plot( TCellsRaw, col = "blue", main = "uncorrected", cex = 0, pch.start = NA, xlim = c(200,400), ylim = c(0,200) )
plot( TCells.drift, col = "red", add = TRUE, cex = 0, pch.start = NA )
plot( TCellsRaw, col = "blue", main = "corrected", cex = 0, pch.start = NA, xlim = c(200,400), ylim = c(0,200) )
plot( TCells.corrected, col = "red", add = TRUE, cex = 0, pch.start = NA )
While the correction is not perfect, at least the tracks resemble their original shape more.
