Spatial summary functions based on point processes
The {mxfda} package provides methods for analyzing
spatial relationships between cell types in single cell imaging data
based on point process theory. The location of cells in image samples
are treated as following a point process, realizations of a point
process are called “point patterns”, and point process models seek to
understand correlations in the spatial distributions of cells. Under the
assumption that the rate of a cell is constant over an entire region of
interest a point pattern will exhibit complete spatial randomness (CSR),
and it is often of interest to model whether cells deviate from CSR
either through clustering or repulsion. When the rate of a cell is not
constant, this CSR assumption is violated and to estimate CSR it is
recommended to use permutations (see C. Wilson et
al. (2022) for explanation).
Spatial summary statistics can be calculated to quantify the
clustering and co-occurrence of cells in a circular region with a
particular radius r. Typically univariate (one cell type) or
bivariate (two cell types) summary statistics are reported, and
inference is obtained by comparing the observed spatial summary
statistic to that obtained under CSR. A popular quantity is Ripley’s
K(r), which studies the number of neighbors to a particular
point within radius r, and has univariate and bivariate
implementations in the {spatstat} package (Baddeley, Rubak, and Turner (2015)). Ripley’s K
is characterized by clustering or repulsion depending on whether it is
above or below the theoretical value of \(\pi
r^2\). Other spatial summary statistics analyze the distance to a
neighbor, and can be interpreted as probabilities of observing a
particular cell type within a radius r. One of these metrics,
G(r), or the nearest neighbor distance distribution, is the
cumulative distribution function of an exponential random variable. More
detailed overviews of spatial summary functions for multiplex imaging
data are provided in C. M. Wilson et al.
(2021) and Wrobel, Harris, and Vandekar
(2023).
Univariate summary functions
Below we calculate univariate Ripley’s K to summarize the spatial
relationship among immune cells in each image. The {mxfda}
package accomplishes this with the function
extract_summary_functions(). Either univariate or bivariate
can be calculated with this function depending on the choice supplied to
the extract_func argument. To calculate a univariate
spatial summary we supply univariate to the
extract_func argument. The summary function that is
calculated depends on the function supplied to summary_fun
which is one of Kest(), Gest(), or
Lest() from the {spatstat.explore} package
(Kcross(), Gcross(), or Lcross()
for bivariate methods). Other options include supplying a vector of
radius values through r_vec, and the a specific edge
correction (see Baddeley, Rubak, and Turner
(2015)). We calculate the K function across a range of radii from
0 to 100 and use the isotropic (“iso”) edge correction. See Baddeley, Rubak, and Turner (2015) for more
details on edge corrections for Ripley’s K and nearest neighbor G.
Options permute_CSR and permutations can be
used if interested in using a measure of the sample-specific CSR instead
of the theoretical, the column specified with markvar just
needs to contain 1 level more than used for summary_func (2
levels for univariate and 3 levels for bivariate) otherwise falls back
to theoretical.
mxFDAobject = extract_summary_functions(mxFDAobject,
extract_func = univariate,
summary_func = Kest,
r_vec = seq(0, 100, by = 1),
edge_correction = "iso",
markvar = "immune",
mark1 = "immune")
#> Using Theoretical Complete Spatial Randomness for Ripley's K
#> ■■■■■■■■■■■■■■■ 46% | ETA: 3s
#> ■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 98% | ETA: 0sRunning this code will calculate univariate Ripley’s K function to
measure spatial clustering of immune cells for each sample, and will
store these spatial summary functions in the
univariate_summaries slot of the mxFDAobject.
To access this slot and view the extracted summary functions, type:
mxFDAobject@univariate_summaries$Kest
#> # A tibble: 24,947 × 6
#> patientImage_id r sumfun csr fundiff `immune cells`
#> <chr> <dbl> <dbl> <dbl> <dbl> <int>
#> 1 2_1 0 0 0 0 6
#> 2 2_1 1 0 3.14 -3.14 6
#> 3 2_1 2 0 12.6 -12.6 6
#> 4 2_1 3 0 28.3 -28.3 6
#> 5 2_1 4 0 50.3 -50.3 6
#> 6 2_1 5 0 78.5 -78.5 6
#> 7 2_1 6 0 113. -113. 6
#> 8 2_1 7 0 154. -154. 6
#> 9 2_1 8 0 201. -201. 6
#> 10 2_1 9 0 254. -254. 6
#> # ℹ 24,937 more rowsNote that the summaries are returned as a dataframe. The variable
sumfun is the estimated summary function value,
csr is the theoretical value under complete spatial
randomness, and fundiff =
sumfun-csr describes the “degree of clustering
beyond what is expected due to chance; in downstream analysis we will
use the fundiff covariate.
Plotting the mxFDA object
{mxfda} has S4 methods for visualization implemented via
the plot() function (see ?plot.mxFDA for
details). The first argument is the mxFDA object followed
by a few options that depend on what plot output is desired. Here, we
want to plot the univariate summary that we just calculated, which was
the K function. By passing in what = 'uni k', the plot
function will extract the univariate K results. We also need to tell
plot() what column is the y-axis which can be
'sumfun' for the observed value, 'csr' for the
theoretical value of complete spatial randomness (CSR), or
'fundiff' which is the difference between the observed K
measure and the theoretical CSR. The output of plot() is a
{ggplot2} object which can then be easily added
to/manipulated as any ggplot plot would.
NOTE: These are the columns when calculating using the
extract_summary_function() of {mxfda} but if
summary data is added from elsewhere with
add_summary_function() then those column names will have to
be used.
plot(mxFDAobject, y = "fundiff", what = "uni k") +
geom_hline(yintercept = 0, color = "red", linetype = 2)
Bivariate summary functions
The extract_summary_functions() function can also be
used to extract bivariate summaries comparing spatial clustering of 2
cell types. We will look at relationship between T-cells and
macrophages. There are a few images that have fewer than 5 T-cells or
macrophages, which makes estimation of spatial summary functions less
stable for those images. To look at T-cells and macrophages, the data
phenotypes and cell locations have to be in long format so we first
create a variable with the cell types ('phenotype') from
the lung_df.
lung_df = lung_df %>%
mutate(phenotype = case_when(phenotype_cd8 == "CD8+" ~ "T-cell",
phenotype_cd14 == "CD14+" ~ "macrophage",
TRUE ~ "other"),
phenotype = factor(phenotype))We then recreate the mxFDAobject
spatial = lung_df %>%
select(-image_id, -gender, -age, -survival_days, -survival_status, -stage)
mxFDAobject = make_mxfda(metadata = clinical,
spatial = spatial,
subject_key = "patient_id",
sample_key = "patientImage_id")Now we calculate the bivariate G function, but can replace
Gcross() with Lcross() or
Kcross() to estimate the L or K bivariate functions
instead. The argument markvar takes the variable that we
created above called 'phenotype', and the 2 cell types that
we are interested in calculating the bivariate G for are
'T-cell' and 'macrophage' so we provide them
to mark1 and mark2, respectively.
mxFDAobject = extract_summary_functions(mxFDAobject,
summary_func = Gcross,
extract_func = bivariate,
r_vec = seq(0, 100, by = 1),
edge_correction = "rs",
markvar = "phenotype",
mark1 = "T-cell",
mark2 = "macrophage")
#> Using Theoretical Complete Spatial Randomness for Bivariate Nearest Neighbor GPlotting bivariate G
Just like with the univariate plots, we can use the
plot() function to plot our mxFDA object
results. The what now is 'bi g',
'bivar g', or 'bivariate g'.
plot(mxFDAobject, y = "fundiff", what = "bi g") +
geom_hline(yintercept = 0, color = "red", linetype = 2)
#> Warning: Removed 712 rows containing missing values or values outside the scale range
#> (`geom_line()`).
