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.
