---
title: "flowPloidy: Flow Cytometry Histograms"
author: "Tyler Smith"
date: "2017-03-15"
output: 
    rmarkdown::html_vignette:
        fig_width: 6
        fig_height: 4
    pdf_document:
        fig_width: 6
        fig_height: 4
        fig_caption: false
        template: pandoc/fpPandoc.latex
vignette: >
  %\VignetteIndexEntry{flowPloidy: FCM Histograms}
  %\VignetteEngine{knitr::rmarkdown}
bibliography: pandoc/flowPloidy.bibtex
---

# Introduction

The application of flow cytometry (FCM) to ploidy assessment is based on a
simple concept. Cell nuclei are stained with DNA-binding fluorochromes, and
the amount of DNA present is determined by measuring the fluorescence
emitted when the nuclei are excited with light of a particular wavelength.
Greater fluorescence means more DNA in the nucleus. That is, a diploid
nucleus should produce half the fluorescence of a tetraploid nucleus, and a
hexaploid nucleus should produce three times more fluorescence than the
diploid.

# Histogram Basics

Applying this concept is complicated by several factors:

* nuclear DNA content varies over the cell cycle
* secondary metabolites may interfere with fluorochrome staining, or may
  themselves fluoresce
* tissue preparations are usually contaminated with non-nuclear cellular
  components, damaged nuclei, and other debris that fluoresce to greater or
  lesser extent
* the accuracy of flow cytometers is limited, and values tend to shift or
  drift over the course of a run
      
We can illustrate these challenges by comparing an ideal FCM DNA histogram
with simulated experimental data. 

```{r ideal histogram, echo = FALSE}
plot(x = 1, type = 'n', xlim = c(0, 256), ylim = c(0, 500),
     ylab = "Nuclei", xlab = "Fluorescence", main = "Ideal Histogram")
abline(h = 0)
lines(x = c(80, 80), y = c(0, 325))
lines(x = c(160, 160), y = c(0, 150))
tmpfun <- function(x) 100 * exp(-x/80)
curve(tmpfun, 80, 160, add = TRUE)
text(x = 80, y = 325, "G1 Peak", pos = 3)
text(x = 160, y = 150, "G2 Peak", pos = 3)
text(x = 120, y = 30, "S Phase", pos = 3)
```

In the ideal case, we have two clear, accurate cell groups, G1 and G2. The
G1 ("gap 1") group represents cells in their 'normal' state. That is, they
have the usual diploid DNA complement. The G2 ("gap 2") group
represents cells that have duplicated their DNA, but haven't yet undergone
mitosis. Consequently, they have twice the diploid DNA complement. In
between the two groups is the S phase ("synthesis") region. This area
of the plot represents cells that are in the process of duplicating their
DNA, in preparation for mitosis.

(N.B. in the flow cytometry literature, clusters of cells, which I've
called "groups" above, are often referred to as cell "populations") 

If our histogram were this clean, we could easily determine the DNA content
of our sample by comparing it to a known standard (see below). In practice,
you'll never see a histogram that looks like this. A more realistic example
follows:

```{r real histogram, echo = FALSE}
set.seed(1234)
plot(x = 1, type = 'n', xlim = c(0, 256), ylim = c(0, 500),
     ylab = "Nuclei", xlab = "Fluorescence", main = "Empirical Histogram")
abline(h = 0)

vals <- rep(0, 256)
vals <- jitter(3500 * dnorm(0:256, mean = 80, sd = 80 * 0.05), amount = 20)
vals <- vals + jitter(200 * exp ((0:-256)/60), amount = 5)
vals <- vals +
  jitter(2500 * dnorm(0:256, mean = 160, sd = 160 * 0.05), amount = 20)
vals[80:160] <- vals[80:160] + jitter(100 * exp((-80:-160)/80),
                                      amount = 5)
vals[vals < 0] <- 0
points(vals, type = 'l')

text(x = 80, y = 450, "G1 Peak", pos = 3)
text(x = 160, y = 200, "G2 Peak", pos = 3)
text(x = 120, y = 100, "S Phase", pos = 3)
text(x = 40, y = 170, "Debris", pos = 3)
text(x = 240, y = 40, "Aggregates", pos = 3)
```

We can still see the G1 and G2 groups, but they aren't exact values.
Instead we see that the fluorescence values for each group are distributed
around a central peak. Similarly, the S phase is not a clean polygon. The
jagged edges of the region are due in part to instrument error: we can't
capture DNA content perfectly. In addition, some of that noise is "real".
That is, the cells in our tissue sample are not proceeding through the cell
cycle in a perfectly regular manner. Even if we could measure the DNA
without error, the S-phase region wouldn't necessarily be a perfectly
smooth region in the histogram.

In addition to the G1, G2, and S phase regions from our ideal histogram,
real histograms also contain debris. This is a mixture of cellular debris,
damaged nuclei, and other contaminants. The dyes we use bind best to DNA,
so most of the debris, which contains relatively little DNA, is found on
the left/lower fluorescence side of the histogram. We may also see
aggregates at the right side of the plot. These are nuclei that have stuck
together in the flow cytometer, and so appear as higher-ploidy nuclei.

# Histogram Analysis

Given that we have a noisy histogram, how do we extract the values of
interest from the data? There are two main approaches. The first, perhaps
more intuitive, approach requires us to manually draw boxes around the
peaks we're interested in, and have the computer calculate the parameters
of interest for us. For the previous example, that would look like this:

```{r manual analysis, echo = FALSE}
set.seed(1234)
plot(x = 1, type = 'n', xlim = c(0, 256), ylim = c(0, 500),
     ylab = "Nuclei", xlab = "Fluorescence",
     main = "Manual Histogram Analysis")
abline(h = 0)

vals <- rep(0, 256)
vals <- jitter(3500 * dnorm(0:256, mean = 80, sd = 80 * 0.05), amount = 20)
vals <- vals + jitter(200 * exp ((0:-256)/60), amount = 5)
vals <- vals +
  jitter(2500 * dnorm(0:256, mean = 160, sd = 160 * 0.05), amount = 20)
vals[80:160] <- vals[80:160] + jitter(100 * exp((-80:-160)/80),
                                      amount = 5)
vals[vals < 0] <- 0
points(vals, type = 'l')
rect(70, 0, 95, 475, border = 2, lwd = 2)
```

There are two main problems with this approach. First, the placement of the
box around a peak is subjective. One of the most important quality-checks
for ploidy analysis is the coefficient of variation (CV) for a peak. We can
easily lower the CV just by drawing a narrower box around our data.
Consequently, our confidence in our results is undermined by not having an
objective and repeatable method for setting the limits around a peak.

The second problem is that we are not accounting for the different types of
cells that contribute to each peak. What we see as the G1 peak also
contains debris and S-phase nuclei:

```{r individual components, echo = FALSE}
set.seed(1234)
plot(x = 1, type = 'n', xlim = c(0, 256), ylim = c(0, 500),
     ylab = "Nuclei", xlab = "Fluorescence",
     main = "Histogram Components")
abline(h = 0)

vals <- rep(0, 257)
sphase <- vals
G1 <- jitter(3500 * dnorm(0:256, mean = 80, sd = 80 * 0.05), amount = 20)
G1[G1<0] <- 0
vals <- G1
debris <- jitter(200 * exp ((0:-256)/60), amount = 5)
debris[debris<0] <- 0
vals <- vals + debris
G2 <- jitter(2500 * dnorm(0:256, mean = 160, sd = 160 * 0.05), amount = 20)
G2[G2<0] <- 0
vals <- vals + G2
sphase[80:160] <- jitter(100 * exp((-80:-160)/80), amount = 5)
sphase[sphase<0] <- 0
vals <- vals + sphase
vals[vals < 0] <- 0
points(vals, type = 'l')
G1[c(1:70, 90:256)] <- 0
points(debris, type = 'l', col = 3, lwd = 2)
points(G1, type = 'l', col = 4, lwd = 2)
sphase[c(1:70, 170:256)] <- 0
points(sphase, type = 'l', col = 2, lwd = 2)
text(x = 25, y = 50, col = 3, "Debris")
text(x = 95, y = 350, col = 4, "G1")
text(x = 120, y = 120, col = 2, "S-Phase")
```

When we draw a box around the G1 peak, we are also including debris and
s-phase cells, which distorts the true distribution of the G1 nuclei DNA
content.

The second approach to FCM histogram analysis uses non-linear regression to
distinguish the different cell populations that contribute to the total
histogram. The G1 and G2 peaks are modeled as normal curves. A variety of
phenomenological and mechanistic functions are available to model the
debris and S-phase components. See @bagwell_1993 and @watson_1992 for a
more detailed presentation of the modeling process.

Using the modeling approach, the analysis of our example data looks like
this:

```{r histogram modeling, echo = FALSE}
set.seed(1234)
plot(x = 1, type = 'n', xlim = c(0, 256), ylim = c(0, 500),
     ylab = "Nuclei", xlab = "Fluorescence",
     main = "Histogram Modeling")
abline(h = 0)

vals <- rep(0, 257)
sphase <- vals
G1 <- jitter(3500 * dnorm(0:256, mean = 80, sd = 80 * 0.05), amount = 20)
G1[G1<0] <- 0
vals <- G1
debris <- jitter(200 * exp ((0:-256)/60), amount = 5)
debris[debris<0] <- 0
vals <- vals + debris
G2 <- jitter(2500 * dnorm(0:256, mean = 160, sd = 160 * 0.05), amount = 20)
G2[G2<0] <- 0
vals <- vals + G2
sphase[80:160] <- jitter(100 * exp((-80:-160)/80), amount = 5)
sphase[sphase<0] <- 0
vals <- vals + sphase
sphase[80:160] <- 100 * exp((-80:-160)/80)
vals[vals < 0] <- 0
points(vals, type = 'l')
points(200 * exp ((0:-256)/60), type = 'l', col = 3, lwd = 2)
points(sphase, type = 'l', col = 2, lwd = 2)
points(2500 * dnorm(0:256, mean = 160, sd = 160 * 0.05), type = 'l',
       lwd = 2, col = 5)
points(3500 * dnorm(0:256, mean = 80, sd = 80 * 0.05), type = 'l',
       col = 4, lwd = 2)
text(x = 25, y = 50, col = 3, "Debris")
text(x = 95, y = 350, col = 4, "G1")
text(x = 120, y = 120, col = 2, "S-Phase")
text(x = 180, y = 100, col = 5, "G2")
```

This is a simulated example, but it illustrates the modeling process very
clearly. We can extract the values of interest directly from the model
components. The regression results include the mean value for the G1 peak
(80), as well as its CV (5%). These values have taken into account the
contribution of the debris curve and the s-phase nuclei, and we'll get the
same results every time - they don't depend on where we decide to draw the
box around each peak.

# Standards

To simplify the presentation of histograms in the previous sections, I
omitted the internal standards used in FCM analysis. This is a critical
component of ploidy analysis, because the fluorescence values for a given
quantity of DNA are not fixed. The measured values vary based on the tissue
preparation and the flow cytometer. Even for the same tissue prep run on
the same flow cytometer, the fluorescence values can vary over the course
of a single day.

If we were to re-run the sample that yielded a mean value of 80 in our
example, we might get a different result. Which means we can't interpret
the raw fluorescence values. Instead, we include a tissue standard in our
preparations, so we can compare our sample peak to the known value. While
the absolute fluorescence values may shift over time, the *ratio* between
our sample mean and the standard mean will stay the same.

Adding a standard to our example data, we get a histogram that looks like
this:

```{r histogram standard, echo = FALSE}
set.seed(1234)
plot(x = 1, type = 'n', xlim = c(0, 256), ylim = c(0, 500),
     ylab = "Nuclei", xlab = "Fluorescence",
     main = "Histogram with Standard")
abline(h = 0)

vals <- rep(0, 257)
sphase <- vals
G1 <- jitter(3500 * dnorm(0:256, mean = 80, sd = 80 * 0.05), amount = 20)
G1[G1<0] <- 0
vals <- G1
debris <- jitter(200 * exp ((0:-256)/60), amount = 5)
debris[debris<0] <- 0
vals <- vals + debris
G2 <- jitter(2500 * dnorm(0:256, mean = 160, sd = 160 * 0.05), amount = 20)
G2[G2<0] <- 0
vals <- vals + G2
sphase[80:160] <- jitter(100 * exp((-80:-160)/80), amount = 5)
sphase[sphase<0] <- 0
standard <- jitter(1500 * dnorm(0:256, mean = 130, sd = 130 * 0.01),
                   amount = 10)
standard[standard<0] <- 0
vals <- vals + sphase
vals <- vals + standard
sphase[80:160] <- 100 * exp((-80:-160)/80)
vals[vals < 0] <- 0
points(vals, type = 'l')
text(x = 160, y = 450, "Standard")
text(x = 95, y = 400, "G1")
```

To analyze this histogram, we add another normal curve to account for the
standard sample. That allows us to compare our sample mean (80) to the
standard mean (130). If the standard was *Solanum lycopersicum*, with a
known 2C value of 1.96pg, then our sample has a 2C value of $\frac{80}{130}
\times 1.96pg = 1.21pg$.

# References