## ----v1, include = FALSE------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
eval = TRUE
)
## ----installation, include = TRUE, eval = FALSE-------------------------------
# if (!requireNamespace("BiocManager")) {
# install.packages("BiocManager")
# }
# BiocManager::install("spatialFDA")
## ----setup, warning = FALSE, message = FALSE----------------------------------
library("spatialFDA")
library("dplyr")
library("ggplot2")
library("tidyr")
library("stringr")
library("dplyr")
library("patchwork")
library("SpatialExperiment")
set.seed(1234)
## ----loading, warning = FALSE, message = FALSE--------------------------------
# retrieve example data from Damond et al. (2019)
spe <- .loadExample(full = TRUE)
spe <- subset(spe, ,patient_id %in% c(6089,6180,6126,6134,6228,6414))
# set cell types as factors
colData(spe)$cell_type <- as.factor(colData(spe)$cell_type)
## ----plotting fovs, warning = FALSE, fig.width=8, fig.height=15---------------
df <- data.frame(spatialCoords(spe), colData(spe))
dfSub <- df %>%
subset(image_name %in% c("E02", "E03", "E04", "E05"))
p <- ggplot(dfSub, aes(x = cell_x, y = cell_y, color = cell_category)) +
geom_point(size= 0.5) +
facet_wrap(~image_name) +
theme(legend.title.size = 20, legend.text.size = 20) +
xlab("x") +
ylab("y") +
labs(color = "cell category")+
coord_equal() +
theme_light()
dfSub <- dfSub %>%
subset(cell_type %in% c("alpha", "beta", "delta", "Th", "Tc"))
q <- ggplot(dfSub, aes(x = cell_x, y = cell_y, color = cell_type)) +
geom_point(size= 0.5) +
facet_wrap(~image_name) +
theme(legend.title.size = 20, legend.text.size = 20) +
xlab("x") +
ylab("y") +
labs(color = "cell type") +
coord_equal() +
theme_light()
wrap_plots(list(p,q), widths = c(1,1), heights = c(1,1), nrow = 2, ncol = 1)
## ----Lfunction, warning = FALSE, message = FALSE------------------------------
metricRes <- calcMetricPerFov(spe = spe, selection = c("alpha", "Tc"),
subsetby = "image_number", fun = "Lcross",
marks = "cell_type",
rSeq = seq(0, 50, length.out = 50),
by = c("patient_stage", "patient_id",
"image_number"),
ncores = 1)
metricRes %>% head(3)
## ----plotLfunction, warning = FALSE, fig.width=8, fig.height=8----------------
# create a unique plotting ID
metricRes$ID <- paste0(
metricRes$patient_stage, "|", metricRes$patient_id
)
# change levels for plotting
metricRes$ID <- factor(metricRes$ID, levels = c("Non-diabetic|6126",
"Non-diabetic|6134",
"Onset|6228","Onset|6414",
"Long-duration|6089",
"Long-duration|6180"))
# plot metrics
plotMetricPerFov(metricRes, correction = "iso", x = "r",
imageId = "image_number", ID = "ID", ncol = 2)
## ----Gfunction, warning = FALSE, message = FALSE------------------------------
metricRes <- calcMetricPerFov(spe = spe, selection = c("alpha", "Tc"),
subsetby = "image_number", fun = "Gcross",
marks = "cell_type",
rSeq = seq(0, 50, length.out = 50),
by = c("patient_stage", "patient_id",
"image_number"),
ncores = 1)
metricRes %>% head(3)
## ----plotGfunction, warning = FALSE, fig.width=8, fig.height=8----------------
# create a unique plotting ID
metricRes$ID <- paste0(
metricRes$patient_stage, "|", metricRes$patient_id
)
# change levels for plotting
metricRes$ID <- factor(metricRes$ID, levels = c("Non-diabetic|6126",
"Non-diabetic|6134",
"Onset|6228","Onset|6414",
"Long-duration|6089",
"Long-duration|6180"))
# plot metrics
plotMetricPerFov(metricRes, correction = "rs", x = "r",
imageId = "image_number", ID = "ID", ncol = 2)
## ----funcBoxPlot, warning = FALSE, results='hide'-----------------------------
# create a unique ID per row in the dataframe
metricRes$ID <- paste0(
metricRes$patient_stage, "x", metricRes$patient_id,
"x", metricRes$image_number
)
#removing field of views that have as a curve only zeros - these are cases where
#there is no cells of one type
metricRes <- metricRes %>% dplyr::group_by(ID) %>% dplyr::filter(sum(rs) >= 1)
collector <- plotFbPlot(metricRes, "r", "rs", "patient_stage")
## ----variancetransform, warning = FALSE---------------------------------------
# can determine with a boxcox transformation what is the ideal parameter
# for the transformation
metricRes$rs <- sqrt(metricRes$rs)
collector <- plotFbPlot(metricRes, 'r', 'rs', 'patient_stage')
## ----fPCA, warning = FALSE----------------------------------------------------
# filter out all rows that have a constant zero part - all r<10
metricRes <- metricRes %>% filter(r > 10)
# prepare dataframe from calcMetricRes to be in the correct format for pffr
dat <- prepData(metricRes, "r", "rs")
# create meta info of the IDs
splitData <- dat$ID %>%
str_replace("-","_") %>%
str_split_fixed("x", 3) %>%
data.frame(stringsAsFactors = TRUE) %>%
setNames(c("condition", "patient_id", "imageId")) %>%
mutate(condition = relevel(condition,"Non_diabetic"))
dat <- cbind(dat, splitData)
# drop rows with NA
dat <- dat |> drop_na()
# calculate the fPCA
pca <- functionalPCA(dat = dat, r = metricRes$r |> unique(), pve = 0.995)
evalues <- pca$evalues
efunctions <- pca$efunctions
# plot the mean curve and the two first eigenfunctions
p_mu <- ggplot(data.frame(r = unique(metricRes$r), mu = pca$mu),
aes(x = r, y = mu)) +
geom_line() +
theme_light() +
xlab("r [µm]")
p_efunction1 <- ggplot(data.frame(r = unique(metricRes$r),
phi1 = pca$efunctions[,1]),
aes(x = r, y = phi1)) +
geom_line() +
theme_light() +
ylim(-0.3,0.3) +
xlab("r [µm]")
p_efunction2 <- ggplot(data.frame(r = unique(metricRes$r),
phi2 = pca$efunctions[,2]),
aes(x = r, y = phi2)) +
geom_line() +
theme_light() +
ylim(-0.3,0.3) +
xlab("r [µm]")
wrap_plots(list(p_mu, p_efunction1, p_efunction2), ncol = 3)
# plot the biplot of the first two PCs
plotFpca(dat = dat, res = pca, colourby = "condition")
# print the eigenvalues
evalues
## ----funcGamG, fig.height=10, warning = FALSE---------------------------------
library('refund')
# create a design matrix
mm <- model.matrix(~condition, data = dat)
colnames(mm)[1] <- "Intercept"
mm %>% head()
r <- metricRes$r |> unique()
# fit the model
mdl <- functionalGam(
data = dat, x = r,
designmat = mm, weights = dat$npoints,
formula = formula(Y ~ 1 + conditionLong_duration +
conditionOnset + s(patient_id, bs = "re")),
family = mgcv::scat(link = "log"),
algorithm = "gam"
)
summary(mdl)
plotLs <- lapply(colnames(mm), plotMdl, mdl = mdl,
shift = mdl$coefficients[["(Intercept)"]])
wrap_plots(plotLs, nrow = 3, axes = 'collect')
## ----contour, warning = FALSE-------------------------------------------------
resid(mdl) |> cor() |> filled.contour(levels = seq(-1, 1, l = 40))
resid(mdl) |> cov() |> filled.contour()
qqnorm(resid(mdl), pch = 16)
qqline(resid(mdl))
## ----intercept, warning = FALSE, eval = TRUE----------------------------------
# look at the smooth random intercepts per patient
data <- coef(mdl)$smterms$`s(patient_id)`$coef
data <- data %>% left_join(splitData %>%
select(patient_id, condition) %>% unique)
p <- ggplot(data, aes(x = x.vec, y = value, colour = condition, group = patient_id)) +
geom_line() +
theme_light() +
geom_smooth(aes(group = 1), col = 'black') +
xlab("r [µm]")
p
## ----sessionInfo--------------------------------------------------------------
sessionInfo()