---
title: "LDAcoop"
output: rmarkdown::html_vignette
vignette: >
%\VignetteIndexEntry{LDAcoop}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
```
# What this package is about
- LDA: limiting dilution assay
- coop: (cellular) cooperation
The LDA is the gold standard for quantification of the clonogenic capacity of
non-adherent cells.
The most common method for calculating this capacity from an LDA experiment
assumes implicitly independence of all cells. Interestingly, this assumption
does not hold for many cell lines (i.e. cooperatively and competitively growing
cells).
### The effect of cellular cooperation
```{r define caption, echo = FALSE}
mycaption <- "Figure: Fitting LDA data;
assuming independence (blue), taking cooperativity into account (grey)"
```
```{r setup,echo=FALSE,fig.width=5, fig.height=3, fig.cap=mycaption}
#install.packages("LDAcoop")
library(LDAcoop)
data(LDAdata)
BT20 <- subset.data.frame(x = LDAdata,
subset = (name == "BT.20") & (Group == 8))
BT20 <- BT20[,c("S-value","# Tested","# Clonal growth","Group","replicate")]
BT20 <- BT20[BT20$`S-value`<1000,]
out <- LDA_prepare_plot(LDA_tab = BT20)
par(
mar = c(2.5, 3.5, 0.5, 0.5),
mgp = c(1.5, 0.5, 0)
)
out[[1]][[1]]$color <- "#545454"
xl <- 750
LDA_plot_activity(LDA_obj = out[[1]],uncertainty.band = T,xlim = c(0,xl))
#statmod::elda(tested = BT20$`# Tested`,
# dose = BT20$`S-value`,
# response = BT20$`# Clonal growth`)
# Lower Estimate Upper
#Group 1 836.5526 645.4182 497.9539
lines(x = c(0,xl),y = c(0,-xl/645.4182),lwd=2,col="blue")
polygon(x = c(0,xl,xl,1),y = c(0,-xl/836.5526,-xl/497.9539,0),border = F,
col = rgb(red = 0,green = 0,blue = 1,alpha = 0.5,maxColorValue = 1))
```
\vspace{0.5cm}
Cellular cooperation induces a non-linearity to the dose-response-curve
(the (log) fraction of non-responding wells over the number of cells seeded)
and thereby biases the estimate of the gold standard analysis method,
which assumes linearity.
In other words: when surrounded by many other cells, the combined
clonogenic activity of
100 cells is higher than the activity of 100 separated single cells.
In terms of single cell activity, 100 cells act as if they were more.
Thus, the gold standard analysis method is syntactically applicable
but its result is biased (i.e. meaningless) for cooperatively growing
cells.
### Robust analysis in presence (and absence) of cellular cooperation
Therefore, this package equips you with the tools you need to robustly quantify
the clonogenic capacity of cell lines as well as to quantify the clonogenic
survival of those cell lines after treatment.
# How to use this package
As known from tools for quantification of the clonogenic capacity from LDA data,
the input is a data frame (table, matrix) with entries for
- cells: the number of cells seeded in a well
- wells: the number of replicate wells
- response: the number of wells with positive response (growth)
- (optional) group: treatment group - e.g. irradiation dose
All commands to get tables with numbers or plots of the activity are explained
below.
## A full experiment
Let's assume, you have conducted an experiment with a set of experimental
conditions and biological replicates.
Your data will look something like this
```{r showData}
#install.packages("LDAcoop")
library(LDAcoop)
data(LDAdata)
head(LDAdata)
```
A table of clonogenic activities, cooperativity coefficients and survival
fractions can be generated from such a data table as follows:
```{r data.table}
BT20 <- subset.data.frame(x = LDAdata,
subset = name == "BT.20")
BT20 <- BT20[,c("S-value","# Tested","# Clonal growth","Group","replicate")]
round(LDA_table(x = BT20),digits = 3)
```
In this table, we find
- `treatment`: the experimental 'Group'
- `act`: the calculated activity for the treatment
- `act.CI.lb` and `act.CI.ub`: its uncertainty (lower and upper bound of the
95%-confidence interval)
- `b`: the cooperativity coefficient (1: no cooperativity)
- `b.pvalue`: the p-value for test hypothesis of b=1
- `SF`: the calculated survival fraction (SF)
- `SF.CI.lb` and `SF.CI.ub`: SF-uncertainty (lb: lower bound, ub: upper bound)
A plot of the experiment and the estimated survival fractions can be generated
as follows:
```{r data.plot,fig.width=6, fig.height=4}
LDA_plot(LDA_tab = BT20,uncertainty.band = T)
```
In case the figure looks overcrowded, `uncertainty.band = FALSE` can switch off
the confidence bands and `xlim` (e.g. `xlim = c(0,100)`) can be used, to adjust
the plotted range of seeded cells.
---
# Mathematical background
## Introduction
Clonogenicity refers to the capacity of single cells to grow into
new colonies.
Since not every single cell will successfully form a new colony,
in a given culture the fraction of cells that do so is often of interest.
This fraction is then called clonogenic activity.
Mathematically it is a probability $p$.
Often the activity is communicated by the reciprocal, which corresponds to a
number of cells (e.g. activity $p = 1/42$ written as $a = 1/p = 42$ for
'one active cell among $a = 42$ non-active cells').
In the LDA, distinct numbers of cells are seeded in wells. The number of
seeded cells per well is often called dose.
Since, we do not want to get confused with the irradiation dose, which is a
frequent treatment in the field of radiation oncology, we denote the number of
cells seeded with $S$.
The readout of each well in the LDA is dichotomous (growth / no growth).
So, from the binomially distributed number of active cells in a well $X$ with
\[P(X=k) = \left( \begin{matrix} s \\ k \end{matrix}\right)
p^{k} (1-p)^{S-k} \]
only $X=0$ and $X \neq 0$ can be observed experimentally.
Thus, it is about the frequency of positive wells as a function of the
number of cells seeded ($k/n = f(S)$).
For dealing with high cell numbers, the binomial distribution is approximated
by the
asymptotically identical (for $S \to \infty$) Poisson distribution (with
$\lambda = S \cdot p$):
\[P_{\lambda}(X=k) = \frac{\lambda^k}{k!} e^{-\lambda}.\]
The basic concept for estimating the clonogenic activity of a cell suspension
through diluting the number of cells in different wells is the following:
- imagine a number of cells in a well, which is so high
that we expect many active cells among them, then this well will be positive.
- imagine a number of cells so low
that the chance of an active cell among them is very low, then this well
will be negative.
- thus, between those two numbers there will be cell numbers, where no-one
can say, whether a colony will grow or not. It is rather a matter of chance.
- from Poisson statistics we know, that a frequency of 37% negative wells will
be observed, when the average number of active cells per well is one
($e^{-1} \approx 0.368$).
So taking the graph 'observed fraction of positive wells over the number
of cells per well', we find the clonogenic activity as the number of cells
(x-axis) just where the curve passes the expectation of 37% negative wells.
## Mathematical model
Let $n$ be the number of wells, and let $\mu$ denote the expected fraction of
negative wells (failures). Following the single-hit Poisson model
(SHPM: $\lambda_{SHPM} = p \cdot S$), the expectation of active cells in a well
is approximately Poisson
distributed and depends linearly on the number of cells seeded $S$ as well as
the probability $p$ of each cell being active (i.e., clonogenic).
We therefore find $P_{\lambda}(k=0) = e^{-p S}$), and thus:
\[\mu = e^{-pS}.\]
The number of failures $r$ is binomially distributed
\[ P(Y=r) = \left( \begin{matrix} n \\ r \end{matrix}\right)
\mu^{r} (1-\mu)^{n-r}. \]
We find (for $\alpha = ln(p)$, $\xi = ln(S)$)
\begin{align*}
\Leftrightarrow \mu &= e^{- p S}\\
\Leftrightarrow ln(\mu) &= - p S \\
\Leftrightarrow -ln(\mu) &= p \cdot S \\
\Leftrightarrow ln(-ln(\mu)) &= \alpha + \xi
\end{align*}
Therefore, technically the identification of the required number of cells is
done via a generalized linear regression model
(binomial family, loglog link function).
The clonogenic activity is $p = e^{\alpha}$.
### Cooperativity
In the presented model, $p$ is independent from $S$.
So implicitly, it is assumed, that 100 single cells in single wells have the
identical chance to result in at least one colony as a single well with those
100 cells would have.
But this is often not the case.
Cells export or release substances into the medium and thereby
communicate and cooperate (or compete).
The combined chance of a group of cells to grow into a clone often exceeds the
sum of the individual chances.
(See CFAcoop, for a number of cells seeded $S$ and a number of resulting
colonies $C$ we find $C = a \cdot S^b$.)
Thus, allowing for the same communication and cooperativity (or competition),
we replace
$\lambda = p \cdot S$ with $\lambda = p \cdot S^b$ analogue to CFA.
The modification reads in both cases as 'a number of $S$ cells shows the same
activity as if they were $S^b$ single cells'.
With this generalization we find
\begin{align*}
P_{\lambda}(k=0) &= e^{-\lambda}\\
\Leftrightarrow \mu &= e^{- p S^b}\\
\Leftrightarrow ln(\mu) &= - p S^b \\
\Leftrightarrow ln(-ln(\mu)) &= \alpha + \xi \cdot b.
\end{align*}
Interestingly, this equation is the same as in the special case of
non-cooperativity, just without the restriction of $b=1$, which corresponds to
the non-cooperative case.
The degree of cooperativity is quantified by the coefficient $b$.
## Activity
When the probability of a cell to grow into a colony is not
independent from the number of surrounding cells, the earlier definition of
clonogenic activity (the chance of an isolated single cell) is pointless.
Thus, the approach to deal with cooperatively growing cells, requires a
generalized definition for clonogenic activity.
For a fixed outcome of e.g. $\lambda = 1$, which is equivalent with
$\mu \approx 0.63$, we still may calculate the required number of cells to be
seeded from the regression coefficients
\[s_{\lambda = 1} = e^{\frac{-\alpha}{b}} = \left(p^{-1}\right)^{\frac{1}{b}}.\]
On average, one out of these cells will grow into a colony.
Therefore, we define this number $s_{\lambda=1}$ as the
(generalized) clonogenic activity.
This is analogue to the CFA approach for cellular cooperation
(which sets the reference at an expectation of $20$ colonies.).
Please note, the non-cooperative case ($b=1$)
results in the identical activity values as from the previous definition.
### Example
Given the data of an LDA experiment, we find the clonogenic activity where
expected. The only change in contrast to conventional LDA data analysis is
the curve, which replaces the linear model.
```{r single.first,fig.width=6, fig.height=3}
cell.line <- unique(LDAdata$name)[2]
AD <- subset.data.frame(LDAdata, subset = (name==cell.line) &
(replicate==1) & (Group == 2))[,4:6]
LDA_plot(LDA_tab = AD,uncertainty = "act", uncertainty.band = T)
LDA_table(x = AD)
full_model_fit <- LDA_activity_single(x = AD)
```
In case, there are biological replicates, the replicate may be indicated by
the fifth column. Thereby the frequency of responding wells can be plotted
separately.
```{r set1, fig.width=6, fig.height=3}
AD <- subset.data.frame(LDAdata, subset = (name==cell.line) &
(Group == 2))[,c(4:6,3,2)]
LDA_plot(LDA_tab = AD,uncertainty = "act", uncertainty.band = T)
LDA_table(x = AD[,1:3],ref_class = 0)
```
Please note that the starting motivation of a concept 'one active cell among
how many?' is still the same. It is just restricted to the special outcome of
37% negative wells.
The choice of 'outcome of 37% negative wells' is quite artificial and the
activity for higher or lower active well ratios would change.
But when it comes to survival fractions, this shift is mostly cut out and
thereby this approach filters (in parts) a cooperativity-bias in those
survival fractions.
## Survival fraction
Even though the concept of clonogenic activity is somewhat fuzzy in the presence
of cellular cooperation, the effect of treatments on this clonogenic capacity
can be investigated quite accurately by the same way it is done for the CFA.
The only requirement is a statement on the outcome of interest. At CFA, a
number of 20 colonies can be agreed upon - and a comparison of the number of
required cells seeded with and without treatment is reasonable.
Analogously, the ratio of cell numbers required to observe the same $\mu$
(e.g. $\mu = e^{-1}$) is of interest, when investigating the effect of
certain treatments.
\[SF_x = \frac{S_0(\lambda_0 = 1)}{S_x(\lambda_x = 1)}\]
can be calculated from the fitted GLM parameters as
\[SF_x = exp\left( \frac{log(p_x)}{b_x} - \frac{log(p_0)}{b_0}\right)\]
Therefore, the interpretation of survival fractions is straight forward.
If the number of cells needed to be seeded with treatment is ten times the
number without treatment, the survival fraction is calculated as $0.1$.
Note that it is implicitly assumed, that those cells that are non-survivors
of the treatment, do not contribute to the cooperative stimulation of the
survivors. Otherwise the treatment effect will be underestimated.
## Uncertainty analysis
Uncertainties for the parameters of the generalized linear model are returned
by the fitting procedure and can be transferred to the clonogenic activity
via the general `predict.glm` function.
Thus, uncertainties of clonogenic activities are calculated the same way
as the activities are (cutting $y=-1$).
Uncertainties of the survival fractions can be generally assessed by two ways
- robust approximation by combination of the 84%-confidence-intervals of
the activity (method (a), see Austin et al. (2002), Payton et al. (2003),
Knol et al. (2011))
- calculation following the law of error propagation through first order Taylor
series expansion (method (b))
In cases of extreme cooperativity, method (b) can be numerically unstable and
therefore we recommend the use of method (a) instead.
Analysis of the calculated uncertainty bounds of the survival fractions
of the 10 non-extreme cell lines under all treatments
(see data delivered with this package),
we find a mean ratio (uncertainty-bound method(a) / uncertainty-bound method(b))
of 0.9969 and a standard deviation of 0.0049.
Thus, the deviation of those two methods is very stable in the order of 1%.
### Example
```{r BT20ep}
LDA_table_act <- LDA_table(x = BT20,uncertainty = "act")
cbind(round(LDA_table_act[,1:5],digits = 1),
round(LDA_table_act[,6:9],digits = 4))
LDA_table_ep <- LDA_table(x = BT20,uncertainty = "ep")
cbind(round(LDA_table_ep[,1:5],digits = 1),
round(LDA_table_ep[,6:9],digits = 4))
```